Based on Numerical Recipes in C++, Sec 11.1

So you want to diagonalize a matrix, do you? Well, if you have a tiny symmetric matrix, you REALLY want to write up the algorithm by hand, and don’t want to spend much time trying to understand the algorithm, then you have come to the right place.

Otherwise, use LAPACK/BLAS to call a highly optimized routine that can work extremely quickly on large matrices. Julia has those libraries built in already. Even if you do call those matrices, you can make them work better by understanding what’s going on underneath the hood, which is why we are going through this now.

Start with a base rotation matrix of the Form $% $

From our starting arbitrary symmetric $A$,

we will run a series of transformations,

where each iteration brings $A$ closer to diagonal form. Thus in our implementing our algorithm, we need to determine two things

• The values of $c$ and $s$
• The pattern of sweeping $p$ and $q$

And in the end we will need to finally determine if this actually converges, and if has any sort of efficiency.

So lets expand one transformation, and we if we can solve for $c$ and $s$.

Determining $s$ and $c$

Given we specifically want $a^{\prime}_{pq}$ to be zero, we re-arrange the last equation,

At first glance, this equation might not look easier to solve for $s$ or $c$. At second glance either. We define a new parameter $t = s/c$, which now makes the equation,

now quite easily solvable by our friendly quadratic formula. Though the book does recommend using a form that pulls out smaller roots through

Then reverse solve back to

Though we could use the expressions above, if we simplify them with our new expressions for $c$ and $s$ analytically, we reduce computational load and round off error. These new expressions are

with the new variable

Convergence

The sum of the squares of the off diagonal elements, choosen in either upper or lower triangles arbitrarily,

Eigenvectors

By forming a product of every rotation matrix, we also come to approximate the matrix $V$ where

and $D$ is the diagonal form of $A$. $V$ is computed through iterative computation

Enough with the talking! LETS COMPUTE STUFF

This makes a matrix easier to look at than when its filled with 1.043848974e-12 everywhere

Now on to the rotations!

We don’t always want to compute the eigenvectors, so those are in the optional entries slot. Both tell the function to compute the vectors with computeV=true and input the V=V after the semicolon.

This function performs one sweep

Just creating some size of matrix

This output is after several sweeps

    5×5 Array{Float64,2}:
2.75547   0.0     0.0      0.0        0.0
0.0      -1.0683  0.0      0.0        0.0
0.0       0.0     1.20611  0.0        0.0
0.0       0.0     0.0      0.734884   0.0
0.0       0.0     0.0      0.0       -1.48581

5×5 Array{Float64,2}:
2.75547       1.03898e-31  6.65473e-46   -1.20864e-52  -4.52913e-60
1.03898e-31  -1.0683       4.25018e-46   -2.9053e-69    1.37672e-80
6.65473e-46   4.25018e-46  1.20611        1.91399e-75   5.09273e-116
-1.20864e-52  -2.9053e-69   1.91399e-75    0.734884      0.0
-4.52913e-60   1.37672e-80  5.09273e-116   0.0          -1.48581

5×5 Array{Float64,2}:
0.702867    0.340672  -0.613045  0.0266825  -0.115695
-0.0750998   0.373435   0.239379  0.773343   -0.446705
0.597584   -0.208098   0.484961  0.308218    0.519041
-0.0642891  -0.754388  -0.459432  0.458861   -0.0716542
-0.37296     0.363431  -0.347289  0.309317    0.715913

1.0794844600380612e-62


Compare the Optimized LAPLACK routine to your results

    Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
eigenvalues:
5-element Array{Float64,1}:
-1.4858101513857376
-1.068298321190819
0.7348837698415409
1.206106933351782
2.755472216052164
eigenvectors:
5×5 Array{Float64,2}:
-0.115695   -0.340672  -0.0266825   0.613045   0.702867
-0.446705   -0.373435  -0.773343   -0.239379  -0.0750998
0.519041    0.208098  -0.308218   -0.484961   0.597584
-0.0716542   0.754388  -0.458861    0.459432  -0.0642891
0.715913   -0.363431  -0.309317    0.347289  -0.37296

      5x5 Array{Float64,2}:
0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0


How long does it take to make a Sweep? How much memory will the computation take? This is dependent on how large the matrix is, and determines whether or not we want to use this algorithm.

    0.000016 seconds (230 allocations: 32.594 KiB)
0.000069 seconds (955 allocations: 196.188 KiB)
0.008129 seconds (24.75 k allocations: 17.372 MiB, 27.48% gc time)


In addition to time per sweep, we need to know how many sweeps we need to run. So again we run it on a 10x10, 20x20, and 100x100. The efficiency of the algorithm would get a lot worse if we have to sweep the 100x100 a bunch of times.

    7×4 Array{Float64,2}:
1.0  1.74923       2.64638      14.7488
2.0  0.0945499     0.422473      2.80609
3.0  0.000314227   0.0162891     0.399226
4.0  6.31792e-10   1.09268e-5    0.0356924
5.0  3.62048e-22   1.18607e-11   0.000598666
6.0  9.38425e-48   7.94096e-26   3.28477e-7
7.0  1.14895e-112  1.23362e-55   6.11775e-13


Well, so we’ve seen how to do one form of exact diagonalization that works, but doesn’t scale very well up to 100x100 matrices. So stay tuned for the Householder method, hopefully coming up soon.

Until then, happy computing :)