Are eigenvectors of real symmetric or hermitian matrix computed from `eigen` all orthogonal

General Usage
,
21

I’m not sure whether eigen will do the check and dispatch appropriate method. Here is a test,

julia> a = rand(200, 200); A = a + transpose(a); issymmetric(A)
true

julia> @code_llvm eigen(A)

;  @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/eigen.jl:234 within `eigen'
define nonnull %jl_value_t* @japi1_eigen_957(%jl_value_t*, %jl_value_t**, i32) #0 {
top:
  %3 = alloca %jl_value_t**, align 8
  store volatile %jl_value_t** %1, %jl_value_t*** %3, align 8
  %4 = load %jl_value_t*, %jl_value_t** %1, align 8
;  @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/eigen.jl:235 within `eigen'
  %5 = call nonnull %jl_value_t* @"j_#eigen#69_958"(i8 1, i8 1, %jl_value_t* nonnull %4)
  ret %jl_value_t* %5
}
22

You have to follow the call chain

23

Agree. I opened the source eigen.jl and relevant lines starting from line 234 are

function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
    AA = copy_oftype(A, eigtype(T))
    isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
    return eigen!(AA; permute=permute, scale=scale, sortby=sortby)
end

It seems that symmetric diagonalization method is not dispatched.

24

Go down the rabbit hole again until the eventual blas call

25

In particular, in eigen! there is:

    issymmetric(A) && return eigen!(Symmetric(A), sortby=sortby)
1 Like
26

Great! I’m convinced.

27

Btw that check is exact, ie any arithmetic error will result in the non-symmetric routine being called

1 Like
28

As far as I can tell, the guts of eigs in matlab are still Arpack.

29

Is there a dedicated way to obtain orthogonal eigenvectors of a normal matrix?
Currently, I use this easy trick.

function eigvecs2(M) 
  r=1.0+im+rand()+im*rand()
  return eigvecs(Hermitian(r'*M + r*M'))
end
1 Like