Increase computation and Plot speed

Question

I have the following functions

EE[n_, m_] := eig[[n]] - eig[[m]];
b[t_, n_, m_] := (1/2)*
Sum[matrix[[n, k]]*
 matrix[[k, m]]*(EE[k, m]*Exp[I*EE[n, k]*t] - 
   EE[n, k]*Exp[I*EE[k, m]*t]), {k, 1, Ncut}];
c[t_, n_] := Sum[b[t, n, m]*Conjugate[b[t, n, m]], {m, 1, Ncut}];
Z[T_] := Sum[Exp[-(1/T)*eig[[n]]], {n, 1, Ncut}];
OTOC[t_, T_] := (1/Z[T])*
Sum[Exp[-(1/T)*eig[[n]]]*c[t, n], {n, 1, Ncut}];

where eig and matrix are known arrays, of size Ncut$\times1$ and Ncut$\times$Ncut respectively. I want to plot $\text{OTOC}(t,T)$ with respect to $t$ for different values of $T$. But even when I try to evaluate a single value, for example OTOC[1,1] it takes ages. This is because I am working with a large matrix, therefore Ncut is very large (it is 150 for my specific case). How can speed this code up? How can I produce plots faster?

Suggestion: You are welcome to run this same code, by generating random eig and matrix arrays, and try to optimize.

---

Edit upon request: If you want to run this code first, you can use the following for eig and matrix and ncut

ncut = 50;
matrix = Table[RandomComplex[], {i, ncut}, {i, ncut}];
matrix = (1/2)*(matrix + ConjugateTranspose[matrix]);
eig = Table[i^2/10, {i, ncut}];

Then run the following

EE[n_, m_] := eig[[n]] - eig[[m]];
b[t_, n_, m_] := (1/2)*
Sum[matrix[[n, k]]*
 matrix[[k, m]]*(EE[k, m]*Exp[I*EE[n, k]*t] - 
   EE[n, k]*Exp[I*EE[k, m]*t]), {k, 1, ncut}];
c[t_, n_] := Sum[b[t, n, m]*Conjugate[b[t, n, m]], {m, 1, ncut}];
Z[T_] := Sum[Exp[-(1/T)*eig[[n]]], {n, 1, ncut}];
OTOC[t_, T_] := (1/Z[T])*
Sum[Exp[-(1/T)*eig[[n]]]*c[t, n], {n, 1, ncut}];

After that you can try to plot the following

Plot[c[t,1],{t,0,10}]

and you can try to optimize the run time for the whole process.

Of course in my case, the arrays eig and matrix are meaningful arrays that I gather from some place, but that should not bother you. The whole point is to optimize this problem for any given eig and matrix.

Answer 1

Initialize

Ncut = 10;
eig = RandomReal[{0, 1}, Ncut];
matrix = RandomReal[{0, 1}, {Ncut, Ncut}];

Testing original code

Timing[Plot[Re[OTOC[t, 100.]], {t, 0., 10.}];]
(*15.3793*)

New code

 OTOCn[ev_, ma_, t_, T_] := Module[{nc, EE, A, B, b, c, z},
  nc = Length[ev];
  EE = Table[ev[[n]] - ev[[m]], {n, nc}, {m, nc}];
  A = ma*Exp[I*EE*t];
  B = ma*EE;
  b = 0.5*(A.B - B.A);
  c = Re[Diagonal[b.ConjugateTranspose[b]]];
  z = Total[Exp[-ev/T]];
  Total[Exp[-(1/T)*ev]*c]
 ]

Testing new code

Timing[Plot[OTOCn[eig, matrix, t, 100.], {t, 0., 10.}];]
(*0.361552*)

Speed up for Ncut=10 is about 42 times. For a large value of Ncut=100 the code needs 36 seconds.