# Improve simulation speed using Compile

I have two square matrices, a hermitian matrix h and an unitary matrix v. I am interested in the following function

step[h_,v_,n_]:=Map[Tr,Dot[NestList[MatrixExp[-I h].#.MatrixExp[I h]&,v,n-1],v]];


which generates a list of n real numbers. The dimensions of the matrices are not so big so MatrixExp is still quite efficient. I want to make the above function faster by using Compile,

stepcp:=Compile[{{m,_Complex,2},{v,_Complex,2},{n,_Integer}}, Map[Tr,Dot[NestList[m.#.ConjugateTranspose[m]&,v,n-1],v]]];


The compiled version is only a few times faster than the uncompiled version. As an example

l = 128;
v = KroneckerProduct[SparseArray[{{1, 1} -> 1., {2, 2} -> -1.}],SparseArray[{i_, i_} -> 1., l]];
h := 1/2 ArrayFlatten[{{PadLeft[1/2 SparseArray[{Band[{1, 2}] -> 1., Band[{2, 1}] -> -1.}, l], {l, 2 l}]}, {PadRight[-(1/2)SparseArray[{Band[{1, 2}] -> 1., Band[{2, 1}] -> -1.}, l], {l, 2 l}]}}] + 1/2 KroneckerProduct[SparseArray[{{1, 1} -> 1., {2, 2} -> -1.}, 2],DiagonalMatrix[SparseArray[Table[RandomReal[{-1., 1.}], {j, l}]]]];

step[m_, u_, n_] := 1/(2 l) Map[Tr,Dot[NestList[MatrixExp[-I m].#.MatrixExp[I m] &, v, n - 1],v]];
stepcp := Compile[{{m, _Complex, 2}, {u, _Complex, 2}, {n, _Integer}},1./(2 l) Map[Tr,Dot[NestList[m.#.ConjugateTranspose[m] &, u, n - 1], u]]];

step[h, v, 64]; // AbsoluteTiming
stepcp[Normal[MatrixExp[-I h]], Normal[v], 64]; // AbsoluteTiming


My questions are:

1. How to compile this function more efficiently?
2. The matrices h and v are SparseArray, why can't I directly pass them into compiled functions, why do I have to use Normal[]?
3. The CompilePrint is shown below, is the MainEvaluate in line 1 and 29 a source of slowdown, if so can it be circumvented?

My eventual goal is to run an efficiently compiled version of step[h,v,n] 10000 times and take the Mean. Any suggestions are greatly appreciated.

• I think you mis-spelled NestList; also the function that returns the complex conjugate is simply Conjugate in Mathematica. Or have you redefined those function? If so you should give us your definition. – MarcoB May 22 '15 at 16:16
• Sorry I meant to say ConjugateTranspose, I have corrected it – user64620 May 22 '15 at 16:20
• I don't get the same CompilePrint output as you do. In particular, you seem to have defined two arguments in your posted function, and yet your CompilePrint output indicates 4 arguments. Are you sure that the code you posted exactly corresponds to your working code? Also, why is $n$ not an argument to your compiled function? – MarcoB May 22 '15 at 16:48
• Look, do yourself a favor: go over your code VERY carefully, and post ONE version of your function that you want to troubleshoot, with the corresponding CompilePrint output. It's very hard to hit a moving target. Your question won't attract many answers if you post incomplete and inconsistent code. Furthermore, the second version of your function, in context, contains a parameter $l$ that the first version of the function does not contain. You should include that parameter in the function definition as well. – MarcoB May 22 '15 at 16:58
• Your question is nice (I didn't see the earlier revisions). But compilation will generally not gain you anything for most linear algebra functions, which are just library calls anyway. Compilation is mainly useful for procedural code. You should probably look at other (algorithmic) ways to optimize the code, such as the approach presented by Simon Woods below. Compilation could then be useful if this causes you to adopt a more procedural approach, but not otherwise. – Oleksandr R. May 23 '15 at 20:25

## 2 Answers

This gives a slight speed up (~40%) for me:

stepcp = Compile[{{m, _Complex, 2}, {u, _Complex, 2}, {n, _Integer}, {l, _Real}},
Block[{
m2 = ConjugateTranspose@m, x = u, ut = Flatten@Transpose@u},
Join[{1.0}, Table[x = m.x.m2; Flatten[x].ut, {i, n - 1}]/(2 l)]]]


I replaced the NestList with a more procedural approach, calculating the trace of the dot product at each step using the relation Tr[a.b] == Flatten[a].Flatten[Transpose[b]]

I also put l in as an argument - there's no real performance increase as a result but I find it neater than having a callback to the main evaluator.

Oleksandr R. points out that the run time is dominated by calls to the math library so compilation is not really helping much. In fact almost all the gain you got from the compiled version versus the original was due to not computing MatrixExp at each step. You get almost the same speed with this:

step[m_, u_, n_, l_] :=
With[{m1 = MatrixExp[-I m], m2 = MatrixExp[I m]},
1/(2 l) Map[Tr, Dot[NestList[m1.#.m2 &, v, n - 1], v]]]

• I find that compilation hardly helps here. With or without it, the run time seems to be completely dominated by MKL calls. – Oleksandr R. May 23 '15 at 20:37
• Thank you for your response. Can you elaborate a bit more on why using a procedural approach is better here? @Simon Woods – user64620 May 24 '15 at 1:08
• @user64620, I needed to use an explicit loop so that I could compute the trace at each step. Since Compile translates NestList into a procedural loop anyway (as you can see from the CompilePrint output), there's nothing to lose from writing it that way. – Simon Woods May 24 '15 at 8:19
• @user64620, More generally, I find it useful to have the Mathematica code closely reflect the compiled instructions, so that I can understand better what the compiler is doing. Sometimes you can spot optimisations more easily that way. Though as Oleksandr pointed out, it makes little difference for this problem since all the time is spent doing the dot products. – Simon Woods May 24 '15 at 8:27
1. Sorry, I don't know. Have you tried compiling to C? One point: you should define stepcp with Set (=) rather than SetDelayed (:=) or you recompile it every time you use it.

2. SparseArray is a specialized format of the high-level Mathematica language. It is not supported by compilation. Many operations are optimized to work on sparse arrays, saving unnecessary computation, which may be why your uncompiled method is proportionately faster than you expect. You would essentially have to reimplement sparse arrays to get the same optimizations in your compiled code.

3. Tr is not a compilable function. (Don't ask me why.) Compiled alternatives that I can think of are not as fast as Tr to begin with so the overhead of a MainEvaluate call seems to balance out. An example alternative to Tr:

Compile[{{m, _Real, 2}}, Sum[CompileGetElement[m, i, i], {i, Length @ m}]]
`
• Thank you for your suggestions. Compiling to C takes about the same amount of time. @Mr.Wizard♦ – user64620 May 23 '15 at 15:35