Worked Analysis of Three-Way Gaussian Data

From a Three-Way Dataset to Two Fitted Structures

A three-way data array can contain two quite different kinds of structure. A systematic signal may be well described by a low-rank Tucker representation, while variation around that signal may have dependence along all three modes. The first structure concerns the mean; the second concerns the covariance. Fitting one does not substitute for fitting the other.

We work with $S$ observed tensors $Y_s\in\mathbb R^{I_1\times I_2\times I_3}$, modeled as $Y_s=M+E_s$. The signal $M$ is approximated by a Tucker tensor, and the errors are modeled as independent, mean-zero tensor-normal observations with separable mode covariances. The worked benchmark below is deliberately balanced to make its covariance calculations exact. It is useful for checking the fitting mechanics, but its sign-pattern residuals are not an iid Gaussian sample; that distinction matters when interpreting diagnostics.

Definition: A three-way tensor-normal model with mean $M$ and positive-definite mode covariance matrices $\Sigma_1,\Sigma_2,\Sigma_3$ has independent errors $E_s=Y_s-M$ satisfying $$ \operatorname{vec}(E_s)\sim N\!\left(0,\ \Sigma_3\otimes\Sigma_2\otimes\Sigma_1\right). $$ A Tucker signal of multilinear rank at most $(r_1,r_2,r_3)$ has the form $M=G\times_1 U_1\times_2 U_2\times_3 U_3$, where $G$ is an $r_1\times r_2\times r_3$ core and each $U_k$ has $r_k$ columns. Here $\times_k$ denotes multiplication along mode $k$.

The separable covariance has a scale nonidentifiability: multiplying one factor by a positive constant and another by its reciprocal leaves the Kronecker product unchanged. We use the determinant-one convention $\det(\Sigma_k)=1$ for each mode. It fixes the factor scales while retaining their shapes. When the mean is estimated, the covariance fitting steps below use residuals $E_s=Y_s-\hat M$; they do not by themselves determine how $\hat M$ was obtained.

A Balanced Synthetic Three-Way Benchmark

Take $I_1=I_2=I_3=2$ and $S=16$. Start with a $16\times16$ Hadamard matrix whose entries are $+1$ or $-1$, whose distinct columns are orthogonal, and whose first column is constant. Let $x_s\in\mathbb R^8$ contain the first eight nonconstant-column entries in row $s$. Each coordinate has zero sample mean, and orthogonality gives

\sum_{s=1}^{16}x_s=0,\qquad \frac{1}{16}\sum_{s=1}^{16}x_sx_s^\mathsf{T}=I_8.

Arrange $x_s$ as a $2\times2\times2$ tensor $X_s$, and set

\Sigma_1= \begin{pmatrix}1/2&0\\0&2\end{pmatrix},\quad \Sigma_2= \begin{pmatrix}1/\sqrt{2}&0\\0&\sqrt{2}\end{pmatrix},\quad \Sigma_3= \begin{pmatrix}1/\sqrt{3}&0\\0&\sqrt{3}\end{pmatrix}.

All three matrices have determinant one. Define $L_k=\Sigma_k^{1/2}$, the symmetric positive-definite square root, and let $R_s=X_s\times_1L_1\times_2L_2\times_3L_3$. Use the rank-one signal $M=3a\otimes b\otimes c$, with $a=c=(1,1)^\mathsf{T}/\sqrt{2}$ and $b=(1,-1)^\mathsf{T}/\sqrt{2}$, and set $Y_s=M+R_s$. This is a transparent covariance-calibrated benchmark: it has the exact second moments specified by the tensor-normal model, but the residual coordinates come from a finite sign design rather than independent Gaussian draws.

Worked Example: Verifying the Benchmark Covariance

