# 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. Oct 15, 2017 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. Oct 15, 2017 at 14:05
• @J.M. I think I found a good way to do it! Oct 15, 2017 at 14:20
• I had some typos in my previous comments. It should say: $a,b$ traverse the range from 1 to $MN$. Oct 16, 2017 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. Apr 4, 2020 at 18:28

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 *)


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$.