FYS415 - Jacobi and Polynomial expansion eigenvalue algorithms - project2 - Solved

Abstract
We will study two algorithms for solving the eigenvalue problem Au = λu; Jacobi’s algorithm and an algorithm by polynomial expansion. We test the algorithms by solving two second order differential equations, first evaluating their correctness with a tridiagonal Toeplitz matrix A, then their overall performance with regards to precision, speed and memory usage on a matrix having non-constant central diagonals. We will see that the Polynomial expansion algorithm exceeds the performance of the Jacobi algorithm on all metrics.

1           Introduction
This work is my submission for the second project in the course FYS4150 given at University of Oslo, autumn 2020 [1],[5]. Parts a)-d) plus g) have been answered. I have implemented the algorithms in C++ using the Armadillo library. You will find the source code at my repository https://github.com/ fonstelien/FYS4150/tree/master

We will study two iterative algorithms for numerical solutions to eigenvalue problems Au = λu and apply them to second order differential equations with Dirichlet boundary conditions on the form

                                                            ,                                                      (1)

and

                                                  .                                          (2)

The first algorithm is the Jacobi eigenvalue algorithm, and the second is a solution based on iterative expansion of the characteristic polynomial p(x) = det(A − xI). I will first introduce the mathematical background, then move over to presenting how the algorithms can be implemented, before I evaluate their performance relative to each other.

1.1         Mathematical basis for the Jacobi eigenvalue method
Given the eigenvalue problem Au = λu where A is symmetrical. The Jacobi algorithm aims to diagonalize A by iterative application of pairwise Givens rotations. A is rotated about an orthogonal matrix S and its transpose S| , such that

S|A = λS|u ⇒ S|AS(S|u) = λS|u.

With S = SmSm−1 ···S1, the Sis will be chosen such that

S|AS = D,

and hence

D(S|u) = λI(S|u) ⇒ Dn×n = diag(λ1,λ2,...,λn).

See [4]. The Givens matrices [6] Si are given by

                                                                                            cosθ      0       ···      −sinθ

                                                                   Ik..                                                                   
Si = Skl 0 G                                                                                , where G.

                     0            0            I                                            .              .         .             . 

                                                                                               sinθ       0       ···       cosθ
After the first dual rotation, we obtain S|1AS1 = B, where Bi is symmetrical with elements

                 , where j 6= k,l.

As outlined in [4], choosing k and l such that the Givens matrix targets A’s largest off-diagonal element akl, will move Bi closer to diagonal form for every iteration. We will therefore choose the angle θ such that bkl = blk = 0, find Bi’s largest off-diagonal element and repeat the process until it converges for some arbitrary ε > maxk6=l Bi.

With v = S|u and Dv = λv we see that the transformation preserves A’s eigenvalues such that D is similar to A, but that D’s eigenvectors v are rotated by S| relative to A’s. However, the transformation preserves the dot product and orthogonality:

v|v = (S|u)|(S|u) = u|(SS|)u = u|u.

Since A is symmetric, its eigenvectors u = [u1 u2 ...] are orthogonal and therefore, assuming normalization, we have that

vi|vj = u|i uj = δij,

where δij is the Kronecker delta.

1.2         Mathematical basis for the Polynomial expansion method
A solution to the eigenvalue problem Au = λu, with A ∈ Rn×n, can always be obtained by solving

                                                        p(x) = det(A − xI) = 0,                                                   (3)

where p(x) is the characteristic polynomial. However, for arbitrary n, an analytic solution to (3) may not always be possible [4]. A numerical approximation is therefore necessary. From Gerschgorin’s circle theorem [6] we know in what area the eigenvalues of matrix A lies in:

 .

The direct approach would therefore be search all rows i and to establish the limits [xmin,xmax], and to iterate over this range with some step length h. However, as shown in [6], we also know that if A ∈ Rn×n is irreducible, tridiagonal and symmetric, the roots of pk(x) = det(Ak − xIk), where Ak is the upper left k × k corner of A, are arranged such that they separate the roots of pk+1(x):

 .

Hence, with an iterative approach where we increment k, we will be able to narrow down the area to search for each of the roots of p(x) (3).

2           Methods
2.1         Jacobi eigenvalue algorithm
As we stated above, if A is symmetrical and S1 is a Givens matrix, S|1AS1 = B Bi, where A is symmetrical. Following [5] we can set s = sinθ, s = cosθ, t = s/c, and rewrite the results from (3) such that

                               bjj =ajj                                                                

