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.

from ngsolve import *
from ngsolve.webgui import Draw

import scipy.linalg

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

/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

We setup a stiffness matrix \(A\) and a mass matrix \(M\), and declare a preconditioner for \(A\):

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

a = BilinearForm(grad(u)*grad(v)*dx)
pre = Preconditioner(a, "multigrid")
a.Assemble()

m = BilinearForm(u*v*dx).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, w, Mu, Au, Mw, Aw = u.vec.CreateVectors(6)

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)
r.SetRandom()
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 / (pi**2))
    # residual
    r.data = Au - lam * Mu
    w.data = pre.mat * r.data
    # w.data = 1/Norm(w) * w
    w *= 1/Norm(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);

21.904404687835495
2.0923100469358484
2.001383689323411
2.0001760769571435
2.0000343152091773
2.0000072299731357
2.0000015502753525
2.000000335620112
2.0000000730586853
2.000000016011766
2.000000003550711
2.0000000008201937
2.00000000021982
2.00000000008751
2.000000000058278
2.0000000000518074
2.000000000050372
2.0000000000500533
2.00000000004998
2.0000000000499645

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 v in u.vec.CreateVectors(2*num):
    vecs.append(v)

for v in u.vecs:
    r.SetRandom()
    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/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 : [12.282227127151252, 81.99992911901202, 91.88901857050824, 104.46884797170395, 118.86745477422326]
1 : [2.053653132433193, 7.283738007133289, 9.223445127135653, 19.786812497260865, 22.139470286251115]
2 : [2.0016408135252726, 5.098538312694389, 5.158130693843754, 9.670403101888818, 11.714654058450128]
3 : [2.0002546080703176, 5.005794119945728, 5.01441914109119, 8.367187346589002, 10.401762259797286]
4 : [2.0000516867144116, 5.000510112486523, 5.001518457988232, 8.111459016213129, 10.13513023490231]
5 : [2.0000112289812426, 5.000072419666186, 5.000181950172341, 8.040423567767577, 10.052179015338654]
6 : [2.0000024840525383, 5.000013283773906, 5.000025493474269, 8.01627943569105, 10.020551253121189]
7 : [2.0000005572115285, 5.000002622076442, 5.00000426217719, 8.006931988818945, 10.008152192727056]
8 : [2.0000001259423077, 5.000000508393762, 5.000000846256888, 8.003035325530684, 10.00323981224326]
9 : [2.000000028656636, 5.000000102185064, 5.000000189311147, 8.0013458970792, 10.001289365422071]
10 : [2.0000000065778663, 5.000000024647198, 5.000000047254059, 8.000600207946556, 10.0005136214664]
11 : [2.000000001545002, 5.000000009500508, 5.000000014869383, 8.000268466870695, 10.00020482434432]
12 : [2.0000000003934337, 5.00000000643932, 5.000000007311151, 8.000120272541873, 10.000081784670416]
13 : [2.0000000001291043, 5.000000005459627, 5.000000005867427, 8.00005394879776, 10.00003273188681]
14 : [2.000000000068246, 5.00000000504763, 5.000000005718079, 8.000024228916963, 10.000013163452858]
15 : [2.0000000000541975, 5.0000000049467905, 5.000000005687488, 8.000010903066325, 10.000005354022317]
16 : [2.0000000000509446, 5.000000004922414, 5.000000005680399, 8.00000492473208, 10.000002236038407]

17 : [2.0000000000501905, 5.0000000049164095, 5.00000000567902, 8.000002241849714, 10.00000099079804]
18 : [2.0000000000500138, 5.00000000491473, 5.000000005678825, 8.000001037461493, 10.000000493320593]
19 : [2.0000000000499734, 5.000000004914405, 5.000000005678364, 8.000000496678462, 10.000000294532889]

The multidim-component select in the Visualization dialog box allows to switch between the components of the multidim-GridFunction (Netgen-gui only).

Implementation using MultiVector

The simultaneous iteration can be optimized by using MultiVectors. These are arrays of vectors of the same format. You can think of a MultiVector with m components of vector size n as an \(n \times m\) matrix.

• a MultiVector consisting of num vectors of the same format as an existing vector vec is created via MultiVector(vec, num).

• we can iterate over the components of a MultiVector, and the bracket operator allows to access a subset of vectors

