# Is there a way to precondition eigensolver?
**URL:** https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790
**Category:** General Usage
**Tags:** question
**Created:** [February 14, 2019, 4:15pm UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790 "2019-02-14T16:15:55Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@Janis\_Erdmanis](https://discourse.julialang.org/u/Janis_Erdmanis)
#### Post date: [February 14, 2019, 4:15pm UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/1 "2019-02-14T16:15:55Z")
Hello everyone. I am solving a nonlinear matrix equation which can be solved be solved iteratively. A single iteration contains construction of a matrix H, its diagonalization and solution of a nonlinear equation for operator eigenvalues. Since eigenvalues and eigenvectors do not change much between iterations, I am looking on the way to precondition the eigensolver `eigen`. However, I don’t see options for that.
Are there other alternative dense eigensolvers for Julia?
# EDIT
To add a little more detail, the matrix equation we are solving is:
2 \epsilon + a \psi + \psi^3 + \frac{1}{2} (\delta a \psi + \psi \delta a)=0
where a is a number, \epsilon, \delta a are matrices and \psi is the unknown. A following code does solve such problem:
```julia-auto
using LinearAlgebra
using IterativeSolvers
using LazyArrays: ⋆
# The equations we are trying to solve
F(a,δam,ϵ,ψ) = 2 .* ϵ + 1/2 .* (δam*ψ + ψ*δam) + a*ψ + ψ^3
function λ(a,Λ)
p = -9*Λ + sqrt(12*a^3 + 81*Λ^2)
return (-2*3^(1/3)*a + 2^(1/3)*p^(2/3))/(6^(2/3)*p^(1/3))
end
### A relaxation parameter
α = 10.
a = 1.
ωD = 10.
N = 100
ϵ = Diagonal(range(-ωD,stop=ωD,length=N))
δa = zeros(Complex{Float64},N)
δa[1] = 2.1
#δa[60] = 1.2
δam = Array{Complex{Float64}}(undef,(N,N))
for m in 1:N
for n in 1:N
if m-n>0
δam[m,n] = δa[m-n]
elseif n-m>0
δam[m,n] = δa[n-m]'
else
δam[m,n] = 0
end
end
end
function solveψ(a,δam,ϵ,N,α)
ψ = zeros(Complex{Float64},(N,N))
U = Matrix{Complex{Float64}}(I,(N,N))
H = zeros(Complex{Float64},(N,N))
δψ = zeros(Complex{Float64},(N,N))
for i in 1:100
H .= 2 .* ϵ + 1/2 .* (δam*ψ + ψ*δam)
# ### Ordinary eigensolver
# h = eigen(H)
# U = h.vectors
# Λ = h.values
# Diagonalizing as much as possible before giving to eigensolver
# Hp = U'*H*U
# @time h = eigen(Hp)
# global U = U*h.vectors
# Λ = h.values
### IterativeSolvers
iter = LOBPCGIterator(H,true,U)
h = lobpcg!(iter)
Λ = h.λ
δψ .= U*Diagonal((x -> λ(a,x)).(Λ))*U' - ψ
if mod(i,10)==0
@show norm(F(a,δam,ϵ,ψ))
end
### Adaptive relaxation
if norm(F(a,δam,ϵ,ψ+α*δψ))>norm(F(a,δam,ϵ,ψ))
α = α/2
continue
end
ψ .= ψ .+ α .* δψ
end
return ψ
end
@time solveψ(a,δam,ϵ,N,α)
```
Previously I thought that preconditioning would increase performance, but even giving an initial guess (which I previously called preconditioning) does not improve the performance.
Are there other changes I could make to improve the performance of `solveψ` function?
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [February 14, 2019, 8:23pm UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/2 "2019-02-14T20:23:38Z")
There’s a Jula implementation of LOBPCG. That might not be what you have in mind when speaking of preconditioning though. If you give more details about your application you might get better answers.
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [February 14, 2019, 9:25pm UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/3 "2019-02-14T21:25:14Z")
I have used `IterativeSolvers.lobpcg` from IterativeSolvers.jl with `Preconditioners.AMGPreconditioner` from Preconditioners.jl before. The latter just wraps AlgebraicMultigrid.jl. The trick is to use the approximate inverse of A - \sigma B as the preconditioner where \sigma is close enough to the eigenvalues satisfying A x = \lambda Bx. That said, if your eigenvectors and eigenvalues are close enough between iterations, you probably won’t need many iterations to converge if you set the starting point properly, so you may not even need a preconditioner. Note that LOBPCG assumes that B is positive definite.
---
### Author:  [@Janis\_Erdmanis](https://discourse.julialang.org/u/Janis_Erdmanis)
#### Post date: [February 14, 2019, 11:24pm UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/4 "2019-02-14T23:24:12Z")
Thanks. I could get `lobpcg` to solve the eigensystem. However, the time it takes is the same as with Julia built in `eigen` solver even if I plug in the initial guess of eigenvectors from previous iterations.
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [February 15, 2019, 5:59am UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/5 "2019-02-15T05:59:10Z")
LOBPCG is tailored towards sparse/matrix-free operators. As far as I know there is no dense eigensolver that is able to take into account starting guesses (or, equivalently, almost diagonal matrices)
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [February 15, 2019, 6:33am UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/6 "2019-02-15T06:33:21Z")
Starting guesses should be helpful in subspace iteration (block inverse power iteration). And that is applicable to both sparse and dense matrices.
---
### Author:  [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)
#### Post date: [February 15, 2019, 6:39am UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/7 "2019-02-15T06:39:04Z")
Yes, but to reliably beat LAPACK with subspace iteration, that I have not seen (but would love to).
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [February 15, 2019, 6:55am UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/8 "2019-02-15T06:55:15Z")
Certainly true without a good guess of the subspace. But with a good guess of the eigenmodes it may be possible.
---
### Author:  [@msmerlak](https://discourse.julialang.org/u/msmerlak)
#### Post date: [May 18, 2022, 10:33am UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/9 "2022-05-18T10:33:58Z")
@Janis_Erdmanis @antoine-levitt I’ve developed a [method](https://github.com/msmerlak/IterativePerturbationTheory.jl) that does just this: efficiently compute (one, few or all) eigenvectors of a near-diagonal matrix. Useful for your purpose?
---
### Author:  [@Janis\_Erdmanis](https://discourse.julialang.org/u/Janis_Erdmanis)
#### Post date: [May 27, 2022, 9:46am UTC](https://discourse.julialang.org/t/is-there-a-way-to-precondition-eigensolver/20790/10 "2022-05-27T09:46:36Z")
It would have been interesting to try at the time when I was still involved in the project. Now though, I have defended my PhD, and I have chosen to follow a different direction where I no longer need to remember the horrors I experienced working on this project/supervisor.
Fortunately, MWE is self-contained. Thus one can easily test whether it improves performance.