Eigenvalue solvers in NGSolve

7.1. Eigenvalue solvers in NGSolve#

The goal of this section is to solve (generalized) matrix eigenvalue problems to find $u\in\mathbb C^N$, $\lambda\in\mathbb C$ such that

\[ A u = \lambda M u \]

for given matrices $A,M\in\mathbb C^{N\times N}$

Shift-and-invert#

Algorithms like the power iteration or the Arnoldi method approximate eigenvectors corresponding to the largest magnitude eigenvalues. However in wave applications typically the smallest magnitude eigenvalues or ones in a specific region are sought.

The following observation transforms the generalized eigenvalue problem into a problem with the same eigenvectors but eigenvalues sorted by closeness to a given shift:

Let $\lambda_0\in \mathbb C$ be not an eigenvalue. Then the problem to find $u\in \mathbb C^N$, $\tilde\lambda\in\mathbb C$ such that

\[ (A-\lambda_0M)^{-1}Mu = \tilde\lambda u \]

has the same eigenvectors as the original problem with the eigenvalues given by $$ \tilde \lambda = \frac{1}{\lambda-\lambda_0}. $$

Thus the largest magnitude eigenvalues $\tilde \lambda$ of the shifted problem correspond to the ones of the original problem closest to the shift $\lambda_0$.

Arnoldi algorithm#

The idea of the arnoldi algorithm is similiar to the one of the power iteration: for a given starting vector $v_0$ the Krylow space of order $k$ is defined by $$V_k:=\mathrm{span}\{v_0,v_1=M^{-1}A v_0,\ldots,v_k=(M^{-1}A)^kv_0\}$$

After a given number of steps the problem (i.e., the large matrices) is projected onto the (small) Krylow space and the small eigenvalue problem is solved. The resulting approximations are approximations to the largest magnitude eigenvalues.

Shift-and-invert Arnoldi algorithm#

It is common to combine these two ideas. Suppose we want to compute the eigenvalues of the Laplacian on the unit_square. Then the matrices in NGSolve are assembled as usual

from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

fes = H1(mesh,order = 4,complex=True)

u,v = fes.TnT()

A = BilinearForm(grad(u)*grad(v)*dx).Assemble().mat
M = BilinearForm(u*v*dx).Assemble().mat

To use ArnoldiSolver we have to supply a multidim GridFunction to store the eigenvectors. The actual Krylow space is twice as large as the dimension of the GridFunction, however half of the approximated eigenvalues are thrown away (rule of thumb)

gfu = GridFunction(fes,multidim=100)
shift = 1

lam = ArnoldiSolver(A,M, shift=shift, vecs = list(gfu.vecs),freedofs=fes.FreeDofs())

import numpy as np
print(np.array(lam)/np.pi**2)

