Flip-Flop Maximum Likelihood for Kronecker Factors

Alternating estimation of Kronecker factors

A Kronecker covariance model replaces one large covariance matrix by factors that describe variation along separate dimensions. This structure can make estimation more manageable, but the two factors enter the likelihood together. Holding one factor fixed turns estimation of the other into a familiar covariance maximization; alternating these conditional maximizations gives the flip-flop algorithm.

The algorithm is an exact block-coordinate maximum-likelihood procedure, not a promise that every dataset has a finite estimate or that every run converges to a unique global maximum. The existence of each update, the scale ambiguity between factors, and the geometry of the likelihood all matter.

For simplicity, the derivation first treats \(M\) as known. If the common mean is unknown, its conditional maximum-likelihood estimate for fixed \(U,V\) is \(\bar Y=n^{-1}\sum_{k=1}^nY_k\); the factor updates then use residuals \(X_k=Y_k-\bar Y\). The formulas below apply to either set of residual matrices when the mean is held fixed during the factor updates.

The likelihood and its conditional maximizers

The matrix-normal density has determinant factor \(|U|^{-q/2}|V|^{-p/2}\) and quadratic form \(\operatorname{tr}(U^{-1}X_kV^{-1}X_k^T)\). Thus, up to a constant independent of \(U,V\), the log-likelihood is

The two factors are coupled in the quadratic term. If \(V\) is fixed, however, define \(A_V=\sum_{k=1}^nX_kV^{-1}X_k^T\). The terms involving \(U\) are \(-nq\log|U|/2-\operatorname{tr}(U^{-1}A_V)/2\). Similarly, with \(U\) fixed, define \(B_U=\sum_{k=1}^nX_k^TU^{-1}X_k\).

To prove the first formula, set \(W=A_V^{1/2}U^{-1}A_V^{1/2}\), which is positive definite. Then \(\log|U|=\log|A_V|-\log|W|\), and \(\operatorname{tr}(U^{-1}A_V)=\operatorname{tr}(W)\). Apart from terms constant in \(W\), the conditional log-likelihood is

If \(\lambda_1,\ldots,\lambda_p\) are the positive eigenvalues of \(W\), this expression is the sum of \(nq\log(\lambda_j)/2-\lambda_j/2\). Each scalar term has derivative \(nq/(2\lambda_j)-1/2\), positive for \(\lambda_j<nq\) and negative for \(\lambda_j>nq\). It has its unique maximum at \(\lambda_j=nq\). Therefore the matrix maximum is \(W=nqI_p\), which implies \(U=A_V/(nq)\). Applying the same argument to \(B_U\), with \(np\) in place of \(nq\), gives the update for \(V\). \(\square\)

The denominators are the numbers of scalar measurements in each direction across the \(n\) matrices: \(nq\) for the \(p\)-dimensional factor, and \(np\) for the \(q\)-dimensional factor. Under the stated positive-definiteness conditions, each update is a genuine conditional maximum, so it cannot decrease the likelihood when the other factor is held fixed.

Worked updates and scale normalization

Even when both updates exist, the individual factors are not identifiable without a convention: for any \(c>0\), \(V\otimes U=(cV)\otimes(U/c)\). This is the scale nonidentifiability of Kronecker factors. A common convention is \(|V|=1\), which selects one representative from each scale-equivalent pair. After an update pair \((U_*,V_*)\), let \(c=|V_*|^{-1/q}\) and replace it by \((U_*/c,cV_*)\). The replacement has determinant-one second factor and leaves the Kronecker covariance unchanged. Such rescaling is not an additional likelihood improvement; it only fixes the scale convention.

When an update does not exist

An update formula produces a positive-definite estimate only when its scatter matrix is positive definite. If \(A_V\) is singular, the conditional likelihood over positive-definite \(U\) has no finite maximizer: a direction with zero scatter can be assigned an eigenvalue of \(U\) approaching zero, increasing the likelihood without bound. The analogous issue occurs for singular \(B_U\). Thus it is not enough to write down the matrix formula; its positive-definiteness must be checked at every step.

What convergence does and does not mean

When every conditional scatter matrix is positive definite, a flip-flop sweep consists of maximizing over \(U\) at fixed \(V\), then maximizing over \(V\) at the updated \(U\). Each step weakly increases the same likelihood. This proves monotonicity of the likelihood values, but it does not by itself prove that the parameter estimates converge, that a maximum exists, or that the limiting fit is globally optimal.

For the likelihood claim, write \(\ell_t\) for the value after sweep \(t\). Exact conditional maximization gives \(\ell_{t+1}\geq\ell_t\). Continuity and compactness imply that \(\ell_t\) is bounded above, so the nondecreasing sequence has a finite limit. It follows that \(\ell_{t+1}-\ell_t\to0\). Let \(z_t\) denote the normalized pair of factors at the start of a sweep, and let \(F\) be the full-sweep update map. Compactness gives convergent subsequences; take one with \(z_{t_j}\to z_*\). Continuity of \(F\) and \(\ell\) yields \(\ell(F(z_*))-\ell(z_*)=\lim_{j\to\infty}(\ell_{t_j+1}-\ell_{t_j})=0\).

The total gain in a sweep is the sum of two nonnegative conditional gains. A zero total gain means both are zero. Because each conditional maximizer is unique, neither factor can be improved in its block at \(z_*\); hence \(F(z_*)=z_*\). At an interior block maximum the derivative with respect to that block vanishes. Both block derivatives therefore vanish at the fixed point, giving joint stationarity in the unrestricted factor coordinates. The normalization only chooses a representative of the scale-equivalent factors. This proves the stated accumulation-point result, not convergence of the entire sequence to one point. \(\square\)

Common pitfall

A common mistake is to say, “Each block update is a maximum, so the algorithm converges to the global maximum.” The first clause only says that each step is optimal with the other factor held fixed. It implies nondecreasing likelihood, not global optimality; a fixed point may be only a local or coordinatewise optimum. Nor does monotonicity ensure a finite endpoint if the likelihood is unbounded or the iterates escape every compact set. Check positive definiteness of both conditional scatters, address the scale ambiguity, and state the boundedness and continuity conditions before interpreting convergence.