What About My Matrix Makes eig Inaccurate?

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(:));

[Paul]

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).

Answer 1

[Steven Lord]

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".

[John D'Errico]

I'd add the condition number of f_inverse is roughly 1e8. Then forming

cov_matrix_P = f_inverse * cov_matrix_G * f_inverse'

would be expected to approximately square the condition number.

In the end though, cov_matrix_P is numerically singular. All computations done with it in double precision will be problematic.

[bil]

I see, the condition number was something I was unaware of! Initially, I saw the determinant of cov_matrix_P was extremely large so I assumed it wasn't a singularity issue. Unfortunately, the f_matrix is from an analytically derived formula and cov_matrix_G is from running a classical circuit and calculating a specific quantity G. I basically wish to propagate the uncertainty forward from G to P so there is not much I can think of adjusting in those regards. I was just wondering if there was some other decomposition method that would retrieve more accurate eigensystems even with ill-conditioned matrices, but perhaps this is just wishful thinking considering how poor my input is.

[John D'Errico]

@bil

WRONG. A large computed determinant does not imply something is non-singular!

For example:

A = magic(6)

A = 6x6

35 1 6 26 19 24 3 32 7 21 23 25 31 9 2 22 27 20 8 28 33 17 10 15 30 5 34 12 14 16 4 36 29 13 18 11

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

det(A)

ans = -3.4498e-09

The problem is, I recall ll magic matrices with an even order will be singular. We can see the matrix is singular, because a linear combination of the rows or columns of A yields the zero vector.

V = double(null(sym(A)))

V = 6x1

2 2 -1 -2 -2 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

A*V

ans = 6x1

0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

But as you see, det(A) did yield a tiny, but non-zero number. Now it is trivial to modify A so that the determinant is as alarge as we wish.

Amod = A*10000

Amod = 6x6

350000 10000 60000 260000 190000 240000 30000 320000 70000 210000 230000 250000 310000 90000 20000 220000 270000 200000 80000 280000 330000 170000 100000 150000 300000 50000 340000 120000 140000 160000 40000 360000 290000 130000 180000 110000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

det(Amod)

ans = -1.3188e+16

As such, Amod must be non-singular, according to your logic, because the determinant was large? Right? But if A was singular because the determinant was tiny, then how can Amod be non-singular?

Amod*V

ans = 6x1

0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Surely you see the lie being unmasked here. Clearly, multiplying A by some constant does nothing to change the singularity of A. I hope this convinces you of the duplicitous nature of the determinant. The determinant is a subtle liar, taught to you by teachers that also very likely fail to recognize the lie themselves.

The source of the problem stems from the fact that you CANNOT use the tools taught to you by your teachers to compute the determinant of a matrix when using floating point arithmetic, especially when the matrices become larger. For example, to compute the determinant of a 100x100 matrix by use of minors would take the lifetime of the universe to perform. Instead, these codes use linear algebra, where that same determinant can be computed in a fraction of a second, even for very large matrices. At the same time, that same linear algebra loses a few bits in the least significant bits of the result. The trick is to use an LU factorization of A, then take the product of the diagonal elements of U.

[L,U,P] = lu(A);
format long g
diag(U)

ans = 6x1

35 35.8857142857143 25.5883757961783 7.64155698683842 5.27176781002638 2.66453525910038e-15

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Do you see what happens now? That last diagonal element is not EXACTLY zero, but by only a tiny amount. And this what your teachers did not teach you, that a determinant is a duplicitous thing. Most likely they are just teaching you what they were themselves taught. And even if your teachers undertand the mathematics of what I have said above, they surely did not appreciate the consequences of doing arithmetic in double precision, using floating point arithmetic.

Back to your original matrix.,,

The matrix f_inverse is not numerically singular. However, it is close enough to the edge, that when you effectively multiply it twice there, the result now IS numerically singular.

Why exactly you are doing what you are doing is uncertain. We would need to untangle what you are doing from the code you wrote. The result is a numerical singularity though.

[Paul]

@Steven Lord or @John D'Errico

cond returns the condition number of a matrix wrt inversion, which is important when solving a system of linear equations.

But the question here is about the accuracy of the eigendecoposition of the matrix. Does cond tell us anything about that?

This link indicates that senstivity wrt inversion is not the same as the sensitivity of the eigenvalues (though it doesn't say anything about the sensitivity of the eigenvectors)

[Steven Lord]

As John said, the determinant is not a good way to decide if a matrix is singular or not. The identity matrix is not singular, so shouldn't a non-zero multiple of the identity matrix also not be singular?

A = 0.1*eye(500);
det(A)

ans = 0

cond(A)

ans = 1

A matrix with two rows that are the same (or almost the same) is singular, right?

A = [1 1; 1 1+eps];
det(A)

ans = 2.2204e-16

cond(A)

ans = 1.6199e+16

But if we multiply A by a large number:

det(1e12*A)

ans = 244140625

cond(1e12*A)

ans = 1.3468e+16

That doesn't make it non-singular, does it?

[John D'Errico]

@Paul

I thought I would get that question. :)

Cond still tells you a lot. It is a good measure of the ratio of the largest and smallest eigenvalues too, even though it is the ratio of singular values.

Consider the very simplest case. That is, note that if your matrix is SPD, then the eigenvalue decomposition and the svd will be identical (with the usual caveat that the vectors will be subject to random sign changes.)

For example:

A = rand(5);
A = A'*A

A = 5x5

1.8025 1.8488 1.6260 2.0895 1.5647 1.8488 2.0629 1.8720 2.2783 1.6881 1.6260 1.8720 2.1674 2.2831 1.4750 2.0895 2.2783 2.2831 2.7133 1.8418 1.5647 1.6881 1.4750 1.8418 1.6651

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[V,D] = eig(A)

V = 5x5

0.5449 0.3833 0.5151 -0.3419 0.4171 0.0057 -0.8615 0.1667 -0.1493 0.4558 0.4138 0.0446 -0.3591 0.7085 0.4425 -0.7243 0.2881 0.2673 0.2122 0.5253 -0.0851 0.1613 -0.7117 -0.5601 0.3827

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 5x5

0.0352 0 0 0 0 0 0.0653 0 0 0 0 0 0.1897 0 0 0 0 0 0.5063 0 0 0 0 0 9.6146

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[u,s,v] = svd(A)

u = 5x5

-0.4171 0.3419 -0.5151 0.3833 -0.5449 -0.4558 0.1493 -0.1667 -0.8615 -0.0057 -0.4425 -0.7085 0.3591 0.0446 -0.4138 -0.5253 -0.2122 -0.2673 0.2881 0.7243 -0.3827 0.5601 0.7117 0.1613 0.0851

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

s = 5x5

9.6146 0 0 0 0 0 0.5063 0 0 0 0 0 0.1897 0 0 0 0 0 0.0653 0 0 0 0 0 0.0352

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

v = 5x5

-0.4171 0.3419 -0.5151 0.3833 -0.5449 -0.4558 0.1493 -0.1667 -0.8615 -0.0057 -0.4425 -0.7085 0.3591 0.0446 -0.4138 -0.5253 -0.2122 -0.2673 0.2881 0.7243 -0.3827 0.5601 0.7117 0.1613 0.0851

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The order of the eigenvalues and singular values has changed. And there are sign flips. But nothing significant is different, and we would expect that to be exactly the case. Of course, this does not mean we should be using svd instead of eig for SPD matrices. eig should be faster. The point is, cond is indeed a good measure of what to expect from eig too.

[Paul]

@Steven Lord

Are you responding to my comment? If so, I'm not sure why you're steering the discussion back to det and cond. I agree that det is not a good way to measure how "close" a matrix is to being singular and that cond is the tool for that job and I don't think I said anything to indicate otherwise.

The question is whether or not cond is the correct measure for getting a feel for whether or not an eigendempoisition is accurate. In other words, is a matrix being poorly conditioned for inversion sufficient to say that that it's poorly conditioned for eigendecomposition.

[Steven Lord]

I think I was writing my comment in support of John's point about the determinant about the same time you were writing your comment but you pressed Submit before I did.

So what exactly leads you to think that the eigenvalues and eigenvectors are inaccurate? Let's go back to your original matrix.

cov_matrix_P = createMatrix();

[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

So the maximum absolute value of an element of the residual is 2. What are the elements of cov_matrix_P*V; what orders of magnitude are they?

format longg
part1 = cov_matrix_P*V

part1 = 30x30

1.0e+00 * -1057277131.90617 63.5487302609211 1.14894262965946 81.3994278654543 0.878619049326739 -70.1941726999822 3.50121914947914 -51.4510066113024 -0.0897804657269083 -28.9812078977131 -15.39058221294 -19.306775290403 0.39550331040483 7.43345902993629 7.26674366332737 5.39183720663283 0.353968453537508 -1.51715028723187 -1.84223835675085 -1.22593032106307 -0.11503243196832 -0.187963700347513 0.454243545513761 0.229951128080962 0.00646238333804196 -0.04646390956877 -0.0290189206185331 -0.0418427666293209 0.0308031264900734 -0.00136956704561005 -15588843482.1225 11.8461823209328 -491.434372698041 -3.57225091904798 531.822461521744 -23.5577800332915 -411.818268483477 0.173313188293465 261.746379941609 -69.899024128943 126.780411977225 -0.0359162423668764 -70.6957533841596 -21.6200759762468 23.8582823055021 0.0274623844359476 -12.9643307181056 3.96596453547921 -3.98382309929421 -0.0168942377325743 1.72481807721642 0.622586864151057 0.338848679155394 0.0067858928645803 0.00808733292054089 0.128883869579374 0.0558031724634705 -0.0195643526457083 0.00146586817995318 -0.00222930155306051 -143071786587.935 -2770.02203257087 -13.7704478591234 -2811.19178965703 -14.8092041958388 1980.03684029865 -102.700850671024 1150.34324923576 2.4281126726588 504.728022806788 268.488539857172 256.408368956956 -5.16647868910719 -72.3446664265808 -70.6320429731391 -37.3800481689404 -2.3707288191646 6.84599843643786 8.2658131370647 3.44409284683839 0.294203423731469 0.243395365891976 -0.550427128268639 -0.0988952203133757 -0.0200608558467471 0.00573501433318845 -0.012268440247298 -0.0203212537095339 0.0196333523433816 -0.00105690425655377 -964884744193.907 -136.270093744543 11446.0141163835 51.6005017372214 -10216.5233365471 355.27601917026 6403.51846727099 -1.93681273105937 -3237.88048658697 673.777168589287 -1223.40342825005 0.319806464322821 514.405859949578 114.729217102376 -126.4712659524 -0.102845633979024 46.4345174662269 -9.01674221568625 9.05187189289521 0.0211622269623427 -1.88048984097324 -0.230155848525556 -0.113708259824242 0.00159210484246806 0.0382941880494067 0.0481267805096201 0.0372285209912731 -0.0136348632714356 0.00118898507498325 -0.00185018639210773 -5021877243514.14 38168.0437094404 115.690360529589 30937.4315746899 112.398914661114 -17255.6723290495 902.627300087253 -7726.10432524119 -17.2642273784695 -2523.10902310471 -1343.36940609746 -910.950943756504 18.1599274807287 167.308497112927 163.162232512038 49.5592111813619 3.11370761255842 -3.11084834975922 -3.68307544138414 0.614114486785927 0.0701460294306108 0.213805592710804 -0.503815563210945 -0.165017130136677 -0.0429351776262561 0.0241246444196488 -0.00308981733909571 -0.0114761994558224 0.0143465760802063 -0.00108634632949318 -20926134494723.8 891.198312071347 -107404.161614559 -326.609219497767 76489.6638880546 -2020.34852758889 -37337.5531677827 6.73277637949605 14213.6651217873 -2122.34662694311 3853.46429305943 -0.668757327879892 -1068.60427145685 -135.794256762994 150.052014064632 0.0393413413978115 -18.1188656720697 -2.20645496617815 2.16021954756216 0.0197492277485079 -1.84541672195809 -0.471739680223795 -0.258222858914922 0.00140438998009945 0.0610527281251111 0.0023345907854466 0.0233319016694293 -0.0105367599068777 0.00297713337661838 -0.00136374390102285 -71740929954878.2 -258220.48290398 -708.493353385245 -163098.229241672 -701.381574672841 68703.0050296606 -3605.03149920282 22000.4108528488 77.8671778268922 4697.70123878098 2502.70051156456 925.943281538208 -17.5742676236017 -41.3964367088218 -40.158773524407 23.4998374215517 1.1357481910857 -5.55661581913696 -6.72820902052624 -1.39149333466875 -0.14745732444429 0.113408028139989 -0.291503689460223 -0.151535301285556 -0.0563192446928492 0.0272338843459716 0.000757919617767898 -0.0108836494329389 0.0173008801252986 -0.00122117461658267 -206253109592342 -2368.11257737395 549125.950006373 1142.13834033091 -298853.565612981 5881.22917059482 105690.568343266 -9.39483983948145 -26793.2644602196 2228.92890227538 -4089.47594562253 -0.00854915065859263 316.058062267808 -67.5705331231684 70.3181772757639 0.0956056677496551 -37.854555291017 3.75476400095284 -4.2373562711851 0.0030491848568789 -0.91324144703865 -0.474460434297481 -0.260828433671896 0.0219798073440372 0.088243121077556 -0.0285247739330405 0.0200197946382071 -0.0227007570511577 0.0219878052724104 -0.00229165784287806 -504179690976176 1030601.29092472 2888.82713485216 477705.140764259 524.187137799591 -136800.776315883 7188.27485055884 -25558.3670288187 66.1152650917443 -1796.10011853948 -956.912144947426 454.986875950735 9.29922743560152 -131.707860820147 -128.635571688558 -18.7546062295537 -3.01994691251745 -3.54363179993859 -4.48430911561573 -2.10400994716064 -0.260037112613648 0.00374432921534008 -0.0457006372161338 -0.0716668774812156 -0.0401932405580333 0.01761788582705 0.0165390198661959 -0.0446139510271653 0.061535751631782 -0.00533432141459828 -1.05877906028574e+15 11512.7243708521 -1737543.40947942 -2336.01379284467 663000.621248506 -6973.23662362695 -145768.96578808 3.64067298525748 15712.3583540608 606.189629018525 -1323.25097086295 0.543715719845421 862.384002631667 60.1315417320251 -82.4869453011908 0.075117510547545 -21.4486167400237 5.60170515052121 -6.81974909558464 0.0119750565935565 0.160050486968813 -0.340877400595533 -0.0996224245338704 0.123809169996761 0.159702008688796 -0.00416972693865268 0.11503171860679 -0.121963018574166 0.121639280856071 -0.0207153504038682

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[minimumValue, maximumValue] = bounds(part1(:))

minimumValue = 
     -5.46629356894126e+15

maximumValue = 
          5221402.33045838

So you're taking the difference of numbers varying between 5 million and -5e15 and the maximum difference is 2. That seems pretty darn good to me!

part2 = V*D;
[maxValue, maxLocation] = max(abs(part1-part2), [], 'all')

maxValue = 
     2

maxLocation = 
    14

part1(maxLocation)

ans = 
     -5.12465103351602e+15

part2(maxLocation)

ans = 
     -5.12465103351603e+15

That maximum difference happens when we're subtracting two extremely large numbers. In fact, let's take a look at a picture of the differences and a picture of one of the parts.

surf(part1-part2)
title('Difference between A*V and V*D')

figure
surf(part1)
title('A*V')

Given the magnitudes of some of the numbers in part1 and part2, that looks pretty flat relatively speaking.

function cov_matrix_P = createMatrix()
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';
end

[Paul]

@Steven Lord

I think that this comment (starting from "So what exactly leads you ...") should be it's own answer, as it gets to directly to the premise of the question (or, more precisely, rejects the premise). If you do move it, feel free to delete this comment.

If you decide to not move it, please clarify who the "you" is in "leads you." As written, it sounds like it refers to me, but I've made no claims in this thread about the accuracy or inaccuracy of the eigendecoposition of cov_matrix_P.

[Paul]

@John D'Errico

Referring to this comment ...

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]

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]

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]

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.