Given a matrix $A$ of size $N \times N$, I want to evaluate the function \begin{equation} F(t) = |\operatorname{Tr}(\exp(-i A t))|^2 \end{equation} as a function of $t$. Here, $A[n]$ is constructed recursively from $A[n-1]$ and $A[n-2]$, so that the size of $A[n]$ is $\operatorname{Fib}[n]$, the $n^{th}$ Fibonacci number. I'd like to be able to push this up to $n \sim 20$, but my current code is taking a long time to evaluate even for $n=12$.

Also, I want to plot this on a log-log scale, but when I use LogLogPlot instead of Plot, it doesn't evaluate even for $n=12$. I'm also confused because when I try to generate a ListPlot instead, it runs out of memory but is still somehow able to generate a Plot (at least for $n=12$). So right now, I first generate the plot and then sample data from that to put it on a log-log scale. I'm certain there has to be a more efficient way of doing this, but cannot pinpoint my error

Here is my current code:

(*Define matrices required for recursion relation*)
Id[n_] := IdentityMatrix[Fibonacci[n - 2]];
Zm[n_] := ConstantArray[0, {Fibonacci[n - 3], Fibonacci[n - 2]}];
J[n_] := ArrayFlatten[{{Id[n]}, {Zm[n]}}];

(*Build Matrix Recursively*)
A[2] = {{0}};
A[3] = {{0, 1}, {1, 0}};
A[n_Integer?Positive] := ArrayFlatten[{{A[n - 1], J[n]}, {Transpose[J[n]], A[n - 2]}}];

(*Define function to be plotted*)
F[n_, t_] := Abs[Tr[N[MatrixExp[-I t A[n]]]]]^2;
plot[n_] := Plot[F[n, t], {t, 0, 10^3}, PlotRange -> All, AxesOrigin -> {0,0}];
data[n_] := Cases[Cases[InputForm[plot[n]], Line[___], Infinity], {_?NumericQ, _?NumericQ}, Infinity];

Here is the output I am currently able to produce:

ListLogLogPlot[data[12], PlotRange -> {{10^-1, 100}, Automatic}]

Output for <span class=$n=12$">

I would appreciate any help for speeding up this code so that it can be pushed up to higher $n$.


1 Answer 1


Matrix construction

First of all, it might be a good idea to use memoization for building the matrices because this avoids recursive calls to A. By using SparseArrays, we can generate the matrices for much higher n:

(*define matrices required for recursion relation*)

Id[n_] := IdentityMatrix[Fibonacci[n - 2], SparseArray];
Zm[n_] := SparseArray[{}, {Fibonacci[n - 3], Fibonacci[n - 2]}, 0];
J[n_] := ArrayFlatten[{{Id[n]}, {Zm[n]}}];

(*build matrix recursively using memoization*)
A[2] = SparseArray[{{0}}];
A[3] = SparseArray[{{0, 1}, {1, 0}}];
A[n_Integer?Positive] := A[n] = ArrayFlatten[{{A[n - 1], J[n]}, {Transpose[J[n]], A[n - 2]}}];


A[30]; // AbsoluteTiming // First


{832040, 832040}

Notice that your code enforces that A[n] has to be recomputed each time you call the function F[n, t]!

Computation of trace of matrix exponential

Since the matrices in question are symmetric and since in the end, only the trace of the matrix exponential is required, we may utilize the functional calculus of self-adjoint operators and the fact that the trace is invariant under conjugation with unitary matrices.

a = N[A[20]];
\[CapitalLambda] = Eigenvalues[Normal[a]]; // AbsoluteTiming // First
f = t \[Function] Abs[Total[Exp[-I t \[CapitalLambda]]]]^2;

tlist = Subdivide[0, 10^3, 10000];
data = Transpose[{tlist, Map[f, tlist]}]; // AbsoluteTiming // First
ListLogLogPlot[data, Joined -> True]



enter image description here

All this is done in less than 12 seconds.


The most expensive part of the computations for larger n is Eigenvalues. Since the whole spectrum is needed here, eigenvalue algorithms for sparse arrays won't work too well in the long run. However, I am pretty sure that the spectrum of A[n] can be determined recursively from the spectra of A[n-2] and A[n-1] without computing A[n]. That would make most of the computations above obsolete.

  • $\begingroup$ Thank you very much for your help! That works much better. I've thought about trying to get the spectrum of $A[n]$ from that of $A[n-1]$ and $A[n-2]$ but couldn't figure out how they were related. It would be great if a recursion relation could be obtained, I would appreciate any ideas you have on that front. $\endgroup$
    – Aegon
    Oct 5, 2019 at 3:49
  • $\begingroup$ Hey @Aegon. Well I thought a bit about it, but this seems to be not as simple as I expected. I tried to figure out a recursive formula for the characteristic polynomials of the matrices A[n], but I finally gave up. Maybe anybody else might find a solution for this. $\endgroup$ Oct 7, 2019 at 20:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.