# How to optimize Mathematica code that depends on eigenvalues of big matrices and big sums?

I've been using Mathematica recently to generate some plots of a few functions. I've been able to get it right after a few questions here.

The resulting code, which works, is this:

M00[nmax_Integer,t_] := SparseArray[Band[{1,1}] -> Table[t^(2n), {n,0,nmax}]]
M11[nmax_Integer,t_] := (1-t^2)*SparseArray[Band[{1,1}] -> Table[n*t^(2(n-1)), {n,0,nmax}]]
M01[nmax_Integer,t_] := Sqrt[1-t^2]*SparseArray[Band[{1,2}] -> Table[Sqrt[n+1]*t^(2n), {n,0,nmax-1}], {nmax+1,nmax+1}]
M10[nmax_Integer,t_] := Transpose[M01[nmax,t]]

rhoplus[nmax_Integer,t_,th_] := (1-t^2)/2*((1+Cos[th])*M00[nmax,t]+(1-Cos[th])*M11[nmax,t]+Sin[th]*(M10[nmax,t]+M01[nmax,t]))
rhominus[nmax_Integer,t_,th_]:= (1-t^2)/2*((1-Cos[th])*M00[nmax,t]+(1+Cos[th])*M11[nmax,t]-Sin[th]*(M10[nmax,t]+M01[nmax,t]))

s[x_] = Piecewise[{{x*Log[2,x], 0<x<1}}]
EntropyPlus[nmax_Integer, t_, th_] := Total[s /@    Eigenvalues[rhoplus[nmax,t,th]]]
EntropyMinus[nmax_Integer, t_, th_] := Total[s /@ Eigenvalues[rhominus[nmax,t,th]]]

SumAuxElement[t_][n_Integer] := ((1 - t^2)/2)*t^(2 n)*(1 + (n + 1) (1 - t^2))
SumAuxList[nmax_Integer, t_] := Array[SumAuxElement[t], nmax + 1, {0, nmax}]
SumAux[t_,nmax_Integer]:=Total[s /@ SumAuxList[nmax, t]]

Discord[nmax_Integer, t_, th_] := 1 + SumAux[t,nmax] - 1/2 (EntropyPlus[nmax, t, th] + EntropyMinus[nmax, t, th])

Plot[Discord[300, Tanh[r], Pi/2], {r, 0, 2.5}]


It works nice, but in my computer it takes forever to run depending on nmax. Indeed for nmax = 300 it takes several minutes to generate the plot. If I need additional plots of things that are written in terms of the same Discord function, soon it becomes unmenageable to work with this.

I wanted to optimize this code, in order to get something faster.

First intuitively speaking, I do notice what perhaps is the biggest two problems:

1. First every time a plot like that is run, the eigenvalues are recomputed for several values of t. This seems to be the greater cause of less of performance.

2. A second issue is that the sums are recomputed everytime for all values of t needed in the plot and the sums do take some time.

The first issue I believe somehow the ideal thing would be to in the beginning of the code compute the eigenvalues for a range of t once and for all and have it already saved to be used later. I don't know how to do it though.

So, using this code as an example, when I have code that needs to perform computations with eigenvalues of big matrices and evaluate big sums like this to generate plots, how can I optmize the code so that it runs faster?

Edit: I've tried out the proposals on the answers, and although both of them give faster results as expected, both in my computer changed a piece of the plot and generated some warnings. The last full code I've tried is this:

range[nmax_] := range[nmax] = Range[0, nmax] // N

