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_: {}] := 
    pos_ :> Subscript[s, StringJoin[ToString /@ pos]], 
    ConstantArray[3, n], is];
(*Identity on vectors*)
I2 = IdentityMatrix@3;
(*Linear map*)

lm[A_, B_] := 
    Join[Dimensions[A][[1 ;; -(1 + TensorRank@B)]], {Times @@ 
(*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[
   # /. Table[
       Variables[#][[i]] -> RandomReal[], {i, Length@Variables@#}] &@
    genH[n, s]
(*Generate instance and check properties*)

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

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

1 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




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





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

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


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

  • $\begingroup$ 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? $\endgroup$ Aug 16, 2018 at 12:36
  • 1
    $\begingroup$ 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. $\endgroup$ Aug 16, 2018 at 13:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.