How to generate computationally a harmonic tensor?

Question

A $n$-th-rank tensor $T_n$ based on $\mathbb{R}^3$ is referred here to as harmonic if it is symmetric in all indices and the linear map of the identity on vectors $I_2$ vanished, denoted as $T_n[I_2] = 0_{n-2}$, where $O_{n-2}$ denotes the $(n-2)$-rank zero tensor. In index notation this mean, for example, for $n=4$ $$ T_{ijkl} = T_{jikl} = T_{ijlk} = \dots \quad \land \quad \sum_{k,l=1}^3 T_{ijkl} \delta_{kl} = \sum_{k=1}^3 T_{ijkk} = 0 \quad \forall i,j \ . $$ I use the term harmonic for these tensors due to their connection to homogeneous polynomials (which are connected to symmetric $n$-th-rank tensors) with vanishing Laplacian. How would you generate a harmonic tensor of rank $n$ computationally?

I know, I can generate a symmetrized array with symbols, solve the linear map condition $T_n[I_2] = 0_{n-2}$, reinsert the solution into $T_n$ and then replace the symbols with some, e.g., random numbers. See upcoming code. But I would like to avoid the symbolic solution since I only want to generate fast some instance of harmonic tensors for high $n$. Symmetrize and RandomReal allow to generate way faster symmetric tensors, but I can not get my head around on how to enforce $T_n[I_2] = 0_{n-2}$ computationally. Any ideas?

(*Generate a tensor of rank n based on symbols s and with index \
symmetries is*)

genT[n_, s_, is_: {}] := 
  Normal@SymmetrizedArray[
    pos_ :> Subscript[s, StringJoin[ToString /@ pos]], 
    ConstantArray[3, n], is];
(*Identity on vectors*)
I2 = IdentityMatrix@3;
(*Linear map*)

lm[A_, B_] := 
  ArrayReshape[A, 
    Join[Dimensions[A][[1 ;; -(1 + TensorRank@B)]], {Times @@ 
       Dimensions@B}]].Flatten[B];
(*Generate harmonic tensor of rank n with symbols s*)

genH[n_, s_] := # /. First@Solve[lm[#, I2] == 0*lm[#, I2]] &@
   genT[n, s, Symmetric@Range@n];
(*Generate a harmonic tensor of rank n*)
genHR[n_] := Module[
   {s},
   # /. Table[
       Variables[#][[i]] -> RandomReal[], {i, Length@Variables@#}] &@
    genH[n, s]
   ];
(*Generate instance and check properties*)

T = genHR[6]; // AbsoluteTiming
TensorSymmetry@T
lm[T, I2]
vars = Variables@T
(*Generate a symmetric tensor*)

genSR[n_] := 
  Normal@Symmetrize[RandomReal[{-1, 1}, ConstantArray[3, n]]];
TensorSymmetry@genSR[6] // AbsoluteTiming

Answer 1

It may help to solve the symbolic equation only once and to generate a (linear) orthoprojektor onto the solution set.

Some preparations:

n = 12;
pat = DeleteDuplicates[Sort /@ Tuples[Range[3], n]];
dim = Length[pat];
totensor = SymmetrizedArray[Thread[pat -> #], ConstantArray[3, n], Symmetric[Range[n]]] &;
dot[A_] := Flatten[A, Join[Table[{i}, {i, 1, n - 2}], {{n - 1, n}}]].Flatten[ IdentityMatrix[3]]

Creating the linear system for the nonzero values of the symmetric tensor:

vals = Array[v, dim];
T = totensor[vals]; // AbsoluteTiming // First
b = Flatten[dot[T]]; // AbsoluteTiming // First
A = Normal@N@CoefficientArrays[b, vals][[2]]; // AbsoluteTiming // First

0.001915

0.600078

0.030182

Computing the orthoprojektor onto the null space of A and checking its conditions:

pr = IdentityMatrix[dim] - PseudoInverse[A].A; //AbsoluteTiming // First
Max[Abs[pr\[Transpose] - pr]]
Max[Abs[pr.pr - pr]]
Max[Abs[A.pr]]

0.322819

1.16864*10^-14

1.14422*10^-14

2.67957*10^-14

Now creating one million value lists for the symmetric tensors at once.

rand = RandomReal[{-1, 1}, {1000000, dim}].pr; // AbsoluteTiming // First

1.24263

The conversion to the actual StructuredArrays with totensor is however still somewhat slow. The reason is of course that these tensors are dense.