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

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!

• I think in this case ConstantArray[0., 500] would be the fasted and most accurate way to compute the result... ;p Commented Jul 1 at 14:30

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];
`