bjk =ajkc − ajls  
                               bjl =ajlc + ajks                       

                                                                                                , where j 6= k,l.                       (4)

bbbkkklll ==0=aakkllcc22θ−+ 22aaklklcscs++aallkks2s2

bkl = 0 is achieved by letting

 .

 √ 

Now, with, we get the roots t = −τ ± τ2 + 1. For numerical stability, we must be careful to pick the right root:

 ,

√ 

and at last we can update the B’s elements in (4) with c = 1/ t2 + 1 and s = tc. Here we also see the algorithm’s greatest weakness, that it will set elements that it earlier had put to zer, to some non-zero value, meaning that the process potentially has to be repeated for the same element again later in the process.

A possible implementation is outlined in Listing 1 below. The Jacobi algorithm runs in O(n3) time, and converges typically after 12n3 to 20n3 operations [4], depending on A and our choice of tolerance ε > maxk6=l Bi.

To verify the algorithm we will apply it on the problem in (1), which has a tridiagonal Toeplitz representation

                                             A                                       (5)

where h = 1.0/n. See [5] for outline. With ε = 1.0 · 10−2, the algorithm performs 113 rotations on a 10 × 10 matrix and gives very good precision even if the tolerance is not too strict, as we see in Table 1 below.

Listing 1: Jacobi eigenvalue algorithm
 

A = Symmetric NxN matrix k,l = row, col indexes eps = tolerance for convergence k,l = MAX(A) \\ Finds max element in A and returns row, col indexes

\\ Rotating

WHILE A(k,l) > eps DO \\ Finding the right rotation a_kk = A(k,k) a_kl = A(k,l) a_ll = A(l,l)

tau = (a_ll - a_kk) / (2*a_kl)

IF (tau > 0) t = 1/(tau + sqrt(tau^2 + 1))

Table 1: Results from running the Jacobi eigenvalue algorithm on a 10 × 10 Toeplitz matrix with tolerance ε = 1.0·10−2. The precision is very high, around 7-8 leading digits, even for relatively high ε. The algorithm ran 113 times before it converged.

λi
exact
numeric
0
8.1014052771e+00
8.1014062027e+00
1
3.1749293434e+01
3.1749318745e+01
2
6.9027853211e+01
6.9027858255e+01
3
1.1691699740e+02
1.1691701506e+02
4
1.7153703235e+02
1.7153703274e+02
5
2.2846296765e+02
2.2846301971e+02
6
2.8308300260e+02
2.8308304077e+02
7
3.3097214679e+02
3.3097216033e+02
8
3.6825070657e+02
3.6825073262e+02
9
3.9189859472e+02
3.9189859742e+02
ELSE t = 1/(tau - sqrt(tau^2 + 1)) END IF

c = 1/SQRT(1 + t^2) s = c*t

\\ Updating A

A(k,k) = a_kk*c^2 - 2*a_kl*s*c + a_ll*s^2 A(k,l) = A(l,k) = 0.

A(l,l) = a_ll*c^2 + 2*a_kl*s*c + a_kk*s^2

FOR i = 0 ... N-1 DO

IF NOT (i == k || i == l) a_ik = A(i,k) a_il = A(i,l)

A(i,k) = A(k,i) = a_ik*c - a_il*s

A(i,l) = A(l,i) = a_il*c + a_ik*s

END IF END DO k,l = MAX(A)

END DO

 

2.2         Polynomial expansion algorithm
We saw above that the roots of p(x) = det(A − xI) all lie in the area

n

|x − aii| ≤ X|aij|.

j=1 j6=i

If we limit our algorithm to positive definite A (symmetric and diagonally dominant [3]), and that every k × k sub-matrix Ak is also positive definite, we can define the limits for where to search for the roots of pk(x) as

                                ; and .                         (6)

Our algorithm will start with k = 2, whose roots are found directly, then incrementally search for the roots of p(k + 1)(x) within the limits in (6) by bisection.

A possible implementation for solving the eigenvalue problem in (1) with the matrix in (5) is outlined in Listing 2 below. As the Jakobi algorithm, the Polynomial expansion algorithm runs in O(n3) time.

In routine p(k, x) we have used the iterative calculation of

                                               pk(x) = (d − x)pk−1(x) − e2pk−2(x),                                         (7)

