Code for Tensor decomposition

Question

I would like to know if there is a package or some MMA code to perform tensor decomposition as e.g. defined in a paper by Robeva "Orthogonal decomposition of symmetric tensors" or some other papers in this field.

Example

Given is the symmetric matrix

$$\textrm{M}=\begin{bmatrix} 0 &0&0&1/2\\ 0&0&-1/2&0\\ 0&-1/2&0&0\\ 1/2&0&0&0 \end{bmatrix}.$$ It is possible to decompose $\textrm{M}$ by tensor powers of orthonormal vectors

$$\textrm{M}=\sum\limits_{i=1}^4 \lambda_i {\bf q}_i^{\otimes 2}\quad\text{with}$$ $${\bf q}_1=\left(\frac{-1}{\sqrt{2}},0,0,\frac{1}{\sqrt{2}}\right), {\bf q}_2=\left(0,\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\right),\\ {\bf q}_3=\left(\frac{1}{\sqrt{2}},0,0,\frac{1}{\sqrt{2}}\right), {\bf q}_4=\left(0,\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}},0\right),\\ \text{where } \langle {\bf q}_i, {\bf q}_j \rangle= \delta_{ij}\text{ and } \lambda_{1,2}=\frac{-1}{2},\lambda_{3,4}=\frac{1}{2}. $$

m={{0,0,0,1/2},{0,0,-1/2,0},{0,-1/2,0,0},{1/2,0,0,0}};
q[1]={-(1/Sqrt[2]),0,0,1/Sqrt[2]};
q[2]={0,1/Sqrt[2],1/Sqrt[2],0};
q[3]={1/Sqrt[2],0,0,1/Sqrt[2]};
q[4]={0,-(1/Sqrt[2]),1/Sqrt[2],0};
\[Lambda][1]=-1/2;
\[Lambda][2]=-1/2;
\[Lambda][3]=1/2;
\[Lambda][4]=1/2;
Sum[\[Lambda][i]TensorProduct[q[i],q[i]],{i,1,4}]==m
(* True *)

The problem

Now is given a concrete symmetric $(n\times n\times n)$-tensor $\mathrm{T}$, and assuming that such a decomposition exists, I want, similarly to the previous example, decompose $\mathrm{T}$ into tensor powers of orthonormal vectors ${\bf q}_i$ with $\langle {\bf q}_i, {\bf q}_j \rangle= \delta_{ij}$

$$\textrm{T}=\sum\limits_{i=1}^p \lambda_i {\bf q}_i^{\otimes 3}$$

with $p\le n$ due to orthogonality. In MMA I would perform

T=Sum[\[Lambda][i]TensorProduct[q[i],q[i],q[i]],{i,1,p}]

Power iteration to find eigenvectors I already tested and it fails presumably as there is no dominant eigenvalue. Therefore I am looking for a code or package that helps to find ${\bf q}_i$ and $\lambda_i$ for symmetric tensors of order 3 and 4.

---

Reference:

Elina Robeva: Orthogonal decomposition of symmetric tensors, SIAM J. Matrix Anal. Appl. 37, p.86, 2016

https://doi.org/10.1137/140989340

https://arxiv.org/abs/1409.6685

Answer 1

If the representation

$$ \mathbf{T}=\sum_{i=1\ldots r} \lambda_i \mathbf{q}_i \otimes \mathbf{q}_i \otimes \mathbf{q}_i, \qquad \mathbf{q}_i\cdot \mathbf{q}_j =\delta_{ij} $$ is possible, you can calculate a double scalar contraction $$ \mathbf{M}=\mathbf{T}\cdot\cdot\mathbf{T}=\sum_{i=1\ldots r} \lambda_i^2 \mathbf{q}_i \otimes \mathbf{q}_i $$ (or $M_{ij}=T_{ikl}T_{klj}$ or any other index permutation in $T_{...}$ due to the symmetry) and then calculate the eigensystem of $\mathbf{M}$, which should give you the $\lambda_i^2$ and the $\mathbf{q}_i$.

Then you have to figure out the original signs of $\lambda_i$, which are lost due to squaring. Some trial and error may lead to the desired result. Additionally, if you have eigenvalue multiplicities, you need to figure out the variable orientation of the eigenvectors in the corresponding eigenspace.

Answer 2

So odeco is short for orthogonal decomposable.

Mathematica has a lot integrated or built in for that purposes.

The best start is to make use of Orthogonalize because You have a matrix that is not orthogonal.

The test for that is OrthogonalMatrixQ.

OrthogonalMatrixQ@m (False)

for the given $m$.

After orthogonalization You have the full arsenal ready at hand.

The full arsenal is given in OrthogonalMatrixQ in the section Applications. For example the JordanDecomposition has access to eigenvalues and eigenvectors.

For example QRDecomposition shows how to use matrix decomposition of all matrices for interpolation as an application. So Mathematica's definition of application is not so standard.

The problem with Your given matrix is that is poses not so general attributes. It is the same as the transposed matrix of it. That stems from being antidiagonal.

ConjugateTranspose[m].m

