Efficiently calculating matrix exponential of a recursively generated matrix

Question

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}]

$n=12$">

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

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*)
ClearAll[A];
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]}}];

Example:

A[30]; // AbsoluteTiming // First
Dimensions[A[30]]

0.363932

{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]

5.38649

6.14281

All this is done in less than 12 seconds.

Remarks

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.