Tensor-Normal Distributions and Mode Covariances

Why a Gaussian Model for a Tensor?

A multiway observation records several indices at once: for example, measurements on people across occasions and outcomes, or an image indexed by row, column, and color channel. Treating all entries as one long vector gives an ordinary multivariate Gaussian model, but its covariance matrix can be large and difficult to interpret. A tensor-normal model imposes a structured alternative: dependence along each mode is described by a separate covariance matrix, and the full covariance is their Kronecker product.

This structure is useful when the modes have meaningful roles and a covariance factor for each mode is a scientifically interpretable description. The representation also makes the convention for ordering the entries essential. We therefore state that convention before giving the distribution.

Definition and Vectorization Convention

Let $X$ be a real $K$-way tensor with dimensions $d_1,\ldots,d_K$, and let $M$ be a fixed tensor of the same dimensions. Write $d=\prod_{k=1}^K d_k$. Throughout, $\operatorname{vec}(X)$ stacks the entries with index 1 varying fastest, then index 2, and so on. Thus for a matrix, $\operatorname{vec}(X)$ stacks its columns.

For each mode $k$, let $\Sigma_k$ be a $d_k$-by-$d_k$ symmetric positive-definite matrix. Positive definiteness ensures that the resulting distribution has a nonsingular density with respect to Lebesgue measure on the $d$ tensor entries.

Definition: The tensor-normal distribution with mean tensor $M$ and mode covariance matrices $\Sigma_1,\ldots,\Sigma_K$, written $X\sim TN(M;\Sigma_1,\ldots,\Sigma_K)$, is the multivariate normal distribution on the tensor entries specified by $$ \operatorname{vec}(X)\sim N\!\left(\operatorname{vec}(M),\ \Sigma_K\otimes\cdots\otimes\Sigma_1\right). $$ The factors are ordered in reverse mode order because index 1 varies fastest in the stated vectorization convention.

Equivalently, this model can be constructed by starting with a tensor $Z$ whose $d$ entries are independent $N(0,1)$ variables, choosing matrices $A_k$ such that $A_kA_k^{\mathsf T}=\Sigma_k$, and setting $X=M+Z\times_1 A_1\times_2\cdots\times_K A_K$. Here $\times_k$ denotes multiplication along mode $k$. This construction gives a direct interpretation: each mode transformation introduces covariance along that mode.

Covariance Between Tensor Entries

The mode covariance matrices determine the covariance of any pair of entries by multiplying one entry from each mode covariance. In particular, for index tuples $i=(i_1,\ldots,i_K)$ and $j=(j_1,\ldots,j_K)$,

\operatorname{Cov}(X_{i_1\cdots i_K},X_{j_1\cdots j_K}) =\prod_{k=1}^K (\Sigma_k)_{i_kj_k}.

This is the entrywise meaning of the separable covariance structure. It does not say that all tensor entries are independent: off-diagonal elements in any mode covariance can induce dependence between entries that differ along that mode.

Worked Example: Reading a Two-Mode Covariance Matrix

Let $X$ be a $2$-by-$2$ random matrix with mean zero and mode covariance matrices $\Sigma_1=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}$ and $\Sigma_2=\begin{pmatrix}2&0\\0&3\end{pmatrix}$. Find the covariance matrix of $\operatorname{vec}(X)$ and interpret two of its entries.

Since vectorization stacks columns, $\operatorname{vec}(X)=(X_{11},X_{21},X_{12},X_{22})^{\mathsf T}$. The covariance matrix is $\Sigma_2\otimes\Sigma_1$, so

\operatorname{Cov}(\operatorname{vec}(X)) = \begin{pmatrix} 2&1&0&0\\ 1&2&0&0\\ 0&0&3&1.5\\ 0&0&1.5&3 \end{pmatrix}.

For example, $\operatorname{Var}(X_{11})=2$, and $\operatorname{Cov}(X_{11},X_{21})=1$. The latter also follows directly from the entry formula: the two entries share mode-2 index 1, contributing $(\Sigma_2)_{11}=2$, and differ across mode 1, contributing $(\Sigma_1)_{12}=0.5$; their product is $2(0.5)=1$. In contrast, $\operatorname{Cov}(X_{11},X_{12})=0$ because the relevant mode-2 covariance entry is zero.

Worked Example: Covariance in a Three-Way Tensor

