Improving Precision of Eigenvectors with Large Eigenvalues

Hi all,

This is a bit of a generic question, but I was hoping someone could provide some insight on how I could improve the precision of eigenvectors using "projection techniques", like in this post, where a matrix will have 1 or 2 very large eigenvalues, and the remaining eigenvalues are much smaller (by several orders of magnitude). They have code written in R, which I am not too familiar with, but is the idea generically that I should subtract out the overlap of smaller eigenvectors with those of larger eigenvectors to improve their precision? When I say precision, what I mean is that applying eig to a matrix M will generate a set of eigenvectors, but once I check M*v - λ*v, the resultant array will deviate from 0, i.e. applying M to the "eigenvector" has resulted in a linear combination of multiple other eigenvectors.

Any guidance is appreciated.

Answer 1

[Christine Tobler]

The linked post is about a symmetric matrix, is this also your case? In that case (if issymmetric returns true for your matrix), the eigenvectors returned by EIG will be orthogonal up to numerical round-off (an exact 0 is not possible in numerical computation, practically speaking).

The proposed solution doesn't improve the eigenvectors, but instead applies a projection to the computed residual vectors, to project out any components along the eigenvectors of the larger eigenvalues. Whether this is useful will depend on your application - that is, do you need to do computations with that residual?

[bil]

Indeed, I am working with a symmetric matrix. Could you elaborate on what you mean by computed residual vectors? I thought the "Analysis" and "Solution" section of the post matched what I wanted pretty well: The eigenvectors that are numerically returned are slightly off/perturbed from the "real" eigenvectors in such a way that applying the matrix on the eigenvector will no longer return only (eigenvalue*eigenvector), but a linear combination of multiple other eigenvector as well, so the issue is much worse when we have a set of eigenvalues where one is significantly larger than the others. Perhaps I'm misunderstanding things though.

[Christine Tobler]

You're right - they don't improve the residual, although it turns out they also don't improve the eigenvectors. I had the R-code autotranslated to MATLAB (don't know much R myself), and here's what it seems to do:

% Matrix copied over from 
% https://stats.stackexchange.com/questions/539436/get-accurate-eigenvectors-when-eigenvalues-are-minuscule
M = [
    0.002307414, -1.318435973, 8.63e-11, -3.33e-11, 0.486238053, 0.41182779, -0.169525454, 0.624275558,
-1.318435973, 753.3425186, -4.93e-08, 1.90e-08, -277.8320732, -235.3147146, 96.86532753, -356.7054684,
8.63e-11, -4.93e-08, 3.23e-18, -1.24e-18, 1.82e-08, 1.54e-08, -6.34e-09, 2.34e-08,
-3.33e-11, 1.90e-08, -1.24e-18, 4.79e-19, -7.01e-09, -5.94e-09, 2.44e-09, -9.00e-09,
0.486238053, -277.8320732, 1.82e-08, -7.01e-09, 102.4642297, 86.78386443, -35.72384951, 131.5526701,
0.41182779, -235.3147146, 1.54e-08, -5.94e-09, 86.78386443, 73.50310589, -30.25693671, 111.4208258,
-0.169525454, 96.86532753, -6.34e-09, 2.44e-09, -35.72384951, -30.25693671, 12.45501408, -45.8654479,
0.624275558, -356.7054684, 2.34e-08, -9.00e-09, 131.5526701, 111.4208258, -45.8654479, 168.8989909
];

% Extract the size, verify M is exactly symmetric, and call eig
n = size(M, 1);
issymmetric(M)

ans = logical
   1

[U, D] = eig(M);

% We can see the residual is accurate to machine precision. Also, the
% eigenvalues are orthogonal to each other up to machine precision.
residual = norm(M*U - U*D)

residual = 4.7680e-13

