# Speed up the computation of the trace of a matrix

I have a symbolic symmetric $6\times6$ matrix:

dim = 6;
SetAttributes[s, Orderless]
S = Array[s, {dim, dim}];


I want to compute the trace of $S^{12}$:

tr = Tr[MatrixPower[S, 12]] // Expand;


but it is taking forever. The calculation seems slightly inefficient because I don't really need $S^{12}$ but only its trace, so I believe there must be some room for improvement. Any suggestion?

For those interested, I was trying to compute a gaussian integral over the space of symmetric matrices. I posted a question on mathoverflow here with the details.

• Do you really need the general symbolic expression for arbitrary input matrices? What's the motivating problem behind this? If you can simplify/specialise the structure of the matrix this probably speeds up the computation of MatrixPower. – Thies Heidecke Feb 5 '18 at 16:49
• @ThiesHeidecke Unfortunately, yes, I need the general symbolic expression. Once I have tr, I am going to integrate it with respect to its matrix elements, $s_{ij}$ (the integration step is very fast -- I have already made it as efficient as possible; the bottleneck is the trace). – AccidentalFourierTransform Feb 5 '18 at 16:50
• I think what actually takes so long is the typesetting of the resulting expression in the frontent, not the actual computation. I'll try to write up an answer to show what i mean. – Thies Heidecke Feb 5 '18 at 16:53
• There will be over 4.5 million terms in your result, which is why Expand takes so long. What kind of integration are you doing? – Carl Woll Feb 5 '18 at 21:51
• I would suggest to do the integration numerically by evaluating the traces for suitable quadrature points (quadrature points are symmetric matrices here)... – Henrik Schumacher Feb 6 '18 at 9:37

(ll = Eigenvalues@S;
trace = Total[ll^12];) // AbsoluteTiming
(*  {2.14433, Null}  *)


Update: Or this?

dim = 6;
power = 12;
SetAttributes[s, Orderless]
S = Array[s, {dim, dim}];

(spS = Most@CoefficientList[CharacteristicPolynomial[S, x], x];
sr = SymmetricReduction[Total[Array[L, dim]^power], Array[L, dim]];
trace = First@sr //. Flatten@
{SymmetricPolynomial[#, Array[L, dim]] -> (-1)^# sp[dim - # + 1] & /@ Range@dim,
c_. (SymmetricPolynomial[#, Array[L, dim]])^2 :>
c (sp[dim - # + 1])^2 & /@ Range@dim};) // AbsoluteTiming

(*  {1.92581, Null}  *)


Checked for power in Range[6]:

Tr@MatrixPower[S, power] - trace // Expand
(*  0  *)


The idea is to use the symmetric reduction to write the sum of the powers of the eigenvalues in terms of the coefficients of the characteristic polynomial. And then substitute.

Here's a nicer way, thanks to @CarlWoll, using the third argument to SymmetricReduction:

trace2 = First@
SymmetricReduction[Total[Array[L, 6]^12], Array[L, 6],
Reverse@Most@CoefficientList[CharacteristicPolynomial[S, x], x] {-1, 1, -1, 1, -1, 1}
]; // AbsoluteTiming

(*  {1.91906, Null}  *)

• Thanks! I am already running another piece of code, so I cannot try your suggestion for now. I'll report back later on. – AccidentalFourierTransform Feb 5 '18 at 16:56
• A problem with this approach is that it generates Root objects, while trace is actually a polynomial in the entries of S. Trying to Simplify it takes even more time than the naïve approach. Any idea how to fix this? – AccidentalFourierTransform Feb 5 '18 at 18:59
• @AccidentalFourierTransform It still takes gigabytes to expand... – Michael E2 Feb 5 '18 at 23:19
• There seems to be a problem with your code, as I get L's in trace[[1]]. Also, simpler is: First @ SymmetricReduction[Total[Array[L, 6]^12], Array[L, 6], Reverse@Most@ CoefficientList[CharacteristicPolynomial[S, x], x] {1, -1, 1, -1, 1, -1}];. Finally, depending on the integration being performed, it probably makes a lot of sense to use parameters in the third argument instead, so that integration can be peformed on 58 relatively simple terms instead of 4.5 million monomials. – Carl Woll Feb 5 '18 at 23:32
• @CarlWoll Thanks. I thought there was a way, but in the beta version the third argument didn't show up in the SymmetricReduction doc page! – Michael E2 Feb 6 '18 at 0:38

The straightforward way is actually quite fast

(S12 = MatrixPower[S, 12];) // AbsoluteTiming
(* {0.009604, Null} *)

(S12tr = Tr[S12];) // AbsoluteTiming
(* {0.061376, Null} *)


but has lots of elements

LeafCount[S12tr]
(* 115747 *)


which makes it slow to work with from there on. You could try to simplify the expression and integrate over that but probably the eigenvalue route @MichaelE2 showed is the better option.

• Yeah, it seems that the slow operation was actually Expand rather than Tr or MatrixPower. I wonder if its possible to speed up Expand though. – AccidentalFourierTransform Feb 5 '18 at 19:06
• @AccidentalFourierTransform I think this method has more promise than mine, actually. – Michael E2 Feb 5 '18 at 19:28

Here is another method to compute the trace of a matrix power. As already mentioned, the trace of a matrix power is equal to the power sum of the matrix's eigenvalues. The first key is to recognize that these power sums can be computed through a linear recurrence relation, where the required coefficients are the coefficients of the matrix's eigenpolynomial. The second key is that the initial conditions for this linear recurrence can be derived from the coefficients as well, by combining the Vieta and Newton-Girard formulae.

To be able to use the Newton-Girard formulae, I'll give an auxiliary routine for generating the required coefficient matrix from the polynomial coefficients:

ngMatrix[p_?VectorQ] := Module[{n = Length[p]},
SparseArray[{Band[{1, 1}] -> 1, {j_, k_} /; j > k :> p[[j - k]]}, {n, n}]]


I'll use a numerical example first for this demo:

n = 4;
mat = Array[Min, {n, n}]
{{1, 1, 1, 1}, {1, 2, 2, 2}, {1, 2, 3, 3}, {1, 2, 3, 4}}

rec =
(-1)^(n - 1) Reverse[Most[CoefficientList[CharacteristicPolynomial[mat, x], x]]]
{10, -15, 7, -1}

init = LinearSolve[ngMatrix[-rec], Range[n] rec]
{10, 70, 571, 4726}

With[{m = 50}, First[LinearRecurrence[rec, init, {m}]]] // AbsoluteTiming
{0.000694122, 8510938110502117856062697655362747468175263710}

Tr[MatrixPower[mat, 50]] // AbsoluteTiming
{0.00263071, 8510938110502117856062697655362747468175263710}


For the OP's symbolic example:

n = 6; SetAttributes[s, Orderless]; S = Array[s, {n, n}];

rec = Simplify[(-1)^(n - 1)
Reverse[Most[CoefficientList[CharacteristicPolynomial[S, x], x]]]];
init = Simplify[LinearSolve[ngMatrix[-rec], Range[n] rec]];


and then evaluate

First[LinearRecurrence[rec, init, {12}]]


However, I've found the timings to be a bit inconsistent; sometimes MatrixPower[] is faster, and sometimes LinearRecurrence[] is faster (not counting the time needed to set up rec and init).

• This is very interesting, never thought about it. I'll play around with your code. Thanks! – AccidentalFourierTransform Mar 24 '18 at 13:59
• My original code had a bug; please try the current version. – J. M. will be back soon Mar 25 '18 at 16:43