Improving Performance of an XY Monte Carlo
You can use the RuntimeTools`Profile
function, shown here, to figure out which parts need optimization the most. It's better to compile everything, if possible, but it is also possible to compile just the functions that need it the most.
Applying this function to 100 iterations of your code, we get
RotationMatrix
really stands out here, so we can replace it by
rotationMatrix = Compile[{th}, {
{Cos[th], -Sin[th]},
{Sin[th], Cos[th]}
},
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];
This simple change brings the total time down from 9.3 seconds to 5.8 seconds.
Predefining those matrices outside the loop using
rotmats = Map[rotationMatrix, α RandomReal[{-pi, pi}, {L, L}], {2}];
and inside the loop u = rotmats[[i, j]].v
brings the time down to 5.0.
Proceeding down the list of things that take time, we find that the condition in the If
statement takes some time.
accept = Compile[{dE, T},
RandomReal[] < Min[1, Exp[-dE/T]],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];
brings the time down to 4.3.
Keep in mind that not all functions are compilable, a list of functions that are compilable is available here. When in doubt whether the code is entirely compilable, use CompiledFunctionTools`CompilePrint
and check for calls to the main evaluator. Sometimes it's worth it to rewrite some code just to make it compilable.
Others will show you more optimizations, I'm sure, I just wanted to share this way of finding opportunities for optimizations.
Here is a compiled version of Sweep
that employs complex numbers instead of 2-vectors and rotation.
cSweep = Compile[{{XY0, _Complex, 1}, {\[Alpha]0, _Real}, {T, _Real}, {n, _Integer},
{top, _Integer, 1}, {bottom, _Integer, 1}, {left, _Integer, 1}, {right, _Integer, 1}
},
Block[{XY, accepts, rand, \[Phi], u, v, z, dE, m, lo, hi, \[Alpha]},
m = Length[XY0];
lo = 0.45 m;
hi = 0.55 m;
XY = XY0;
\[Alpha] = \[Alpha]0;
Do[
accepts = 0.;
rand = RandomReal[{0., 1.}, m];
\[Phi] = RandomReal[\[Alpha] {-Pi, Pi}, m];
Do[
v = Compile`GetElement[XY, i];
u = v Exp[I Compile`GetElement[\[Phi], i]];
z = Plus[
Compile`GetElement[XY, Compile`GetElement[top, i]],
Compile`GetElement[XY, Compile`GetElement[bottom, i]],
Compile`GetElement[XY, Compile`GetElement[left, i]],
Compile`GetElement[XY, Compile`GetElement[right, i]]
];
dE = Re[Conjugate[z] (u - v)];
If[Compile`GetElement[rand, i] < Exp[Min[Max[-700., -dE/T], 0.]],
XY[[i]] = u;
accepts++;
];
, {i, 1, m}];
If[accepts < lo,
\[Alpha] *= .95;
,
If[accepts > hi, \[Alpha] *= 1.05];
];
, {n}];
XY
],
CompilationTarget -> "C",
RuntimeOptions -> "Speed"
];
Usage example:
L = 20;
T = 1.0;
\[Alpha] = 1.;
With[{idx = Partition[Range[L^2], L]},
bottom = Flatten[RotateLeft[idx]];
top = Flatten[RotateRight[idx]];
right = Flatten[idx[[All, RotateLeft[Range[L]]]]];
left = Flatten[idx[[All, RotateRight[Range[L]]]]];
];
XY0 = Exp[I RandomReal[{-Pi, Pi}, L^2]];
XY = cSweep[XY0, \[Alpha], T, 100000, top, bottom, left, right]; // AbsoluteTiming // First
3.3154
Please double-check whether the implementation is correct.
Probably it would be more efficient to adjust the algorithm such that all entries of XY
are modified at once. But of course, the rejection step has to be adjusted accordingly in order to guarantee the same equilibrium probability distribution for the stochastic process.