Flip-Flop Maximum Likelihood for Separable Covariance

Why Covariance Factors Are Updated in Turn

A separable covariance model replaces one large covariance matrix for a tensor with a collection of smaller mode-specific matrices. This makes estimation more manageable, but the likelihood still couples the factors: the update for one mode depends on the current values of the others. The flip-flop algorithm handles that coupling by maximizing the likelihood with respect to one factor at a time, then cycling through the modes.

The key calculation is a conditional maximum-likelihood update. It has the form of a covariance estimate based on mode-wise unfoldings, with the other modes’ covariance factors used to weight the columns. Because factor scales can be shifted between modes without changing their Kronecker product, an implementation must also specify how it normalizes the factors.

Model and Notation for the Update

Suppose \(X_1,\ldots,X_N\) are independent \(K\)-way tensor observations, each of dimensions \(p_1\times\cdots\times p_K\). Write \(P=\prod_{j=1}^K p_j\), and let \(\Sigma_j\) be a positive-definite \(p_j\times p_j\) covariance factor for mode \(j\). We use the tensor-normal convention that the covariance of the vectorized tensor is the Kronecker product of the mode covariances, in the order corresponding to the chosen vectorization.

If the common mean tensor is unknown, its maximum-likelihood estimate is the sample mean \(\bar X=N^{-1}\sum_{r=1}^N X_r\). For covariance estimation, use residual tensors \(Y_r=X_r-\bar X\). If the mean is specified, use \(Y_r=X_r-\mu\) instead. The derivation below applies to either set of residuals; in the unknown-mean case it describes the covariance part of the profile likelihood.

Definition: The mode-\(k\) unfolding \(Y_{r(k)}\) arranges the mode-\(k\) fibers of \(Y_r\) as columns of a \(p_k\times q_k\) matrix, where \(q_k=P/p_k\). Let \(B_k\) be the \(q_k\times q_k\) covariance induced on those columns by all modes other than \(k\). It is the Kronecker product of the other mode covariance factors, ordered to match the chosen unfolding. Thus \(B_k\) is positive definite when all \(\Sigma_j\), \(j\ne k\), are positive definite.

The trace term in the tensor-normal log-likelihood can be written, for any mode \(k\), as \(\operatorname{tr}\{\Sigma_k^{-1}Y_{r(k)}B_k^{-1}Y_{r(k)}^{T}\}\). Also, the determinant contribution of \(\Sigma_k\) is \(q_k\log|\Sigma_k|\). These identities follow from the Kronecker covariance and the standard determinant and trace identities for Kronecker products and unfoldings.

Deriving the Conditional Update

Hold every factor except \(\Sigma_k\) fixed. Terms in the log-likelihood that depend on \(\Sigma_k\) are

\ell_k(\Sigma_k) =-\frac{Nq_k}{2}\log|\Sigma_k| -\frac{1}{2}\operatorname{tr}\!\left(\Sigma_k^{-1}S_k\right) +\text{constant}, \qquad S_k=\sum_{r=1}^N Y_{r(k)}B_k^{-1}Y_{r(k)}^T .

Theorem: If \(S_k\) is positive definite, the unique conditional maximum-likelihood estimate of \(\Sigma_k\), with the other factors fixed, is \(\widehat{\Sigma}_k=S_k/(Nq_k)\).

Proof. Set \(C=S_k/(Nq_k)\), which is positive definite under the hypothesis, and write \(\Sigma_k=C^{1/2}MC^{1/2}\) for a positive-definite matrix \(M\). Apart from constants independent of \(M\), maximizing \(\ell_k\) is equivalent to minimizing

Nq_k\log|M|+Nq_k\operatorname{tr}(M^{-1}).

If \(\lambda_1,\ldots,\lambda_{p_k}\) are the positive eigenvalues of \(M\), the expression divided by \(Nq_k\) is \(\sum_{a=1}^{p_k}(\log\lambda_a+1/\lambda_a)\). For \(\lambda>0\), the derivative of \(\log\lambda+1/\lambda\) is \((\lambda-1)/\lambda^2\), so this function decreases up to \(1\), increases after \(1\), and has its unique minimum at \(1\). Consequently every \(\lambda_a=1\) at the unique minimum; hence \(M=I\), \(\Sigma_k=C\), and the stated update follows. \(\square\)

