Matrix exponentiation using gsl library (Computational Quantum Mechanics)

In summary: H*ewIn summary, the author is trying to solve a two dimensional coupled oscillator problem numerically. He found that using the matrix representation of the creation and annihilation operator was the easiest way. He then found the required hamiltonian by multiplying the state ψ with the unitary operator U=exp(-(i/h)HΔt). The author did the problem with MATLAB using expm function. However, he found that getting better computational result he has 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.
  • #1
sijokjoseph
6
0
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 ?
 
Physics news on Phys.org
  • #2
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.
 
  • #3
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
 
  • #4
sijokjoseph said:
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):
Code:
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).
 
  • #5
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
 
  • #6
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);
}
 

1. What is matrix exponentiation and how is it used in computational quantum mechanics?

Matrix exponentiation is a mathematical operation used to raise a square matrix to a given power. In computational quantum mechanics, it is used to simulate the time evolution of quantum systems, where the Hamiltonian of the system is represented by a matrix. By exponentiating this matrix, we can obtain the time evolution operator, which can then be applied to the initial quantum state to calculate its state at any given time.

2. What is the gsl library and how does it help with matrix exponentiation?

The gsl library, or GNU Scientific Library, is a collection of mathematical functions and algorithms commonly used in scientific computing. It includes a function called gsl_matrix_exponential which can be used to exponentiate a matrix efficiently and accurately. This function is especially useful in computational quantum mechanics as it allows for faster and more accurate simulations of quantum systems.

3. How do I use the gsl library for matrix exponentiation in my code?

To use the gsl library for matrix exponentiation, you first need to include the gsl/gsl_matrix.h header file in your code. Then, you can use the gsl_matrix_exponential function to exponentiate a matrix. Make sure to provide the matrix as the input and store the resulting exponentiated matrix in a new variable.

4. Are there any limitations to using the gsl library for matrix exponentiation?

The gsl library is a powerful tool for matrix exponentiation, but it does have some limitations. It can only be used for square matrices and may not be suitable for very large matrices. Additionally, it requires some knowledge of linear algebra to properly use the function and interpret the results.

5. Can I use other libraries or methods for matrix exponentiation in computational quantum mechanics?

Yes, there are other libraries and methods available for matrix exponentiation in computational quantum mechanics. Some popular alternatives include the Eigen library and the Krylov subspace method. It is important to choose the method that best suits your specific needs and to properly understand the limitations and assumptions of each method.

Similar threads

  • Quantum Physics
Replies
1
Views
858
  • Quantum Physics
Replies
2
Views
907
Replies
2
Views
679
  • Quantum Physics
Replies
1
Views
657
Replies
3
Views
939
Replies
41
Views
8K
Replies
4
Views
3K
  • Quantum Physics
Replies
31
Views
3K
Replies
1
Views
980
Replies
2
Views
2K
Back
Top