# Insufficient orthogonality of eigenvectors

In `eigen.jl` function `eigen()` calls `eigen!()` which appears to call
`LAPACK.geevx!()` for all types of matrices (symmetric and non-symmetric).
For symmetric matrices this results in insufficiently orthogonal eigenvectors

• the orthogonality is worse than `norm(A)*eps()`.
For symmetric matrix the wrapper should call `LAPACK.syev!`.
This is illustrated by the following example in Julia 1.0.2:
``````# Example
using LinearAlgebra
import Random
Random.seed!(123)
n=1600
D=rand(n)
a=rand(n)
A=diagm(0=>D)+a*a'
λ,X=eigen(A)
μ,Y=LAPACK.syev!('V','U',deepcopy(A))
# Relative residuals and orthogonality
opnorm(A*X-X*diagm(0=>λ))/opnorm(A), opnorm(X'*X-I),
opnorm(A*Y-Y*diagm(0=>μ))/opnorm(A), opnorm(Y'*Y-I)
``````
1 Like

We do call a symmetric solver when the matrix is symmetric but we don’t call `xsyev`. Instead, we call Dhillon’s `xsyevr` because it’s faster and I thought it had comparable accuracy. You can see the method here and confirm by calling something like

``````julia> λ,X=LAPACK.syevr!('V', 'A', 'L', copy(A), 0.0, 0.0, 0, 0, -1.0);

julia> opnorm(A*X-X*diagm(0=>λ))/opnorm(A), opnorm(X'*X-I)
(8.825242876889729e-16, 2.70267075939732e-12)
``````

Eventually, we should make it possible to select the algorithm. We are already in the process of adding this functionality to `svd`, see https://github.com/JuliaLang/julia/pull/31057.

5 Likes