Speed up the computation of the trace of a matrix

How about this?:

(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}  *)

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.

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).