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 \begin{equation} P_{pq} = \begin{pmatrix} 1& & & && && 0\\ & \ddots &&&& & \\ && c & \cdots & s && \\ &&\vdots& 1 & \vdots &&\\ && -s & \cdots & c && \\ &&&&& \ddots & \\ 0&&& & && 1 \end{pmatrix} \end{equation}
From our starting arbitrary symmetric A, \begin{equation} A^{T} = A \end{equation} we will run a series of transformations, \begin{equation} A^{\prime}= P^{T}_{pq} \cdot A \cdot P_{pq} \end{equation} where each iteration brings A closer to diagonal form. Thus in implementing our algorithm, we need to determine two things
So lets expand one transformation, and we if we can solve for $c$ and $s$.
\begin{align} a^{\prime}_{rp} & = c a_{rp} - s a_{rq} \\ a^{\prime}_{rq} & = c a_{rq} + s a_{rp} \\ a^{\prime}_{pp} & = c^2 a_{pp} + s^2 a_{qq} -2 sc a_{pq} \\ a^{\prime}_{qq} & = s^2 a_{qq} + c^2 a_{qq} + 2sc a_{pq} \\ a^{\prime}_{pq} & = \left( c^2-s^2 \right) a_{pq} + sc \left(a_{pq} - a_{qq} \right) \end{align}Given we specifically want $a^{\prime}_{pq}$ to be zero, we re-arrange the last equation, \begin{equation} \frac{c^2-s^2}{2 sc} = \frac{a_{pq}-a_{qq}}{2 a_{pq}} =\theta \end{equation} At first glance, this equation might not look easier to solve for $s$ or $c$. Second either. We define a new parameter $t = s/c$, which now makes the equation, \begin{equation} \frac{1-t^2}{2 t} = \theta \;\;\;\; \implies \;\;\; t^2 -2 \theta t -1=0, \end{equation} now quite easily solvable by our friendly quadratic formula. Though the book does recommend using a form that pulls out smaller roots through \begin{equation} t=\frac{\text{sgn}(\theta)}{|\theta| + \sqrt{\theta^2 + 1} }. \end{equation} Then reverse solve back to \begin{align} c&=\frac{1}{\sqrt{t^2+1}}\\ s&=tc \end{align}
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 \begin{align} a^{\prime}_{pq} & = 0\\ a^{\prime}_{qq} & = a_{qq} + t a_{qp} \\ a^{\prime}_{pp} &= a_{pp} - t a_{pq} \\ a^{\prime}_{rp} &= a_{rp} - s \left( a_{rq} +\tau a_{rp} \right) \\ a^{\prime}_{rq} &= a_{rq} + s \left( a_{rp} -\tau a_{rq} \right) \end{align} with the new variable \begin{equation} \tau = \frac{s}{1+c} \end{equation}
The sum of the squares of the off diagonal elements (choosen in either upper or lower triangles arbitrarily) \begin{equation} S=\sum\limits_{r < s} |a_{rs}|^2 \end{equation}
By forming a product of every rotation matrix, we also come to approximate the matrix $V$ where \begin{equation} D = V^{T} \cdot A \cdot V \end{equation} and $D$ is the diagonal form of $A$. $V$ is computed through iterative computation \begin{align} V^{\prime} & = V \cdot P_i \\ v^{\prime}_{rs} &= v_{rs} \\ v^{\prime}_{rp} &= c v_{rp} - s v_{rq} \\ v^{\prime}_{rq} &= s v_{rp} + c v_{rq} \end{align}
In [1]:
using LinearAlgebra
In [2]:
# First, Lets make our nice, helpful functions
## A function to look at the convergence
function convergence(A::Array)
num=0.0
l=size(A)[1]
for ii in 1:(l-1)
for jj in (ii+1):l ## just looking at the lower triangle
num+=A[ii,jj]^2
#println(ii,' ',jj,' ',num,' ',A[ii,jj])
end
end
return num
end
Out[2]:
In [3]:
# This makes a matrix easier to look at when its filled
# with 1.043848974e-12 everywhere
function roundmatrix(A::Array,rtol::Real)
Ap=copy(A)
for ii in 1:length(A)
if abs(Ap[ii])<rtol
Ap[ii]=0
end
end
return Ap;
end
Out[3]:
In [4]:
## Here we create a random symmetric matrix
function makeA(n)
A=randn(n,n);
for ii in 1:n
A[ii,1:ii]=transpose(A[1:ii,ii])
end
V=Matrix{Float64}(I,n,n) #initializing the orthogonal transformation
return A,copy(A),V
end
## One A returned will be stored to compare initial and final
Out[4]:
In [5]:
#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.
function Rotate(A::Array,p::Int,q::Int; computeV=false, V::Array=Matrix{Float64}(I,1,1) )
θ=(A[q,q]-A[p,p])/(2*A[p,q]);
t=sign(θ)/(abs(θ)+sqrt(θ^2+1));
c=1/sqrt(t^2+1)
s=t*c
τ=s/(1+c)
l=size(A)[1]
Ap=copy(A[:,p])
Aq=copy(A[:,q])
for r in 1:l
A[r,p]=Ap[r]-s*(Aq[r]+τ*Ap[r])
A[r,q]=Aq[r]+s*(Ap[r]-τ*Aq[r])
A[p,r]=A[r,p]
A[q,r]=A[r,q]
end
A[p,q]=0
A[q,p]=0
A[p,p]=Ap[p]-t*Aq[p]
A[q,q]=Aq[q]+t*Aq[p]
if computeV==true
Vp=copy(V[:,p])
Vq=copy(V[:,q])
for r in 1:l
V[r,p]=c*Vp[r]-s*Vq[r]
V[r,q]=s*Vp[r]+c*Vq[r]
end
return A,V
else
return A;
end
end
Out[5]:
In [6]:
# This function performs one sweep
function Sweep(A;compV=false,V=Matrix{Float64}(I,1,1))
n=size(A)[1]
for ii in 2:n
for jj in 1:(ii-1) ## Just over one triangle
if compV==false
A=Rotate(A,ii,jj)
else
A,V=Rotate(A,ii,jj;computeV=true,V=V);
end
end
end
if compV==false
return A
else
return A,V
end
end
Out[6]:
In [7]:
A,A0,V=makeA(5);
In [8]:
## keep evaluating for a couple iterations
## watch how it changes
A,V=Sweep(A;compV=true,V=V);
display(roundmatrix(A,1e-10))
display(A)
display(V)
display(convergence(A))
In [9]:
## Compare the Optimized LAPLACK routine to your results
eigen(A0)
Out[9]:
In [10]:
## A good check to make sure V is an orthonomal transformation
roundmatrix(V*A*transpose(V)-A0,1e-12)
Out[10]:
In [11]:
# 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
A,A0,V=makeA(10);
@time Sweep(A);
A,A0,V=makeA(20);
@time Sweep(A);
A,A0,V=makeA(100);
@time Sweep(A);
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.
In [12]:
A10,Ap10,V=makeA(10);
A20,Ap20,V=makeA(20);
A100,Ap100,V=makeA(100);
nsweep=collect(1:7);
conv10=zeros(7)
conv20=zeros(7)
conv100=zeros(7)
for i in nsweep
A10=Sweep(A10)
A20=Sweep(A20)
A100=Sweep(A100)
conv10[i]=convergence(A10)
conv20[i]=convergence(A20)
conv100[i]=convergence(A100)
end
[nsweep conv10/10 conv20/20 conv100/100]
Out[12]:
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 :)
In [ ]: