This page was generated from unit-2.2-eigenvalues/pinvit.ipynb.

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.

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\):

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

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
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.

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:

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.

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)

26.34675597290611
2.090269418996443
2.0020128284539975
2.00021345673707
2.0000425921423406
2.0000107252282793
2.0000027865736625
2.000000734794761
2.0000001955368636
2.000000052415031
2.000000014141689
2.000000003853284
2.000000001075648
2.0000000003233405
2.000000000119015
2.0000000000634066
2.0000000000482387
2.0000000000440985
2.0000000000429616
2.000000000042654

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:

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:

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:

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 : [11.962000545545939, 81.00148263192348, 86.174969300157969, 93.248627919198157, 107.36842729454897]
1 : [2.0297798958306665, 7.6530515256426561, 11.936113980365072, 14.251301733990736, 25.19652815778473]
2 : [2.0013614139265532, 5.0837890117240434, 5.5954700798537624, 8.6605528154326201, 12.293023053875208]
3 : [2.0001345995706448, 5.0066489593755872, 5.0240330699133633, 8.1567892520788714, 10.359917249924854]
4 : [2.0000218040290521, 5.0007085579047752, 5.0018502770853868, 8.0551559229546985, 10.106580757936934]
5 : [2.0000044304533584, 5.0001033903743775, 5.0002364704622204, 8.0227837763085645, 10.04000065548658]
6 : [2.0000009358127913, 5.0000185081833042, 5.0000395025462199, 8.0100212497455203, 10.015412064704607]
7 : [2.0000002029488404, 5.0000035117426664, 5.0000080105058062, 8.0044757448363537, 10.005986972653934]
8 : [2.0000000444568888, 5.0000006982581446, 5.0000017969055808, 8.0020116888482704, 10.002328549516426]
9 : [2.0000000098570712, 5.0000001490533839, 5.0000004240423745, 8.0009071959972768, 10.000906471383379]
10 : [2.0000000022263777, 5.0000000362095998, 5.0000001048248146, 8.0004098057309481, 10.000353209610019]
11 : [2.0000000005319776, 5.0000000120163355, 5.000000028843262, 8.0001853439587016, 10.000137757931329]
12 : [2.0000000001529461, 5.0000000066606756, 5.0000000105246611, 8.0000838945403157, 10.000053820396143]
13 : [2.0000000000676144, 5.0000000054485012, 5.0000000060775092, 8.0000380098991482, 10.000021095330789]
14 : [2.0000000000482707, 5.0000000049771511, 5.0000000051834785, 8.0000172415909869, 10.00000833332431]
15 : [2.0000000000438565, 5.0000000047157398, 5.0000000051143552, 8.0000078381062831, 10.000003353733582]
16 : [2.0000000000428422, 5.0000000046506736, 5.0000000050989337, 8.00000357843013, 10.000001410432802]
17 : [2.0000000000426064, 5.0000000046354254, 5.0000000050955409, 8.000001648399472, 10.00000065173904]
18 : [2.0000000000425531, 5.0000000046314153, 5.0000000050939697, 8.0000007736211209, 10.000000355499722]
19 : [2.0000000000425415, 5.0000000046293387, 5.0000000050941091, 8.0000003770649251, 10.00000023979284]

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