# How to generate computationally a harmonic tensor?

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; // 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 // AbsoluteTiming


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, 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]


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][]; // 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.

• Hey Henrik! Cool approach using the orthoprojector and the kernel of A! Is there any standard algorithm for this kind of things well known in the math literature? – Mauricio Fernández Aug 16 '18 at 12:36
• Hey Mauricio! What I used is one of the many useful properties of the Moore-Penrose pseudoinverse. Alternatively, you can obtain the projector by using NullSpace: pr = Transpose[#].# &[NullSpace[A]] -- provided that A is a floating point matrix (otherwise, NullSpace[A] might not be orthonormalized). However, using Pseudoinverse turned out to be faster. – Henrik Schumacher Aug 16 '18 at 13:28