3
$\begingroup$

I am writing a script to simulate the spins. I have used "For" loop, but the execution takes long time as I have a lot of dataset. Has there any way to make it faster ? Here is the code:

Ix = 1/2 {{0, 1}, {1, 0}}; Iy = 1/2 {{0, -I}, {I, 0}}; Iz = 
 1/2 {{1, 0}, {0, -1}}; Io = IdentityMatrix[2];
I1z = KroneckerProduct[Iz, Io];
I2z = KroneckerProduct[Io, Iz];
I1x = KroneckerProduct[Ix, Io];
I2x = KroneckerProduct[Io, Ix];
I1y = KroneckerProduct[Iy, Io];
I2y = KroneckerProduct[Io, Iy];
 \[Rho]1={{1/2, 0, I/2, 0}, {0, -(1/2), 0, I/2}, {-(I/2), 0, 1/2, 
     0}, {0, -(I/2), 0, -(1/2)}}

  \[Omega]D = 2 \[Pi] 20; 
  \[Omega]CS = 2 \[Pi] 5 ; 
  \[Omega]r = 2 \[Pi] 15;
  \[Gamma] = \[Beta] = \[Pi]/4;

 aq = 100(*ms*); DwL = 0.002(*Dwell time in ms*);
 time = Table[i, {i, 0, aq, DwL}];
 
 U0t = ConstantArray[0, {Length@time - 1, 2, 2}]; \[Rho]t = 
  ConstantArray[0, {Length@time, 2, 2}];
 Mxt = Myt = ConstantArray[0, Length@time - 1];
    
  H0 = \[Omega]D (\[Sqrt]2 Sin[
         2 \[Beta]] Cos[\[Omega]r t + \[Gamma]] - (Sin[\[Beta]])^2 Cos[
         2 \[Omega]r t + 2 \[Gamma]]) (2 I1z.I2z - I1x.I2x - I1y.I2y) +
         \[Omega]CS (-\[Sqrt]2 Sin[
         2 \[Beta]] Cos[\[Omega]r t + \[Gamma]] + (Sin[\[Beta]])^2 Cos[
         2 \[Omega]r t + 2 \[Gamma]]) (I1z + I2z);  

   integrals = Table[Integrate[H0, {t, time[[i]], time[[i + 1]]}],{i, 1, Length@time - 1}];
   
    For[i = 1, i < Length@time, i++, dt = time[[i + 1]] - time[[i]];
        U0t[[i]] = MatrixExp[-I integrals[[i]]];
        \[Rho]t[[i]] = U0t[[i]].\[Rho]1.U0t[[i]]\[ConjugateTranspose];
         Mxt[[i]] = Tr[\[Rho]t[[i]].I1x];
         Myt[[i]] = Tr[\[Rho]t[[i]].I1y];
    ]

Any suggestion will be really helpful!

$\endgroup$
6
  • 1
    $\begingroup$ Start here: mathematica.stackexchange.com/q/134609/12 Do not use a For loop. Once it's rewritten in more structured form, it'll almost certainly be faster, easier to parallelize, and much easier to see where speed could be improved. $\endgroup$
    – Szabolcs
    Nov 21, 2023 at 12:09
  • $\begingroup$ Also, avoid unnecessary symbolic computation. Do you really need these integrals symbolically? I didn't wait for them to finish computing. I expect the result is huge, unproductive to use, and takes ages to evaluate numerically. $\endgroup$
    – Szabolcs
    Nov 21, 2023 at 12:12
  • $\begingroup$ yes, I need the integration, there is no way to avoid.. I tried to integrate first before use it in the loop. but that does not improve much. $\endgroup$
    – P Pyne
    Nov 21, 2023 at 12:15
  • 2
    $\begingroup$ @Szabolcs Integrate uses the numeric branch here, because the interval endpoints are machine precision numbers. The problem is just that this are so many integrals. So the key here is actually to do the integral symbolically once. ;) $\endgroup$ Nov 21, 2023 at 12:15
  • 1
    $\begingroup$ [Rho]1={{1/2, 0, I/2, 0}, {0, -(1/2), 0, I/2}, {-(I/2), 0, 1/2, 0}, {0, -(I/2), 0, -(1/2)}}. $\endgroup$
    – P Pyne
    Nov 21, 2023 at 12:24

2 Answers 2

5
$\begingroup$

The time consuming step is not the For loop, but the many integrals. Therefore, calculate the indefinite integral first:

