# Accelerating sum over permutations of matrix elements

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

• Hi. In LaTeX you have a $3 \times 3$ matrix, but your code produces a $4 \times 4$ matrix, which I find confusing. Also, you write $RXR-X=0$, are you missing a transpose somewhere, and why do you write $=0$? Also, can you explain in words what the meaning of the $(i,j)$-th entry is supposed to be? Commented Aug 13, 2022 at 5:41
• Apologies for a couple of typos: my understanding is that CoefficientArrays generates the coefficients of the parameters $x$ in the linear equation $g(x)=0$: I guess we could drop the 0 if that is implied? The resulting matrix $A$ is expected to mean that $A.\mathrm{vec}(X)=0$ (where $\mathrm{vec}(X)$ is Flatten@X in Mathematica-language) has the same set of solutions as the original $R\cdot X \cdot R^\mathrm{T}-X=0$, as you'd get from CoefficientArrays. Does that make sense? Commented Aug 13, 2022 at 9:47

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

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

• This is helpful, thank you - in my timings on small examples, with dim=15 and power=2 my method took 0.19s and 1.6MB memory; this approach took 0.11s and 8.7MB memory. With dim=15 and power=3 my method took 19s and 47MB; this method took 8.6s and 1.6GB. So there's definitely a memory cost: I'm hoping to get to dim=15 and power=6 which I think my workstation will struggle with. Do you have thoughts on how to reduce the memory cost? As symL and symR are very sparse it feels like it should be possible to only calculate the entries which will eventually be required? Commented Aug 15, 2022 at 19:23
• Please see the edit. Your dim == 15 and power == 6 are much bigger than I anticipated. (Sure you need that? In symbolic form?) The code I given here will not directly work. One could try to split things into pieces and avoid large intermediate results. For example, use the power == 3 result to get the power == 6 result: Calculate aux = symmetricTensorPower[mat,3], then aux2 = KroneckerProduct[aux,aux], then auxL.aux2.auxR where auxL is "defined by" symL[15,6] = auxL.symL[15,3] or so. I have not tried, and even this will generate huge intermediate results. Commented Aug 15, 2022 at 20:13
• Thanks for the suggestions, I will try them out - the R matrix I'm using is not symbolic, but does need to be to high precision (working precision around 30): I'm finding the nullspace of the result of this calculation, and the precision is needed to avoid errors in that step. If a non-machine-precision method is out there to help with memory that would also be great (I tried Compile but that doesn't like non-machine-precision) Commented Aug 16, 2022 at 13:27
• One can calculate the nullspace without constructing the matrix.I think you mean the nullspace of symmetricTensorPower[mat,power] minus identity, which is the eigenspace for the eigenvalue $1$ of symmetricTensorPower[mat,power]. Suppose mat is diagonalizable with eigenvectors $e_i$ and eigenvalues $\lambda_i$ (possibly complex). Then a basis of eigenvectors of symmetricTensorPower[mat,power] is given by symmetrizing $e_{i_1}\otimes\cdots \otimes e_{i_{\text{power}}}$ with $i_1\leq i_2\leq \ldots\leq i_{\text{power}}$ which has eigenvalue $\lambda_{i_1}\cdots\lambda_{i_{\text{power}}}$. Commented Aug 16, 2022 at 13:56
• Thank you, I will look into symmetrizing the tensor product of the eigenvectors of the underlying mat - accepting the answer given as solving the problem as stated! Commented Aug 17, 2022 at 13:56