Hi bil,
While it doesn't seem to help in this case wrt to the max_error figure-of-merit, you may want to consider revising the code to ensure that f_inverse and cov_matrix_P exactly symmetric (as it appears they should be based on the code).
What About My Matrix Makes eig Inaccurate?
6 views (last 30 days)
Show older comments
Hi all,
I want to understand what about my matrix is causing eig to return inaccurate eigensystem calculations. The code to reproduce the matrix (called 'cov_matrix_P') is attached below. I am wondering if it has to with the fact that there are several orders of magnitude between the largest and smallest elements in cov_matrix_P. Most importantly, I would like to know what about my matrix is causing these inaccuracies, and if possible, how I could improve the result.
N = 30;
f_matrix = zeros(N, N);
even_sum_matrix = zeros(N, N);
odd_sum_matrix = zeros(N, N);
for t = 1:N
for s = 1:N
even_sum = 0;
odd_sum = 0;
% Calculate even_sum
for r = 0:2:min(t, s)
if N-r < s-r || N-s < t-r
even_sum = 0;
continue
end
even_sum = even_sum + nchoosek(N, r) * nchoosek(N-r, s-r) * nchoosek(N-s, t-r) / (nchoosek(N, s) * nchoosek(N, t));
end
% Calculate odd_sum
for r = 1:2:min(t, s)
if N-r < s-r || N-s < t-r
odd_sum = 0;
continue
end
odd_sum = odd_sum + nchoosek(N, r) * nchoosek(N-r, s-r) * nchoosek(N-s, t-r) / (nchoosek(N, s) * nchoosek(N, t));
end
f_matrix(t, s) = even_sum - odd_sum;
even_sum_matrix(t, s) = even_sum;
odd_sum_matrix(t, s) = odd_sum;
end
end
[f_eigvector,f_eigvalues] = eig(f_matrix);
diag_elements = diag(f_eigvalues);
inverted_elements = 1./diag_elements;
f_inverse = f_eigvector * diag(inverted_elements) * f_eigvector';
variances = [0.9236396396396394, 0.9925285285285287, 0.9966406406406404, 0.9997037037037036, 1.0001001001001, 1.0008568568568565, 1.0008568568568565, 0.9999759759759757, 1.0006766766766766, 0.9999759759759757, 1.0008568568568565, 0.9998438438438437, 1.0008568568568565, 0.992892892892893, 0.9995555555555556, 1.000676676676677, 1.001001001001001, 1.0001001001001, 0.9982942942942948, 1.0005165165165162, 0.9997037037037036, 0.9982942942942948, 0.9992352352352354, 1.0006766766766766, 0.9995555555555556, 1.000424424424425, 0.9978618618618617, 0.9984984984984983, 0.9980820820820822, 1.0001001001001];
cov_matrix_G = diag(variances);
cov_matrix_P = f_inverse * cov_matrix_G * f_inverse';
[V,D] = eig(cov_matrix_P);
eig_check = cov_matrix_P * V - V * D;
max_error = abs(eig_check);
max_error = max(max_error(:));
Accepted Answer
Steven Lord
on 13 Jul 2024
N = 30;
f_matrix = zeros(N, N);
even_sum_matrix = zeros(N, N);
odd_sum_matrix = zeros(N, N);
for t = 1:N
for s = 1:N
even_sum = 0;
odd_sum = 0;
% Calculate even_sum
for r = 0:2:min(t, s)
if N-r < s-r || N-s < t-r
even_sum = 0;
continue
end
even_sum = even_sum + nchoosek(N, r) * nchoosek(N-r, s-r) * nchoosek(N-s, t-r) / (nchoosek(N, s) * nchoosek(N, t));
end
% Calculate odd_sum
for r = 1:2:min(t, s)
if N-r < s-r || N-s < t-r
odd_sum = 0;
continue
end
odd_sum = odd_sum + nchoosek(N, r) * nchoosek(N-r, s-r) * nchoosek(N-s, t-r) / (nchoosek(N, s) * nchoosek(N, t));
end
f_matrix(t, s) = even_sum - odd_sum;
even_sum_matrix(t, s) = even_sum;
odd_sum_matrix(t, s) = odd_sum;
end
end
[f_eigvector,f_eigvalues] = eig(f_matrix);
diag_elements = diag(f_eigvalues);
inverted_elements = 1./diag_elements;
f_inverse = f_eigvector * diag(inverted_elements) * f_eigvector';
variances = [0.9236396396396394, 0.9925285285285287, 0.9966406406406404, 0.9997037037037036, 1.0001001001001, 1.0008568568568565, 1.0008568568568565, 0.9999759759759757, 1.0006766766766766, 0.9999759759759757, 1.0008568568568565, 0.9998438438438437, 1.0008568568568565, 0.992892892892893, 0.9995555555555556, 1.000676676676677, 1.001001001001001, 1.0001001001001, 0.9982942942942948, 1.0005165165165162, 0.9997037037037036, 0.9982942942942948, 0.9992352352352354, 1.0006766766766766, 0.9995555555555556, 1.000424424424425, 0.9978618618618617, 0.9984984984984983, 0.9980820820820822, 1.0001001001001];
cov_matrix_G = diag(variances);
Let's look at the elements of the matrix whose eigenvalues you're computing.
format longg
cov_matrix_P = f_inverse * cov_matrix_G * f_inverse'
[minimumValue, maximumValue] = bounds(cov_matrix_P(:))
That's a wide range of elements. How well or ill conditioned is your matrix?
cond(cov_matrix_P)
That's pretty bad, considering the rule of thumb given in the first section of the Wikipedia page on condition numbers.
[V,D] = eig(cov_matrix_P);
eig_check = cov_matrix_P * V - V * D;
max_error = abs(eig_check);
max_error = max(max_error(:))
You asked "how I could improve the result." A flippant answer would be to not compute eigenvalues of a matrix with elements that cover such a wide range of orders of magnitude. But without knowing what exactly you're trying to do, what problem you're trying to solve, where did the elements of the variances matrix come from, why are you taking the sums of products of binomial coefficients, etc. it's likely going to be difficult to offer more targeted guidance for potential ways to solve the problem "better".
13 Comments
Paul
on 22 Jul 2024
Edited: Paul
on 23 Jul 2024
Hi Christine,
Thanks for responding. A few followups if you don't mind.
The doc page condeig says that "Large condition numbers imply that A is near a matrix with multiple eigenvalues."
a) should that really say "repeated eigenvalues." ?
b) What if the matrix actually has repeated eigenvalues? Is condeig meaningful in that case? Does it matter whether or not the repeated eigenvalues have geometric multiplicity > 1 (edit: should say: geometric multiplicity < algebraic multiplicity)? For example
A = eye(2), eig(A), condeig(A)
A(1,2) = 0.01, eig(A),condeig(A)
c) if A is "near a matrix with multiple eigenvalues" does that necessarily mean that the computed eigenvalues of A suffer from inaccuracy? What about the computed eigenvectors?
d) if yes to c), does condeig provide a quantitative measure of that inaccuracy? I see that, in the above case, the condition is 4.5e13 but the eigenvalues are exactly correct, though maybe that's not a good example becasue eig uses a special solver for a triangular matrix?
e) What does condeig == 1 for a symmetric matrix mean? It seems that a symmetric matrix could be arbitrarily close to being a matrix with multiple (repeated?) eigenvalues.
A = [1 eps;eps 1];
format long e
eig(A)
condeig(A)
Here's the case from the link I posted above
G = gallery('grcar',100);
cond(G)
figure
semilogy(condeig(G))
As I understand it, the cond means that solving linear systems of the form G*x = b should be o.k.
But what do the large values of condeig (for all eigenvalues) indicate? Don't trust eig(G) (or other functions that might use eig, like mpower, ^ in some cases)?
Christine Tobler
on 23 Jul 2024
Edited: Christine Tobler
on 23 Jul 2024
Hi Paul,
Quick aside to start: Most of my reply is based on sections 7.2 of Matrix Computations by Golub and Van Loan (third edition).
a) yes, repeated eigenvalues is meant (although I'm not sure of the exact difference).
b) If the matrix has repeated eigenvalues, the result of condeig for these repeated eigenvalues is not meaningful (to the best of my knowledge). Citing Section 7.2.3 "Sensitivity of Repeated Eigenvalues": "In general, if lambda is a defective eigenvalue of A, then O(eps) perturbations in A can result in O(eps^(1/p)) perturbations in lambda if lambda is associated with a p-dimensional Jordan block."
This also points to geometric multiplicity < algebraic multiplicity being a problem, since that would result in a Jordan block.
For a simple eigenvalue, a perturbation A+dA where ||dA||_2 = eps will result in a perturbation of this eigenvalue of order c*eps, where c is the element of condeig(A) associated with the eigenvalue.
The statement about being close to a repeated eigenvalue is due to another bit of perturbation theory: The norm of the smallest perturbation A+dA for which a given simple eigenvalue of A becomes repeated is proportional to 1/sqrt(c^2-1), where c>1 is the output of condeig for that eigenvalue.
c) This only applies to those specific eigenvalues that have a high condition number, not to all eigenvalues of the matrix. Yes, the computed eigenvectors will also suffer. As mentioned above, if there are repeated eigenvalues this analysis doesn't apply anymore, it's only relevant for simple eigenvalues.
d) For a simple eigenvalue, condeig provides a quantitative measure of how much the computed eigenvalue might be off from the original one. An example of a triangular matrix will have eig just pick the diagonals of that matrix as eigenvalues, so no perturbations are applied. This is similar to a linear system: a diagonal matrix might have a high condition number, but in many cases this won't have an immediate effect, since no factorization is necessary and so less round-off error is introduced there.
e) The symmetric case is different because these matrices have orthogonal eigenvalue bases. This means a small perturbation in A will never cause a larger perturbation of any simple eigenvalue of A. The theory I'm reading doesn't directly address multiple eigenvalues in a symmetric matrix, but since these will never form a Jordan block, I would say that the same applies there.
For the matrix from gallery, yes since it has high condition numbers for most eigenvalues, you should not trust the eigenvalues computed by EIG too much. Here's an illustration:
G = gallery('grcar',100);
[U, D] = eig(G);
F = randn(size(G)); F = F / norm(F);
[U2, D2] = eig(G + eps*F);
plot(diag(D), 'o'); hold on; plot(diag(D2), 'x')
norm(G*U - U*D)
norm(G*U2 - U2*D2)
Both of the outputs here, the original and the perturbed matrix, have quite low residuals, so they satisfy G*x = lambda*x up to round-off error. However, a small change in G can lead to a much larger change in the eigenvalues - this is what the condeig result alerts about.
Whether this matters to an algorithm that uses an eigenvalue decomposition will depend on the algorithm. For mpower, there is an alternative algorithm using the Schur decomposition that gets used for non-symmetric matrices in more ill-conditioned cases.
More Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)