Insufficient orthogonality of eigenvectors

Specific Domains Numerics
1

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
2

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