Hi! How to decrease the computation time of calculating eigenvalues and eigenvectors of large sparse matrices like 10000X10000(I need all the eigenvalues so I can't use eigs)?

I need to find all the eigenvalues of a large matrix which is 10000X10000, and I have to loop it for some 1000 times.My computer has sufficient memory, the problem is it's taking too much time. It is taking 1min or so to calculate the eigenvalues for one time. I cannot use parallel computing because I will go to the next loop based on the eigenvalues I got in the previous loop. But I know the nature of the eigenvalues that I get. The eigenvalues are symmetric(This may not be the proper term) what I mean is for example 2, -2, 1, -1, 2+3i, -2-3i, 3.5, -3.5. So is there any way I could make use of this information to run my code faster?

Thanks in advance.

[Matt J]

It might be possible, if we knew what kind of matrix structure guarantees the kind of eigenvalue symmetry that you describe. How do you know mathematically that the matrix has this property?

[Yaswanth Sai]

Hi Matt, I have already run the code several times but it took more time, as I have already said, and the result is what I put was just that result. Only special property of the matrix that I know is the half of the matrix contains zeros and the their positions are upper right rectangle and lower left rectangle. Something like the one below:

A = [0 0 2 3;
     0 0 5 8;
     4 5 0 0;
     8 9 0 0]; % A = [O X;Y O]

[Matt J]

How does that guarantee symmetric eigenvalues?

[Yaswanth Sai]

I don't know about that. I have mentioned that just because that's the only property that the matrix contains, as far as I know. So I wanted to know whether that can be used somehow to reduce my computation time.

Answer 1

[John D'Errico]

The simple answer is, I doubt it. But that might not be the complete answer, or even the correct one. And see that any answer that requires any significant amount of computation on the matrix will probably slow down things enough that just calling eig on the complete matrix will be faster. Regardless, it is worth a bit of thought.

But your question reduces to computing the roots of the characteristic polynomial, IF you knew that all of the roots were duplicates, but with a sign difference. So if lambda is a root, then so is -lambda.

In that case, you could write your polynomial in the form

(x - lambda1)*(x + lambda1)*(x - lambda2)*(x + lambda2)* ...

Of course, this reduces to

(x^1 - lambda1^2)*(x^2 - lambda2^2)* ...

Essentially your polynomial will have ALL even powers of x, and ONLY even powers of x. If there are any odd powers at all with non-zero coefficients, then your claim about the roots fails completely.

But if you do have only those even powers, then there is a simple transformation available. Just use y=x^2. So create a new polynomial where you divide all the even powers by 2. This new polynomial will be of lower degree than the original, just as if you were able to turn a 10kx10k matrix into a 5kx5k matrix. Compute the roots, and take the square root, apply a +/- to that root, and you are done.

Note this is not equivalent to computing the eigenvalues of A*A. See that A*A will not only be far less sparse than is A, but A*A will still be of the same size as is A. If A has the desired property of signed symmetric eigenvalues, then A*A will have eigenvalues of multiplicity 2.

The problem is, knowing only that the characteristic polynomial must be of the form that it only has even powers of lambda is not sufficient to determine the patterning of your matrix, or for me to know how to reduce your matrix to one of half the size. This must be the case, since there are lots of matrices which will result in the same characteristic polynomial. However, it may be the case that you know some factor of importance here that you are not telling us. For example, how is it that you know your matrix has eigenvalues with the property that you claim it does?

[the cyclist]

I don't know the theory here, but empirically I can verify that "block antidiagonal" (not sure if that is the correct term) of the form

Z = zeros(2);

B1 = rand(2);
B2 = rand(2);

A = [Z, B1; B2, Z];

"always" have opposite-signed paired eigenvalues, as OP claims.

Here "always" means "I did it enough times that I am convinced it is always true".

[John D'Errico]

Ok. I was just coming to the same conclusion. Matrices with that pattern will indeed have the property that the characteristic polynomial has only even powers. And that in turn will result in eigenvalues that are paired with opposite signs, since we can write the characteristic polynomial as I did. So now the pertinent question is, can one find a reduction for matrices with that shape? I.e., one that is equivalent to finding the eigenvalues of an n/2xn/2 matrix? And, one needs to do so without taking more composite time than finding the eigenvalues of the nxn matrix in the first place?

[Yaswanth Sai]

Hi John, the symmetric eigenvalues is a result which I found after I have run my code through many cases. This result matches with the physical problem that I am dealing with. I am working on waves and the eigenvalues in my case are wave-numbers, the values should be the way they are as the waves can propagate in either direction. Apart from this I don't know much about the matrix. And coming to the method of reducing the matrix, how exactly am I going to do the transformation using 'eig' as I am pretty sure it would take more time than before if I code to find the roots of the characteristic polynomial instead of using eig. One more thing, is there any way I can calculate only the real eigenvalues(not real parts) using eigs? Thank you.

[John D'Errico]

Oh. Of course. Actually, the answer seems simple in hindsight, once I started to do some reading about the determinants of blocked matrices. Note that until I saw your comment that the matrix was of anti-block diagonal form, I would probably not have come to these conclusions as I did.

Given a matrix with anti-block diagonal 2x2 form, so a matrix of the form

A = [0,B;C,0]

where 0 is interpreted as a matrix of all zeros of compatible size with the square matrices B and C. The characteristic polynomial of A will have only even powers. (Proof deferred, but not difficult I think.)

This clearly implies the eigenvalues will have the claimed property, of being signed duplicates.

As well, the eigenvalues of A can be computed directly from B and C from eig(B*C). Again, proof deferred. I'm being lazy here, and it is not my thesis to write. ;-) For example...

A1 = rand(5);
A2 = rand(5);
Z = zeros(5);
A = [Z,A1;A2,Z];

eig(A)
ans =
          -2.34529174173102 +                     0i
           2.34529174173102 +                     0i
       -1.2490009027033e-16 +      0.53779419573051i
       -1.2490009027033e-16 -      0.53779419573051i
         -0.479030555460411 +                     0i
          0.479030555460411 +                     0i
         -0.207603454338438 +    0.0921646900938415i
         -0.207603454338438 -    0.0921646900938415i
          0.207603454338439 +    0.0921646900938414i
          0.207603454338439 -    0.0921646900938414i

eig(A1*A2)
ans =
           5.50039335383173 +                     0i
         -0.289222596961427 +                     0i
           0.22947027306471 +                     0i
         0.0346048641531584 +    0.0382674160630264i
         0.0346048641531584 -    0.0382674160630264i

sqrt(eig(A1*A2))
ans =
           2.34529174173102 +                     0i
                          0 +     0.537794195730511i
          0.479030555460411 +                     0i
          0.207603454338439 +    0.0921646900938416i
          0.207603454338439 -    0.0921646900938416i

As you can see, the square roots of the eigenvalues of A1*A2 all appear as the eigenvalues of A, with both signs observed.

Now with some thought, the above claims I have left unproved can be shown to follow.

So, IF your goal is to find only the eigenvalues of the anti-block-diagonal matrix in question, your problem does have a solution, in the sense that eig will be significantly faster on a 5kx5k matrix than it will be on a 10kx10k matrix) even with the cost of a matrix multiply thrown in.

If you also need the eigenvectors, I need some aspirin.

[Yaswanth Sai]

Thanks a lot John, this method would actually reduces the computation time a lot. By the way I actually need eigenvectors, but I think I would be able to figure it out someway using 'eigs'. Thanks again.

[John D'Errico]

As a starter on proving that the eigenvalues of an ant-blockdiagonal matrix has sign paired eigenvalues, we need only one simple result. That is, the determinant of the 2x2 block matrix (with square blocks) can be computed as

det([A,B;C,D]) = det(A*D - B*C)

which applies IF the product C*D commutes, thus C*D=D*C. (Of course, proof of this claim also would be required. Yeah, I know, I'm a lazy guy, but then, I finished school a million years ago and I'm retired, so I can afford it.)

How does this help? First, what is the characteristic polynomial of the anti-block diagonal matrix? For nxn sub-matrices, it can be written as

A = [0,B;C,0]
det(A-lambda*I(2*n)) = det(lambda^2*I(n) - C*B)

here I(n) is the nxn identity matrix, and I(2*n) is the 2nx2n identity matrix.

Again, by the claim above, this applies as long as C commutes with an identity matrix, which clearly it does. That final form for the determinant clearly implies the characteristic polynomial will have only even powers of lambda, which in turn implies the signed pairs of eigenvalues. It also tells us how to compute the desired eigenvalues, as sqrt(eig(C*B)).

I'll let you dot the i's, but this should get you close enough.

I don't think there is anything more to be done with this though. However much the speed boost is, you will need to accept it. Even parallel processing will not get you more of a gain in speed, since eig will already use multi-threaded computation for a problem that large. The only way to get more of a gain will then be to use a faster computer.

[Christine Tobler]

Nice solution, John! By the way, we can get the eigenvectors of A from the eigenvectors of B*C, too. Think of the eigenvalue problem in A as two coupled problems:

 B*y = lambda*x
 C*x = lambda*y

By using the substitution y = C*x*(1/lambda), you get the eigenvalue problem B*C*x = lambda^2*x, and can compute eigenvector y directly from eigenvector x. Recombining them gives the eigenvectors of A:

 [X, D2] = eig(B*C);
 D = diag(sqrt(diag(D2)));
 Y = C*X/D;
 % C*X == Y*D
 % B*Y == X*D

 U = [X X; Y -Y];
 DA = blkdiag(D, -D);
 % A*U == U*DA

[John D'Errico]

Thanks Christine. Neat.

[Yaswanth Sai]

This is a kind of follow-up to the previous question, but this time the matrix is of a different form but still gives eigenvalues that are symmetric. It's hard to explain the pattern. So I will give you an example.

a =

     0     0     0     0     1     1     0     1
     0     0     0     0     1     1     0     1
     0     0     0     0     0     0     1     0
     0     0     0     0     1     1     0     1
     1     1     0     1     0     0     1     0
     1     1     0     1     0     0     1     0
     0     0     1     0     1     1     0     1
     1     1     0     1     0     0     1     0

The ones(1s) just denote the positions of the elements with a non zero value. While observing the matrices please look at each block separately i.e, the top left block should always be a zero matrix and if we take the top left block you can observe a pattern(observe the columns and their repeatability). Let's take a random matrix k of the above form

k = rand(8).*a

k =

         0         0         0         0    0.2954    0.1930         0    0.4316
         0         0         0         0    0.3661    0.9196         0    0.7179
         0         0         0         0         0         0    0.2182         0
         0         0         0         0    0.6302    0.5509         0    0.8900
    0.0685    0.5280         0    0.3843         0         0    0.4340         0
    0.2052    0.1851         0    0.3996         0         0    0.7848         0
         0         0    0.1693         0    0.9444    0.2577         0    0.8935
    0.5751    0.4641         0    0.5554         0         0    0.3313         0
eig(k)

ans =

   1.7579 + 0.0000i
  -1.7579 + 0.0000i
   0.0000 + 0.4743i
   0.0000 - 0.4743i
   0.2006 + 0.0000i
  -0.2006 + 0.0000i
   0.0304 + 0.0000i
  -0.0304 + 0.0000i

So the eigenvalues behavior is similar to that of anti diagonal block matrices. Is there any way I could manipulate the matrices, reduce their size and get the required eigenvalues?

[John D'Errico]

I'd do some reading here and here to start.

[Yaswanth Sai]

Thank you for the two sources, John, I will go through them.

[John D'Errico]

The answer is not obvious to me at the moment, but I have a funny feeling it will be clear at some point. Perhaps the approach Christine took would get you there? Thus, given the matrix

A = [0,B;C,D]

Assume a solution exists of the form

A*[x;y] = lambda*[x;y]

Then we have two equations:

B*y = lambda*x
C*x + D*y = lambda*y

For any non-zero scalar lambda, we can now eliminate x simply as

x = B*y/lambda

So for non-zero lambda, the problem reduces to an expression that looks vaguely like a generalized eigenvalue problem.

(C*B)*y = (lambda^2*I - lambda*D)*y

where I is an nxn identity of the same size as D. Of course, when D was the zero matrix, things got simple before.

While this looks promising, it has one major problem. A standard trick is to complete the square on the right hand side. So, add D*D/4*y to both sides...

(C*B + D*D/4)*y = (lambda^2*I - lambda*D + D*D/4)*y

(C*B + D*D/4)*y = (lambda*I - D/2)^2 * y

But at the moment I don't see this yielding a solution yet, as it still is not either a generalized eigenvalue problem, nor a standard eigenvalue problem.