and to avoid overflow, we pull the 1/h2 factor in A (5) out of the calculations and instead multiply the resulting eigenvals by this factor to obtain the right answers. However, for the problem in (2), where the diagonal elements of A are on the form 2/h2 + (ih)2, avoiding numerical imprecision and overflow is not possible with the characteristic polynomial on the form that it has in (7), so we have to divide pk(x) by pk−1(x) to achieve the numerically better form

                                                    qk(x) = (d − x) − e2/qk−1(x).                                              (8)

See [2] for further details.

Verification of the algorithm is again done on a 10 × 10 matrix, and we see in Table 2 that it produces good results, with 4-6 leading digits precision for a root search ε = 1.0 · 10−6. For the root search, maxiter is necessary in case the range does not converge due to numerical imprecision, but does not have to be set very high since the search range is divided by 2 for every iteration and quickly approaches 0. In Table 2 maxiter was set to 50.

Listing 2: Polynomial expansion eigenvalue algorithm for Toeplitz matrix
 

N = rows in the A matrix eigenvals = array of length N

\\ Finding roots by polynomial expansion

Table 2: Results from running the Polynomial expansion eigenvalue algorithm on a 10 × 10 Toeplitz matrix with tolerance ε = 1.0 · 10−6 for the root search. Precision is around 4-6 leading digits.

λi
exact
numeric
0
8.1014053e+00
8.1014062e+00
1
3.1749293e+01
3.1749319e+01
2
6.9027853e+01
6.9027858e+01
3
1.1691700e+02
1.1691702e+02
4
1.7153703e+02
1.7153703e+02
5
2.2846297e+02
2.2846302e+02
6
2.8308300e+02
2.8308304e+02
7
3.3097215e+02
3.3097216e+02
8
3.6825071e+02
3.6825073e+02
9
3.9189859e+02
3.9189860e+02
eigenvals[0] = d - sqrt(-e); \\ roots for k=2 eigenvals[1] = d + sqrt(-e);

\\ Finding roots by polynomial expansion

FOR k = 3 ... N-1 DO

\\ Find the first k-1 roots x_min = 0. i = 0

WHILE i < k-1 DO x_max = eigenvals[i] eigenvals[i] = bisection_root_search(k, x_min, x_max) x_min = x_max i++

END DO

\\ find the kth root x_min = x_max x_max = d + 2*ABS(e)

eigenvals[i] = bisection_root_search(k, x_min, x_max)

END DO

\\ Avoiding overflow for N > approx. 90 h = 1./N \\ step length eigenvals = eigenvals / h^2

\\ Finds the root of pk(x) in the range [x_min, x_max]

ROUTINE bisection_root_search(k, x_min, x_max) eps = convergence tolerance maxiter = maximum number of iterations

x_left = x_min x_right = x_max

i = 0

WHILE (x_right-x_left)/(x_max-x_min) > eps && i < maxiter DO x_mid = (x_right + x_left) / 2

IF p(k, x_mid) == 0 DO RETURN x_mid

END DO

IF p(k, x_left)*p(k, x_mid) < 0 DO x_right = x_mid

ELSE x_left = x_mid END DO

END DO

RETURN x_right

END ROUTINE

\\ Calculates the characteristic polynomial pk(x)

ROUTINE p(k, x)

\\ NOTE: We have taken out the h^2 factor from d,e to avoid overflow for N > approx. 90 d = the center diagonal element of A (constant) e = the first off-diagonal elements of A (constant)

pk_2 = d - x

IF k == 1 DO

RETURN pk_2 END DO

pk_1 = (d - x)*pk_2 - e*e

IF k == 2 DO

RETURN pk_1

END DO

FOR i = 3 ... k DO pk = (d - x)*pk_1 - e^2*pk_2 pk_2 = pk_1 pk_1 = pk i++

END DO

RETURN pk

END ROUTINE

 

3           Results
Figure 1 shows how the number of rotations before convergence increases for increasing size of the problem. We see that it increases exponentially. In Figure 2 we see how the eigenvector has been rotated in the process.

 

Figure 1: Number of rotations needed for the Jacobi algorithm to converge relative to the size of the problem. ε = 1.0 · 10−4.

If we test the Jacobi algorithm on the problem in 2, where the diagonal elements are given on the form [5]