The condition \(S_k>0\) matters: it ensures the conditional maximizer is an interior, positive-definite covariance matrix. If \(S_k\) is singular, the displayed update is not positive definite, and there is no positive-definite conditional maximizer given by this formula. The algorithm therefore requires checking that its mode-wise scatters are positive definite (or using an explicitly specified regularized method, which is a different estimation procedure).

Worked Updates

Worked Example: A Two-Way Update with Identity Weighting

Take \(N=2\) residual matrices of size \(2\times2\), and hold the second-mode covariance at \(\Sigma_2=I_2\). Let \(Y_1=\begin{pmatrix}1&0\\0&1\end{pmatrix}\) and \(Y_2=\begin{pmatrix}1&1\\0&1\end{pmatrix}\). For the mode-1 update, \(q_1=2\) and \(B_1=\Sigma_2=I_2\), so the denominator is \(Nq_1=4\).

The two contributions to \(S_1\) are \(Y_1Y_1^T= \begin{pmatrix}1&0\\0&1\end{pmatrix}\) and \(Y_2Y_2^T=\begin{pmatrix}2&1\\1&1\end{pmatrix}\). Therefore

S_1=\begin{pmatrix}3&1\\1&2\end{pmatrix}, \qquad \widehat{\Sigma}_1=\frac{1}{4} \begin{pmatrix}3&1\\1&2\end{pmatrix} = \begin{pmatrix}0.75&0.25\\0.25&0.50\end{pmatrix}.

This estimate is positive definite: its determinant is \(0.75(0.50)-0.25^2=0.3125>0\), and its first diagonal entry is positive. The calculation is a conditional update only; a later mode-2 update will use this new mode-1 factor.

Worked Example: A Mode-One Update from a Three-Way Sample

Consider two residual tensors of dimensions \(2\times2\times2\), with all three covariance factors initially equal to \(I_2\). For each tensor, display its two mode-1 slices as matrices indexed by \(i=1,2\):

Y_1(:,:,1)=\begin{pmatrix}1&0\\0&0\end{pmatrix}, \quad Y_1(:,:,2)=\begin{pmatrix}0&0\\1&0\end{pmatrix}; \qquad Y_2(:,:,1)=\begin{pmatrix}0&0\\2&0\end{pmatrix}, \quad Y_2(:,:,2)=\begin{pmatrix}0&1\\0&0\end{pmatrix}.

Using columns in the order \((j,k)=(1,1),(1,2),(2,1),(2,2)\), the mode-1 unfoldings are \(Y_{1(1)}=\begin{pmatrix}1&0&0&0\\0&1&0&0\end{pmatrix}\) and \(Y_{2(1)}=\begin{pmatrix}0&0&0&1\\2&0&0&0\end{pmatrix}\). Since the other factors are identities, \(B_1=I_4\). The denominator is \(Nq_1=2\cdot4=8\), and direct multiplication gives

S_1=Y_{1(1)}Y_{1(1)}^T+Y_{2(1)}Y_{2(1)}^T =\begin{pmatrix}2&0\\0&5\end{pmatrix}, \qquad \widehat{\Sigma}_1=\begin{pmatrix}1/4&0\\0&5/8\end{pmatrix}.

To keep this update in a determinant-one normalization for modes 1 and 2, compute \(c=|\widehat{\Sigma}_1|^{1/2}=\sqrt{(1/4)(5/8)}=\sqrt{5/32}\). Replace \(\Sigma_1\) by \(\widehat{\Sigma}_1/c\) and multiply \(\Sigma_3\) by \(c\). The normalized mode-1 factor is \(\operatorname{diag}(\sqrt{2/5},\sqrt{10}/2)\), whose determinant is \(1\). The mode-3 factor becomes \(cI_2\). This rescaling does not change the full Kronecker covariance: dividing one factor by \(c\) and multiplying another by \(c\) leaves their Kronecker product unchanged.

Worked Example: Updating with a Nonidentity Other-Mode Factor

Return to the two residual matrices in the first example, but now hold \(\Sigma_2=\operatorname{diag}(1,4)\). Then \(B_1^{-1}=\operatorname{diag}(1,1/4)\): the second column is downweighted because its covariance is four times as large. The weighted contributions are

