# Cannot plot a function that calls Eigensystem on a large matrix

Is there a way to plot an expensive function without resorting to making a list of values and using ListPlot?

I have a function called gap[h,size] which computes the difference between the lowest two eigenvalues in a large sparse matrix. h is a parameter in the matrix that varies from 0 to 1 while size controls the size of the matrix. The matrix contains 2^size * 2^size entries.

Computing gap[h,size] gets expensive as size gets larger, but it's reasonable to get all the way up to size=18 at which point each evaluation takes around 8 seconds.

Here is my problem: Although I can compute a list of values of gap[h,18] at different values of h, I cannot plot the function gap[h,18] at all. Even at a more modest value of size=12,

Plot[
gap[h, 12],
{h, 0, 1},
PlotPoints -> 10, MaxRecursion -> 0
]

will not run, even though computing

Table[gap[h, 12], {h, 0, 1, 0.1}]

takes only half a second. During evaluation Mathematica runs a single CPU core at 100% but doesn't hog any memory.

Here is some sample code to try for yourself. The matrix is already diagonal, but you cannot plot gap for values above size=10.

largeMatrix[h_, size_] :=
SparseArray[
{i_, i_} -> h,
{2^size, 2^size}
]
sparseIdentity[length_] :=
SparseArray[{i_, i_} -> 1, {2^length, 2^length}];
sparseEigensystem[matrix_, length_, valuesKept_] :=
Block[{values, vectors},
{values, vectors} =
Eigensystem[matrix - length*sparseIdentity[length], valuesKept];
{values + length, vectors}
]
ClearAll[gap]
gap[h_, length_] /; NumericQ[h] :=
gap[h, length] = Block[{values, vectors},
{values, vectors} =
sparseEigensystem[largeMatrix[h, length], length, 2];
EuclideanDistance @@ values
]
• By "will not run" do you mean "doesn't terminate in the time I've waited"? Oct 11, 2022 at 22:15
• You could try Chebyshev interpolation. There's StatisticsLibraryBarycentricInterpolation to help with this. Search the site or see mathematica.stackexchange.com/questions/191070/… Oct 11, 2022 at 22:39
• Try wrapping your Plot with Monitor to see where it gets stuck. Oct 11, 2022 at 22:43
• Have you tried using Evaluate? Plot[Evaluate@gap[h, 12], ...] Oct 11, 2022 at 23:04
• Thanks everyone, these are useful suggestions and I will try them. Sorry to not put the matrix in the question, but it's too big of course. Oct 11, 2022 at 23:12