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?

enter image description here

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.

  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    May 22, 2015 at 16:16
  • $\begingroup$ Sorry I meant to say ConjugateTranspose, I have corrected it $\endgroup$
    – user64620
    May 22, 2015 at 16:20
  • $\begingroup$ 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? $\endgroup$
    – MarcoB
    May 22, 2015 at 16:48
  • 2
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    May 22, 2015 at 16:58
  • 2
    $\begingroup$ 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. $\endgroup$ May 23, 2015 at 20:25

2 Answers 2


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

stepcp = Compile[{{m, _Complex, 2}, {u, _Complex, 2}, {n, _Integer}, {l, _Real}},
    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]]]
  • 2
    $\begingroup$ I find that compilation hardly helps here. With or without it, the run time seems to be completely dominated by MKL calls. $\endgroup$ May 23, 2015 at 20:37
  • $\begingroup$ Thank you for your response. Can you elaborate a bit more on why using a procedural approach is better here? @Simon Woods $\endgroup$
    – user64620
    May 24, 2015 at 1:08
  • $\begingroup$ @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. $\endgroup$ May 24, 2015 at 8:19
  • $\begingroup$ @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. $\endgroup$ May 24, 2015 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[Compile`GetElement[m, i, i], {i, Length @ m}]]
  • $\begingroup$ Thank you for your suggestions. Compiling to C takes about the same amount of time. @Mr.Wizard♦ $\endgroup$
    – user64620
    May 23, 2015 at 15:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.