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 - 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[2, 1], 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)?