# Derivative of eigenvalues and eigenvectors of Hermitian Matrix by automatic differentiation

**URL:** <https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563>\
**Category:** General Usage\
**Tags:** question\
**Created:** [June 10, 2018, 5:41am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563 "2018-06-10T05:41:32Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Chong\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chong_wang/32/20307_2.png) [@Chong\_Wang](https://discourse.julialang.org/u/Chong_Wang)\
**Post date:** [June 10, 2018, 5:41am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/1 "2018-06-10T05:41:32Z")

</div>

I am recently looking into automatic differentiation, as implemented by `ForwardDiff.jl`. The manual says functions that call blas are not supported, therefore any functions containing `eig` are not supported. I think I partly understand the reason behind this statement.

However, I noticed that `eigh` (eigen decomposition of Hermitian matrix) is supported by `AlgoPy` ([https://github.com/b45ch1/algopy](https://github.com/b45ch1/algopy)), a package providing automatic differentiation in python. This is quite unexpected since derivatives of eigenvectors are not well defined since eigenvectors can be multiplied by arbitrary number. I wonder how `AlgoPy` supports eigen decomposition and whether it’s doable in julia.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 10, 2018, 6:44am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/2 "2018-06-10T06:44:02Z")

</div>

> [@Chong\_Wang](#):
>
> This is quite unexpected since derivatives of eigenvectors are not well defined since eigenvectors can be multiplied by arbitrary number.

Once you assume a normalization, you can define derivatives. The simplest way is to use implicit differentiation and add the normalization to the constraints. This can be tedious, but I find the [Matrix Cookbook](http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf) and [Old and new matrix algebra useful for statistics](https://tminka.github.io/papers/matrix/) helpful.

I have worked out a subcase of SVD, I am planning to make a PR when I make it work generally.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [June 10, 2018, 7:23am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/3 "2018-06-10T07:23:13Z")

</div>

Hermitian matrices have orthogonal eigenvectors, so the only sensible normalisation is that they each have norm 1, giving an orthogonal matrix Q that diagonalises the matrix. This is unique up to sign / phase.

---

<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:** [June 10, 2018, 7:55am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/4 "2018-06-10T07:55:51Z")

</div>

That’s in the case of numerically well-separated eigenvalues. For repeated/almost-repeated eigenvalues, the naive definition will fail (although perturbation of a repeated eigenvector makes sense if you only use the _subspace_ and not the individual eigenvectors), and you have to project out the components related to the other eigenvectors or you will get NaN/noise. I don’t know if there’s any good, generic way to do this.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [June 10, 2018, 8:14am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/5 "2018-06-10T08:14:28Z")

</div>

It’s still well-defined if you run auto-differentiation on the QR algorithm itself. The only place of discontinuity is choice of when to deflate, but this discontinuity won’t introduce NaNs.

---

<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:** [June 10, 2018, 8:25am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/6 "2018-06-10T08:25:57Z")

</div>

Sure, but still the eigenvector perturbation goes like 1/delta, if delta is the separation between the eigenvalue of interest and the rest of the spectrum. This means that if you’re interested in the subspace and want to compute the derivative of O(1) quantities, eg \<x1, Bx1\> + \<x2, Bx2\> where x1 and x2 are the eigenvectors associated with the smallest eigenvalues of A (this is a typical use case in quantum mechanics) and do it naively on a matrix with |lambda2 - lambda1| = delta, and those two eigenvalues well-separated from the rest of the spectrum, you will be relying on cancellation of 1/delta terms, and the relative accuracy on your output will go as like eps/delta, where eps is the machine epsilon. The correct way to do it is to compute the perturbation of x1 orthogonal to x2, and reciprocally, but this of course assumes that you are not interested in separating x1 and x2, so it cannot be done in a completely generic way.

---

<div class="post-metadata">

**Author:** ![Chong\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chong_wang/32/20307_2.png) [@Chong\_Wang](https://discourse.julialang.org/u/Chong_Wang)\
**Post date:** [June 10, 2018, 8:36am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/7 "2018-06-10T08:36:33Z")

</div>

> [@Tamas\_Papp](#):
>
> This can be tedious, but I find the [Matrix Cookbook](http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf) and [Old and new matrix algebra useful for statistics](https://tminka.github.io/papers/matrix/) helpful.

Thank you. I am aware of the methods in the two links. However, I don’t think it has anything to do with auto-differentiation. Maybe auto-differentiation of QR mentioned by @dlfivefifty is what `Algopy` does. However, if QR is not implemented in julia, as the current situation for `qrfact`, is it still possible to do auto-differentiation easily?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 10, 2018, 8:46am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/8 "2018-06-10T08:46:31Z")

</div>

Maybe I was not clear: AFACT AD for `eig` is currently missing and needs to be _implemented_, hence the suggestions.

---

<div class="post-metadata">

**Author:** ![Chong\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chong_wang/32/20307_2.png) [@Chong\_Wang](https://discourse.julialang.org/u/Chong_Wang)\
**Post date:** [June 10, 2018, 8:48am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/9 "2018-06-10T08:48:16Z")

</div>

I see. Thank you.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [June 10, 2018, 8:56am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/10 "2018-06-10T08:56:50Z")

</div>

> [@Chong\_Wang](#):
>
> However, if QR is not implemented in julia, as the current situation for `qrfact`

It is implemented.

```julia
julia> using ForwardDiff

julia> A = ForwardDiff.Dual{Float64}.(rand(5, 2), ones(5,2))
5×2 Array{ForwardDiff.Dual{Float64,Float64,1},2}:
 Dual{Float64}(0.444375,1.0) Dual{Float64}(0.870728,1.0)
 Dual{Float64}(0.322077,1.0) Dual{Float64}(0.349103,1.0)
 Dual{Float64}(0.225048,1.0) Dual{Float64}(0.287789,1.0)
 Dual{Float64}(0.746621,1.0) Dual{Float64}(0.0858256,1.0)
 Dual{Float64}(0.523616,1.0) Dual{Float64}(0.324125,1.0)

julia> qrfact(A)
Base.LinAlg.QR{ForwardDiff.Dual{Float64,Float64,1},Array{ForwardDiff.Dual{Float64,Float64,1},2}} with factors Q and R:
ForwardDiff.Dual{Float64,Float64,1}[Dual{Float64}(-0.408481,-0.138573) Dual{Float64}(-0.779144,0.319615) … Dual{Float64}(-0.235038,-0.187226) Dual{Float64}(-0.336001,-0.139441); Dual{Float64}(-0.296062,-0.353418) Dual{Float64}(-0.18001,0.0177585) … Dual{Float64}(0.89992,-0.111596) Dual{Float64}(0.235075,-0.0145282); … ; Dual{Float64}(-0.686313,0.392397) Dual{Float64}(0.569669,0.392032) … Dual{Float64}(-0.00134686,-0.094033) Dual{Float64}(-0.449298,-0.0753754); Dual{Float64}(-0.481321,0.0006333) Dual{Float64}(0.0394445,0.227912) … Dual{Float64}(-0.365894,-0.173119) Dual{Float64}(0.793217,-0.103401)]
ForwardDiff.Dual{Float64,Float64,1}[Dual{Float64}(-1.08787,-2.07905) Dual{Float64}(-0.733478,-2.43997); Dual{Float64}(0.0,0.0) Dual{Float64}(-0.733004,-0.174499)]

```

---

<div class="post-metadata">

**Author:** ![Chong\_Wang](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chong_wang/32/20307_2.png) [@Chong\_Wang](https://discourse.julialang.org/u/Chong_Wang)\
**Post date:** [June 10, 2018, 8:57am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/11 "2018-06-10T08:57:57Z")

</div>

Incredible. Thank you.

---

<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:** [June 10, 2018, 9:02am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/12 "2018-06-10T09:02:30Z")

</div>

QR for orthogonal matrix factorization != QR for eigenvalues: [QR algorithm - Wikipedia](https://en.wikipedia.org/wiki/QR_algorithm)

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [June 10, 2018, 11:30am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/13 "2018-06-10T11:30:57Z")

</div>

[https://github.com/andreasnoack/LinearAlgebra.jl](https://github.com/andreasnoack/LinearAlgebra.jl)

Contains generic methods for eigenSelfAdjoint.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [June 10, 2018, 12:04pm UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/14 "2018-06-10T12:04:35Z")

</div>

> [@antoine-levitt](#):
>
> The correct way to do it is to compute the perturbation of x1 orthogonal to x2, and reciprocally, but this of course assumes that you are not interested in separating x1 and x2,

Not sure if this is the same thing as you are suggesting, but if I wanted to know the eigenspace corresponding to an eigenvalue `λ`, I would call `Q = nullspace(A - λ*I)`. I believe if `A` is Hermitian then one can first tridiagonalize and then use a single QR decomposition to determine the nullspace, this should be fine with auto-differentiation.

Here’s a quick implementation:

```julia
julia> λ = [1,1,2,3,4,5,6]; Q = qr(randn(length(λ))).Q; A = Q*diagm(0 => λ)*Q'; # Random matrix with redundant eigenvalues

julia> Q₁, T = hessenberg(A-I);

julia> Q₂ = qr(T).Q;

julia> q = Q₁*Q₂[:,end-1:end];

julia> norm(A*q - λ[1]*q)
8.740803670116966e-15

```

This algorithm will work fine with dual numbers. The apparent ill-posedness of the dimension of the nullspace to perturbation is resolved because it’s really an “epsilon” null-space, as we choose to treat small numbers in the QR decomposition as zero.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 10, 2018, 12:22pm UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/15 "2018-06-10T12:22:44Z")

</div>

> [@dlfivefifty](#):
>
> This algorithm will work fine with dual numbers.

Yes, but I think one would also need to obtain d\lambda.

---

<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:** [June 10, 2018, 12:37pm UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/16 "2018-06-10T12:37:44Z")

</div>

To be more clear:

```julia
B = randn(3,3)
pert = randn(3,3)
pert = pert+pert'
function E(y, δ)
    A = Symmetric(full(Diagonal([1.0, 1.0+δ, 2.0]))+y*pert)
    X = eig(A)[2]
    trace(B*X[:,1:2]*X[:,1:2]')
end

```

Differentiate that wrt y around 0, for representative values of delta = 0.0, 1e-10, 0.1. The method you suggest will give an inaccurate result for delta = 1e-10, because the perturbation to X will be of order 1e10, even though the function is smooth wrt y and the perturbation to E will be of order 1.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [June 10, 2018, 7:48pm UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/17 "2018-06-10T19:48:24Z")

</div>

Ah, in your example you need pivoting to get it to reveal the rank:

```julia
julia> δ = 0.000001; A = full(Diagonal([1.0, 1.0+δ, 2.0]))+y*pert

julia> Q₁, T = hessenberg(A-I);

julia> Q₂ = qr(T, Val(true)).Q;

julia> q = Q₁*Q₂[:,end-1:end];

julia> norm(A*q - q)
8.941524032091245e-5

```

Looking at `T`, the choice of pivots in rank-revealing QR should be robust to perturbation:

```julia
julia> T
3×3 Array{Float64,2}:
 1.08233e-5 3.53086e-5 1.69407e-21
 3.53086e-5 0.0651445 0.246914   
 0.0 0.246914 0.934715 

```

---

<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:** [June 10, 2018, 8:54pm UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/18 "2018-06-10T20:54:14Z")

</div>

> [@antoine-levitt](#):
>
> Sure, but still the eigenvector perturbation goes like 1/delta, if delta is the separation between the eigenvalue of interest and the rest of the spectrum.

This is not generally true. I think you’re thinking of the usual expression for [first-order perturbation theory of the eigenvectors](https://en.wikipedia.org/wiki/Eigenvalue_perturbation), which has a δ in the denominator. However, there is also a matrix element in the numerator, which typically also vanishes with δ if there is an eigenvalue crossing as a function of some parameter, in which case there is no divergence in the derivative.

(There _is_ a √δ singularity in the derivative at points where two eigen_vectors_ merge, i.e. where you have a defective matrix, a 2×2 Jordan block, sometimes called an “exceptional point”; see e.g. Ref. 39 of [this paper](https://www.osapublishing.org/oe/abstract.cfm?uri=oe-25-11-12325). But this can’t happen in a Hermitian problem.)

In any case, I don’t find that textbook eigenvector perturbation-theory expression very useful in practice (for one thing, it requires a sum over _all_ eigenvectors, which is impractical in large sparse/structured problems). Often, you are really trying to compute the derivative of one _scalar function_ (or a few functions) of an eigenvector/eigenvalue, in which case there is a more efficient expression that can be derived by the adjoint method. In particular, if you can compute the derivatives of the _matrix_ that you are computing eigenvalues/eigenvectors of (e.g. by ForwardDiff), then there is a relatively simple formula for derivatives of any function of the eigenvalues and/or eigenvectors, reviewed in section 4 of my course notes:

> **[adjoint.pdf](https://math.mit.edu/~stevenj/18.336/adjoint.pdf)**
>
> 213.98 KB

A complication is that, right at the point of an eigenvalue crossing, both eigenvalues and eigenvectors cease to be differentiable in the usual sense. One option is to use [generalized gradients](https://www.sciencedirect.com/science/article/pii/S0022123685711172) e.g. in [eigenvalue optimization](https://link.springer.com/article/10.1007/BF01742705) (essentially what is called “degenerate perturbation theory” in QM textbooks). Alternatively, in some cases you can formulate eigenvalue optimization problems [as a sequence of SDPs](https://www.osapublishing.org/oe/abstract.cfm?uri=oe-22-19-22632).

---

<div class="post-metadata">

**Author:** ![JaredCrean2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jaredcrean2/32/3574_2.png) [@JaredCrean2](https://discourse.julialang.org/u/JaredCrean2)\
**Post date:** [June 11, 2018, 1:37am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/19 "2018-06-11T01:37:23Z")

</div>

A useful collection of matrix AD results was published [here](https://link.springer.com/chapter/10.1007/978-3-540-68942-3_4), and the author has the pdf on his website [here](https://www.cs.ox.ac.uk/files/723/NA-08-01.pdf).

---

<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:** [June 11, 2018, 11:31am UTC](https://discourse.julialang.org/t/derivative-of-eigenvalues-and-eigenvectors-of-hermitian-matrix-by-automatic-differentiation/11563/20 "2018-06-11T11:31:25Z")

</div>

Yes, and that’s the correct way to differentiate this kind of problem. ForwardDiff, however, does not know that and has to differentiate the full matrix of eigenvectors. I’m not claiming that the problem is intractable when you know what you want to do, I’m just saying that I don’t see how to do it in full generality using ForwardDiff.

Your argument about the numerator does not apply in the case of avoided crossings. If it’s more clear, use the same example as above but with the (essentially equivalent but simpler) matrix

```julia
[y δ 0.0;
 δ -y 0.0;
 0.0 0.0 1.0]

```

The first eigenvector switches from e2 when y \<\< -delta to e1 when y \>\> delta, and its derivative is O(1/delta). Accordingly, if you attempt to differentiate naively (ie differentiate the matrix X of eigenvectors) as a function of y, you will rely on cancellation of O(1/delta) terms and lose `log10(delta)` digits in the result, even though the function itself is perfectly well-defined (indeed, it is constant in this simplified case).