# What About My Matrix Makes eig Inaccurate?

73 ビュー (過去 30 日間)
bil 2024 年 7 月 13 日 23:44

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 件のコメント-1 件の古いコメントを表示-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'
cov_matrix_P = 30x30
1.0e+00 * 93.8043019442907 1360.01645173048 12470.8340160598 84179.1781741391 438157.9425107 1825653.67855217 6258809.24139487 17994088.9920376 43986107.3929892 92370800.6143717 167946514.421767 265915351.256035 368190699.188861 447088665.786112 476894484.530979 447088612.098099 368190659.069842 265915456.673998 167946587.698871 92370630.9603515 43986022.3205832 17994280.1056806 6258877.0729833 1825505.83054011 438122.244866562 84254.3962813616 12481.980735042 1337.32202058047 92.2419557900037 3.07795146321052 1360.01645173048 20064.2245109128 184038.364026197 1241085.34106718 6459821.56497993 26918287.1188646 92282939.8092452 265310352.522754 648544483.784898 1361946137.06402 2476259958.91919 3920743404.68677 5428724746.10349 6592023308.15625 7031490151.59952 6592022161.82391 5428724717.26371 3920745421.31514 2476259994.04333 1361943102.01769 648544461.678827 265313616.464764 92282947.5614691 26915860.5886842 6459820.39168697 1242274.84927721 184038.27898586 19717.8741007064 1360.05610025603 45.3880988580989 12470.8340160598 184038.364026197 1689210.48911282 11391213.0443175 59286520.5905216 247049247.555479 846958923.547858 2434977936.59785 5952229783.51951 12499708021.6831 22726700580.0167 35983928188.9904 49823921738.1568 60500480271.5677 64533836363.4026 60500472842.3599 49823923058.9349 35983942553.0357 22726698329.7226 12499685035.8867 5952232089.14133 2435003793.75894 846957395.026471 247029247.102205 59287145.5202459 11401389.8406955 1689072.41883416 180967.349218876 12482.3772188182 416.56264734115 84179.1781741391 1241085.34106718 11391213.0443175 76823788.4429797 399836148.723596 1666111976.76203 5711931187.17745 16421644636.0418 40142213098.0415 84298773850.503 153270219333.44 242677787964.115 336015541456.766 408018875825.114 435220071950.064 408018834969.383 336015539573.813 242677874450.585 153270221626.985 84298628813.2141 40142211654.1901 16421812167.4412 5711931693.57481 1665980127.63903 399836072.325106 76891657.452467 11391207.2518728 1220455.08169295 84181.8317427187 2809.28984241382 438157.9425107 6459821.56497993 59286520.5905216 399836148.723596 2081006950.75264 8671527976.99572 29728529985.8135 85468704662.568 208925772087.638 438744778211.662 797716204958.43 1263050349647.58 1748839783605.31 2123591186642.71 2265163557101.76 2123590926094.23 1748839753535.99 1263050853697.4 797716246160.195 438743971442.962 208925740419.019 85469612258.2422 29728545199.2036 8670825951.45108 2081002617.25021 400193356.793071 59287127.0906648 6352028.78227955 438136.025267943 14621.3443807851 1825653.67855217 26918287.1188646 247049247.555479 1666111976.76203 8671527976.99572 36134223846.1113 123878671379.473 356147574694.086 870592426187.198 1828247117595.76 3324079111450.69 5263123688856.91 7287405522702.58 8848992706790.22 9438923989228.59 8848991523490.58 7287405482198.13 5263125883186.39 3324079160787.13 1828243681334.72 870592395125.552 356151398679.768 123878682278.484 36131283346.8626 8671526329.41648 1667603500.99252 247049124.100474 26468866.7886876 1825710.6906639 60927.1829807804 6258809.24139487 92282939.8092452 846958923.547858 5711931187.17745 29728529985.8135 123878671379.473 424692488388.889 1220978542218.92 2984646201023.91 6267767385486.64 11395918842924.4 18043532781741 24983364488501.8 30336943575195.8 32359401473277 30336939853269.9 24983364468758.3 18043539982332.7 11395918855324.6 6267755860245.87 2984646207155.91 1220991507894.76 424692476427.326 123868642420.706 29728536929.2182 5717034167.02116 846956896.134621 90743042.182896 6259072.19488388 208876.290267877 17994088.9920376 265310352.522754 2434977936.59785 16421644636.0418 85468704662.5679 356147574694.086 1220978542218.92 3510278317533.23 8580772293471.75 18019650627775.5 32762938300374.7 51874637996525.8 71826454632443 87217839631167.2 93032347880417.5 87217829476846.3 71826454232628.9 51874658237141.6 32762938787342 18019617797001.3 8580771986914.56 3510315451728.17 1220978649764.42 356118783800.639 85468688412.9102 16436308688.3832 2434976716.44988 260883641.037636 17994651.6008862 600513.018508992 43986107.3929892 648544483.784898 5952229783.51951 40142213098.0415 208925772087.638 870592426187.198 2984646201023.91 8580772293471.75 20975447103416 44048509345549.6 80088043383802 126806033704593 175577667224040 213201457612854 227414851512584 213201431454561 175577665869713 126806084309552 80088045036257.1 44048428348581.4 20975446056873.6 8580863413224.65 2984646576626.31 870521944998.366 208925708162.328 40178075677.328 5952229319.87583 637722427.826094 43987399.0385379 1467936.33253365 92370800.6143717 1361946137.06402 12499708021.6831 84298773850.503 438744778211.662 1828247117595.76 6267767385486.64 18019650627775.5 44048509345549.6 92502019617861.8 168185163817793 266293099638757 368713694978084 447723784838796 477571960281201 447723728063657 368713692926104 266293207311508 168185166317132 92501848738990.8 44048507772100.7 18019842283692.2 6267767937504.05 1828099032300.64 438744694799.833 84374095395.6352 12499701758.7745 1339219245.74966 92373688.6075544 3082671.84963881
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
[minimumValue, maximumValue] = bounds(cov_matrix_P(:))
minimumValue =
0.205203174074088
maximumValue =
2.46562161696071e+15
That's a wide range of elements. How well or ill conditioned is your matrix?
cond(cov_matrix_P)
ans =
3.99843362237122e+17
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(:))
max_error =
2
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 件のコメント11 件の古いコメントを表示11 件の古いコメントを非表示
Paul 2024 年 7 月 22 日 21:57

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 = 2x2
1 0 0 1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
ans = 2x1
1 1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
ans = 2x1
1 1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
A(1,2) = 0.01, eig(A),condeig(A)
A = 2x2
1.0000 0.0100 0 1.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
ans = 2x1
1 1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
ans = 2x1
1.0e+13 * 4.5036 4.5036
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
V = 2x2
1.0000e+00 -1.0000e+00 0 2.2204e-14
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
D = 2x2
1 0 0 1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
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)
ans = 2x1
9.999999999999998e-01 1.000000000000000e+00
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
condeig(A)
ans = 2x1
1 1
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Here's the case from the link I posted above
G = gallery('grcar',100);
cond(G)
ans =
3.591110462980591e+00
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 2024 年 7 月 23 日 13:27

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)
ans = 8.0384e-14
norm(G*U2 - U2*D2)
ans = 7.9470e-14
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.

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeLinear Algebra についてさらに検索

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by