0
$\begingroup$

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
   ]
$\endgroup$
12
  • $\begingroup$ By "will not run" do you mean "doesn't terminate in the time I've waited"? $\endgroup$
    – lericr
    Commented Oct 11, 2022 at 22:15
  • 1
    $\begingroup$ You could try Chebyshev interpolation. There's Statistics`Library`BarycentricInterpolation to help with this. Search the site or see mathematica.stackexchange.com/questions/191070/… $\endgroup$
    – Michael E2
    Commented Oct 11, 2022 at 22:39
  • $\begingroup$ Try wrapping your Plot with Monitor to see where it gets stuck. $\endgroup$
    – Chris K
    Commented Oct 11, 2022 at 22:43
  • $\begingroup$ Have you tried using Evaluate? Plot[Evaluate@gap[h, 12], ...] $\endgroup$
    – Bob Hanlon
    Commented Oct 11, 2022 at 23:04
  • $\begingroup$ 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. $\endgroup$
    – Diffycue
    Commented Oct 11, 2022 at 23:12

1 Answer 1

4
$\begingroup$

When I ran

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

I noticed it got stuck on 0. Avoiding integers fixes it:

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

runs in 0.08s.

$\endgroup$

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.