Eig function for symmetric matrices
显示 更早的评论
Two of the properties of symmetric matrices are that their eigenvalues are always real, and that they are always orthogonally diagonalizable. Given any two distinct eigenvalues, the corresponding eigenvectors are orthonormal. Yet if some eigenvalue L with multiplicity greater than 1 appear, it is necessary to use, for example, Gram Schmidt for ensuring orthogonality between the vectors representing the span of Ker(A - Id L). I was so glad to see that computing the MatLab function eig over symmetric matrices having multiple eigenvalues would output an orthogonal eigenvector matrix, meaning that MatLab doesn't only normalize the vectors, but also it make them orthogonal! I tried it on a few matrices, I just need to know if I was "lucky" with them, and if usually eig doesn't do that, or, elsewise, what procedure is implemented in eig for dealing with this situation. Thank you for your help.
回答（2 个）
Matt J
2016-10-4
0 个投票
I doubt you were just lucky. I've never seen EIG fail to give orthogonalized eigenvectors for symmetric matrices if it can recognize the matrix as symmetric. If there's a chance you'll have numerical error in the creation of the matrix (such that it is not perfectly symmetric), it would be safest to use SVD instead.
1 个评论
Walter Roberson
2016-10-4
Or force them to be symmetric by using
A = (A + A.')/2;
Steven Lord
2016-10-4
编辑：Steven Lord
2016-10-4
According to the documentation page for eig, specifically the section describing the output argument V:
[V,D] = eig(A) returns matrix V, whose columns are the right eigenvectors of A such
that A*V = V*D. The eigenvectors in V are normalized so that the 2-norm
of each is 1.
If A is real symmetric, then the right eigenvectors, V, are orthonormal.
[Edited to fix formatting of quoted text.]
13 个评论
Lorenzo
2018-9-20
Hi Steven,
I noticed that if A is complex and hermitian with repeated eigenvalues, then the right eigenvectors, V, fail to be orthonormal (V is not unitary). Do you know a way to fix this?
Bruno Luong
2018-9-20
编辑：Bruno Luong
2018-9-20
Lorenzo, make sure you input matrix is exactly Hermitian so that MATLAB engine can correctly detect it form. This can accomplish by "symmetrizing" it.
>> [Q,~]=qr(randn(3)+1i*rand(3));
>> A=Q'*diag([1 1 100])*Q
A =
15.3593 + 0.0000i -7.9491 - 8.0340i 23.9826 +22.6380i
-7.9491 + 8.0340i 9.8955 - 0.0000i -25.9424 + 0.8862i
23.9826 -22.6380i -25.9424 - 0.8862i 76.7452 + 0.0000i
>> A'-A
ans =
1.0e-14 *
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.1776i -0.3553 - 0.3553i
0.0000 + 0.0000i 0.3553 - 0.3553i 0.0000 + 0.0000i
>> [V D]=eig(A)
V =
0.9246 + 0.0000i 0.2770 + 0.2614i -0.4084 + 0.0100i
0.0868 - 0.0878i -0.2996 + 0.0102i 0.8103 + 0.0000i
-0.2620 + 0.2473i 0.8747 + 0.0000i 0.4039 - 0.1157i
D =
1.0e+02 *
0.0100 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 1.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0100 - 0.0000i
>> V'*V % not orthonormal
ans =
1.0000 + 0.0000i 0.0000 - 0.0000i -0.4417 + 0.0108i
0.0000 + 0.0000i 1.0000 + 0.0000i -0.0000 + 0.0000i
-0.4417 - 0.0108i -0.0000 - 0.0000i 1.0000 + 0.0000i
>> A=0.5*(A'+A) % symmetrize it
A =
15.3593 + 0.0000i -7.9491 - 8.0340i 23.9826 +22.6380i
-7.9491 + 8.0340i 9.8955 + 0.0000i -25.9424 + 0.8862i
23.9826 -22.6380i -25.9424 - 0.8862i 76.7452 + 0.0000i
>> A'-A
ans =
0 0 0
0 0 0
0 0 0
>> [V D]=eig(A)
V =
0.3099 + 0.2813i -0.5848 - 0.5812i 0.2770 + 0.2614i
-0.7875 + 0.0165i -0.5382 - 0.0087i -0.2996 + 0.0102i
-0.4521 + 0.0000i 0.1746 + 0.0000i 0.8747 + 0.0000i
D =
1.0000 0 0
0 1.0000 0
0 0 100.0000
>> V'*V % now V is orthonormal
ans =
1.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 - 0.0000i
-0.0000 + 0.0000i 1.0000 + 0.0000i -0.0000 + 0.0000i
0.0000 + 0.0000i -0.0000 + 0.0000i 1.0000 + 0.0000i
>>
Bruno Luong
2018-9-20
The second way is calling SVD which return orthonormal basis in any case, and happens to be eigen vectors as well for Hermitian matrix, even slightly off by numerical error:
[Q,~]=qr(randn(3)+1i*rand(3));
A=Q'*diag([1 1 100])*Q
A =
33.3677 + 0.0000i -30.6395 - 1.8052i 19.3512 +28.9867i
-30.6395 + 1.8052i 30.1043 - 0.0000i -19.9347 -26.3597i
19.3512 -28.9867i -19.9347 +26.3597i 38.5280 + 0.0000i
[V D]=eig(A)
V =
0.8204 + 0.0000i 0.3175 + 0.4756i 0.6104 - 0.1457i
0.3772 - 0.0222i -0.3270 - 0.4325i 0.7529 + 0.0000i
-0.2383 + 0.3569i 0.6157 + 0.0000i 0.1977 + 0.0178i
D =
1.0e+02 *
0.0100 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 1.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0100 + 0.0000i
V'*V % V is NOT orthonormal
ans =
1.0000 + 0.0000i -0.0000 - 0.0000i 0.7440 - 0.1776i
-0.0000 + 0.0000i 1.0000 + 0.0000i -0.0000 + 0.0000i
0.7440 + 0.1776i -0.0000 - 0.0000i 1.0000 + 0.0000i
[U D V]=svd(A)
U =
-0.5718 - 0.0000i -0.0000 + 0.0000i -0.8204 + 0.0000i
0.5413 - 0.0319i 0.4742 - 0.5817i -0.3772 + 0.0222i
-0.3418 + 0.5121i -0.1567 - 0.6421i 0.2383 - 0.3569i
D =
100.0000 0 0
0 1.0000 0
0 0 1.0000
V =
-0.5718 + 0.0000i 0.0000 + 0.0000i -0.8204 + 0.0000i
0.5413 - 0.0319i 0.4742 - 0.5817i -0.3772 + 0.0222i
-0.3418 + 0.5121i -0.1567 - 0.6421i 0.2383 - 0.3569i
V'*A*V
ans =
100.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
-0.0000 - 0.0000i 1.0000 - 0.0000i -0.0000 + 0.0000i
0.0000 - 0.0000i -0.0000 - 0.0000i 1.0000 - 0.0000i
V'*V % V is orthonormal
ans =
1.0000 + 0.0000i -0.0000 + 0.0000i 0.0000 + 0.0000i
-0.0000 - 0.0000i 1.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 - 0.0000i 0.0000 - 0.0000i 1.0000 + 0.0000i
The second way is calling SVD which return orthonormal basis in any case, and happens to be eigen vectors as well for Hermitian matrix
That is true as long as U=V' and a Hermitian matrix always has one SVD satisfying this. However, here the SVD is non-unique. Can we be sure Matlab will always return a Hermitian symmetric U,V pair even in non-unique cases?
Bruno Luong
2018-9-20
编辑：Bruno Luong
2018-9-20
If H is input hermitian and G:=sqrtm(H) (G is hermitian as well)
If [U,D,V] is (any) svd(H), then [U,sqrt(D),V] is svd of G. Therefore by simple calculation V'*D*V is svd of G'*G = G*G = H. QED.
BTW you might ask the same question when you propose using ORTH() to fix V.
Bruno Luong
2018-9-20
编辑：Bruno Luong
2018-9-20
Assuming first H is full-rank, and let define R := V'*U.
From H = H' we have
U*D*V' = V*D*U'
or in other term R*D*R = D;
That means
R(i,:)*D(i,i)*R(:,j) = D(i,i)*delta(i,j)
with delta(.) is kronecker's symbol.
Since D(i,i) ~= 0, we conclude
R(i,:)*R(:,j) = delta(i,j)
Meaning: R^2 = R*R = I;
Thus R = I (let's recall R is orthonormal); meaning
U=V;
If H is not full-rank, U=V is generally not true but still this is true V*D*V' = H. I won't show that.
Christine Tobler
2018-9-20
编辑：Matt J
2018-9-20
If A is symmetric (or complex hermitian), both U and V will be orthogonal matrices. However, they need not be identical if the matrix A is not positive definite.
Example:
>> A = diag([-3, 3]);
>> [U, D] = eig(A)
U =
1 0
0 1
D =
-3 0
0 3
>> [U, S, V] = svd(A)
U =
1 0
0 1
S =
3 0
0 3
V =
-1 0
0 1
This is because S must have nonnegative values on its diagonal, while D can have any real numbers.
This comes down to the step saying that R*R = I implies R = I; in fact, R can be any diagonal matrix with 1 and -1 on the diagonal.
Matt J
2018-9-20
This comes down to the step saying that R*R = I implies R = I; in fact, R can be any diagonal matrix with 1 and -1 on the diagonal.
It can be other things as well, e.g., R=[1,1;1,-1]/sqrt(2)
Bruno Luong
2018-9-20
编辑：Bruno Luong
2018-9-20
Yeap Christine and you catch the mistake.
Actually what the proof need to work is one assumption
R (=V'*U) be in SO(n) and not only on U(n).
So in case H definite positive MATLAB SVD always return R in SO(n).
Lorenzo
2018-9-20
Thank you Bruno and to all in this discussion. In the meantime I'm using schur to diagonalize an hermitian (or hermitian up to numerical errors) matrix as it seems to always return a unitary matrix. Also it seems not slower (actually a tiny bit faster) than eig. Moreover it does not have the problems of svd and I can directly get the eigenvalues (which I need).
This begs the question: shouldn't schur be the method of choice to diagonalize hermitian matrices?
Bruno Luong
2018-9-20
Great find of using SCHUR
Christine Tobler
2018-9-21
Calling SCHUR for a hermitian matrix should be slower than calling EIG. I still think the problem is that your matrix is not exactly hermitian.
Here's an example:
>> n = 1000;
>> A = randn(n) + 1i*randn(n);
>> A = A*diag(rand(n, 1))*A'; % nearly hermitian
>> norm(A - A')
ans =
7.3104e-13
>> tic; [U, D] = eig(A); toc
Elapsed time is 4.432694 seconds.
>> tic; [U, T] = schur(A); toc
Elapsed time is 3.521355 seconds.
>> tic; [U, D] = eig((A+A')/2); toc
Elapsed time is 0.637474 seconds.
Because A is not hermitian, EIG uses the non-hermitian algorithm, which is based on calling SCHUR, and then recovering the eigenvectors from the Schur vectors. Since A is nearly hermitian, its eigenvectors are nearly identical to the Schur vectors.
However, for a hermitian matrix, the Schur decomposition's matrix T is diagonal. There is no need to compute the upper triangle of T, as is done in SCHUR, and the hermitian algorithm is faster because it doesn't do this.
类别
在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息
2016-10-4
2018-9-21
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
