Improve speed of Nested Matrix exponential

Question

I have a short Code example and I want to know whether it is possible to improve the speed. The L matrix in my code is a SparseMatrix, but for this question I just use something similiar.

n = 20;
m = 5;

L[r_] := L[r] = 
  Table[RandomComplex[{-1 + -r*I, 1 + r*I}], {i, 1, n^2}, {j, 1, n^2}]

v = Table[SparseArray[{i, i} -> 1, {n, n}], {i, 1, n}];

timestep = 1;

t1 = AbsoluteTime[];
matrix = Table[
   NestList[
    Transpose@
      Partition[MatrixExp[timestep*L[j], Flatten[Transpose@#1]], n] &,
     v[[i]], m], {j, 1, 2}, {i, 1, 2}];
AbsoluteTime[] - t1

This should run with large n and m, and i should be equal to n in the end (the i in the Table for matrix).

I also tried to use ParallelTable, which is quite faster but sometimes it stops calculating... and I dont understand how to know when it happens..maybe it depends on the size of the memory that is available?

I would be thankful for any help!

EDIT: Correct definition for L:
Prelude:

n = 20;
m = 5;
timestep = 1;
startingconditions = 1;
hamiltonian = 
  DiagonalMatrix[Table[-1, {i, 1, n - 1}], 1] + 
    DiagonalMatrix[Table[-1, {i, 1, n - 1}], -1] + 
    SparseArray[{{1, n} -> -1, {n, 1} -> -1}] // SparseArray;
A = SparseArray[{#}, Dimensions[hamiltonian]] & /@ 
   Most[ArrayRules[Sqrt[Abs[hamiltonian]]]];
id = IdentityMatrix[n];
v = Table[SparseArray[{i, i} -> 1, {n, n}], {i, 1, n}];

Clear[L]
L[\[Alpha]_] := 
  L[\[Alpha]] = -(1 - \[Alpha]) I (KroneckerProduct[id, hamiltonian] -
        KroneckerProduct[Transpose[hamiltonian], id]) + \[Alpha] Sum[
      KroneckerProduct[Conjugate[L], 
        L] - (1/2) (KroneckerProduct[id, ConjugateTranspose[L].L] + 
          KroneckerProduct[Transpose[L].Conjugate[L], id]), {L, A}];

This gives me the coefficientmatrix for a system of ODE. The system of ODE's I want to solve by using something like this but, if possible, with improved speed (similiar formula as above)

t1 = AbsoluteTime[];
matrix = Table[
   NestList[
    Transpose@Partition[MatrixExp[L[j], Flatten[Transpose@#1]], n] &, 
    v[[i]], m]
   , {j, 0.1, 1, 0.1}, {i, 1, n}];
AbsoluteTime[] - t1

Maybe it would better to delete the upper part? but you gave a nice answer to it so I will leave it there okay? Maybe the new definition of L will change something. In the end, m is about 200, n as high as possible, and timestep something between 1 and 5. I hope you can optimize something:)

Edit2:
I've played little bit with your comments and my program looks now like this

n = 20;
m = 30;
timestep = 1;
hamiltonian = 
  DiagonalMatrix[Table[-1, {i, 1, n - 1}], 1] + 
    DiagonalMatrix[Table[-1, {i, 1, n - 1}], -1] + 
    SparseArray[{{1, n} -> -1, {n, 1} -> -1}] // SparseArray;
A = SparseArray[{#}, Dimensions[hamiltonian]] & /@ 
   Most[ArrayRules[Sqrt[Abs[hamiltonian]]]];
id = IdentityMatrix[n];
v = Table[SparseArray[{i, i} -> 1, {n, n}], {i, 1, n}];
rev = Table[SparseArray[{i} -> 1, {n^2}], {i, 1, n^2}];
Clear[L]
L[\[Alpha]_] := 
  L[\[Alpha]] = -(1 - \[Alpha]) I (KroneckerProduct[id, hamiltonian] -
        KroneckerProduct[Transpose[hamiltonian], id]) + \[Alpha] Sum[
      KroneckerProduct[Conjugate[L], 
        L] - (1/2) (KroneckerProduct[id, ConjugateTranspose[L].L] + 
          KroneckerProduct[Transpose[L].Conjugate[L], id]), {L, A}];

and now

 t1 = AbsoluteTime[];
    vReshape = Flatten[Transpose@#1] & /@ v;
    matrix2 = 
      Transpose[
       Map[Partition[#, n]\[Transpose] &, 
        ParallelTable[
         matExp = 
          Table[MatrixExp[timestep*L[j], rev[[i]]], {i, 1, 
             n^2}]\[Transpose];
         NestList[#.matExp &, vReshape, m], {j, 0.1, 1, 0.1}], {3}], {1,
         3, 2}];
    AbsoluteTime[] - t1

this now has the advantage that in the NestList command the MatrixExponential has not to be calculated so many times and with

In[54]:= t1 = AbsoluteTime[];
MatrixExp[timestep*L[0.1], rev[[1]]];
AbsoluteTime[] - t1

Out[56]= 0.013133

that is quite faster than

In[57]:= t1 = AbsoluteTime[];
MatrixExp[timestep*L[0.1]];
AbsoluteTime[] - t1

Out[59]= 1.504880

the code becomes about a factor 1.5 faster for n=150, m=200, and {j,0.1,1,0.1}. But this is only true as long m is larger than n as far as I realised. But now there a still some questions open:
1) Is it possible to further increase the speed of this code for large n ( n> 150) and m( m about 200 to 400)
3) Is there a possibility to achieve more free memory i.e. by compression or something? The underlying target is to choose n as high as possible
4) Also I want to use all Kernels, so Parallelmap or ParallelTable should be included
5) Is there a further possibility to gain speed by programming this in another language like C or Fortran etc. ?

thanks in advance

Answer 1

• Use RandomComplex[{-1 + -r*I, 1 + r*I}, {n^2, n^2}] instead of
  Table, generate a random n*n matrix is much faster than generate a random number n^2 times.
• MatrixExp[mat, v] equal to MatrixExp[mat].v, and MatrixExp[mat] will be use for many times. So you can store it to avoid repeating calculation. (because of the floating error, the result of MatrixExp[mat].v may different from MatrixExp[mat, v], you can decide which to use )

With these two item the code can be written like this, which is about four times faster than the code you gave.

In[219]:= n = 20;
m = 5;
Clear[L];
L[r_] := L[r] = RandomComplex[{-1 + -r*I, 1 + r*I}, {n^2, n^2}];

v = Table[SparseArray[{i, i} -> 1, {n, n}], {i, 1, n}];

timestep = 1;

t1 = AbsoluteTime[];
matrix = Table[
   Table[
    matExp = MatrixExp[timestep*L[j]]; 
    NestList[Transpose@Partition[matExp.Flatten[Transpose@#1], n] &, 
     v[[i]], m], {i, 1, 2}]
   , {j, 1, 2}];
AbsoluteTime[] - t1

Out[227]= 0.5463890

• It can be prove that matExp.Flatten[Transpose@#1] & /@ v and (Flatten[Transpose@#1] & /@ v).matExp\[Transpose] are always same, so you can use matrix dot matrix instead calculating matrix dot vector for several times. That can significant speed up your code.

`

In[228]:= t1 = AbsoluteTime[];
vReshape = Flatten[Transpose@#1] & /@ v[[1 ;; 2]];
matrix2 =
  Transpose[Map[Partition[#, n]\[Transpose] &,
    Table[
     matExp = MatrixExp[timestep*L[j]]\[Transpose];
     NestList[#.matExp &, vReshape, m],
     {j, 1, 2}], {3}], {1, 3, 2}];
AbsoluteTime[] - t1

Out[231]= 0.2831995

This code is a bit hard to understand, but it really has a high effiency.