tab1[nmax_Integer, t_] := tab1[nmax, t] = t^(2. range[nmax])
tab2[nmax_Integer, t_] := tab2[nmax, t] = # t^(2. (# - 1.)) &@range[nmax]
tab3[nmax_Integer, t_] := tab3[nmax, t] = Sqrt[# + 1.] t^(2. #) &@Most[range[nmax]]

xM00[nmax_Integer, t_] := SparseArray[Band[{1, 1}] -> tab1[nmax, t]]
xM11[nmax_Integer, t_] := SparseArray[Band[{1, 1}] -> (1. - t^2.) tab2[nmax, t]]
xM01[nmax_Integer, t_] := SparseArray[Band[{1, 2}] -> Sqrt[1. - t^2.] tab3[nmax, t], {nmax + 1, nmax + 1}]
xM10[nmax_Integer, t_] := Transpose[xM01[nmax, t]]

xrhoplus[nmax_Integer, t_, th_] := (1. - t^2.)/2. ((1. + Cos[th]) xM00[nmax, t] + (1. - Cos[th]) xM11[nmax, t] + Sin[th] (xM10[nmax, t] + xM01[nmax, t]))
xrhominus[nmax_Integer, t_, th_] := (1. - t^2.)/2. ((1. - Cos[th]) xM00[nmax, t] + (1. + Cos[th]) xM11[nmax, t] - Sin[th] (xM10[nmax, t] + xM01[nmax, t]))

cs = Compile[{{x, _Real}}, If[0. < x < 1., x Log[2, x], 0.],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True]

xEntropyPlus[nmax_Integer, t_, th_] := Total[cs@Eigenvalues[xrhoplus[nmax, t, th], Method -> "Banded"]]
xEntropyMinus[nmax_Integer, t_, th_] := Total[cs@Eigenvalues[xrhominus[nmax, t, th], Method -> "Banded"]]

xSumAuxElement = Compile[{{t, _Real}, {n, _Real,1}}, ((1. - t^2.)/2.) t^(2. n) (1. + (n + 1.) (1. - t^2.)), CompilationTarget -> "C", Parallelization -> True]

xSumAuxList[nmax_Integer, t_] := xSumAuxElement[t, range[nmax]]

xSumAux[t_, nmax_Integer] := Total[cs@xSumAuxList[nmax, t]]

xDiscord[nmax_Integer, t_, th_] := 1 + xSumAux[t, nmax] - 1/2 (xEntropyPlus[nmax, t, th] + xEntropyMinus[nmax, t, th])

Plot[xDiscord[300, Tanh[r], Pi/2.], {r, 0, 2.5}]


The output then is the following: So there are some warning messages and the piece of the plot to the left of 0.5 vanished. Is this perhaps something about my computer or am I missing something of the code?

I have further evaluated CompilePrint@cs and got:

    1 argument
1 Boolean register
1 Integer register
6 Real registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {Listable}

R0 = A1
R2 = 1.
I0 = 2
R3 = 7.
R1 = 0.
Result = R4

1   B0 = R1 < R0 < R2 (tol R3)
2   if[ !B0] goto 7
3   R4 = Log[ I0, R0]
4   R5 = R0 * R4
5   R4 = R5
6   goto 8
7   R4 = R1
8   Return


# The eigensolver

The default eigensolver method for SparseArrays is "Arnoldi"; it is not good at computing the whole spectrum of an operator. Luckily, you matrices are tridiagonal, so you can use a banded solver. So, better use

EntropyPlus[nmax_Integer, t_, th_] :=
Total[s /@ Eigenvalues[rhoplus[nmax,t,th], Method -> "Banded"]]
EntropyMinus[nmax_Integer, t_, th_] :=
Total[s /@ Eigenvalues[rhominus[nmax,t,th], Method -> "Banded"]]


This should get you rid of quite some message handling (precision issues). Moreover, it should be more accurate and faster, especially for larger matrices.

# Compiling some helper function

Using the following should also help. In particular, cf is implement such that underflow of machine precision numbers cannot happen.

cs = Compile[{{x, _Real}},
If[0. < x < 1., x Log2[x], 0.],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
];
EntropyPlus[nmax_Integer, t_, th_] :=
Total[cs@Eigenvalues[rhoplus[nmax, t, th], Method -> "Banded"]]
EntropyMinus[nmax_Integer, t_, th_] :=
Total[cs@Eigenvalues[rhominus[nmax, t, th], Method -> "Banded"]]
cf = With[{threshold = 0.5 Log[\$MinMachineNumber]},
Compile[{{t, _Real}, {n, _Integer}},
Block[{x},
If[Log[t] (2 n) < threshold,
0.,
x = ((1. - t^2) 0.5) t^(2 n) (1. + (n + 1) (1. - t^2));
If[0. < x < 1., x Log2[x], 0.]
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]
];
SumAux[t_, nmax_Integer] := Total[cf[t, Range[0, nmax]]]

• Thanks @HendrikSchumacher ! I have two questions: the first of those is that with this specific method I got the message "The method "Banded" accepts only sparse matrices with elements that are machine-real or machine-complex numbers." and the plot now misses a part (from r = 0 to r = 0.9 the graph disappeared). Is something more that needs to be changed to use this method? The second question is: and in a more general context (not tridiagonal matrices) is there some general approach to get a faster code? For instance, my idea of "precomputing" the eigenvalues even makes some sense? Thanks! – user1620696 Jan 28 '19 at 14:44
• "The method "Banded" accepts only sparse matrices with elements that are machine-real or machine-complex numbers." -- This indicates that you have tried to feed Eigenvalues with symbolic matrices. Maybe that is because Plot first tries to evaluate its first argument symbolically. You can prevent that by defining Discord  such that enforces numeric arguments to execute with Discord[nmax_Integer, t_?NumberQ, th_?NumberQ] := ... – Henrik Schumacher Jan 28 '19 at 14:59
• Towards the second question: You have a two-parameter family of matrices. Maybe that can be exploited in a way such that eigenvalue information can recycled or inexpensively compute from precomputed data. But that really depends on the family oft matrices. I don't know of a general way to do that. – Henrik Schumacher Jan 28 '19 at 15:02
• Well, I'm giving to Eigenvalues the matrix rhoplus[nmax,t,th]. I think that somehow by standard Mathematica evaluates this a symbolic matrix and then later feeds the values. I tried redefining Discord as you said and further redefined all the functions definitions to use t_?NumberQ and th_?NumberQ but the message persists. Do you have any idea what may be wrong? Thanks again! – user1620696 Jan 28 '19 at 16:19
• Are nmax, t, th really numeric? Notice that I meant to add Method -> "Banded" into your preexisting code. Executing Eigenvalues[rhoplus[nmax, t, th], Method -> "Banded"] as standalone will (of course) produce nonsense... – Henrik Schumacher Jan 28 '19 at 16:27
1. the Mxx functions waste time with Table; instead use an appropriate Range for all of them; additionally, use memoization (probably helpful...); use the modified Mxx's to define modified versions of rhosign functions:

(* the utility of N will become apparent later *)
range[nmax_] := range[nmax] = Range[0, nmax] // N

(* tables needed to define the Mxx functions *)
tab1[nmax_Integer, t_] := tab1[nmax, t] = t^(2. range[nmax])
tab2[nmax_Integer, t_] := tab2[nmax, t] = # t^(2. (# - 1.)) &@range[nmax]
tab3[nmax_Integer, t_] := tab3[nmax, t] = Sqrt[# + 1.] t^(2. #) &@Most[range[nmax]]

(* the modified Mxx functions using the tables above, that in turn use range *)
xM00[nmax_Integer, t_] := SparseArray[Band[{1, 1}] -> tab1[nmax, t]]
xM11[nmax_Integer, t_] := SparseArray[Band[{1, 1}] -> (1. - t^2.) tab2[nmax, t]]
xM01[nmax_Integer, t_] := SparseArray[Band[{1, 2}] -> Sqrt[1. - t^2.] tab3[nmax, t], {nmax + 1, nmax + 1}]
xM10[nmax_Integer, t_] := Transpose[xM01[nmax, t]]

(* the rhosign modified functions *)
xrhoplus[nmax_Integer, t_, th_] := (1. - t^2.)/ 2. ((1. + Cos[th]) xM00[nmax, t] + (1. - Cos[th]) xM11[nmax, t] + Sin[th] (xM10[nmax, t] + xM01[nmax, t]))
xrhominus[nmax_Integer, t_, th_] := (1. - t^2.)/ 2. ((1. - Cos[th]) xM00[nmax, t] + (1. + Cos[th]) xM11[nmax, t] - Sin[th] (xM10[nmax, t] + xM01[nmax, t]))

2. use Compile to define the Entropysign functions, as already suggested in the answer by Henrik Schumacher:

cs = Compile[{{x, _Real}},

If[0. < x < 1., x Log[2, x], 0.],

CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]

(* Use Method -> "Banded" for Eigenvalues as suggested by **Henrik Schumacher** *)
xEntropyPlus[nmax_Integer, t_, th_] := Total[cs@Eigenvalues[xrhoplus[nmax, t, th], Method -> "Banded"]]
xEntropyMinus[nmax_Integer, t_, th_] := Total[cs@Eigenvalues[xrhominus[nmax, t, th], Method -> "Banded"]]

3. use Compile to assist the SumAux functions:

(* the second argument for this function makes use of the fact that range is N - see below *)
xSumAuxElement = Compile[{{t, _Real}, {n, _Real, 1}},

((1. - t^2.)/2.) t^(2. n) (1. + (n + 1.) (1. - t^2.)),

CompilationTarget -> "C",
Parallelization -> True
]

xSumAuxList[nmax_Integer, t_] := xSumAuxElement[t, range[nmax]]

(* another point where Compile helps... *)
xSumAux[t_, nmax_Integer] := Total[cs@xSumAuxList[nmax, t]]


Most of the work is done; all that's left is defining an appropriately modified main function (Discord) ie:

 xDiscord[nmax_Integer, t_, th_] := 1 + xSumAux[t, nmax] - 1/2 (xEntropyPlus[nmax, t, th] + xEntropyMinus[nmax, t, th])


Finally, evaluating Plot[xDiscord[300, Tanh[r], Pi/2.], {r, 0, 2.5}] returns almost instantly: • Thanks @user42582. It is much faster but I have two issues with it. The first is that I don't get theright plot. It misses a part, in other words, to the left of something around 0.5 in the x axis, the plot vanishes. A second issue is that a get a few warning messages: "Eigenvalues: The method "Banded" accepts only sparse matrices with elements that are machine-real or machine-complex numbers" and "CompiledFunction: Argument Eigenvalues[SparseArray[Automatic,<<3>>],Method->Banded] at position 1 should be a machine-size real number". Any idea why it doesn't fully work in my computer? – user1620696 Jan 29 '19 at 12:19
• @user1620696 I'm not sure what's wrong; I have just rechecked the code in my answer and everything evaluates without generating messages; please give a detailed account of the code you used and relevant input – user42582 Jan 29 '19 at 14:06
• I've added one edit with the code I've ran and the output. The code I've tried is basically yours, since it worked perfectly for you. There are some warnings and the plot misses one piece. Is this perhaps something related to my computer or am I missing something about the code proposal? I think this is not about the code, since I've tried exactly the one which works for you. – user1620696 Jan 29 '19 at 14:25
• @user1620696 I've run the code you provided in the edit-which actually seems identical to the code in my answer, like you said-and I encounter no message or plot problems; if the problem persists even after a clean run (quitting the kernel and evaluating everything on a fresh kernel) then the most probable culprit I can think of is some issue with Compile; please evaluate the following and report back : load package << "CompiledFunctionTools`" and then evaluate eg CompilePrint@cs – user42582 Jan 29 '19 at 14:44
• I have evaluated the code you mention and posted the output. By the way the compiler I'm using is the Visual Studio compiler. I have it in my machine because I work with Visual Studio as well, and thought it would be enough for Mathematica as well. Is perhaps the compiler the problem ? – user1620696 Jan 29 '19 at 14:58