aii = 2/h2 + (ih)2
we see that the Jacobi algorithm quickly starts to be too time- and potentially memory-consuming demanding, since the A matrix must be extended to a large n before we get satisfactory results for the lower eigenvalues. And in contrast to the simple problem with constant diagonal elements, we must also demand a smaller tolerance ε. See Table 3, where a 500×500 matrix has been solved with tolerance ε = 1.0 · 10−6, and some values for ρmax. We see that the precision decreases for growing ρmax, as could be expected since the step length h grows. CPU time is more or less constant, depending only on n and ρmax.

The Polynomial expansion algorithm performs better for this type of problem. In Table 4 we see that for a similarly conditioned input, we get better precision and a 10-fold increase in performance. However, also here we see that with higher ρmax, precision start to be reduced, even if it is less than for

 

Figure 2: The Jacobi algorithm only preserves the eigenvalues. Here we see the original eigenvector u and the eigenvector after the Givens rotations S|u.

the Jacobi method. Here, to avoid numerical overflow, we must calculate the characteristic polynomials using (8).

Table 3: Results from running the Jacobi eigenvalue algorithm on a 500 × 500 matrix with different ρmax. Tolerance ε = 1.0·10−6. Precision is only 3 leading digits for the lower eigenvalues with, and not usable for th rest. Precision decreases with increasing ρmax. CPU time around 70 seconds for all calculations.

λi
exact
ρmax = 4.0
ρmax = 10.0
ρmax = 20.0
0
3
2.982e+00
2.955e+00
2.910e+00
1
7
6.976e+00
6.932e+00
6.863e+00
2
11
1.104e+01
1.091e+01
1.083e+01
3
15
1.554e+01
1.490e+01
1.479e+01
4
19
2.100e+01
1.888e+01
1.876e+01
5
23
2.767e+01
2.287e+01
2.273e+01
6
27
3.560e+01
2.686e+01
2.670e+01
7
31
4.477e+01
3.085e+01
3.067e+01
8
35
5.518e+01
3.483e+01
3.464e+01
Table 4: Results from running the Polynomial expansion eigenvalue algorithm on a 500 × 500 matrix with various ρmax. Tolerance ε = 1.0 · 10−6. Precision is up to 5 leading digits for the lower eigenvalues and lower ρmax, but, as for Jacobi, the results are not usable for the rest. CPU time about 8 seconds for all calculations.

 

λi              exact       ρmax = 4.0               ρmax = 10.0            ρmax = 20.0            ρmax = 100.0

 

 

4           Conclusion
Our main findings are that the Jacobi eigenvalue algorithm performs well on solving the eigenvalue problem for Toeplitz matrices, since we can accept a large convergence constant ε for the off-diagonal elements, but not so well for the more general case where we must increase n in order to give a proper representation of the problem. The value of the diagonal elements increase by (ih)2, meaning that the addition and subtraction of diagonal elements in (3) is a source for numerical imprecision, which is accumulated due to the high number of rotations needed, which grow by O(n2) [4].

The Polynomial expansion algorithm performs better, both with regards to time and memory use, and since over/underflow problems can be largely mitigated by calculating the characteristic polynomial with qk(x) = (d − x) − e2/qk−1(x), this algorithm is the overall better between the two, and may be the only feasible for problems of size over n = 104 to 105, due to its memory economy.

References
[1]  Fys4150 computational physics. https://www.uio.no/studier/emner/matnat/fys/FYS4150/, Autumn semester 2020.

[2]  W Barth, RS Martin, and JH Wilkinson. Calculation of the eigenvalues of a symmetric tridiagonal matrix by the method of bisection. Numerische Mathematik, 9(5):386–393, 1967.

[3]  Keith               Briggs.   Diagonally            dominant              matrix. https://mathworld.wolfram.com/DiagonallyDominantMatrix.html. Visited on 29/09/20.

[4]  Morten Hjort-Jensen.  Computational physics, lecture notes fall 2015.

https://github.com/CompPhysics/ComputationalPhysics/blob/master/doc/Lectures/lectures2015.pdf, University of Oslo, August 2015.

[5]  Morten            Hjort-Jensen.        Project   2,             computational     physics  fys4150.

http://compphysics.github.io/ComputationalPhysics/doc/Projects/2020/Project2/html/Project2.html, University of Oslo, Autumn semester 2020.

[6]  Tom Lyche.    Lecture  notes      for           mat-inf   4130,      2017. https://www.uio.no/studier/emner/matnat/math/MATINF4130/h17/book2017.pdf, University of Oslo, 2017.