Taking a matrix to a very large power.

So I'm interested in taking a matrix (actually an adjacency matrix of 1s and 0s, typically of size ~ 1600 x 1600 and quite sparse) to as large a power as I possibly can. I was actually very surprised how well simply mpower(A,n) works to do this - I've been able to get results for a 1600 x 1600 matrix to the power of 200 within a second (don't know how accurate this is likely to be but I'm really only interested in the relative values of elements and the results seem reasonable within the context of what I'm doing). Basically the problem I'm running into is that I would like to go to higher powers, but the numbers being produced are so large that MATLAB is returning Inf for all elements.

As I've mentioned I'm only actually interested in the relative sizes of elements, so I'm wondering if there's some trick (I suspect involving taking logarithms somewhere) that would let me go to even higher matrix powers?

Answer 1

[John D'Errico]

Actually, taking the log of your elements will not work, since raising a matrix to some power requires addition too. If you were raising the individual elements, then the log would apply.

The problem with raising a matrix to some large power is you will quickly overflow the dynamic range of double precision. mpower does as well as you can, but it will fail at some point. What matters are the eigenvalues of the matrix, as it is those eigenvalues that will be effectively raised to a power.

You do have some options, if you are looking for a REALLY high power. If that power is high enough, then all that matters is the very largest eigenvalue, as long as it is significantly distinct from the remaining eigenvalues.

First, you want to compute the eigenvalues and eigenvectors. Only the largest few such eigenvalues matter here.

I'll create a moderately sparse matrix of 0's and 1's as an example.

A = double(rand(1600) > 0.01);

[V,D] = eig(A);
values = diag(D);
values(1:10)
ans =
         1583.96955854307 +                     0i
         4.00155305448781 +     0.559803131102432i
         4.00155305448781 -     0.559803131102432i
         4.01385193091935 +                     0i
         3.60025510613369 +      1.66146362049061i
         3.60025510613369 -      1.66146362049061i
         3.25068698344655 +      2.25060776388857i
         3.25068698344655 -      2.25060776388857i
         3.10943173157108 +      2.44959623206898i
         3.10943173157108 -      2.44959623206898i

Now, if we know that A can be written as

A = V*D*V'

where V is orthogonal, then

A^2 = V*D*V' * V*D*V'

So

A^2 = V*D^2*V'

Larger powers are as easy.

A^N = V*D^N*V'

But, remember that D^N will be dominated by that largest eigenvalue. In this case, that largest eigenvalue is nearly 400 times as large as the remainder of the eigenvalues.

r = abs(D(1,1))/abs(D(2,2))
r =
       392.02

Since raising a diagonal matrix (D) to a power merely raises the diagonal elements to that power, all we care about is the point where all other terms are less than eps times the first.

log(eps)/log(r)
ans =
      -6.0361

So if we raise this matrix to the 7th power or larger, then all other terms are irrelevant! We can ignore them completely, as long as we work in double precision.

D1 = values(1);
V1 = V(:,1);

So now, as long as that power is sufficiently large (here, N >= 7) then A^N will be given as (remember that D1 is a scalar)

V*V1'*D1^N

Again, D1 is a SCALAR! So, if N is so large that you would get overflow, than all that matters is you multiply the entire matrix by some HUGE constant! If N is too large here,

log(realmax)/log(D1)
ans =
       96.337

So in this case, if N is larger than roughly 96, we will see overflows. Everything turns into inf around that value of N.

Ap = A^96;
Ap(1,1)
ans =
    9.3846e+303

Ap = A^97;
Ap(1,1)
ans =
    1.4865e+307

Ap = A^98;
Ap(1,1)
ans =
    Inf

We cannot go past that point. However, you don't really care in double precision. Of course, IF you recognize that

A^N = (V1*V1')*D1^N

then you should see that larger values of N just multiply the simple rank 1 matrix V1*V1' by some huge scalar value.

All of this is valid as long as you raise your matrix to a sufficiently high power. How high depends on the ratio of the magnitudes of the two largest eigenvalues.

[Henry Charlesworth]

Thanks this is a really useful answer. I think I have another similar but simpler solution which is just to just take the matrix to some smaller power where you don't get overflow, then divide this by some big number to make all the elements relatively small before multiplying it by another matrix power - e.g. to do A^N I would do (A^n / somenumber)*A^(N-n) (or maybe renormalize multiple times), and because I only care about the relative sizes of elements of A^N this should be OK.

[John D'Errico]

It all depends on what you want out the end. The continuously renormalizing approach is viable, but it has its own issues you will need to deal with.

Answer 2

[Star Strider]

Providing, of course, that all the matrix elements are greater than zero, the log would likely work. See the funm function for a list of all the matrix functions. It should make your coding them easier.

I suspect mpower and the others are based on the Cayley-Hamilton theorem and its corollaries. If you have access to T-S Chen Linear System Theory and Design, 3d Ed., Oxford 1999, see P 63-5 for a concise, thorough discussion.