# 2.2 Programming an eigenvalue solver¶

We solve the generalized eigenvalue problem

$A u = \lambda M u,$

where $$A$$ comes from $$\int \nabla u \nabla v$$, and $$M$$ from $$\int u v$$, on the space $$H_0^1$$.

This tutorial shows how to implement linear algebra algorithms.

[1]:

import netgen.gui
from netgen.geom2d import unit_square
from ngsolve import *
import math
import scipy.linalg
from scipy import random

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


We setup a stiffness matrix $$A$$ and a mass matrix $$M$$, and declare a preconditioner for $$A$$:

[2]:

fes = H1(mesh, order=4, dirichlet=".*")
u = fes.TrialFunction()
v = fes.TestFunction()

a = BilinearForm(fes)
pre = Preconditioner(a, "multigrid")

m = BilinearForm(fes)
m += u*v*dx

a.Assemble()
m.Assemble()

u = GridFunction(fes)


The inverse iteration is

$u_{n+1} = A^{-1} M u_n,$

where the Rayleigh quotient

$\rho_n = \frac{\left \langle A u_n, u_n\right \rangle}{\left \langle M u_n, u_n\right \rangle}$

converges to the smallest eigenvalue $$\lambda_1$$, with rate of convergence $$\lambda_1 / \lambda_2,$$ where $$\lambda_2$$ is the next smallest eigenvalue.

The preconditioned inverse iteration (PINVIT), see [Knyazef+Neymeyr], replaces $$A^{-1}$$ by an approximate inverse $$C^{-1}$$:

$\begin{split}\rho_n = \frac{\left \langle A u_n, u_n\right \rangle}{\left \langle M u_n, u_n\right \rangle} \\ w_n = C^{-1} (A u_n - \rho_n M u_n) \\ u_{n+1} = u_n + \alpha w_n\end{split}$

The optimal step-size $$\alpha$$ is found by minimizing the Rayleigh-quotient on a two-dimensional space:

$u_{n+1} = \operatorname{arg} \min_{v \in \{ u_n, w_n\}} \frac{\left \langle A v, v\right \rangle}{\left \langle M v, v\right \rangle}$

This minimization problem can be solved by a small eigenvalue problem

$a y = \lambda m y$

with matrices

$\begin{split}a = \left( \begin{array}{cc} \left \langle A u_n, u_n \right \rangle & \left \langle A u_n, w_n \right \rangle \\ \left \langle A w_n, u_n \right \rangle & \left \langle A w_n, w_n \right \rangle \end{array} \right), \quad m = \left( \begin{array}{cc} \left \langle M u_n, u_n \right \rangle & \left \langle M u_n, w_n \right \rangle \\ \left \langle M w_n, u_n \right \rangle & \left \langle M w_n, w_n \right \rangle \end{array} \right).\end{split}$

Then, the new iterate is

$u_{n+1} = y_1 u_n + y_2 w_n$

where $$y$$ is the eigenvector corresponding to the smaller eigenvalue.

## Implementation in NGSolve¶

First, we create some help vectors. CreateVector creates new vectors of the same format as the existing vector, i.e., same dimension, same real/ complex type, same entry-size, and same MPI-parallel distribution if any.

[3]:

r = u.vec.CreateVector()
w = u.vec.CreateVector()
Mu = u.vec.CreateVector()
Au = u.vec.CreateVector()
Mw = u.vec.CreateVector()
Aw = u.vec.CreateVector()


Next, we pick a random initial vector, which is zeroed on the Dirichlet boundary.

Below, the FV method (short for FlatVector) lets us access the abstract vector’s linear memory chunk, which in turn provides a "numpy" view of the memory. The projector clears the entries at the Dirichlet boundary:

[4]:

r.FV().NumPy()[:] = random.rand(fes.ndof)
u.vec.data = Projector(fes.FreeDofs(), True) * r


Finally, we run the PINVIT algorithm. Note that the small matrices $$a$$ and $$m$$ defined above are called asmall and msmall below. They are of type Matrix, a class provided by NGSolve for dense matrices.

[5]:

for i in range(20):
Au.data = a.mat * u.vec
Mu.data = m.mat * u.vec
auu = InnerProduct(Au, u.vec)
muu = InnerProduct(Mu, u.vec)
# Rayleigh quotient
lam = auu/muu
print (lam / (math.pi**2))
# residual
r.data = Au - lam * Mu
w.data = pre.mat * r.data
w.data = 1/Norm(w) * w
Aw.data = a.mat * w
Mw.data = m.mat * w

# setup and solve 2x2 small eigenvalue problem
asmall = Matrix(2,2)
asmall[0,0] = auu
asmall[0,1] = asmall[1,0] = InnerProduct(Au, w)
asmall[1,1] = InnerProduct(Aw, w)
msmall = Matrix(2,2)
msmall[0,0] = muu
msmall[0,1] = msmall[1,0] = InnerProduct(Mu, w)
msmall[1,1] = InnerProduct(Mw, w)
# print ("asmall =", asmall, ", msmall = ", msmall)

