# Insufficient orthogonality of eigenvectors

**URL:** <https://discourse.julialang.org/t/insufficient-orthogonality-of-eigenvectors/22193>\
**Category:** Numerics\
**Created:** [March 22, 2019, 12:10pm UTC](https://discourse.julialang.org/t/insufficient-orthogonality-of-eigenvectors/22193 "2019-03-22T12:10:01Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![ivanslapnicar](https://avatars.discourse-cdn.com/v4/letter/i/ecb155/32.png) [@ivanslapnicar](https://discourse.julialang.org/u/ivanslapnicar)\
**Post date:** [March 22, 2019, 12:10pm UTC](https://discourse.julialang.org/t/insufficient-orthogonality-of-eigenvectors/22193/1 "2019-03-22T12:10:02Z")

</div>

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:

```julia
# 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)

```

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [March 22, 2019, 1:12pm UTC](https://discourse.julialang.org/t/insufficient-orthogonality-of-eigenvectors/22193/2 "2019-03-22T13:12:28Z")

</div>

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](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/symmetric.jl#L514) and confirm by calling something like

```julia
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](https://github.com/JuliaLang/julia/pull/31057).