def eig_jacobi(A):
tol=10**(-4)
n = np.shape(A)[0]
P = np.eye(n)
while off(A) > tol:
r, s = offabsmaks(A)
theta = np.arctan( 2*A[r-1,s-1]/(A[s-1,s-1]-A[r-1,r-1]) )/2
Pk = P_k(r, s, theta, n)
A = Pk.T @ A @ Pk
P = P @ Pk
return A, P2 Project: Jacobi’s method for finding eigenvalues
We will study the Jacobi method, which can be used to find all eigenvectors and eigenvalues of a symmetric matrix numerically. It does so by constructing orthogonal matrices P_1,P_2,..., and successively applying similarity transformations \begin{eqnarray*} A_1 &=& A \\ A_2 &=& P_1^TA_1P_1 \\ A_3 &=& P_2^TA_2P_2 \\ \vdots & & \vdots \\ A_{k+1} &=& P_k^TA_kP_k \end{eqnarray*}
Exercise 2.1 Set B=P^TAP with A symmetric and P orthogonal. Show that B is symmetric.
Exercise 2.2 Show that the matrices A_1,A_2,A_3,... all have the same eigenvalues. Show also that, if \mathbf{v} is an eigenvector of A_k, then (P_k)^T\textbf{v} is an eigenvector for A_{k+1} and that, if \mathbf{v} is an eigenvector of A, then (P_1\cdots P_k)^T\textbf{v} is an eigenvector of A_{k+1}.
Thus, if P is a matrix consisting of orthonormal eigenvectors for A (so that A=PDP^T), (P_1\cdots P_k)^TP is a matrix consisting of orthonormal eigenvectors of A_{k+1}. It follows that A_{k+1}=(P_1\cdots P_k)^TPDP^T(P_1\cdots P_k), so that A=PDP^T=(P_1\cdots P_k)A_{k+1}(P_1\cdots P_k)^T \tag{2.1} The idea is now to choose the P_i so that the matrices A_i get closer to zero outside the diagonal. We will later see that we can achieve this.
Exercise 2.3 One can show that the eigenvalues of a matrix can be written as a continuous function of the matrix entries. Assume that the off-diagonal entries of A_{k+1} in Equation 2.1 all converge to 0. Explain that the diagonal entries of A_{k+1} converge to the eigenvalues of A.
We will not go into similar arguments for the convergence of P_1\cdots P_k to a matrix with eigenvectors as columns.
After each transformation A_{k+1} = P_k^TA_kP_k we will therefore compute \text{off}(A) = \sqrt{\sum_{i\neq j}^n a_{ij}^2}, and if this value is less than a given tolerance, A is considered ‘diagonal enough’, and the algorithm terminates.
The matrices P_k will have the following form P_k = \begin{pmatrix} 1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 0 \\ 0 & 1 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cos\theta & \cdots & \sin\theta & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & -\sin\theta & \cdots & \cos\theta & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1 \end{pmatrix} \tag{2.2} where (P_k)_{rr} = (P_k)_{ss} = \cos\theta, (P_k)_{rs} = -(P_k)_{sr} = \sin\theta, (P_k)_{ii} = 1 (for i\neq r,s) and 0 else.
Exercise 2.4 Show that the matrices P_k defined by Equation 2.2 are orthogonal, so that (P_k)^{-1} = (P_k)^T = \begin{pmatrix} 1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 0 \\ 0 & 1 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cos\theta & \cdots & -\sin\theta & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \sin\theta & \cdots & \cos\theta & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & 0 & 1 \end{pmatrix} \tag{2.3}
Exercise 2.5 Why is P_k often called a rotation matrix?
Exercise 2.6 Write a function which takes the position (r,s), the angle \theta, and the dimension n as arguments, and returns the matrix Equation 2.2.
Exercise 2.7 Write a function which takes a real n\times n-matrix A as input, and returns \text{off}(A).
Exercise 2.8 Write a function which takes a real n\times n-matrix as input, and returns the position (r,s) of the element outside the diagonal with biggest absolute value.
Exercise 2.9 The quantity \sqrt{\sum_{i,j} a_{ij}^2} is called the Frobenius norm of A, and is written \|A\|_F. Show that, if P is orthogonal and B=PA, A and B have the same Frobenius norm. We say that the Frobenius norm is invariant under orthogonal transformations.
Exercise 2.10 Show that a_{rr}^2 + a_{ss}^2 - b_{rr}^2 - b_{ss}^2=2b_{rs}^1-2a_{rs}^2. \tag{2.4} by restricting B=P^TAP to row/columns r and s,
If i,j are both different from r and s, we have that b_{ij}=a_{ij}. We also need to compute b_{rs}.
Exercise 2.11 Show that \begin{equation} b_{rs} = b_{sr} = \sin\theta\cos\theta (a_{rr}-a_{ss}) + (\cos^2\theta-\sin^2\theta)a_{rs}. \end{equation}
From the first point above it follows that \sum_{i\neq r,s} a_{ii}^2 = \sum_{i\neq r,s} b_{ii}^2. Adding a_{rr}^2+a_{ss}^2+b_{rr}^2+b_{ss}^2 to both sides gives \sum_{i=1}^n a_{ii}^2 + b_{rr}^2 + b_{ss}^2 = \sum_{i=1}^n b_{ii}^2 + a_{rr}^2 + a_{ss}^2, which can be written as \|A\|_F^2 - \text{off}(A)^2 + b_{rr}^2 + b_{ss}^2 = \|B\|_F^2 - \text{off}(B)^2 + a_{rr}^2 + a_{ss}^2. Since \|A\|_F=\|B\|_F this can be written as \text{off}(B)^2 = \text{off}(A)^2 + a_{rr}^2 + a_{ss}^2 - b_{rr}^2 - b_{ss}^2 = \text{off}(A)^2 + 2b_{rs}^1 - 2a_{rs}^2, \tag{2.5} where we used Equation 2.4. Thus, the off-diagonal contribution \text{off}(B) only changes due to the (r,s) and (s,r) components,
Exercise 2.12 Explain that r,s,\theta will give a B with smallest possible \text{off}(B) when r,s are chosen so that |a_{rs}| becomes biggest possible, and b_{rs}=0, i.e., when \sin\theta\cos\theta (a_{rr}-a_{ss}) + (\cos^2\theta-\sin^2\theta)a_{rs} = 0 \tag{2.6} The new value for \text{off}(B) is then \text{off}(B)^2 = \text{off}(A)^2 - 2a_{rs}^2
Note that @los-denne0 also can be written as \tan(2\theta)=\frac{2a_{rs}}{a_{ss}-a_{rr}}. \tag{2.7}
Exercise 2.13 Show that when r,s,\theta are chosen as above, we have that \text{off}(B)^2 = \left(1- \frac{2}{n(n-1)}\right)\text{off}(A)^2. Explain from this why \text{off}(B) converges to 0.
Exercise 2.14 Study the following code.
Explain all parts of this code, and which exercise over that are used. What will the code print?
To test this we create a random 6\times 6 matrix A. A+A^T is then a symmetric matrix:
A = np.random.rand(6,6)
A = A+A.TOur method gives the following approximation to the eigenvalues:
D, P = eig_jacobi(A)
np.diag(D)array([ 5.97750786, -0.14043933, 0.52028617, -0.54697347, -1.38928858,
1.17552043])We first verify that the columns in P are eigenvectors. This is the same as testing that PDP^T is an orthogonal diagonalisation of A:
A - P @ D @ P.Tarray([[-6.66133815e-16, -5.55111512e-16, -4.99600361e-16,
2.22044605e-16, 4.44089210e-16, 0.00000000e+00],
[-1.11022302e-16, -2.22044605e-16, -2.22044605e-16,
-2.22044605e-16, 4.44089210e-16, 2.22044605e-16],
[-4.44089210e-16, -1.11022302e-16, 0.00000000e+00,
0.00000000e+00, 2.22044605e-16, -2.22044605e-16],
[ 4.44089210e-16, -4.44089210e-16, 2.22044605e-16,
0.00000000e+00, -1.11022302e-16, 0.00000000e+00],
[ 4.44089210e-16, 3.33066907e-16, 4.44089210e-16,
1.11022302e-16, 7.77156117e-16, 2.22044605e-16],
[ 0.00000000e+00, 0.00000000e+00, -2.22044605e-16,
1.11022302e-16, 2.22044605e-16, -4.44089210e-16]])It is seen that A-PDP^T is not exactly zero, but very close to zero. This is due to roundoff errors on the computer. We then test that we get the same eigenvalues as returned by np.linalg.eig:
D, P = np.linalg.eig(A)
Darray([ 5.97750786, -1.38928858, -0.54697347, -0.14043933, 0.52028617,
1.17552043])We see that eig_jacobi and np.linalg.eig give approximately the same values, but they are listed in a different order. On a computer we have no control over the order eigenvalues are listed - the algorithm decides this along the way.
The matrix A may be very large, and the function P_k allocates a matrix of the same size.
Write a function which implements the transformation A\to AP_k and avoids this problem