[-6.74934666e-16-8.78821179e-20j  1.00000000e+00-6.08403476e-16j
00000000e+00+8.05443239e-16j  2.00000000e+00+5.20489798e-16j
00000000e+00+2.01681759e-15j  4.00000000e+00-1.30117264e-15j
00000001e+00+8.01484341e-16j  5.00000001e+00+6.46410877e-16j
00000006e+00-2.47620708e-14j  9.00000016e+00+1.52618897e-15j
00000009e+00-1.10193340e-16j  1.00000002e+01+1.15520533e-14j
00000002e+01-1.08967891e-15j  1.30000007e+01+1.11715711e-14j
30000009e+01+1.88158637e-14j  1.60000028e+01-4.32291758e-13j
60000015e+01+5.12239111e-14j  1.80000042e+01+3.73278413e-13j
70000035e+01-1.92334169e-13j  1.70000021e+01-1.24099531e-15j
00000073e+01-3.71454602e-15j  2.00000056e+01-7.93094869e-14j
60000162e+01-1.65079680e-14j  2.60000288e+01+1.27398572e-14j
50000277e+01+4.62916153e-13j  2.50000142e+01-8.14723057e-15j
50000220e+01+3.73688515e-13j  2.50000187e+01+4.26481558e-14j
90000474e+01-3.28919139e-13j  2.90000339e+01-9.75418779e-14j
20000788e+01+8.60500533e-14j  3.40000783e+01+6.17382516e-14j
40001126e+01+4.86709431e-14j  3.60000803e+01-1.52990145e-13j
60001502e+01+1.10338670e-14j  3.70001838e+01+2.72452232e-13j
70000935e+01+1.20941288e-13j  4.00001506e+01+6.88834080e-14j
00002180e+01-9.75347516e-14j  4.10002660e+01+7.79875926e-15j
10002006e+01+2.63619619e-15j  4.50003008e+01+1.82009881e-13j
50003802e+01-3.21046569e-14j  4.90004482e+01-3.34645685e-14j
90008512e+01-2.61089528e-14j  5.00004637e+01+2.45789066e-13j
00009064e+01+1.23895874e-13j  5.00006775e+01-1.46127832e-13j
20006175e+01+3.89780986e-14j  5.20009507e+01+1.78267542e-13j
30006263e+01-3.36937613e-13j  5.30011080e+01+1.19616276e-13j
80013329e+01-4.42517799e-13j  5.80009120e+01+5.05880085e-13j
10014716e+01-4.78070899e-13j  6.10016770e+01+1.30385657e-13j
40027880e+01+7.43520634e-14j  6.40016106e+01+2.03544840e-13j
80032402e+01+2.86235584e-14j  6.80020988e+01-3.85101131e-13j
50033885e+01-6.86418175e-15j  6.50020873e+01+1.79386381e-13j
50017260e+01+2.21141052e-13j  6.50018360e+01-1.08542399e-13j
20040556e+01-4.11383105e-13j  7.30042300e+01-4.10773212e-15j
30028447e+01-4.15033689e-13j  7.40035331e+01+9.78191276e-14j
40052227e+01-2.27287059e-13j  8.00064645e+01+6.73536110e-13j
00044111e+01+3.78685690e-13j  8.10092967e+01-1.18899667e-12j
10046709e+01-7.65123230e-14j  8.20053210e+01+2.10845626e-13j
20099275e+01+7.33332118e-13j  8.90099731e+01+8.68987104e-13j
90068408e+01-1.73751869e-13j  9.00077004e+01+1.63789643e-12j
00133902e+01-5.97674248e-13j  8.50110027e+01+6.10089734e-13j
50060808e+01+2.03110687e-12j  8.50077161e+01-3.80393249e-11j
50077806e+01+3.08080810e-12j  9.70109077e+01-1.05184791e-09j
70154319e+01-2.74004793e-09j  9.80150627e+01+6.89073404e-11j
01025450e+02-3.13008097e-06j  1.01014447e+02+5.64661220e-06j
04016863e+02+7.27162850e-06j  1.04027639e+02+4.83998586e-06j
00025792e+02+4.20671505e-04j  1.00020378e+02+2.03286777e-04j
00014558e+02-2.74822442e-04j  1.06016864e+02+9.98479517e-05j
06024000e+02+1.04940924e-04j  1.09034108e+02-3.36169721e-05j
09019233e+02+1.02890581e-04j  1.13030593e+02-2.67955914e-03j
13035965e+02+3.42096879e-04j  1.37066144e+02-5.94184599e-03j]

Draw(gfu)

BaseWebGuiScene

Explicit methods#

Explicit eigenvalue solvers in NGSolve include PINVIT and LOBPCG

from ngsolve.solvers import *

help(PINVIT)
help(LOBPCG)

Help on function PINVIT in module ngsolve.eigenvalues:

PINVIT(mata, matm, pre, num=1, maxit=20, printrates=True, GramSchmidt=True)
    preconditioned inverse iteration

Help on function LOBPCG in module ngsolve.eigenvalues:

LOBPCG(mata, matm, pre, num=1, maxit=20, initial=None, printrates=True, largest=False)
    Knyazev's cg-like extension of PINVIT

lam,vecs = PINVIT(A,M,M.Inverse(),num=4,maxit=10000,printrates=False)
print(np.array(lam)/np.pi**2)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 lam,vecs = PINVIT(A,M,M.Inverse(),num=4,maxit=10000,printrates=False)
      2 print(np.array(lam)/np.pi**2)

File ~/bin/ngsolve/lib/python3.10/dist-packages/ngsolve/eigenvalues.py:104, in PINVIT(mata, matm, pre, num, maxit, printrates, GramSchmidt)
    101 asmall = InnerProduct (vecs, mata * vecs)
    102 msmall = InnerProduct (vecs, matm * vecs)
--> 104 ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)
    105 lams = Vector(ev[0:num])
    106 if printrates:

File ~/.local/lib/python3.10/site-packages/scipy/linalg/_decomp.py:576, in eigh(a, b, lower, eigvals_only, overwrite_a, overwrite_b, turbo, eigvals, type, check_finite, subset_by_index, subset_by_value, driver)
    572         lwork_args = {'lwork': lw}
    574     drv_args.update({'uplo': uplo, 'jobz': _job})
--> 576     w, v, *other_args, info = drv(a=a1, b=b1, **drv_args, **lwork_args)
    578 # m is always the first extra argument
    579 w = w[:other_args[0]] if subset else w

KeyboardInterrupt: 

lam,vecs = LOBPCG(A,M,M.Inverse(),num=4,maxit=1000,printrates=False)
print(np.array(lam)/np.pi**2)

[-1.58926725e-16  1.00000000e+00  1.00000000e+00  2.00000000e+00]