Exercise "matrix diagonalization"

1. (6 points) Jacobi diagonalization with cyclic sweeps

   1. Implement the Jacobi eigenvalue method for real symmetric matrices using the cyclic sweep algorithm: eliminate off-diagonal elements in a predefined sequence which spans all off-diagonal elements, for example row after row, and repeat the sweeps until converged.

   2. Make some interesting examples proving your implementation works as intended.

   Typically one stores the diagonal elements of the matrix in a separate vector and only updates the strict upper triangular part of the matrix. The user can then always restore their symmetric matrix from the remaining lower triangle.

2. (3 points) Jacobi diagonalization eigenvalue-by-eigenvalue

   1. Implement the following variant of Jacobi diagonalization algorithm:

      1. Eliminate the off-diagonal elements in the first row only (by running the sweeps not over the whole matrix, but only over the first row until the off-diagonal elements of the first row are all zeroed). Argue, that the corresponding diagonal element is the largest (or lowest) eigenvalue. Find out how to change the formulas for the rotation angle to obtain the lowest (largest) eigenvalue.
      2. Eliminate the off-diagonal elements in the second row. Argue, that the eliminated elements in the first row will not be affected. Argue that the corresponding diagonal element is the second largest (lowest) eigenvalue.
      3. Continue elimination of the subsequent rows until the whole matrix is diagonalized.

      Compare the number of Jacobi rotations for the cyclic strategy and for the consecutive-eigenvalues strategy.

   2. Now, instead of eliminating the off-diagonal elements of a row consecutively one after another, do find the largest off-diagonal element of the row and eliminate it until all off-diagonal elements of the row are eliminated.

      Find out whether this gives any improvements in i) the number of rotations, ii) the run-time.

3. (1 point)

   Implement the following strategy in the Jacobi diagonalization method:

   1. Eliminate the upper right quadrand of the matrix. You obtain the middle eigenvalue and two decoupled matrices of half-size: the upper left quadrand and the lower right quadrand.
   2. Apply the algorithm reqursively to the two half-size matrices until you get 1-by-1 matrix which is already diagonal.

   Find out whether this algoritm gives any speed improvement compared to your other algorithms.

4. (0 point) Comparison with library routines

   Compare the speed of your implementations with library implementations (e.g. GSL, Armadillo, and NumPy).