eval,evec = scipy.linalg.eigh(a=asmall, b=msmall)
# print (eval, evec)
u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w

Draw (u)

24.673435460042366
2.085481530496192
2.001556570708932
2.0002004276552072
2.000043255889899
2.0000108184520875
2.0000027553853923
2.000000711578449
2.0000001852459497
2.0000000485803833
2.00000001282033
2.000000003420145
2.0000000009379275
2.0000000002805387
2.000000000105925
2.000000000059452
2.000000000047055
2.0000000000437455
2.000000000042863
2.0000000000426237


### Simultaneous iteration for several eigenvalues¶

Here are the steps for extending the above to num vectors.

Declare a GridFunction with multiple components to store several eigenvectors:

[6]:

num = 5
u = GridFunction(fes, multidim=num)


Create list of help vectors, and a set of random initial vectors in u, with zero boundary conditions:

[7]:

r = u.vec.CreateVector()
Av = u.vec.CreateVector()
Mv = u.vec.CreateVector()

vecs = []
for i in range(2*num):
vecs.append (u.vec.CreateVector())

for v in u.vecs:
r.FV().NumPy()[:] = random.rand(fes.ndof)
v.data = Projector(fes.FreeDofs(), True) * r


Compute num residuals, and solve a small eigenvalue problem on a 2 $$\times$$ num dimensional space:

[8]:

asmall = Matrix(2*num, 2*num)
msmall = Matrix(2*num, 2*num)
lams = num * [1]

for i in range(20):

for j in range(num):
vecs[j].data = u.vecs[j]
r.data = a.mat * vecs[j] - lams[j] * m.mat * vecs[j]
vecs[num+j].data = pre.mat * r

for j in range(2*num):
Av.data = a.mat * vecs[j]
Mv.data = m.mat * vecs[j]
for k in range(2*num):
asmall[j,k] = InnerProduct(Av, vecs[k])
msmall[j,k] = InnerProduct(Mv, vecs[k])

ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)
lams[:] = ev[0:num]
print (i, ":", [lam/math.pi**2 for lam in lams])

for j in range(num):
u.vecs[j][:] = 0.0
for k in range(2*num):
u.vecs[j].data += float(evec[k,j]) * vecs[k]

Draw (u)

0 : [13.34372915709268, 73.111393440121432, 86.80910160699402, 100.49611089084627, 109.70970770933609]
1 : [2.0427515735268815, 8.3910124497260625, 13.149594193752973, 20.063722678649157, 28.147570168086272]
2 : [2.0012236523371114, 5.1671574763539052, 7.8706384541163965, 12.311964144125524, 12.758372946430553]
3 : [2.0001026494647833, 5.0053850934540209, 5.3924700875536562, 8.8163910630083642, 10.160731142117031]
4 : [2.0000172984479208, 5.0004810651809253, 5.0350609655998761, 8.1520590840006566, 10.024407683473084]
5 : [2.0000037952077121, 5.0000896264650327, 5.0034580754336888, 8.0311002305530668, 10.006524423130674]
6 : [2.0000008648763741, 5.00001950961222, 5.000408148741716, 8.006648032793434, 10.002173040202404]
7 : [2.0000002033013922, 5.0000044279161253, 5.0000582692600037, 8.0015013228991876, 10.00081422339221]
8 : [2.000000048224535, 5.0000010240633346, 5.0000098819079595, 8.0003714380158932, 10.00031708780746]
9 : [2.0000000115711067, 5.0000002425112671, 5.0000019030905101, 8.0001034997344469, 10.000125645935084]
10 : [2.0000000028143097, 5.0000000606915105, 5.000000400473291, 8.0000338953652932, 10.000050214677737]
11 : [2.0000000007122534, 5.000000018084517, 5.0000000907963065, 8.0000126926737725, 10.000020212594483]
12 : [2.0000000002048997, 5.0000000080364293, 5.0000000240791458, 8.0000052525130645, 10.000008214346929]
13 : [2.0000000000820326, 5.0000000056523337, 5.0000000092173984, 8.0000022950450944, 10.000003401601523]
14 : [2.0000000000521689, 5.0000000050453206, 5.000000005876224, 8.0000010402211217, 10.000001467133856]
15 : [2.000000000044893, 5.0000000047834776, 5.0000000052218132, 8.0000004904075386, 10.000000688574564]
16 : [2.0000000000431135, 5.0000000046706381, 5.0000000051197135, 8.0000002461093676, 10.00000037496177]
17 : [2.000000000042677, 5.000000004639392, 5.0000000050995759, 8.0000001368935489, 10.000000248552418]
18 : [2.0000000000425722, 5.0000000046320769, 5.0000000050953508, 8.0000000879268836, 10.000000197578768]
19 : [2.0000000000425469, 5.0000000046296869, 5.0000000050940852, 8.0000000659422206, 10.000000177016078]


The multidim-component select in the Visualization dialog box allows to switch between the components of the multidim-GridFunction.