Let $E_{s,(1)}$ be the mode-1 matricization of $X_s$, a $2\times4$ matrix. For any fixed $4\times4$ matrix $A$, the $(i,i')$ entry of $\sum_s E_{s,(1)}AE_{s,(1)}^\mathsf{T}$ is a sum over the two remaining coordinate indices $q,r$. Orthogonality of the eight columns gives

\left(\sum_{s=1}^{16}E_{s,(1)}AE_{s,(1)}^\mathsf{T}\right)_{ii'} =16\,\delta_{ii'}\sum_{q=1}^{4}A_{qq}.

Thus the full matrix is $16\operatorname{tr}(A)I_2$, not $\operatorname{tr}(A)I_2$. For example, with $A=I_4$, the sum is $16\cdot4 I_2=64I_2$. This factor of 16 is essential in the mode-wise covariance update.

Applying $L_1,L_2,L_3$ to the three modes transforms the identity empirical covariance into the Kronecker covariance $\Sigma_3\otimes\Sigma_2\otimes\Sigma_1$. Therefore the centered residuals $R_s$ have exactly that empirical covariance when the divisor is 16. The result follows directly from the linear transformation rule for covariance.

Fitting Separable Covariance by Flip-Flop Updates

For fixed residuals, write $E_{s,(k)}$ for the mode-$k$ matricization, with $I_k$ rows and $\prod_{\ell\ne k}I_\ell$ columns. Holding the other two factors fixed, the unnormalized maximum-likelihood update is

\widetilde\Sigma_k= \frac{1}{S\prod_{\ell\ne k}I_\ell} \sum_{s=1}^{S}E_{s,(k)} \left(\bigotimes_{\ell\ne k}\Sigma_\ell^{-1}\right)_{\!k} E_{s,(k)}^\mathsf{T}.

The subscript $k$ indicates that the inverse factors are arranged in the order matching the columns of the mode-$k$ matricization. Under the determinant-one convention, normalize each update as $\Sigma_k^{\text{new}}=\widetilde\Sigma_k/\det(\widetilde\Sigma_k)^{1/I_k}$. The flip-flop algorithm cycles through the modes, updating one factor at a time while holding the others fixed.

Result: For the balanced benchmark, the determinant-one flip-flop updates recover all three specified covariance factors after one complete sweep from identity initialization. Indeed, if the empirical residual covariance is $\Sigma_3\otimes\Sigma_2\otimes\Sigma_1$, then the mode-$k$ update with the other factors set to their true values equals $\Sigma_k$. If some other factors are identity, the update is a positive scalar multiple of $\Sigma_k$; determinant-one normalization removes that scalar.

To see the scalar claim for mode 1 when both other factors start at identity, let $A=I_4$ in the preceding orthogonality calculation. Before the $L_k$ transformations, the mode-1 sum is $64I_2$; after transforming the other two modes and mode 1, the sum is $16\operatorname{tr}(\Sigma_2)\operatorname{tr}(\Sigma_3)\Sigma_1$. The update denominator is $S I_2 I_3=16\cdot4=64$. Hence

\widetilde\Sigma_1= \frac{\operatorname{tr}(\Sigma_2)\operatorname{tr}(\Sigma_3)}{4}\Sigma_1.

The multiplier is positive, and $\det(\Sigma_1)=1$, so determinant-one normalization returns $\Sigma_1$ exactly. The same argument applies to modes 2 and 3. More generally, when the other factors are correct, their inverse weighting cancels their covariance contributions and the update is exactly the remaining factor.

Worked Example: Checking One Flip-Flop Sweep

Initialize all three factors at $I_2$. Since the true factors have traces $5/2$, $\sqrt{2}+1/\sqrt{2}=3/\sqrt{2}$, and $\sqrt{3}+1/\sqrt{3}=4/\sqrt{3}$, the unnormalized mode-1 update is

\widetilde\Sigma_1= \frac{(3/\sqrt{2})(4/\sqrt{3})}{4} \begin{pmatrix}1/2&0\\0&2\end{pmatrix} =\sqrt{\frac{3}{2}}\Sigma_1.

Its determinant is $3/2$, so dividing by $(3/2)^{1/2}$ gives the determinant-one update $\Sigma_1$. With this factor in place and the third factor still $I_2$, the mode-2 update is $\widetilde\Sigma_2=(\operatorname{tr}\Sigma_3/2)\Sigma_2=(2/\sqrt{3})\Sigma_2$. Its determinant is $4/3$, and normalization again gives $\Sigma_2$. Once modes 1 and 2 are correct, mode 3's update is $\Sigma_3$. Thus the sweep recovers the specified shapes exactly.

A numerical implementation would monitor both the factor changes, for example $\|\Sigma_k^{(t+1)}-\Sigma_k^{(t)}\|_F/\|\Sigma_k^{(t)}\|_F$, and the observed-data objective. In this exact benchmark those changes are zero after the sweep, up to rounding. In ordinary finite Gaussian samples they are generally nonzero, so a stopping rule should use a declared tolerance rather than expect exact equality.

Likelihood and Residual Checks

With the mean fixed, the negative log-likelihood, up to an additive constant, is

\frac{S}{2}\sum_{k=1}^{3} \left(\prod_{\ell\ne k}I_\ell\right)\log\det(\Sigma_k) +\frac12\sum_{s=1}^{S} \operatorname{vec}(E_s)^\mathsf{T} (\Sigma_3^{-1}\otimes\Sigma_2^{-1}\otimes\Sigma_1^{-1}) \operatorname{vec}(E_s).

Under the determinant-one convention the log-determinant terms vanish. In the benchmark, whitening each residual by the true factors returns $X_s$; the total whitened squared norm is therefore $\sum_s\|x_s\|^2=16\cdot8=128$. The objective is $64$, excluding the constant. This supplies a direct check on the fitted likelihood calculation. It is not a universal target value for other data sets.

A covariance diagnostic asks whether the whitened residuals have approximately zero mean and identity empirical covariance. A Gaussian diagnostic asks more: whether their distribution is plausibly Gaussian, including its tail behavior and joint shape. A low-rank signal diagnostic asks yet another question: how much of the signal remains outside the selected Tucker subspace. These checks should not be collapsed into one score.

Worked Example: Tucker Reconstruction and a Residual Pitfall

The chosen signal is $M=3a\otimes b\otimes c$, with all three vectors of unit norm. It has Tucker rank $(1,1,1)$, core equal to the scalar $3$, and squared Frobenius norm $3^2=9$. Because $\sum_sR_s=0$, the sample mean is $\bar Y=M$. Projecting $\bar Y$ onto the one-dimensional mode spaces spanned by $a,b,c$ gives the core $3$, so its Tucker reconstruction is exactly $M$. The signal reconstruction error $\|\bar Y-\hat M\|_F^2$ is zero.

After subtracting this fitted signal and whitening by the true $L_k^{-1}$, the residual tensors are $X_s$. Their empirical mean is zero and their empirical covariance is $I_8$, as required by the covariance check. Yet every whitened coordinate is $+1$ or $-1$, so its empirical fourth moment is $1$, not the Gaussian value $3$. Each whitened vector also has squared norm exactly $8$. These patterns reflect the benchmark's balanced sign construction, not evidence that a genuine iid Gaussian sample must have constant norm. They demonstrate why matching covariance does not establish Gaussianity.

For a general Tucker fit, one can report the relative reconstruction error $\|\bar Y-\hat M\|_F^2/\|\bar Y\|_F^2$, specifying the ranks and fitting procedure. Truncating singular vectors of mode unfoldings is a useful computational construction, but it should not automatically be described as the globally optimal Tucker approximation. Here the signal has exactly the stated rank, so the zero error follows from the explicit rank-one form, not from a general optimality claim.

Interpreting a Three-Way Fit

A useful report separates the fitted objects and their checks: give the Tucker ranks and reconstruction error for the signal; give the determinant-normalized covariance factors and the flip-flop stopping criterion for the residual dependence; and report whitened residual covariance and distributional diagnostics separately. For data generated by an actual Gaussian simulation, finite-sample residuals will not generally whiten to an exact identity matrix, and likelihood convergence does not certify model adequacy.

Scale: state the factor normalization, since raw factor scales are not separately identifiable.
Convergence: report factor changes and objective changes; a stable objective alone can conceal a slowly changing factorization.
Reconstruction: identify whether the error is for the fitted mean signal or for the full data including residual noise.
Diagnostics: distinguish a second-moment check from a check of Gaussian shape, and interpret each in light of sample size and data construction.

Common pitfall: The identity $\frac{1}{16}\sum_s x_sx_s^\mathsf{T}=I_8$ does not imply that a sum over 16 observations has the same scale as an average. In the mode-1 calculation, $\sum_s E_{s,(1)}AE_{s,(1)}^\mathsf{T}=16\operatorname{tr}(A)I_2$. Omitting the factor 16 gives an incorrectly scaled update. Determinant normalization may hide that error in a special shape-only calculation, but the unnormalized update and likelihood arithmetic would still be wrong.

Check Your Understanding

Use the definitions and calculations above to check the distinctions among signal, covariance, and residual diagnostics.

Why does the mode-1 sum with $A=I_4$ equal $64I_2$, rather than $4I_2$?
What does determinant-one normalization resolve, and what feature of each covariance factor does it retain?
For the benchmark, why does the sample mean recover $M$ exactly?
Why is identity empirical covariance of the whitened residuals not enough to conclude that those residuals are Gaussian?
In a non-balanced Gaussian sample, which convergence quantities would you monitor during flip-flop fitting?