Accelerating sum over permutations of matrix elements

Question

I am trying to short-cut the use of a CoefficientArrays call by manually calculating the resulting matrix of coefficients myself (this avoids using symbolic arrays and is therefore quicker). The calculation that I'd like to shortcut looks like $R\cdot X \cdot R^\mathrm{T} - X=0$ (where $X$ is symmetric): as an example, given a matrix $$ X=\left(\begin{array}{c,c} x[1,1] & x[1,2]\\ x[1,2] & x[2,2] \end{array}\right) $$ and automorphism matrix $$ R=\left(\begin{array}{c,c} a & b\\ c & d \end{array}\right) $$ the CoefficientArrays entry looks like $$ \left(\begin{array}{c,c,c} a^2-1 & 2 a b & b^2 \\ a c & b c + a d - 1 & b d \\ c^2 & 2 c d & d^2 - 1 \end{array}\right) $$

The code I'm using to calculate this in Mathematica looks like

coeff[automorph_, xarray_] := Module[{subfunc, indexfn, permind},
   (*Get the unique matrix element IDs*)
   indexfn = (Level[#, {1}] & /@ (DeleteDuplicates@Flatten@xarray));
   (*Find all permutations of these elements*)
   permind = Map[Permutations[#] &, indexfn];
   (*Define a function that works for each row of the coefficient array matrix*)
   subfunc[indices_] := 
    Total /@ 
     Map[Map[(automorph[[indices[[1]], #[[1]]]] automorph[[indices[[2]], #[[2]]]]
           - DiscreteDelta[indices[[1]] - #[[1]], indices[[2]] - #[[2]]]) &, #] &, permind];
   (*Map over rows*)
   subfunc /@ indexfn
   ];

and is called like

coeff[{{a, b}, {c, d}}, {{x2[1, 1], x2[1, 2]}, {x2[1, 2], x2[2, 2]}}]

This works for small examples as here, but I am attempting to scale this up to higher-dimensional tensor $R$ and $X$ (e.g. $X_{ijk}R^i_{i'}R^j_{j'}R^k_{k'}-X_{i'j'k'}$) and it does not scale particularly well. Could anyone help with how to speed up this function (or come up with another way to calculate the result of CoefficientArrays more rapidly)?

Answer 1

Intro. Given a matrix mat, the matrices of the 2nd, 3rd, 4th, ... tensor power of mat are

KroneckerProduct[mat,mat]
KroneckerProduct[mat,mat,mat]
KroneckerProduct[mat,mat,mat,mat]
...

See tensorPower below.

Closely related to this, the powers of mat also act on the symmetric tensor product spaces. To obtain a matrix for this, we just have to take tensorPower and pre- and post-multiply with appropriate matrices symL and symR.

Code.

(* rectangular, sparse matrices that satisfy symL.symR = 1 *)
symR[dim_,power_]:=symR[dim,power]=With[{X=Map[Sort,Tuples[Range[1,dim],power]]},
    With[{Y=MapIndexed[(#1->#2[[1]])&,DeleteDuplicates[X]]},
        SparseArray[MapIndexed[({#2[[1]],#1/.Y}->1)&,X]]]];
symL[dim_,power_]:=symL[dim,power]=Transpose[symR[dim,power]]//#/Map[Total,#]&;

(* tensor power *)
tensorPower[mat_?MatrixQ,power_]:=If[power==1,mat,KroneckerProduct@@ConstantArray[mat,power]];
symmetricTensorPower[mat_?MatrixQ,power_]:=With[{dims=Dimensions[mat]},
    symL[dims[[1]],power].tensorPower[mat,power].symR[dims[[2]],power]];

Examples.

symmetricTensorPower[{{a,b},{c,d}},2]

gives

To get the matrix OP has given, one just has to subtract the identity matrix, which is easy enough.

Similarly

symmetricTensorPower[{{a,b},{c,d}},3]

gives

Comments.

• There is some ambiguity in the choice of symL and symR, corresponding to the choice of a basis on the symmetric tensor product spaces. Different choices give similar matrices. The choice I made is consistent with the example OP has provided.
• I have not done any timing, since I do not know which values of dim and power are relevant for OP. But the generation of the sparse matrices symL and symR should be negligible, and apart from that, the code relies on a single KroneckerProduct call.

---

Memory considerations. The dimension of the tensor product and the symmetric tensor products are, respectively,

tdim[dim_,power_]:=dim^power;
sdim[dim_,power_]:=Binomial[dim+power-1,power];

This means for example that symR is a sparse matrix of size tdim x sdim.

The code given in this answer first constructs a matrix of size tdim x tdim using KroneckerProduct, and then reduces it to size sdim x sdim by multiplying with symL and symR. This means that if tdim is much bigger than sdim, then this code will produce an intermediate object that is much bigger than the final object.

Here are two examples:

tdim[15,3]/sdim[15,3]
(* about 5 *)

tdim[15,6]/sdim[15,6]
(* about 300 *)

In the second case, dim == 15 and power == 6, the overhead is very big.

Comment. The case dim == 15 and power == 6 mentioned by OP will be challenging regardless of the method. In fact

sdim[15,6]
(* 38760 *)

This means that in this case, the final matrix will have size 38760 x 38760 with entries that are order 6 polynomials in various symbols. That will take up a lot (!) of memory. Even if the entries of this matrix were simply floating point numbers, this matrix would require about 12 gigabytes. In the case of polynomials of order 6, it will be much more.