# Strategy for finding eigenvalues of an infinite banded matrix

**URL:** https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513
**Category:** Numerics
**Created:** [November 21, 2023, 8:20am UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513 "2023-11-21T08:20:32Z")
**Posts on this page:** 10
**Page:** 1

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [November 21, 2023, 8:20am UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/1 "2023-11-21T08:20:32Z")

</div>

So I’ve got a (Hermitian) like the following, and I want to find its first few eigenvalues:

\begin{pmatrix} \ddots & \ddots &\ddots\\ \ddots & T\_{i,i} &T\_{i+1,i} &T\_{i+2,i} & 0 & 0 & 0 &\\ \ddots & T\_{i,i+1} &T\_{i,i} &T\_{i+1,i} &0 & 0 & 0 &\\ & T\_{i,i+2} &T\_{i,i+1} &T\_{A} &T\_{B} & 0 & 0 &\\ & 0 &0 &T\_{B} &T\_{A} &T\_{i+1,i} &T\_{i+2, i} &\\ & 0 &0 &0 &T\_{i,i+1} &T\_{i,i} &T\_{i+1,i} & \ddots\\ & 0 &0 &0 &T\_{i,i+2} &T\_{i,i+1} &T\_{i,i} &\ddots\\ &&&&& \ddots &\ddots &\ddots \end{pmatrix}

The idea is that it stretches off to infinity in both directions. For the bulk of the matrix, the main diagonal band just consists of T\_{i,i}, and the first diagonal bands above and below that are always T\_{i+1,i} = T\_{i,i+1}, and the next diagonal diagonal bands above and below are T\_{i+2,i} = T\_{i,i+2}.

The only exception is that in the centre of the matrix, there’s this block shown with the T\_{A} and T\_{B}. All the entries in the matrix are negative, and the block around T\_{A} and T\_{B} should be the dominant entries for the first few eigenvalues.

Are there any clever strategies (ideally julia packages!) for finding the eigenvalues of this matrix? My approach thus far has just been brute force making a big banded matrix and adjust the size until things converge, but I strongly suspect there’s a better way to go about this. I know about InfiniteLinearAlgebra.jl but that unfortunately doesn’t seem to support my problem.

Any advice would be appreciated!

---

<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: [November 21, 2023, 9:08am UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/2 "2023-11-21T09:08:56Z")

</div>

