Matrix exponentiation using gsl library (Computational Quantum Mechanics)

[sijokjoseph]

Dear Friends,
I have to solve a two dimensional coupled oscillator problem numerically.
Hence the easiest way is to use the matrix representation of the creation and annihilation operator and then find the required hamiltonian. To find the time evolution of the state I have to multiply the state ψ with the unitary operator U=exp(-(i/h)HΔt). Where H is the coupled nonlinear harmonic oscillator hamiltonian.
I did the problem with MATLAB using expm function. Getting better computational result I have to choose the size of the matrix H as 13225 by 13225. In MATLAB doing this integration for 1000 points will take more than 24 hours. Hence this is not a nice option.
Hence I started programming using C gsl_libraries. Now the problem comes, how will I exponentiate a gsl_complex matrix ? that is the first difficulty. The second one is how will assign complex numbers to a gsl_matrix_complex. I am trying to diagonalize the hamiltonian then find the eigen values (D) and eigen vectors (V)then calculate V*Diag(i*dt*λ)*inv(V) . which will give me the expm result. But I don't know how to multiply a gsl_matrix with a complex number, do you have any ideas, regarding this problem ?

 

[cgk]

I don't know about gsl, but there might be simpler solutions to your problem.

Diagonalizing large matrices takes time. On a quad-core i7 I'd estimate something on the order of 5--10 minutes for a 13500^2 matrix. That cannot be avoided. Of course, if you have a time-independent Hamiltonian, then as you correctly note, you can determine all evolution points from the single diagonalization; however, you could do the very same thing in Matlab or scipy, too (i.e., using expm if you need multiple expm(t *M) with the same M and different t is very inefficient). Note also that expm(2*M) = expm(M) * expm(M), it is thus sufficient to calculate the expm a single time to get all your time points: Just keep on powering the matrix you get the first time for the time step. You would thus only need one diagonalization or expm, and 999 matrix-matrix multiplications.

Note also that in order to propagate a wave function, it is not actually required to diagonalize the Hamiltonian. It is enough to be able to calculate contractions of the Hamiltonian with the current wave function, because this will give you the time derivative of the wave function at the current point in time. You could thus simply integrate the wave function as an ordinary differential equation in time, without ever doing anything except matrix * vector multiplications, so only O(N^2) operations (unlike above, which are all O(N^3)). This would allow for going for much larger matrix dimensions.

 

[sijokjoseph]

Dear cgk,
Thank you for the response, It is enough to be able to calculate contractions of the Hamiltonian with the current wave function, because this will give you the time derivative of the wave function at the current point in time. You could thus simply integrate the wave function as an ordinary differential equation in time, without ever doing anything except matrix * vector multiplications, so only O(N^2) operations (unlike above, which are all O(N^3)). This would allow for going for much larger matrix dimensions.
What you said is right, but exponentiating a matrix is doing the same with higher accuracy. For example, exp(-i/h*HΔt)=1-i/h*HΔt+(h*HΔt)^2/2+...
What you told is truncating the last terms. But MATLAB uses pade approximation in expm function.
My actual problem is how to multiply a gsl_matrix H with a complex number i ?
Thanking you cgk
Sijo K Joseph

 

[cgk]

My actual problem is how to multiply a gsl_matrix H with a complex number i ?

This is where we do not agree. I think your main problem[1] is that you are using a bad algorithm for propagating your wave function (that is, a bad factorization). If you do it in a slightly different way, then you can easily adjust your working Matlab code to do the calculation in minutes, for hundreds of thousands of time points.

Consider the folloing approach (written in python):

from scipy import *
from scipy.linalg import *

def PropagateWf(v0, dt, ew,ev):
   """return normalize(exp(-i dt H) |v0>).
      Input: ew,ev are the eigenvalues (ew) and eigenvectors (ev) of H."""
   # transform v into eigenbasis of H [O(N^2)] and cast to complex
   vt = (1.+0.j) * dot(transpose(conj(ev)), v0)
   # apply exp(-i dt H) in diagonal H basis [O(N^1)]
   for i in range(vt.shape[0]):
      vt[i] = exp(-1j*dt*ew[i]) * vt[i]
   # transform back to original basis [O(N^2)). Normalize is for complex
   # arguments dt, where exp(-i dt H) is not unitary.
   return dot(ev, vt/norm(vt))

# make an example: set H as random symmetric matrix and
# initial wave function v0 as random vector.
from numpy.random import random
N = 10
H = random((N,N))
H = H + H.T
v0 = random(N)
v0 /= norm(v0)

# compute eigenvalues (ew) and eigenvectors (ev) of H [O(N^3)]
ew,ev = eigh(H)

# propagate in real time -> energy should stay the same
vt = PropagateWf(v0, +10., ew,ev)
print "\nreal-time propagation:"
print "<v0|H|v0> (initial):        %12.6f" % dot(conj(v0), dot(H,v0))
print "<vt|H|vt> (after dt=+10.):  %12.6f" % dot(conj(vt), dot(H,vt)).real

