What is multilinear principal components analysis? I've gotten a lot of usage out of principal component analysis, and after recently learning the basics of performing canonical polyadic decomposition I was intrigued to learn that there exists a multilinear principal components analysis (MPCA).
The trouble is, the wikipedia page doesn't really clarify what it is mathematically. What it says is:

MPCA computes a set of orthonormal matrices associated with each mode of the data tensor which are analogous to the orthonormal row and column space of a matrix computed by the matrix SVD. This transformation aims to capture as high a variance as possible, accounting for as much of the variability in the data associated with each data tensor mode(axis).

On the face of it this usage of orthogonal matrices along modes of the tensor sounds similar to higher-order singular value decomposition which is a special case of Tucker decomposition.
What is multilinear principal components analysis?

I just tracked down this lengthy paper which might hold the answer.

Here are some exerpts from Lu, Plataniotis and Venetsanopoulos 2006.

Let $\{ \mathcal{A}_m, m=1,\cdots, M \}$ be a set of $M$ tensor samples in $\mathbb{R}^{I_1} \otimes \mathbb{R}^{I_2} \otimes \cdots \mathbb{R}^{I_N}$. The total scatter of these tensors is defined as: $\Psi_{\mathcal{A}} = \sum_{m=1}^M \vert| \mathcal{A}_m - \bar{\mathcal{A}} \vert|_F^2$, where $\bar{\mathcal{A}}$ is the mean tensor calculated as $\bar{\mathcal{A}} = \frac{1}{M} \sum_{m=1}^M \mathcal{A}_m$. The $n$-mode total scatter matrix of these samples is then defined as: $\mathbf{S}_{T_{\mathcal{A}}}^{(n)} = \sum_{m=1}^M (\mathbf{A}_{m(n)} - \bar{\mathbf{A}_{(n)}}) (\mathbf{A}_{m(n)} - \bar{\mathbf{A}_{(n)}})^T$, where $\mathbf{A}_{m(n)}$ is the $n$-node unfolded matrix of $\mathcal{A}_m$. The statement above leads to the following formal definition of the problem to be solved:


A set of $M$ tensor objects $\{ \chi_1, \chi_2, \cdots, \chi_M \}$ is available for training. Each tensor object $\chi_m \in \mathbb{R}^{I_1 \times I_2 \times \cdots \times I_N}$ assumes values in a tensor space $\mathbb{R}^{I_1} \otimes \mathbb{R}^{I_2} \otimes \cdots \mathbb{R}^{I_N}$, where $I_n$ is the $n$-mode dimension of the tensor. The MPCA objective is to define a multilinear transformation $\{ \tilde{\mathbf{U}}^{(n)} \in \mathbb{R}^{I_n \times P_n, n = 1, \cdots, N} \}$ that maps the original tensor space $\mathbb{R}^{I_1} \otimes \mathbb{R}^{I_2} \otimes \cdots \mathbb{R}^{I_N}$ into a tensor subspace $\mathbb{R}^{P_1} \otimes \mathbb{R}^{P_2} \otimes \cdots \mathbb{R}^{P_N}$ (with $P_n < I_n$, for $n=1,\cdots,N$): $\mathcal{Y}_m = \chi_m \times _1\tilde{\mathbf{U}}^{(1)^T} \times _2\tilde{\mathbf{U}}^{(2)^T} \cdots \times _N \tilde{\mathbf{U}}^{(N)^T}$, $m=1, \cdots, M$, such that $\{ \mathcal{Y}_m \in \mathbb{R}^{P_1} \otimes \mathbb{R}^{P_2} \otimes \cdots \mathbb{R}^{P_N}, m=1, \cdots, M \}$ captures most of the variation observed in the original tensor objects, assuming that these variation are measured by the total scatter.


In other words, the MPCA objective is the determination of the $N$ projection matrices $\{ \tilde{\mathbf{U}}^{(n)} \in \mathbb{R}^{I_n \times P_n}, n=1,\cdots,N \}$ that maximizes the total tensor scatter $\Psi_{\mathcal{Y}}$: $$\{ \tilde{\mathbf{U}}^{(n)}, n=1, \cdots, N \} = \arg\max_{\tilde{\mathbf{U}^{(1)}}, \tilde{\mathbf{U}^{(2)}}, \cdots, \tilde{\mathbf{U}^{(N)}}} \Psi_{\mathcal{Y}}.$$

I guess this technically answers the question, but I am going to have to ruminate on this some more.
