# Smallest magnitude sparse generalized eigenvalues

**URL:** <https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095>\
**Category:** General Usage\
**Created:** [April 15, 2025, 12:51pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095 "2025-04-15T12:51:50Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [April 15, 2025, 12:51pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/1 "2025-04-15T12:51:50Z")

</div>

Hi,

I am trying to find the `neig=5` smallest eigenvalues (in magnitude) of a generalised problem

```julia
(K-val*M)*vec = 0

```

where `K` is real symmetric positive _indefinite_ and `M` is real symmetric positive definite. This is close to a question in a [previous thread](https://discourse.julialang.org/t/smallest-magnitude-eigenvalues-of-the-generalized-eigenvalue-equation-for-a-large-sparse-matrix/75485/11) but here I try to use `KrylovKit.jl`.

My toy problem below could use better forms of `K` and `M` than generalised sparse storage, but my real world problems starts with sparse matrices.

The problem is parametrised by `N` so that the eigenfreqency (I do mechanics) `f`, computed from the eigenvalue is independent of `N` (asymptoticaly, at least). The game: scale up. But moderately: let’s stay on a personal computer for now.

```julia
using SparseArrays,KrylovKit

N = 1000

neig = 5

K = spdiagm(N,N,-1=>range(-N,-N,N-1),0=>range(2N,2N,N), 1=>range(-N,-N,N-1))
M = spdiagm(N,N,0=>range(1/N,1/N,N))

vals, vecs, info = KrylovKit.geneigsolve((K/N,M*N),randn(N),
                        neig,
                        EigSorter(abs; rev = false);                                
                        maxiter = 200,    
                        krylovdim = 100,
                        tol = 1e-4,
                        verbosity = 3,
                        orth = KrylovKit.ModifiedGramSchmidtIR(),
                        issymmetric = true,
                        isposdef = true)

f=sqrt.(vals[1:neig])./(2π)*N

```

`N=1000` works, but I hit the wall before `N=10000` ☹ . I get some gain from scaling the matrices `(K/N,M*N)` so that they have “kind of one” terms… does this mean I need a preconditionner? “Krylov” would suggest that, but the manual states “\*KrylovKit.jl does currently not provide a general interface for including preconditioners”.

I also gain from a lax `tol` - all eigenvalues but the last are returned with better precision.

Any ideas on how to tweak the solver?

---

<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:** [April 15, 2025, 1:01pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/2 "2025-04-15T13:01:40Z")

</div>

For the smallest-magnitude eigenvalues of an indefinite matrix, you usually want to apply Krylov methods to the _inverse_ problem

Kx = \lambda Mx \Longleftrightarrow K^{-1} y = \lambda^{-1} M^{-1} y

where y = Mx, and then find the _largest_ magnitude of \lambda^{-1}. You don’t compute K^{-1} explicitly, but instead form (if possible) its sparse factorization (e.g. sparse LU or sparse LDLᵀ) and use that as a linear operator (e.g. with `InverseMap` in LinearMaps.jl). (Your M matrix is apparently diagonal, which would make it easy to invert explicitly.

---

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [April 15, 2025, 1:42pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/3 "2025-04-15T13:42:17Z")

</div>

Hi Steven,

Thank you. I noted that KrylovKit (and I reckon, iterative eigensolvers in general) would rather deal with the eigenvalues of largest magnitude. And you seem to say that my asking for the _smallest_ eigenvalues is what is giving the solver a hard time. Very important input.

The catch is that my `K` (in general, as opposed to my toy example), is not necessarily invertible (and may have negative eigenvalues). But OK, I’ll ponder what algebraic trick I can pull to transform my problem into a more palatable form (i.e. find largest eigenvalues). I will have to think very carefully about what properties I can assume `K` and `M` to have in general… As you point out, I do not need to _invert_ any matrix, just to provide a matvec for any linear expression.

I’ll take it from there!

😃

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 15, 2025, 1:44pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/4 "2025-04-15T13:44:02Z")

</div>

Would this help?  
[PetrKryslUCSD/VibrationGEPHelpers.jl: Generalized Eigenvalue Problem helper functions](https://github.com/PetrKryslUCSD/VibrationGEPHelpers.jl)

---

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [April 15, 2025, 3:46pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/5 "2025-04-15T15:46:05Z")

</div>

Would this? 😃  
A blueprint on how to use KrylovKit, Arpack and more, to compute eigenmodes in structures!  
Thanks!

---

<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:** [April 15, 2025, 11:21pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/6 "2025-04-15T23:21:22Z")

</div>

> [@Philippe\_Maincon1](#):
>
> The catch is that my `K` (in general, as opposed to my toy example), is not necessarily invertible (and may have negative eigenvalues).

If it’s not invertible, then the problem is even easier: you know that there is an eigenvalue of zero, and just need to solve a linear system to find it.

(Whereas negative eigenvalues are no problem.)

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 15, 2025, 11:54pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/7 "2025-04-15T23:54:55Z")

</div>

Not necessarily: in mechanics the fundamental frequency is always non-zero (ie. the first non-zero after the zeros).

---

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [April 16, 2025, 5:13am UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/8 "2025-04-16T05:13:10Z")

</div>

If an object floats in space, it will have 6 zero eigenfrequencies, and 6 “stiff body motion” eigenvectors. Then I am interested in “a few” more eigenfrequencies - “at what notes does it sound if I hit it”. These will be positive eigenvalues.

Then comes the fact that _if_ I apply boundary conditions, I do not do it in the usual manner of “suppressing” or “condensing” some degrees of freedom, which would (typically) lead to a symmetric-positive-definite K. Rather, I introduce Lagrange multipliers (new degrees of freedom, augmented problem), resulting in an indefinite K (negative eigenvalues, but typically no zero eigenvalue). I have not thought through what story the eigenvectors with negative eigenvalues will tell, and I look forward to discover that.

I’ll explore the ideas in VibrationGEPHelpers now.

---

<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:** [April 16, 2025, 12:50pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/9 "2025-04-16T12:50:33Z")

</div>

> [@Philippe\_Maincon1](#):
>
> f an object floats in space, it will have 6 zero eigenfrequencies, and 6 “stiff body motion” eigenvectors. Then I am interested in “a few” more eigenfrequencies - “at what notes does it sound if I hit it”. These will be positive eigenvalues.

One simple option is to add a tiny shift to your matrix, A + (\text{small}) I to your matrix that makes it invertible. Even better, if you have an order-of-magnitude estimate for the lowest nonzero eigenfrequencies, e.g. from dimensional analysis, you could subtract that estimate as a shift. That would not only make the matrix invertible, it would accelerate convergence of the desired eigenfrequencies.

(For example, a shift strategy to deal with rigid-body modes is mentioned in [this article](https://www.mdpi.com/2075-1702/11/2/2180).)

---

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [April 16, 2025, 1:12pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/10 "2025-04-16T13:12:20Z")

</div>

Whether I use the algorithm corresponding to the option `KrylovKit` in VibrationGEPHelpers, or a modification using an LU factorisation

```julia
luK = lu(K)
L,U,s,p,q = luK.L,luK.U,luK.Rs,luK.p,luK.q
vals, vecs, info = KrylovKit.eigsolve(x->L\((s.*(M*((U\x)[q])))[p]),...
for vecsᵢ ∈ vecs  
    vecsᵢ .= (U\vecsᵢ)[q]
end

```

my toy problem is solved with excellent stability, tested up to `N = 1000000`. Which is an improvement, to say the least. Next step: the real world!

In both cases: the problem is transformed in to an “inverse” problem (with inverted eigen values as suggested by Steven), that is an ordinary (as opposed to generalized) eigenproblem (as suggested by Petr).

Thank you Petr and Steven!

---

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [April 16, 2025, 1:17pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/11 "2025-04-16T13:17:14Z")

</div>

Thank you, I’ll study this.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [April 16, 2025, 1:23pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/12 "2025-04-16T13:23:15Z")

</div>

You can certainly use the GEP for the same purpose.  
Shifting is mentioned in the docs

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/9/09481eb9e378313f9a1b3c778b8801d70399dabb.png)

---

<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:** [April 16, 2025, 11:03pm UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/13 "2025-04-16T23:03:59Z")

</div>

> [@Philippe\_Maincon1](#):
>
> In both cases: the problem is transformed in to an “inverse” problem (with inverted eigen values as suggested by Steven), that is an ordinary (as opposed to generalized) eigenproblem (as suggested by Petr).

The problem with changing to an ordinary eigenproblem is that, if you are not careful, you lose the Hermitian/real-symmetric property. Here, since M is diagonal, you can easily scale by M^{1/2} factors to get a Hermitian ordinary eigenvalue problem if you want.

---

<div class="post-metadata">

**Author:** ![Philippe\_Maincon1](https://avatars.discourse-cdn.com/v4/letter/p/ec9cab/32.png) [@Philippe\_Maincon1](https://discourse.julialang.org/u/Philippe_Maincon1)\
**Post date:** [April 17, 2025, 6:13am UTC](https://discourse.julialang.org/t/smallest-magnitude-sparse-generalized-eigenvalues/128095/14 "2025-04-17T06:13:13Z")

</div>

Oh… the operator `x->L\((s.*(M*((U\x)[q])))[p])` of my inverted problem will generally not be exactly symmetric, even though K and M are.

But KrylovKit, if told `ishermitian=true` stays `Float64`, and gives good results. If not, I could always rewrite the above operator to something like `x->U\(L\(M*x))` (ignoring the details here), which must be symmetric (baring finite precision).

The fact that M is banded, is just my toy problem, not my general problem…

Petr: GEP can do it - no doubt. I am infamous for… learning by doing 😃 . And my linear algebra got upgraded in the past couple of days, thanks to you both.