Y_1B_1^{-1}Y_1^T=\begin{pmatrix}1&0\\0&1/4\end{pmatrix}, \qquad Y_2B_1^{-1}Y_2^T=\begin{pmatrix}5/4&1/4\\1/4&1/4\end{pmatrix}.

Adding them and dividing by \(Nq_1=4\) gives \(\widehat{\Sigma}_1=\begin{pmatrix}9/16&1/16\\1/16&1/8\end{pmatrix}\). Its determinant is \(9/128-1/256=17/256>0\), so this conditional update is positive definite. It differs from the identity-weighted result because the mode-2 covariance changes how the columns of each unfolding contribute to the mode-1 scatter.

Flip-Flop Cycle and Scale Normalization

A flip-flop cycle applies the conditional formula successively to modes \(1,\ldots,K\), using the latest available values of the other factors each time. Repeat cycles until a chosen stopping rule is met, such as a small relative change in log-likelihood or in the factors under a fixed normalization. The formula specifies an exact conditional maximization, not a guarantee that a particular stopping rule has found the global maximum.

For definiteness, one convenient scale convention is to require \(|\Sigma_j|=1\) for \(j=1,\ldots,K-1\), leaving the overall scale in \(\Sigma_K\). After a raw update \(A_k=S_k/(Nq_k)\) for \(k<K\), calculate \(c=|A_k|^{1/p_k}\), replace \(\Sigma_k\) by \(A_k/c\), and multiply \(\Sigma_K\) by \(c\). After an update of mode \(K\), leave that factor unrestricted; the other factors already have determinant one. This convention is imposed after every update, and the transfer preserves the represented covariance, so it does not change the likelihood. Other consistent scale conventions are possible, but an algorithm must state which one it uses.

Result: Each raw factor update, when its scatter is positive definite, cannot decrease the likelihood with the other factors held fixed. Normalizing by a scale transfer leaves the full covariance unchanged, so it also leaves the likelihood unchanged. Therefore the likelihood is nondecreasing across valid updates and cycles.

Justification. The theorem proves that the raw update maximizes the conditional log-likelihood over positive-definite choices of the updated factor. The factor value before the update is one such choice, so the conditional likelihood after the update is at least as large. A normalization transfer changes neither the Kronecker covariance nor the mean, and therefore changes no probability density. Applying these two facts at each step proves the claim. This argument establishes monotonicity only: it does not establish that the factor sequence converges, that a limit is unique, or that the attained point is a global maximum. \(\square\)

Common Pitfalls and Practical Interpretation

Omitting the column count. The denominator is \(Nq_k\), not just \(N\). Each tensor contributes \(q_k\) mode-\(k\) fibers.
Using an unweighted unfolding scatter. Unless all other factors are identities, the update uses \(B_k^{-1}\) between the unfolding and its transpose.
Normalizing a factor without compensation. Dividing a factor by \(c\) alone changes the covariance model and generally changes the likelihood. The scale must be transferred to another factor.
Interpreting monotonicity as global optimization. Conditional maximization ensures nondecreasing likelihood, not a global solution or convergence of the factors. Initialization and numerical checks remain important.
Ignoring a singular scatter. If \(S_k\) is not positive definite, the stated interior update is unavailable. The issue cannot be fixed by merely writing down its inverse or treating a singular matrix as a valid covariance factor.

The mode-wise formula is useful because it turns a joint covariance problem into a sequence of familiar positive-definite matrix updates. Its practical value depends on keeping the unfolding order, the induced column covariance, the \(Nq_k\) denominator, and the scale convention consistent throughout the computation.

Check Your Understanding

Use the definitions and update derivation above to answer the following questions.

For dimensions \(3\times2\times4\), what are \(q_2\) and the denominator in the mode-2 update for \(N=5\) observations?
Why does the mode-\(k\) scatter include \(B_k^{-1}\) between the unfolding and its transpose?
In the three-way example, verify that the normalized mode-1 factor has determinant one.
What condition on \(S_k\) makes the stated conditional update a positive-definite maximizer?
Why does transferring a scale between two factors preserve the likelihood, and what does likelihood monotonicity fail to prove?