What About My Matrix Makes eig Inaccurate?
59 次查看(过去 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".
11 个评论
Paul
2024-7-14,20:13
Yes, if A is SPD (note that cov_matrix_P as computed isn't, though it's clearly intended to be) then its eigenalues are its singular values, in theory.
However, I still don't see how, in practice, it follows that cond tells us anything about how accurately the eigenvalues and eigenvectors can be computed, and the link I provided previously claims that cond, in fact, does not.
Christine Tobler
about 7 hours 前
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!