Suppose $X$ has dimensions $2$-by-$2$-by-$2$, with mode covariance matrices $\Sigma_1=\begin{pmatrix}1&0.2\\0.2&1\end{pmatrix}$, $\Sigma_2=\begin{pmatrix}2&0.5\\0.5&1\end{pmatrix}$, and $\Sigma_3=\begin{pmatrix}4&1\\1&3\end{pmatrix}$. Compute the covariance between $X_{121}$ and $X_{222}$, and their variances.

The index pairs are $(1,2)$ in mode 1, $(2,2)$ in mode 2, and $(1,2)$ in mode 3. Therefore

\operatorname{Cov}(X_{121},X_{222}) =(\Sigma_1)_{12}(\Sigma_2)_{22}(\Sigma_3)_{12} =(0.2)(1)(1)=0.2.

For an entry's variance, use the diagonal element from each mode matrix: $\operatorname{Var}(X_{121})=(\Sigma_1)_{11}(\Sigma_2)_{22}(\Sigma_3)_{11}=(1)(1)(4)=4$. Similarly, $\operatorname{Var}(X_{222})=(1)(1)(3)=3$. The covariance is not itself a correlation; the correlation here is $0.2/\sqrt{4\cdot3}=0.2/\sqrt{12}$.

Deriving the Vectorized Normal Representation

The vectorized representation follows from the action of mode products under the chosen ordering. If $A_k$ is a $d_k$-by-$d_k$ matrix, applying $A_k$ along mode $k$ transforms the vectorized tensor by a Kronecker product with $A_k$ in the $k$th position. Applying the transformations in all modes gives $\operatorname{vec}(X-M)=(A_K\otimes\cdots\otimes A_1)\operatorname{vec}(Z)$. This is the standard mode-product and Kronecker identity.

Theorem: If $Z$ has independent $N(0,1)$ entries and $X=M+Z\times_1 A_1\times_2\cdots\times_K A_K$, where $A_kA_k^{\mathsf T}=\Sigma_k$ and each $\Sigma_k$ is positive definite, then $$ \operatorname{vec}(X)\sim N\!\left(\operatorname{vec}(M),\Sigma_K\otimes\cdots\otimes\Sigma_1\right). $$

Proof. Put $B=A_K\otimes\cdots\otimes A_1$. The mode-product identity gives $\operatorname{vec}(X-M)=B\operatorname{vec}(Z)$. Since $\operatorname{vec}(Z)$ has mean zero and identity covariance, the transformed vector has mean zero and covariance $BB^{\mathsf T}$. The Kronecker product rule gives

BB^{\mathsf T} =(A_KA_K^{\mathsf T})\otimes\cdots\otimes(A_1A_1^{\mathsf T}) =\Sigma_K\otimes\cdots\otimes\Sigma_1.

A linear transformation of a multivariate normal vector is multivariate normal. Adding $\operatorname{vec}(M)$ shifts its mean to $\operatorname{vec}(M)$, proving the claim. $\square$

The Tensor-Normal Density

For a density in the original tensor entries, write $\Delta=X-M$. The determinant of the full covariance is $\prod_{k=1}^K|\Sigma_k|^{d/d_k}$. To see the exponent, a Kronecker factor $\Sigma_k$ appears once for each combination of indices in all other modes, namely $d/d_k$ times. The quadratic form can be written either in vectorized form or as the squared Frobenius norm after applying the inverse square root of each mode covariance.

Theorem: If each $\Sigma_k$ is symmetric positive definite, the density of $X\sim TN(M;\Sigma_1,\ldots,\Sigma_K)$ is $$ f(X)= (2\pi)^{-d/2} \left\{\prod_{k=1}^K|\Sigma_k|^{-d/(2d_k)}\right\} \exp\!\left[ -\frac12\operatorname{vec}(\Delta)^{\mathsf T} (\Sigma_K^{-1}\otimes\cdots\otimes\Sigma_1^{-1}) \operatorname{vec}(\Delta) \right]. $$ Equivalently, the quadratic form in the exponent is $\|\Delta\times_1\Sigma_1^{-1/2}\times_2\cdots\times_K\Sigma_K^{-1/2}\|_F^2$.