(* {{1/4, 0, 0, 0}, {0, 1/4, 0, 0}, {0, 0, 1/4, 0}, {0, 0, 0, 1/4}} *)

shows that there is a factor of $4$ from keeping it being diagonal after multiplication of the transposed matrix from the left.

Check primarily:

Det@m

(* 1/16 *)

if the matrix has potential for such a linear power development.

Eigenvalues@m

(* {-(1/2), -(1/2), 1/2, 1/2} *)

gives the eigenvalues of the matrix. This confirms the development has chances to exist.

Eigenvectors@m

(* {{-1, 0, 0, 1}, {0, 1, 1, 0}, {1, 0, 0, 1}, {0, -1, 1, 0}} *)

offers You the basis of the two subspaces that the $m$ decomposes into. Mind these eigenvectors are not normalized! Up to that normalization factor which is $\frac{1}{\sqrt{2}}$ the basis is already in the order You have written down. That is universally the case.

Name the eigenvalues and eigenvectors in the manner You already did and verify the formula. That is plainly all that is to do.

$\lambda=Eigenvalues@m$


$q=Normalize@Eigenvector@m$

with Normalize.

Mind to make use of the Part built-in instead of $[..]$ in mathematical notation. Then You are done. Use have already done so much of the work.

A simple hind. In Mathematica fundamental built-ins are implemented only once. So the built in decomposition make use of the built-ins I have named.

You have given TensorProduct and $\otimes$, \otimes already.

Mind that orthogonality is a special case of unitarity reduced for matrices with real not anymore with complex elements in the matrix slots. So everything can be done with complex numbers as well.

T=Sum lambda[[i]] MatrixPower[TensorProduct[q[[i]],q[[i]]],3],{i,1,p}]

With TensorProduct You make matrices with the required component. This is not the Power replacement for matrices. This does MatrixPower.

So the given symbol is decomposed into TensorProduct and MatrixPower!

I suggest to write it better as:

$(q_{i}\otimes q_{i})^{n}$ written as (q_{i}\otimes q_{i})^{n}

to avoid confusion.

From this question higher order tensor of variable rank a notion that higher order tensors can be dealt with in Wolfram Language and therefore this works. The guy from the answer uses the built-in Array with is again an algorithmic approach to create, generate a list. Stephen Wolfram's first book did serve the same purposes.

So what to do is stick to the advice given in my answer. The question imposes by the references confirmed restriction on the matrices. This works for rotations on the unit circle. The MatrixPower is that used for rotation multiplications as well. So this are roots one the unity matrix. That are two constraints in one for the $nxnxn$-matrices targeted at. The are expected to multiply to the unit matrix in the $nxnxn$-matrices space.

The eigenvalue decomposition is valid in the $nxnxn$-matrices space as well if the is handled as a constraint. A simple generative path to understand that is extent the well understandable $3x3$-matrices spaces to the next dimension by using only 1 dimension different from zero. Then in the subspace of $3x3$-matrices the above sketches methods are valid again. The more general cases can be understood by making the higher spaces separable as stems from the eigenvalue - eigenvector - method. So make the 4th dimension of the matrices diagonal or antidiagonal, as in the given example matrix.

As that remain visually unclear and incomprehensible. But this works. The concepts can be drawn from the much richer in literature presented visualization of the unit-sphere. And this whole deals about rotation on the unit sphere. Some author extent than to mirrorings and shearings and so one. This is a well known set of established mathematical procedures.

I recommend the article of Barabara Edwards and Michael B. Ward about "The role of Mathematical definitions in Mathematics and in Undergraduate Mathematics Courses". They name as necessary features the criterion of hierarchy, criterion of existences, criterion of equivalence, criterion of acclimatization. This can be found by using google.de

You need to make use of the criterion of acclimatization to accept that MatrixPower is dealing with the $nxnxn$-matrices brilliantly. This might not be the only method doing this but it is only that poses the attributes of having a unit matrice and that has divisors in the sense of rotations. How these may rotate is what ever it will be. Important is there are some that behave like the $nxn$-matrices subspace matrices. That can be understood for $n=3$.

The method addressed is the Laplace expansion. That is well defined for $nxnxn$-matrices. And implemented in Wolfram language. Nothing changes other than the $nxnxn$-matrices. Since the are more eigenvalues to be taken into account with more eigenvectors and corresponding eigenspaces. Still the method does not change. This is only summing over more linearly envolved eigenvalues.

The sum here is simpler because this is a Linear combination of matrices in matrices spaces representing eigenspaces but not in the vectorspace but in the matrice(-sub-)spaces. Look at this to get a better idea in 3D and general what matrix rotations are: matrix+rotations.

This shows the roots of unity: Root of unity. Of interest in the understanding is section about periodicity for this question. Some ideas from the following sections are nice for understanding too.

What is used in Wolfram language is mainly the property Listable for generating and treating the next step $nx$ that is needed. This corresponds to the expansion of Laplace. Whether a built-in has this can be checked or explored by Attributes.

Here is some example and insight on how this looks in 4D: Rotations in 4-dimensional Euclidean space. But formally there is no change. So this involves already many different cases.