Use of Compile[] to optimize the code to run faster

Question

I am working on a code that takes long time to evaluate, So I come to know that use of Compile[] function can speed up significantly. Here is my original code.

    Ix = PauliMatrix[1]; Iy = PauliMatrix[2]; Iz = PauliMatrix[3];
\[Rho]sample[\[Beta]_, \[Gamma]_] := Cos[\[Beta]] Cos[\[Gamma]] Iz;
Ud[\[Beta]_, \[Gamma]_] := 
  Table[i Cos[\[Beta]] Cos[\[Gamma]] (Ix-Iy+Iz), {i, 1, 500}];

 \[Rho]dc[\[Beta]_, \[Gamma]_] := 
  Dot[#, \[Rho]sample[\[Beta], \[Gamma]], #\[ConjugateTranspose]] & /@
    Ud[\[Beta], \[Gamma]];
FreeEvolve[\[Beta]_, \[Gamma]_] := 
  Tr[(Ix - I Iy)\[ConjugateTranspose].#] & /@ \[Rho]dc[\[Beta], \
\[Gamma]];
Selectedangle = 
  Flatten[Outer[{#1, #2} &, Range[0, \[Pi], \[Pi]/200], 
    Range[0, 2 \[Pi], \[Pi]/30]], 1];
FID = Total[
    ParallelMap[FreeEvolve[#[[1]], #[[2]]] &, 
     Selectedangle]]; // AbsoluteTiming

Following the code, I tried with Compile[].

compiledRhoDc = 
  Compile[{{\[Beta], _Real}, {\[Gamma], _Real}}, 
   Block[{ud, rhoSample, rhoDcList, i}, ud = Ud[\[Beta], \[Gamma]];
    rhoSample = \[Rho]sample[\[Beta], \[Gamma]];
    rhoDcList = 
     Table[Dot[ud[[i]], rhoSample, ConjugateTranspose[ud[[i]]]], {i, 
       Length[ud]}];
    rhoDcList], RuntimeAttributes -> {Listable}, 
   Parallelization -> True];
compiledFreeEvolve = 
  Compile[{{\[Beta], _Real}, {\[Gamma], _Real}}, 
   Block[{rhoDcList, CxConj, CyConj, result, i}, 
    CxConj = ConjugateTranspose[Ix];
    CyConj = ConjugateTranspose[Iy];
    rhoDcList = compiledRhoDc[\[Beta], \[Gamma]];
    result = 
     Table[Tr[(CxConj + I CyConj).rhoDcList[[i]]], {i, 
       Length[rhoDcList]}];
    result], RuntimeAttributes -> {Listable}, Parallelization -> True];
compFID = 
   Total[ParallelMap[compiledFreeEvolve[#[[1]], #[[2]]] &, 
     Selectedangle]]; // AbsoluteTiming

When I try to evaluate, it gives some error although it returns same output as that of original code. The errors appear as follows:

*Compile::part: Part specification ud[[i]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.

Compile::cplist: Compile`$15 should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function.*

Any lead will be really helpful to improve the efficiency using Compile[] or other ways. Thanks in advance!

Answer 1

The strength of Compile is that it accelerates loop constructs (like Table,Sum,Map). So it is better to put the loops within Compile. Also it is a good idea to simplify the symbolic expressions before compilation. Here the intermediate values are matrices, but the result for each {i, \[Beta], \[Gamma]} is just a number. Moreover, the numbers for all {\[Beta], \[Gamma]} are summed. So it makes some sense to add them together directly; this saves many memory operations (which are way more expensive than actual flops). Here is a short implementation of a function cf that should return the result of FID.

Block[{A, B, whatWeWantToCompute, Ix, Iy, Iz, i, \[Beta], \[Gamma]},
 Ix = PauliMatrix[1];
 Iy = PauliMatrix[2];
 Iz = PauliMatrix[3];
 A = i Cos[\[Beta]] Cos[\[Gamma]] (Ix - Iy + Iz);
 B = Cos[\[Beta]] Cos[\[Gamma]] Iz;
 whatWeWantToCompute = 
  ComplexExpand[
   Tr[Dot[(Ix - I Iy)\[ConjugateTranspose], 
     Dot[A, B, A\[ConjugateTranspose]]]]];
 
 cf = With[{code = N@whatWeWantToCompute},
   Compile[{{n, _Integer}, {\[Beta]list, _Real, 1}, {\[Gamma]list, _Real, 1}},
    Table[ Sum[code, {\[Beta], \[Beta]list}, {\[Gamma], \[Gamma]list}], {i, 1, n}],
    CompilationTarget -> "C",
    RuntimeOptions -> "Speed"
    ]
   ];
 ]

Now I get FID in 0.026 milliseconds like this:

\[Beta]list = Subdivide[0., \[Pi], 200];
\[Gamma]list = Subdivide[0., 2 \[Pi], 30];
result = cf[500, \[Beta]list, \[Gamma]list];

Moreover, when you inspect whatWeWantToCompute, then you will see that it looks like this:

(4 - 4 I) i^2 Cos[\[Beta]]^3 Cos[\[Gamma]]^3

So we can get the result in roughly 1/500 of the time with this:

result2 = Range[1, 500]^2 cf[1, \[Beta]list, \[Gamma]list][[1]];

Or entirely without Compile like this:

result3 = (4 - 4 I) Range[1, 500]^2 Total[TensorProduct[Cos[\[Beta]list]^3, Cos[\[Gamma]list]^3], 2];

Or by observing that it is not neccessary to form the tensor product first:

result4 = (4 - 4 I) Range[1, 500]^2 Total[Cos[\[Beta]list]^3] Total[
 Cos[\[Gamma]list]^3]; // AbsoluteTiming

Or by observing that Total[Cos[\[Beta]list]^3] should actually be 0 because $\cos^3$ is an odd function on the interval $[0,\pi]$ (and because you are using a symmetric quadrature rule):

result5 = ConstantArray[0, 500];