orthogonality = norm(U'*U - eye(n))

orthogonality = 9.7243e-16

% Extract eigenvalues into a vector, logical map for the two "large"
% eigenvalues.
lambda = diag(D);
largeEigenvalue = lambda >= max(lambda) * 1e-12;

for ii=1:500
    % Make a random vector dx that should only be comprised of components
    % along the small eigenvectors:
    randVec = randn(8, 1);
    randVec(largeEigenvalue) = 0; % set to zero for large eigenvalues
    dx = U*randVec;
    dx = dx / norm(dx);
    normOrig(ii) = norm(M*dx);

    % Correct this random vector dx by projecting away components along the
    % large eigenvectors that have been introduced by round-off error in
    % U*randVec operation.
    dxCorrected = dx - U(:, largeEigenvalue) * (U(:, largeEigenvalue)' * dx);
    dxCorrected = dxCorrected / norm(dxCorrected);
    normCorrected(ii) = norm(M*dxCorrected);
end

semilogy(normOrig, 'x'); hold on; semilogy(normCorrected, 'o')
legend('orig', 'corrected')

So in this case, the accuracy of the original and the corrected version are about the same, and about match the accuracy of the corrected result from the blog post you linked.

Perhaps the eigenvalues and eigenvectors computed there are less accurate? Or I may have mixed something up - as I said, I don't know R very well and am relying on some automatic translation.

[bil]

Thanks a lot for the translated code. Indeed it seems that R might have lower precision than Matlab with calculating eigenvectors. To confirm my understanding: while the code in the post is dealing with an arbitrary vector living in the space spanned by the smaller eigenvectors, if I were to apply this technique directly one of the eigenvectors, I would subtract the projection onto all of the other eigenvectors right (not just the biggest eigenvectors)?

[Christine Tobler]

Very generally, you want to orthogonalize (meaning subtract the projection) against those eigenvectors in which you have higher trust. Orthogonalizing against eigenvectors who you don't trust to be accurate would be just spreading the noise around. Think about 2 vectors x and y. If you make xnew perpendicular to y and ynew perpendicular to x, the new vectors (xnew, ynew) won't be more perpendicular to each other than the old ones.

For example, in an Arnoldi algorithm, you would orthogonalize a candidate for a new eigenvector against all eigenvectors that have already been computed accurately.

Can you say more about your problem? Where are your eigenvectors coming from, and what problem are you aiming to solve by using projection?

[bil]

I actually initially encountered this problem while working in Python and realized that Matlab had significantly better accuracy than any library I was aware of. I think much of the inaccuracies is related to the fact that I am interested in finding the eigensystem of the inverse of a matrix, and that for larger size, the matrix slowly becomes singular, thus resulting in inaccuracies in subsequent calculations. This has already been very helpful for me though. Much appreciated.

[Torsten]

I think much of the inaccuracies is related to the fact that I am interested in finding the eigensystem of the inverse of a matrix,

If the matrix is invertible, shouldn't both A and inv(A) have the same eigensystem ?

If lambda ~= 0 is eigenvalue of A, then 1/lambda is eigenvalue for inv(A) to the same eigenvectors.

[bil]

It is, but for some reason, calling something like scipy.linalg.inv in Python versus taking the reciprocals of the eigenvalues for the pre-inverted matrix give different results, and the differences get quite significant as I push the matrix size.

[Torsten]

I don't understand what you are doing. Determining eigenvalues and eigenvectors for a given matrix will automatically give you eigenvalues and eigenvector of its inverse. No need to use inversion or pre-inversion or whatever.

[bil]

Sorry, I should have been more explicit in my previous responses/original question. But I am actually interested in finding the eigensystem of the inverse of a certain matrix M. So I initially did that by using a numpy library to invert M and then used the scipy function to find its eigensystem. But I suspect this might not be the most accurate and the most straightforward method is, as you said, to just take the reciprocal of the eigenvalues of M to get the eigenvalues of , and the eigenvectors are the same for both matrices.

[Torsten]

But I suspect this might not be the most accurate and the most straightforward method is, as you said, to just take the reciprocal of the eigenvalues of M to get the eigenvalues of M^(-1), and the eigenvectors are the same for both matrices.

Yes, at least I cannot think of any advantage to work with the inverse for your case.

[Christine Tobler]

Closing the loop, I agree with Torsten that computing the eigenvalues and eigenvectors of the original matrix and then inverting the eigenvalues will be more accurate than calling eig on the inverse.