• linear operator application is optimized for MultiVector

• vector operations are optimized and called as mv * densematrix: \(x = y * mat\) results in x[i] = sum_j y[j] * mat[j,i] (where x and y are 'MultiVector’s, and mat is a dense matrix)

• pair-wise inner products of two MultiVectors is available, the result is a dense matrix: InnerProduct(x,y)[i,j] = InnerProduct(x_j, y_i). Before ngsolve-version 6.2.2301-242, the result was the transposed matrix.

• mv.Orthogonalize() uses modified Gram-Schmidt to orthogonalize the vectors. Optionally, a matrix defining the inner product can be provided.

• with mv.Append(vec) we can add another vector to the array of vectos. A new vector is created, and the values are copied

• mv.AppendOrthogonalize(vec) appends a new vector, and orthogonalizes it against the existing vectors, which are assumed to be orthogonal.

uvecs = MultiVector(u.vec, num)
vecs = MultiVector(u.vec, 2*num)

for v in vecs[0:num]:
    v.SetRandom()
uvecs[:] = pre * vecs[0:num]
lams = Vector(num * [1])

for i in range(20):
    vecs[0:num] = a.mat * uvecs - (m.mat * uvecs).Scale (lams)
    vecs[num:2*num] = pre * vecs[0:num]
    vecs[0:num] = uvecs

    vecs.Orthogonalize() # m.mat)

    asmall = InnerProduct (vecs, a.mat*vecs)
    msmall = InnerProduct (vecs, m.mat*vecs)

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

    uvecs[:] = vecs * Matrix(evec[:,0:num])

0 : [282.18581484482945, 888.3473212212376, 932.0000741821314, 994.9016375381786, 1094.325839330091]
1 : [2.008741149057017, 8.124820085794749, 8.649649692520288, 17.0759075512237, 21.169320225202185]
2 : [2.0001326550027776, 5.104323161656279, 5.199614785284877, 8.427570263744133, 10.913848743411808]
3 : [2.000015944572217, 5.00569779393365, 5.019855046909897, 8.044264225073231, 10.247421690695443]
4 : [2.0000030502909336, 5.000454283190818, 5.002190420599743, 8.009314561171152, 10.085061309507664]
5 : [2.0000006301643105, 5.000055900102042, 5.000268920081747, 8.00286975501685, 10.031607338130547]
6 : [2.000000134440842, 5.000009866443959, 5.000037995188132, 8.001062702160999, 10.012080387964293]
7 : [2.0000000292802897, 5.000002066773891, 5.000006198389284, 8.000429781056999, 10.004665446917265]
8 : [2.0000000064909114, 5.000000465286217, 5.000001153081171, 8.000182390867032, 10.001810171756057]
9 : [2.0000000014855384, 5.000000110814313, 5.000000238647958, 8.000078907791982, 10.000704743905999]
10 : [2.0000000003726153, 5.000000030037269, 5.000000055531016, 8.000034444463514, 10.000274873135252]
11 : [2.0000000001230527, 5.0000000113311325, 5.000000016342957, 8.000015109601318, 10.000107512651228]
12 : [2.000000000066622, 5.000000006891378, 5.000000007664294, 8.000006656365239, 10.000042167339084]
13 : [2.0000000000537796, 5.000000005477917, 5.000000006044898, 8.000002954253397, 10.000016636426409]
14 : [2.0000000000508398, 5.000000005048776, 5.0000000057645595, 8.00000132941281, 10.000006640165635]
15 : [2.0000000000501634, 5.000000004945979, 5.0000000056993885, 8.000000615951656, 10.000002724377104]
16 : [2.0000000000500067, 5.000000004922144, 5.000000005683639, 8.000000302241785, 10.000001187959652]
17 : [2.00000000004997, 5.000000004916688, 5.000000005680261, 8.000000164269744, 10.000000584899265]
18 : [2.0000000000499623, 5.0000000049149005, 5.000000005678951, 8.000000103524545, 10.000000347887639]
19 : [2.000000000049961, 5.0000000049145354, 5.000000005678877, 8.000000076775622, 10.00000025471206]

The operations are implemented using late evaluation. The operations return expressions, and computation is done within the assignment operator. The advantage is to avoid dynamic allocation. An exception is InnerProduct, which allows an expression in the second argument (and then needs vector allocation in every call).