Derivation. In the construction above, the transformation from $\operatorname{vec}(Z)$ to $\operatorname{vec}(X)$ has matrix $B=A_K\otimes\cdots\otimes A_1$. Its absolute determinant is $|B|=\prod_{k=1}^K|A_k|^{d/d_k}=\prod_{k=1}^K|\Sigma_k|^{d/(2d_k)}$. The standard normal density of $\operatorname{vec}(Z)$ is $(2\pi)^{-d/2}\exp(-\operatorname{vec}(Z)^{\mathsf T}\operatorname{vec}(Z)/2)$. Changing variables divides this density by $|B|$. Also, $\operatorname{vec}(Z)=B^{-1}\operatorname{vec}(\Delta)$, so its squared Euclidean norm becomes $\operatorname{vec}(\Delta)^{\mathsf T}(B^{-\mathsf T}B^{-1})\operatorname{vec}(\Delta)$. Because $BB^{\mathsf T}=\Sigma_K\otimes\cdots\otimes\Sigma_1$, the matrix $B^{-\mathsf T}B^{-1}$ is the inverse of that covariance. Substitution gives the stated density. $\square$

Worked Example: Evaluating a Tensor-Normal Density

Let $X$ be a $2$-by-$2$ tensor with mean zero and mode covariances $\Sigma_1=\begin{pmatrix}1&0\\0&4\end{pmatrix}$ and $\Sigma_2=\begin{pmatrix}9&0\\0&16\end{pmatrix}$. Evaluate its density at $X=\begin{pmatrix}1&2\\2&-2\end{pmatrix}$.

Column-stacking gives $\operatorname{vec}(X)=(1,2,2,-2)^{\mathsf T}$, and the covariance is $\Sigma_2\otimes\Sigma_1=\operatorname{diag}(9,36,16,64)$. Thus the quadratic form is

q=\frac{1^2}{9}+\frac{2^2}{36}+\frac{2^2}{16}+\frac{(-2)^2}{64} =\frac19+\frac19+\frac14+\frac1{16} =\frac{77}{144}.

The square root of the determinant of the covariance is $\sqrt{9\cdot36\cdot16\cdot64}=576$. Equivalently, the density's mode-factor term is $|\Sigma_1|^{-1}|\Sigma_2|^{-1}=4^{-1}144^{-1}=1/576$, since $d/(2d_1)=d/(2d_2)=1$. Therefore

f(X)=\frac{1}{576(2\pi)^2}\exp\!\left(-\frac{77}{288}\right).

This value is a density at the specified four-dimensional point, not the probability that a continuous random tensor equals that point.

Interpretation and a Common Pitfall

The tensor-normal assumption replaces one unrestricted $d$-by-$d$ covariance matrix with $K$ smaller mode covariance matrices and a prescribed separable pattern. This can make a Gaussian model interpretable and computationally manageable when dependence is organized by modes. It is also a substantive restriction: the covariance between any two entries must factor into one contribution per mode.

Pitfall: treating the entries as independent. The model is not an independent-entry model unless the resulting full covariance is diagonal. Nonzero off-diagonal entries in a mode covariance can create dependence between tensor entries.
Pitfall: reversing the Kronecker factors without changing vectorization. Under column-stacking, a matrix uses $\Sigma_2\otimes\Sigma_1$, not $\Sigma_1\otimes\Sigma_2$. A different stacking convention changes the displayed order and the correspondence between vector positions and tensor entries.
Pitfall: interpreting each $\Sigma_k$ as the full covariance of the tensor. Each matrix describes the covariance factor for one mode; the covariance of all entries is their Kronecker product.

Key takeaway: With index 1 varying fastest, a tensor-normal distribution has vectorized covariance $\Sigma_K\otimes\cdots\otimes\Sigma_1$. Its entrywise covariances are products of mode-specific covariances, and its density is the multivariate-normal density expressed through this structured covariance.

Check Your Understanding

Use the stated convention that index 1 varies fastest.

For a three-way tensor, write the Kronecker-product order in the covariance of its vectorization.
For two tensor entries, explain how their covariance is obtained from the mode covariance matrices.
What condition on each mode covariance matrix ensures the tensor-normal density is nonsingular?
Why does the covariance matrix of a vectorized $2$-by-$2$ tensor use $\Sigma_2\otimes\Sigma_1$ under column-stacking?
Explain why a tensor-normal model does not generally imply independence of all entries.