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

**URL:** <https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893>\
**Category:** General Usage\
**Tags:** package, linearalgebra\
**Created:** [May 10, 2021, 2:49pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893 "2021-05-10T14:49:21Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 10, 2021, 2:49pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/1 "2021-05-10T14:49:21Z")

</div>

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,

```julia
eigvals, eigvecs = eigen(A)

```

Are all calculated eigenvectors `eigvecs` normalized and mutually orthogonal, meaning

```julia
norm(eigvecs*eigvecs' - Matrix(I, size(A, 1), size(A,1))) < tol

```

where `tol` is a very small number, say `1E-15` ?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 10, 2021, 3:15pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/2 "2021-05-10T15:15:55Z")

</div>

Short answer is yes. This `eigen(A)` eventually calls `LAPACK.ggev!` which gives both the eigenvalues and eigenvectors.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 10, 2021, 3:24pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/3 "2021-05-10T15:24:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 10, 2021, 3:45pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/4 "2021-05-10T15:45:48Z")

</div>

The answer is probably yes. But I do not see why ` LAPACK.ggev!` call guarantees the eigenvectors are orthonormal.

---

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 10, 2021, 3:49pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/5 "2021-05-10T15:49:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 10, 2021, 3:55pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/6 "2021-05-10T15:55:34Z")

</div>

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!`

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 10, 2021, 3:57pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/7 "2021-05-10T15:57:43Z")

</div>

> [@Oscar\_Smith](#):
>
> Short answer is yes. This `eigen(A)` eventually calls `LAPACK.ggev!` which gives both the eigenvalues and eigenvectors.

Not for Hermitian matrices. For the Hermitian case [it calls `syevr`](https://github.com/JuliaLang/julia/blob/e4f79b7167312dab649218b53b84c0829ea0c103/stdlib/LinearAlgebra/src/symmetriceigen.jl#L4) (i.e. `dsyevr` for real-symmetric and `zheevr` for complex-Hermitian), which guarantees that the eigenvectors are orthonormal (up to floating-point accuracy).

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [May 10, 2021, 5:12pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/8 "2021-05-10T17:12:45Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 10, 2021, 5:43pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/9 "2021-05-10T17:43:51Z")

</div>

> [@antoine-levitt](#):
>
> Only eigenvectors associated with different eigenvalues of normal matrices are orthogonal. For repeated eigenvalues, they can be arbitrary.

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 Hermitian-eigensolver routines.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [May 10, 2021, 5:53pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/10 "2021-05-10T17:53:48Z")

</div>

Right, I should have specified.

Btw I’m not sure whether it actually happens that ggev does give (very) non-orthonormal eigenvectors for (close to) hermitian matrices or if the non-hermitian solver tries to keep eigenvectors as orthogonal as they can be anyway for stability.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 10, 2021, 6:06pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/11 "2021-05-10T18:06:18Z")

</div>

> [@antoine-levitt](#):
>
> Btw I’m not sure whether it actually happens that ggev does give (very) non-orthonormal eigenvectors for (close to) hermitian matrices or if the non-hermitian 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?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [May 10, 2021, 6:17pm UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/12 "2021-05-10T18:17:15Z")

</div>

Turns out no, arbitrary small perturbations of matrices with degenerate eigenvalues have non-orthogonal eigenvectors:

```julia
julia> A = eigen(I+1e-14*randn(2,2)).vectors; norm(A'A-I)
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 floating-point arithmetic level. Anyway, the short answer is: wrap in `Hermitian`.

---

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 11, 2021, 2:01am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/13 "2021-05-11T02:01:12Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 11, 2021, 2:30am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/14 "2021-05-11T02:30:53Z")

</div>

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()`.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 11, 2021, 2:34am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/15 "2021-05-11T02:34:09Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 11, 2021, 3:19am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/16 "2021-05-11T03:19:22Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 11, 2021, 3:25am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/17 "2021-05-11T03:25:58Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 11, 2021, 4:24am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/18 "2021-05-11T04:24:02Z")

</div>

`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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 11, 2021, 4:28am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/19 "2021-05-11T04:28:10Z")

</div>

from what I can tell, `eigs` in MATLAB uses ARPACK as well, and also only gives you the largest eigenvalues.

---

<div class="post-metadata">

**Author:** ![ZhiyuanYao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zhiyuanyao/32/15488_2.png) [@ZhiyuanYao](https://discourse.julialang.org/u/ZhiyuanYao)\
**Post date:** [May 11, 2021, 5:34am UTC](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893/20 "2021-05-11T05:34:09Z")

</div>

For old version of `MATLAB `,Yes. But later MATLAB adopted a new method as far as I know.

[Next page](https://discourse.julialang.org/t/are-eigenvectors-of-real-symmetric-or-hermitian-matrix-computed-from-eigen-all-orthogonal/60893.md?page=2)