You can represent bi-infinite operators in InfiniteLinearAlgebra.jl by interlacing the positive and negative coefficients as a 2x2 block matrix, see eg. an example of a [periodic Schrödinger operator](https://github.com/JuliaLinearAlgebra/InfiniteLinearAlgebra.jl/blob/master/examples/periodicschrodinger.jl).

Probably your best bet then is inverse iteration with Rayleigh shifts. `\` should work with the infinite block matrix and adaptively find the best truncation, so the complexity would be O(n).

There are alternative algorithms for finding eigenvalues using contour integrals. There’s a recent reference for infinite nonlinear eigenvalue problems by [Colbrook and Townsend](https://arxiv.org/abs/2305.01691) that should give some pointers for the linear case.

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [November 21, 2023, 11:17am UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/3 "2023-11-21T11:17:15Z")

</div>

> [@dlfivefifty](#):
>
> You can represent bi-infinite operators in [InfiniteLinearAlgebra.jl](https://juliahub.com/ui/Packages/InfiniteLinearAlgebra) by interlacing the positive and negative coefficients as a 2x2 block matrix, see eg. an example of a [periodic Schrödinger operator](https://github.com/JuliaLinearAlgebra/InfiniteLinearAlgebra.jl/blob/master/examples/periodicschrodinger.jl).

Oh of course, yeah that transformation makes sense, thanks for the help. For those reading, the idea would be that I should re-arrange rows and columns such that my matrix is of the form

\begin{pmatrix} \hat{T}\_{ab} & \hat{T}\_{1} & \hat{T}\_{2} & \hat{0} & \hat{0} \\ \hat{T}\_{1} & \hat{T}\_{0} & \hat{T}\_{1} & \hat{T}\_{2} & \hat{0} \\ \hat{T}\_{2} & \hat{T}\_{1} & \hat{T}\_{0} & \hat{T}\_{1} & \ddots \\ \hat{0} & \hat{T}\_{2} & \hat{T}\_{1} & \hat{T}\_{0} & \ddots \\ \hat{0} & \hat{0} & \ddots & \ddots & \ddots \end{pmatrix}

where

\hat{T}\_{ab} = \begin{pmatrix} Ta & Tb\\ Tb & Ta \end{pmatrix}

\hat{T}\_{0} = \begin{pmatrix} T\_{i,i} & 0 \\ 0 & T\_{ii} \end{pmatrix}

\hat{T}\_{1} = \begin{pmatrix} T\_{i+1,i} & 0 \\ 0 & T\_{i+1,i} \end{pmatrix}

\hat{T}\_{2} = \begin{pmatrix} T\_{i+2,i} & 0 \\ 0 & T\_{i+2,i} \end{pmatrix}

and this is a form that InfiniteLinearAlgebra.jl is able to handle.

* * *

However, the the example Schrödinger example you quote errors for me. I’ve raised an issue here: [Periodic Schrödinger example errors · Issue #153 · JuliaLinearAlgebra/InfiniteLinearAlgebra.jl · GitHub](https://github.com/JuliaLinearAlgebra/InfiniteLinearAlgebra.jl/issues/153)

* * *

> [@dlfivefifty](#):
>
> Probably your best bet then is inverse iteration with Rayleigh shifts. `\` should work with the infinite block matrix and adaptively find the best truncation, so the complexity would be O(n).
> 
> There are alternative algorithms for finding eigenvalues using contour integrals. There’s a recent reference for infinite nonlinear eigenvalue problems by [Colbrook and Townsend](https://arxiv.org/abs/2305.01691) that should give some pointers for the linear case.

Thank you for these suggestions, I’ll investigate them!

---

<div class="post-metadata">

### Author: ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [November 21, 2023, 12:16pm UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/4 "2023-11-21T12:16:20Z")

</div>

Am I reading this right that your problem has 7 parameters, i.e. all your stuff with and `i` in it is independent of `i`?

In that case, you should first figure out the infinite dimensional part by hand, physicist-style, so that you end up with a nice small problem that you then either solve by hand or via numerics.

You do that by writing down the linear recursion for large indices; you solve that via the standard way (i.e. you get a matrix for the recursion, then compute its eigenvalues and eigenvectors).

You next decide how happy you are with solutions (i.e. vectors) that are exponentially (or polynomially) growing. You probably are unhappy with these guys? (the spectrum of an infinite matrix depends on the space you choose! Your question is meaningless without specifying the growth conditions on your vectors)

That gives you conditions: The exploding modes must be vanish. This gives you a bunch of matching conditions in the finite part, where T\_A and T\_B play a role. And ta-da, you have finite number of equations you need to deal with.

Once you have this finite number of equations for us, somebody here can probably recommend you a julia package to deal with it (e.g. to get a full bifurcation analysis or something like that).

Or are you asking for further help with reducing that to a finite dimensional problem? (that’s a math question, not a julia question, but people here will be happy to oblige)

This looks like it would fit well as a quantum mechanics or functional analysis homework exercise.

edit: The same applies if your parameters have known asymptotics. Then the keyword is “stable manifolds”; i.e. you consider this as a dynamical systems problem for `i -> Inf` and `i -> -Inf`.

---

<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: [November 21, 2023, 1:09pm UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/5 "2023-11-21T13:09:34Z")

</div>

> [@Mason](#):
>
> For the bulk of the matrix, the main diagonal band just consists of […] The only exception is that in the centre of the matrix, there’s this block shown […]

In other words, you have an infinite periodic system with a “defect” in the centre. This kind of thing is quite common in physics, e.g. in solid-state physics.

> [@foobar\_lv2](#):
>
> You do that by writing down the linear recursion for large indices; you solve that via the standard way (i.e. you get a matrix for the recursion, then compute its eigenvalues and eigenvectors).

This particular problem has lots of structure that you can exploit.

I would start by analyzing the periodic system without the defect — this has a continuous spectrum which can be characterized via Bloch/Floquet theory. You look for quasi-periodic solutions (periodic solutions x\_n = x\_{n+3} times e^{ikn}), which gives you a small 3x3 eigenproblem for each k, giving you eigenvalues \lambda\_{1:3}(k) describing 3 “bands” for k \in [0,2\pi/3).

Then, to introduce the defect (T\_A and T\_B) there are a variety of strategies. If the desired eigenvalue will appear in a gap of the continuous spectrum (which is often true for extremal eigenvalues in perturbations of 1d-periodic systems under [quite mild assumptions](https://math.mit.edu/~stevenj/papers/ParzygnatLe10.pdf)), then it will be exponentially localized and you can simply truncate the matrix with exponentially small error. More generally, you can match boundary conditions with your quasi-periodic solutions on either side of the defect, obtaining a equations for the solutions (and a nonlinear equation for the eigenvalue) similar to the undergraduate problem of transmission through a potential well or RCWA numerical methods in scattering theory.

(The advantage of starting with the continuous spectrum of the periodic system, i.e. the “essential spectrum”, is that it gives a lot of understanding of what happens when you perturb it, so that you aren’t simply blindly relying on numerical solutions.)

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [November 21, 2023, 2:14pm UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/6 "2023-11-21T14:14:22Z")

</div>

Thank you very much @foobar_lv2 and @stevengj, these are some great hints. I’m used this sort of problem without the defect, but never learned these methods of dealing with it when the defect is included.

---

<div class="post-metadata">

### Author: ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 21, 2023, 7:37pm UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/7 "2023-11-21T19:37:27Z")

</div>

> [@Mason](#):
>
> first few eigenvalues

Does this refer to the smallest or to the largest eigenvalues by magnitude?

---

<div class="post-metadata">

### Author: ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 21, 2023, 8:59pm UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/8 "2023-11-21T20:59:12Z")

</div>

> [@Mason](#):
>
> All the entries in the matrix are negative

Some of the entries are obviously zero, so not negative. So, did you actually mean “all five real parameters are negative” or did you mean “all five real parameters are nonpositive”?

---

<div class="post-metadata">

### Author: ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 23, 2023, 9:02am UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/9 "2023-11-23T09:02:37Z")

</div>

I’m playing with solving truncated finite versions of this algebraically using Groebner.jl. Are there any unstated relations between the magnitudes of the five parameters? For example, a bound on the magnitude of some parameter, or a polynomial inequality that would connect the magnitudes of two different parameters?

FTR, I don’t think this will lead anywhere because the algorithm may be too slow to solve systems corresponding to a large matrix, but I’m curious.

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [November 23, 2023, 10:42am UTC](https://discourse.julialang.org/t/strategy-for-finding-eigenvalues-of-an-infinite-banded-matrix/106513/10 "2023-11-23T10:42:54Z")

</div>

> [@nsajko](#):
>
> Does this refer to the smallest or to the largest eigenvalues by magnitude?

Sorry for the late reply. I mean I want the most negative eigenvalues.

> [@nsajko](#):
>
> did you actually mean “all five real parameters are negative” or did you mean “all five real parameters are nonpositive”?

I mean that all five parameters are negative.

> [@nsajko](#):
>
> Are there any unstated relations between the magnitudes of the five parameters? For example, a bound on the magnitude of some parameter, or a polynomial inequality that would connect the magnitudes of two different parameters?

There are, but their relationships are pretty complicated. In general though, I’d assume that T\_{i,i+2} are similar in magnitude to T\_{i,i+1}, whereas the T\_{i,i} are significantly larger. The region of the parameter space I’m most interested in is the crossover region where the most-negative eigenvalue solutions favour being trapped by T\_{a} and T\_{b}, versus when they are unbound and only see the bulk T values. T\_{a} is slightly more negative than T\_{b} always.
