# Improve speed of Nested Matrix exponential

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. ?

• Welcome! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Jul 25 '16 at 13:08
• Let me simplify your construction of L a bit: n = 20; hamiltonian = SparseArray[{Band[{1, 2}] -> -1, Band[{2, 1}] -> -1, {1, n} -> -1, {n, 1} -> -1}, {n, n}]; A = SparseArray[MapIndexed[MapAt[Prepend[First[#2]], #, 1] &, Most[ArrayRules[Sqrt[Abs[hamiltonian]]]]], Prepend[Dimensions[hamiltonian], 2 n]]; id = IdentityMatrix[n, SparseArray]; v = SparseArray[{i_, i_, i_} -> 1, {n, n, n}]; – J. M.'s discontentment Jul 26 '16 at 10:32
• and then L[α_] := L[α] = α Sum[KroneckerProduct[Conjugate[L], L] - (KroneckerProduct[id, ConjugateTranspose[L].L] + KroneckerProduct[Transpose[L].Conjugate[L], id])/2, {L, A}] - (1 - α) I (KroneckerProduct[id, hamiltonian] - KroneckerProduct[Transpose[hamiltonian], id]) – J. M.'s discontentment Jul 26 '16 at 10:32
• @J.M. thank you for your modification :) it looks very nice ! but there is an argument missing in the first Prepend command to define A. Do you think there is a possibility to improve the speed of the nested list including the MatrixExponential with this definition of L? – Max Jul 27 '16 at 9:03
• @Max, it should work in the current version; what version are you using? – J. M.'s discontentment Jul 27 '16 at 10:06

• 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.

• The reason for using the action form of MatrixExp[] is that it can be prohibitive to store the entire matrix if one only wants the resulting vector for his application. Yet again the trade-off between speed and storage… – J. M.'s discontentment Jul 25 '16 at 13:31
• @J.M. Yes, but matExp has the same size with L[r]. Since he has stored L[1] and L[2], I think storing matExp is acceptable. – wuyingddg Jul 25 '16 at 14:10
• okay thats extremely nice :) Is it further useful to distribute these calculations to several kernels with parallelmap or paralleltable ( I think paralleltable would be the more appropriate choice isnt it? or is it useful to use both?). Also it would be interesting to know how much memory has to be available to distribute these calculations to 2, 4 or 16 kernels, do you maybe also know this? :)or can you tell me how to calculate this on my own?:) – Max Jul 25 '16 at 19:12
• @Max I'm afraid not. You can just replace Table into ParallelTable and will found the code become slower. That is because your j is too small. If the range in Table is {j, 1, 5} or larger, ParallelTable can have higher effiency. – wuyingddg Jul 26 '16 at 1:23
• ah nice :) yes, in the end I will have {j,1,10}. Further I tried some examples and maybe it doesnt become clear but the code for L was just an example... The expression I use gives me a L which is a SparseArray. Further I tried to calculate MatrixExp[L] and found that it is extremely slow for Sparsearrays ( it becomes faster by using Normal but for large n it is nevertheless slower then MatrixExp[L,v]`. I will edit my question with the right definition for L. Maybe it was little unclear, sorry for that. – Max Jul 26 '16 at 9:41