# propagate in imaginary time -> should go to ground state
vt = PropagateWf(v0, -50.j, ew,ev)
print "\nimaginary time propagation:"
print "<v0|H|v0> (initial):        %12.6f" % dot(conj(v0), dot(H,v0))
print "<vt|H|vt> (after dt=-50j):  %12.6f" % dot(conj(vt), dot(H,vt)).real
print "Ground state energy:        %12.6f" % ew[0]

Here you do a single O(N^3) operation[2] (the spectral decomposition of H in ew,ev=eigh(H)) and then can propagate any wave function to any point in time in O(N^2), via the PropagateWf function. Note that this function does *not* calculate expm(-i dt H), but rather the contraction of expm(-i dt H) |v0> with a given trial vector, and doing the latter is orders of magnitude cheaper and just as exact.

I would recommend to think a second about this point before going on with the gsl problem.


[1] Of course you should also learn how to solve the complex number problem gsl if you apply it in practice... but that is not a physics problem. Note that using C/gsl, python/scipy, or Matlab will make next to no difference in this particular case; in fact Matlab will probably be the fastest because it calls a good actual BLAS like MKL or ACML to do the spectral decomposition, followed by scipy (calls a bad BLAS), and then gsl (if it uses its own spectral decomposition routine instead of calling BLAS).

[2] Another side note: This O(N^3) operation can also be avoided by integrating the system as an ordinary differential equation in time as I mentioned. This does not imply an O(dt) error in the time step--you can also use a better integrator and get a smaller error (e.g., the 4:4 Runge Kutta is simple to implement and will yield an O(dt^4) error).

 

[sijokjoseph]

Dear cgk,
Thank you very much for your effort. Even though, I don't know python, it really
helped me. You are right, it is more wiser to implement faster algorithms in MATLAB itself.
Again thanking you.
Sijo K Joseph

 

[sijokjoseph]

GSL Complex Matrix Expoential C code

By the way, I have solved my problem by writing this code in gsl. This is ok, but the lack of sparse matrices makes the program slower. From my personal experience, It is better to start with ARPACK to have the fast quantum codes in C.#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
void my_gsl_complex_matrix_exponential(gsl_matrix_complex *eA, gsl_matrix_complex *A, int dimx)
{
int j,k=0;
gsl_complex temp;
gsl_matrix *matreal =gsl_matrix_alloc(2*dimx,2*dimx);
gsl_matrix *expmatreal =gsl_matrix_alloc(2*dimx,2*dimx);
//Converting the complex matrix into real one using A=[Areal, Aimag;-Aimag,Areal]
for (j = 0; j < dimx;j++)
for (k = 0; k < dimx;k++)
{
temp=gsl_matrix_complex_get(A,j,k);
gsl_matrix_set(matreal,j,k,GSL_REAL(temp));
gsl_matrix_set(matreal,dimx+j,dimx+k,GSL_REAL(temp));
gsl_matrix_set(matreal,j,dimx+k,GSL_IMAG(temp));
gsl_matrix_set(matreal,dimx+j,k,-GSL_IMAG(temp));
}

gsl_linalg_exponential_ss(matreal,expmatreal,.01);

double realp;
double imagp;
for (j = 0; j < dimx;j++)
for (k = 0; k < dimx;k++)
{
realp=gsl_matrix_get(expmatreal,j,k);
imagp=gsl_matrix_get(expmatreal,j,dimx+k);
gsl_matrix_complex_set(eA,j,k,gsl_complex_rect(realp,imagp));
}
gsl_matrix_free(matreal);
gsl_matrix_free(expmatreal);
}

int main()
{
int dimx=4;
int i, j ;
gsl_matrix_complex *A = gsl_matrix_complex_alloc (dimx, dimx);
gsl_matrix_complex *eA = gsl_matrix_complex_alloc (dimx, dimx);

for (i = 0; i < dimx;i++)
{
for (j = 0; j < dimx;j++)
{
gsl_matrix_complex_set(A,i,j,gsl_complex_rect(i+j,i-j));
if ((i-j)>=0)
printf("%d+%di ",i+j,i-j);
else
printf("%d%di ",i+j,i-j);

}
printf(";\n");
}
my_gsl_complex_matrix_exponential(eA,A,dimx);
printf("\n Printing the complex matrix exponential\n");
gsl_complex compnum;
for (i = 0; i < dimx;i++)
{
for (j = 0; j < dimx;j++)
{
compnum=gsl_matrix_complex_get(eA,i,j);
if (GSL_IMAG(compnum)>=0)
printf("%f+%fi\t ",GSL_REAL(compnum),GSL_IMAG(compnum));
else
printf("%f%fi\t ",GSL_REAL(compnum),GSL_IMAG(compnum));
}
printf("\n");
}
return(0);
}