# Eigenvalues and eigenvectors of tensors

I have a linear equation where the unknown is a matrix:

$$\sum_{pq} T_{ijpq}X_{pq} = \lambda X_{ij}$$

This is an eigenvalue problem for the tensor $T_{ijpq}$. Mathematica Eigenvalues and Eigenvectors allows me to compute eigenvalues and eigenvectors of a matrix, which means I have to flatten $T_{ijpq}$ before I can apply these functions. But then the eigenvectors will be flattened too, and it is not obvious to me how I can recover the $X_{pq}$.

I can think of a contrived way of doing this, but I feel it should not be that complicated in Mathematica.

• "I can think of a contrived way of doing this" - would you mind posting the "contrived" method as a self-answer? Then people can propose improvements. – J. M.'s technical difficulties Oct 15 '17 at 13:55
• @J.M. The idea is to transform the pair of indices $i,j$ to a single index $a$, and the pair of indices $p,q$ to a single index $b$. In my application, $i$ and $p$ traverse the same range, from $1$ to $M$, while $j$ and $q$ traverse a range from $1$ to $N$. Then $a,b$ traverses a ranges from 1 to $M^2$ and from 1 to $N^2$, respectively. But it gets messy when one tries to do the index mapping in any specific way. I am trying do it like this now. Then to read the eigenvector one needs inverse index mappings. – becko Oct 15 '17 at 14:05
• @J.M. I think I found a good way to do it! – becko Oct 15 '17 at 14:20
• I had some typos in my previous comments. It should say: $a,b$ traverse the range from 1 to $MN$. – becko Oct 16 '17 at 3:35
• Does 'flattening' leave the eigenvalues intact ? I can't mathematically understand why that's the case. Can someone enlighten me regarding this perhaps. – Lelouch Apr 4 at 18:28

I think I found a simple way.

mat = Flatten[T, {{1,2},{3,4}}]
eig = Eigenvalues[mat]
eiv = Eigenvectors[mat]


Then reshape the eigenvectors using ArrayReshape:

ArrayReshape[#, {m,n}] & /@ eiv


where $m,n$ are the ranges of the indices $p,q$.

I have only two improvements.

First, use Eigensystem instead of Eigenvalues and Eigenvectors. Otherwise, you compute the eigenvalues twice.

Second, instead of Mapping ArrayReshape onto eiv, you can apply it directly, which also improves performance. (This is less an issue because Eigensystem consumes most of the ressources.

n = 10;
m = 12;
T = RandomReal[{-1, 1}, {m, n, m, n}];
mat = Flatten[T, {{1, 2}, {3, 4}}];
{eig, eiv} = Eigensystem[mat];
eiv = ArrayReshape[eiv, {m n, m, n}];


And here is a test

Norm[
Table[
Norm[
Sum[T[[All, All, k, l]] eiv[[i, k, l]], {k, 1, m}, {l, 1, n}] -
eig[[i]] eiv[[i]],
Infinity],
{i, 1, m n}],
Infinity]

(* 5.12285*10^-14 *)