What About My Matrix Makes eig Inaccurate?
61 ビュー (過去 30 日間)
古いコメントを表示
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(:));
1 件のコメント
Paul
2024 年 7 月 14 日 14:44
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).
採用された回答
Steven Lord
2024 年 7 月 13 日 23:53
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".
12 件のコメント
Christine Tobler
約9時間 前
Hi Paul,
I agree that the condition of the eigenvalues doesn't map to the cond function call on A. Instead, condeig(A) will give the condition of each eigenvalue, as described in the link you posted.
In the case of a symmetric matrix, the condition number of each eigenvalue will be 1.
In practice, that isn't the full story, though. Eigenvalues smaller than eps * (the largest eigenvalue) are likely to just be noise, unless the matrix has some special structure (diagonal, tridiagonal, block, ...), because the initial orthogonal transformations inside of EIG will introduce noise of the same level as these eigenvalues.
So cond can be useful as an indication of the largest eigenvalue divided by the smallest. A large condition number doesn't mean the whole result of EIG is inaccurate - just that all residuals are bounded w.r.t. norm(A), and that might just leave the smallest of the eigenvalues as noise.
Paul
約1時間 前
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? 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)?
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!