# Shift-Inverse diagonalization in Julia

**URL:** <https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994>\
**Category:** Performance\
**Created:** [September 29, 2022, 3:41pm UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994 "2022-09-29T15:41:39Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [September 29, 2022, 3:41pm UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/1 "2022-09-29T15:41:39Z")

</div>

What is the best technique to use the shift-inverse diagonalisation technique in Julia? I have used `Arpack.jl` but it doesn’t seem to be much faster; I am not sure how it handles the inverse part. So I would like to understand a better technique to get some eigenpairs of a large sparse hermitian matrix around a specified point in the spectrum, say the middle.

Furthermore, I need this in my research related to a random field spin chain model, where there is a random parameter in my Hamiltonian, the matrix I want to diagonalise. And because of the randomness, I need to diagonalise several instances of the Hamiltonian, and finally calculate some averages. So I was wondering if there is a nicer way to do this parallelly?

Thanks.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [September 29, 2022, 3:44pm UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/2 "2022-09-29T15:44:23Z")

</div>

Perhaps you should look into [GitHub - JuliaLinearAlgebra/ArnoldiMethod.jl: Implicitly Restarted Arnoldi Method, natively in Julia](https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl), they include [a section](https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable/usage/02_spectral_transformations.html) on shift-inverse diagonalization

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [September 29, 2022, 4:06pm UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/3 "2022-09-29T16:06:03Z")

</div>

Thanks for letting me know about this library. It seems it’ll work for me. However, I was also wondering if you could tell me whether this technique is threadsafe. I mean if I have to do this diagonalization many times with a slightly different matrix and each time I calculate some quantity out of the obtained eigenvectors.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [September 29, 2022, 4:37pm UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/4 "2022-09-29T16:37:43Z")

</div>

Yes, I think this package was developed as a thread-safe alternative to ARPACK

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [September 30, 2022, 4:18am UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/5 "2022-09-30T04:18:22Z")

</div>

Ok thanks!

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [October 7, 2022, 9:42am UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/6 "2022-10-07T09:42:02Z")

</div>

Hi, sorry to bother you again. The method described in the reference your provided fails while using `sparse matrices`. Since I am dealing with large matrices and many of them, it saves me a lot of memory allocation if I make use of the sparsity in the matrix. So is there any way to make it work? I tried to convert the matrix into dense form before factorization, which, if fact, helps. However, I was wondering if there was a better way.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [October 7, 2022, 9:47am UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/7 "2022-10-07T09:47:28Z")

</div>

This is strange, sounds like a bug. Could you post the error message? It should definitely work for sparse matrices (barring numerical issues like invertibility)

---

<div class="post-metadata">

**Author:** ![devanshu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/devanshu/32/37104_2.png) [@devanshu](https://discourse.julialang.org/u/devanshu)\
**Post date:** [October 7, 2022, 9:57am UTC](https://discourse.julialang.org/t/shift-inverse-diagonalization-in-julia/87994/8 "2022-10-07T09:57:29Z")

</div>

I get this error while doing

```julia
F = factorize(A - s * I)
A_lmap = LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, x), size(A, 1), ismutating=true)

decomp, = partialschur(A_lmap, nev=n, tol=1e-5, restarts=100, which=LM())
λs_inv, X = partialeigen(decomp)

```

> ERROR: LoadError: MethodError: no method matching ldiv!(::SuiteSparse.CHOLMOD.Factor{Float64}, ::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})

However, the error’s gone if I use `lu()` instead of `factorize()`, nevertheless the memory allocation is very large, compared to converting first to a dense matrix!