integ[t_] = Integrate[H0 // FullSimplify, t]

enter image description here

Now you can calculate the table of integrals:

integs =  Table[integ[time[[i + 1]]] - integ[time[[i]]], {i, 1,    Length@time - 1}]

With this we can get U0t;

U0t =  MatrixExp[-I #] & /@ integs;

And from this we can get the rest:

ρ1 = {{1/2, 0, I/2, 0}, {0, -(1/2), 0, I/2}, {-(I/2), 0, 1/2, 
    0}, {0, -(I/2), 0, -(1/2)}};
ρt = # . ρ1 . #\[ConjugateTranspose] & /@ U0t;
Mxt = Tr[# . I1x] & /@ ρt;
Myt = Tr[# . I1y] & /@ ρt;
$\endgroup$
5
$\begingroup$

It takes already extremely long to compute integral. Note that your integrands are basically a functions times a constant $4 \times 4$ matrix. So just integrating the function and then multiplying by the matrix should help already. We can accelerate this further by employing the Fundamental Theorem of Calculus: We use Integrate to compute anti-derivatives symbolically and then compute differences of their function values to get the integrals. Finally we use TensorProduct to multiply the scala integrals with the matrix:

F[t_] = ComplexExpand[ Re[Integrate[\[Omega]D (\[Sqrt]2 Sin[2 \[Beta]] Cos[\[Omega]r t + \[Gamma]] - (Sin[\[Beta]])^2 Cos[2 \[Omega]r t + 2 \[Gamma]]), t]]];
G[t_] = ComplexExpand[ Re[Integrate[\[Omega]CS (-\[Sqrt]2 Sin[2 \[Beta]] Cos[\[Omega]r t + \[Gamma]] + (Sin[\[Beta]])^2 Cos[2 \[Omega]r t + 2 \[Gamma]]), t]]];
\[CapitalDelta]F = Differences[F[time]];
\[CapitalDelta]G = Differences[G[time]];
integrals =  TensorProduct[\[CapitalDelta]F, (2 I1z . I2z - I1x . I2x - I1y . I2y)] + TensorProduct[\[CapitalDelta]G, (I1z + I2z)];

All entries of integrals compute in a fraction of a second.

And here some further performance improvements. Mathematica is an interpreted language, and those do typically not well with for loops. While you code would be nearly optimimal in C or C++, there are better ways to achive your goals in Mathematica. For example, you do not need to preallocate arrays for your results if you can use command Map (/@). Btw., Map is one of the reasons why I adopted Mathematica in the first place; it is so much easier to understand, less error prone and often faster at runtime than Do, For, or Table.

Next thing we can use is to let Dot operation thread over the list of U0t. Here we can use Compile with the option RuntimeAttributes -> {Listable}.

Finally, Tr[A.B] is always a bad idea, because A.B is an $O(n^3)$ operation while Dot[Flatten[A].Flatten[B]] does the same and only requires $O(n^2)$ time. By clever reshaping of \[Rho]t with Flatten, we can even compute all Dots with a single call to a matrix-vector Dot. The backend BLAS library is highly optimized and performs the task 35 times faster on my machine. (And I think the Flatten is actually for free here.)

\[Rho]1 = {{1/2, 0, I/2, 0}, {0, -(1/2), 0, I/2}, {-(I/2), 0, 1/2, 0}, {0, -(I/2), 0, -(1/2)}};

cf = With[{\[Rho]1 = \[Rho]1},
   Compile[{{U, _Complex, 2}}, U . \[Rho]1 . U\[ConjugateTranspose], 
    RuntimeAttributes -> {Listable}]
   ];


U0t = MatrixExp /@ (-I integrals);
\[Rho]t = cf[U0t];

With[{\[Rho] = Flatten[\[Rho]t, {{1}, {2, 3}}]},
  Mxt = \[Rho] . Flatten[I1x];
  Myt = \[Rho] . Flatten[I1y];
  ];
$\endgroup$
5
  • $\begingroup$ This tensor product thing reduces the time significantly (the whole script runs within 5seconds and before it took 5hours). many thanks $\endgroup$
    – P Pyne
    Nov 21, 2023 at 12:24
  • $\begingroup$ You're welcome! $\endgroup$ Nov 21, 2023 at 12:27
  • $\begingroup$ In the script I actually have to vary the [Beta] and [Gamma] values as well. The code works fine even when I vary the values but a bit slower again to perform all the integration. Has there a way to get rid of it? Here is the code: [Beta] = [Gamma] = Table[i, {i, [Pi]/36, 179 [Pi]/180, [Pi]/36}]; [Omega]Dt = [Omega]D ([Sqrt]2 Sin[ 2 [Beta]] Cos[[Omega]r t + [Gamma]] - (Sin[[Beta]])^2 Cos[ 2 [Omega]r t + 2 [Gamma]]) ; Integrate[[Omega]Dt, t]; $\endgroup$
    – P Pyne
    Nov 21, 2023 at 17:15
  • 1
    $\begingroup$ Delete the definition of \[Beta] and \[Gamma], and replace F[t_] = ... by F[t_, \[Beta]_, \[Gamma]_] = ... and G[t_] = ... by G[t_, \[Beta]_, \[Gamma]_] = ... as well. This does the symbolic integration once, for arbitrary values of \[Beta] and \[Gamma]. Then in the loop over \[Beta] and \[Gamma] do just \[CapitalDelta]F = F[time, \[Beta], \[Gamma]] and \[CapitalDelta]G = G[time, \[Beta], \[Gamma]] and proceed from there as before. $\endgroup$ Nov 21, 2023 at 18:32
  • $\begingroup$ Thank you for your suggestion. It helps. $\endgroup$
    – P Pyne
    Nov 22, 2023 at 16:05

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.