In Fortran/C, for a real symmetric or complex Hermitian matrix A
, I can call Lapack dsyev
/zheev
and the resulting eigenvectors are normalized and mutually orthogonal. What method is Julia using for diagonalizing such matrix? I search the source of eigen.jl
but did not find such dsyev
or zheev
function calls.
So when one use the method eigen()
from LinearAlgebra
to calculate the eigensystem,
eigvals, eigvecs = eigen(A)
Are all calculated eigenvectors eigvecs
normalized and mutually orthogonal, meaning
norm(eigvecs*eigvecs'  Matrix(I, size(A, 1), size(A,1))) < tol
where tol
is a very small number, say 1E15
?
Short answer is yes. This eigen(A)
eventually calls LAPACK.ggev!
which gives both the eigenvalues and eigenvectors.
1 Like
Also, it might be worth noting that the specific behavior of eigen
will depend on the structure of your matrix, and what low level work it actually does are an implementation detail, and won’t be the same for all types of matrices.
The answer is probably yes. But I do not see why LAPACK.ggev!
call guarantees the eigenvectors are orthonormal.
Yes, but what matters here is the implementation detail. I need to know what basic Lapack or other function calls the command eigen()
invoked to be convinced that the eigenvectors are mutually orthogonal.
The reason the eigenvectors are orthonormal is that eigenvectors are always orthonormal. Otherwise, they aren’t eigenvectors.
Also for symmetric/hermetian matrices, it will actually call LAPACK.sygvd!
Not for Hermitian matrices. For the Hermitian case it calls syevr
(i.e. dsyevr
for realsymmetric and zheevr
for complexHermitian), which guarantees that the eigenvectors are orthonormal (up to floatingpoint accuracy).
2 Likes
Well, no. Only eigenvectors associated with different eigenvalues of normal matrices are orthogonal. For repeated eigenvalues, they can be arbitrary. As far as I know eigen
checks if the matrix is hermitian up to machine precision and dispatches to the appropriate routine. You can bypass this by wrapping the matrix in Hermitian
.
5 Likes
To clarify, (the basis of) eigenvectors of normal matrices (e.g. Hermitian) can always be chosen orthonormal even for repeated (multiplicity > 1) eigenvalues, and LAPACK guarantees an orthonormal choice for its Hermitianeigensolver routines.
3 Likes
Right, I should have specified.
Btw I’m not sure whether it actually happens that ggev does give (very) nonorthonormal eigenvectors for (close to) hermitian matrices or if the nonhermitian solver tries to keep eigenvectors as orthogonal as they can be anyway for stability.
Since the algorithm starts by finding the Schur form, I think it follows that it should stay close to orthonormal for nearly Hermitian matrices (where the Schur form is nearly diagonal), even for repeated eigenvalues?
Turns out no, arbitrary small perturbations of matrices with degenerate eigenvalues have nonorthogonal eigenvectors:
julia> A = eigen(I+1e14*randn(2,2)).vectors; norm(A'AI)
1.2919032002219082
That’s just because those eigenvectors are those of randn(2,2) which have no reason to be orthonormal. For orthogonality to hold I think you need the perturbation to be much smaller than the eigenvalue gap. In practice however this kind of trouble probably doesn’t happen very often since to get repeated eigenvalues of hermitian matrices you usually need some symmetry, so the only way you could run into trouble is if that symmetry is not reflected at the floatingpoint arithmetic level. Anyway, the short answer is: wrap in Hermitian
.
3 Likes
I had the same problem of diagonalizing a matrix with many 0 eigenvalues, which is the story behind this question. So to solve the problem of realistic numerical error, It seems that I just need to wrap Symmetric()
for real and Hermitian()
for complex to make sure eigenvectors are orthonormal.
Let me to see if I get everything correct. To battle numeric error of real symmetric matrix or complex hermitian matrix and make sure the eigenvectors are orthonormal, the options are

eigen(Symmetric(A))
or eigen(Hermitian(A))
for real case

eigen(Hermitian(A))
for complex case
If S
is a sparse matrix, to achieve the same goal, one just need to use Matrix
to convert the sparse matrix to a normal matrix type, wrap it using corresponding Symmetric
/Hermitian
method and then use eigen()
.
If S
is sparse, you probably don’t want to be using eigen
since it will be really inefficient. You can take Symmetric
of a sparse matrix without problems, (but eigen(that)
doesn’t work because it would be really inefficient to do so).
To get the eigensystem of a sparse matrix S
, besides using eigen(Symmetric(Matrix(S)))
, do you have better (more efficient) solutions? Apparently eigen(Symmetric(S))
does not work since no method matching.
I think you will want to use Arpack.jl
for this. I’m not 100% sure why we don’t have these methods in Base
.
ARPACK
uses iterative method to perform sparse diagonalization and is suitable to get only a small subset of eigensystem. The performance of ARPACK
as of now seems not so good, it consumes way too memory and is also very slow compared to eigs
in MATLAB.
from what I can tell, eigs
in MATLAB uses ARPACK as well, and also only gives you the largest eigenvalues.
For old version of MATLAB
,Yes. But later MATLAB adopted a new method as far as I know.