I have a matrix of the form
Jy = [[ 0 -im*N1 0 0 ....... ] [im*N1 0 -im*N2 0 ..... ] [0 im*N2 0 -im*N3 ....] ........................................... ]
It is tridiagonal, with identically zero diagonal. It is also hermitian as I tried to draw here. It’s size is
2j+1 (this is angular momentum operator Jy in Jz basis in quantum mechanics).
Its eigenvalues are also known to me:
j, j-1, j-2,.....0, -1,.....-j. Using the fact that its tridiagonal, I was able to write a routine to find its eigenvectors which is as follows:
function renormalize!(o, e, n, d) for k ∈ 1:2:n o[k]/=d; e[k+1]/=d end end function eigvecJy(v, λ) len = length(v) o = zeros(len+1); e = zeros(len+1) #odd and even η = √(1+λ^2/v^2) o=1/η; e = λ/(v*η); for n ∈ 3:2:len-1 o[n] = (-e[n-1]*λ + o[n-2]*v[n-2])/v[n-1] renormalize!(o, e, n, √(1 + o[n]^2)) e[n+1] = (o[n]*λ + e[n-1]*v[n-1])/v[n] renormalize!(o, e, n, √(1 + e[n+1]^2)) end o[end] = (-e[len]*λ + o[len-1]*v[len-1])/v[len] d = √(1+o[end]^2); renormalize!(o, e, len-1, d); o[end]/=d; return o, e #eigenvector=o+i e end
which is slower than
eigen in extracting all the eigenvectors.
In fact, in my custom method I only calculate
j+1 vectors and others are obtained by a symmetry. Still eigen is faster (and it also gives eigenvalues).
I am very surprised to see this. I don’t understand why … maybe my function isn’t efficiently written. I guess
eigen is from
LAPACK library, but still I have used so many physical inputs. Is this what I must accept?
It is important to point out why I did this: I need the rotation matrix
At first I tried using
Jy is hermitian, so as per Julia docs its exponential is calculated using
eigen. Or so I thought. But
exp clearly does not use
eigen because its WAYY slower than my function, and furthermore, it gives NaN all over the place for
j>1000. (I am using Julia 1.4.1.) Even
scipy was faster than that - so I sticked with it till I found it lacking.
That lead me to write my own function. But when I curiously compared it with
eigen directly, I found it faster than my custom implementation. If I knew eigen was so damn fast I wouldn’t even have bothered spending time on this. But now that I have, I can’t accept it!