Worked Kronecker Covariance Estimation

From Kronecker Structure to Numerical Estimation

A Kronecker covariance model replaces one large covariance matrix with two smaller factors. Its practical appeal is also a computational challenge: estimating the factors requires alternating between them, and their scales are not separately identified. A worked calculation makes it possible to check both the matrix arithmetic and what the resulting likelihood does—and does not—tell us.

Suppose $X_1,\ldots,X_n$ are independent, mean-zero $p\times q$ random matrices with a common matrix-normal distribution. We use the convention that the row covariance is $U$, the column covariance is $V$, and column-stacking gives $\operatorname{Cov}(\operatorname{vec}(X_i))=V\otimes U$. Both factors are required to be symmetric positive definite.

Model and flip-flop updates: The density of each observation is $$ f(X_i)=(2\pi)^{-pq/2}|U|^{-q/2}|V|^{-p/2} \exp\left\{-\frac12\operatorname{tr}(U^{-1}X_iV^{-1}X_i^{T})\right\}. $$ For fixed $V$, the row-factor update is $U_{\mathrm{new}}=(nq)^{-1}\sum_{i=1}^n X_iV^{-1}X_i^T$. For fixed $U$, the column-factor update is $V_{\mathrm{new}}=(np)^{-1}\sum_{i=1}^n X_i^TU^{-1}X_i$. Each update is a conditional maximum likelihood update when its displayed scatter matrix is positive definite.

Why the Conditional Updates Have These Forms

The updates are the flip-flop maximum likelihood equations. The following derivation records the conditional calculation needed to interpret the numerical example. With $V$ fixed, terms in the log likelihood that depend on $U$ are

\ell(U\mid V)=-\frac{nq}{2}\log|U|-\frac12\operatorname{tr}(U^{-1}S_U)+C, \qquad S_U=\sum_{i=1}^n X_iV^{-1}X_i^T,

where $C$ does not depend on $U$. If $S_U$ is positive definite, put $B=S_U^{1/2}U^{-1}S_U^{1/2}$, also positive definite. Since $|U|=|S_U|/|B|$ and $\operatorname{tr}(U^{-1}S_U)=\operatorname{tr}(B)$, maximizing over $U$ is equivalent to maximizing

\frac{nq}{2}\log|B|-\frac12\operatorname{tr}(B) =\sum_{j=1}^{p}\left(\frac{nq}{2}\log\lambda_j-\frac{\lambda_j}{2}\right),

where $\lambda_j$ are the positive eigenvalues of $B$. Each summand has derivative $nq/(2\lambda_j)-1/2$, which is positive for $\lambda_j<nq$, zero at $nq$, and negative for $\lambda_j>nq$. Thus the unique maximum is attained when $B=nqI_p$, giving $U=S_U/(nq)$. The same argument with $U$ fixed gives the column update: use $S_V=\sum_iX_i^TU^{-1}X_i$ and the $q$ eigenvalues of $S_V^{1/2}V^{-1}S_V^{1/2}$, whose maximizing value is $np$. This proves the stated conditional updates when the corresponding scatters are positive definite. $\square$

Worked Alternating Updates

Worked Example: Two Flip-Flop Sweeps for Three Matrices

Take $n=3$, $p=q=2$, and suppose the three observed residual matrices are

X_1=\begin{pmatrix}2&0\\0&1\end{pmatrix},\qquad X_2=\begin{pmatrix}0&1\\1&0\end{pmatrix},\qquad X_3=\begin{pmatrix}1&1\\1&-1\end{pmatrix}.

Start with $V_0=I_2$. The first row scatter is

\sum_{i=1}^3X_iX_i^T =\begin{pmatrix}4&0\\0&1\end{pmatrix} +\begin{pmatrix}1&0\\0&1\end{pmatrix} +\begin{pmatrix}2&0\\0&2\end{pmatrix} =\begin{pmatrix}7&0\\0&4\end{pmatrix}.

Because $nq=6$, the first row update is $U_1=\operatorname{diag}(7/6,\,2/3)$, with inverse $\operatorname{diag}(6/7,\,3/2)$. The column scatter, using this updated $U_1$, is

\begin{aligned} \sum_{i=1}^3X_i^TU_1^{-1}X_i &=\begin{pmatrix}24/7&0\\0&3/2\end{pmatrix} +\begin{pmatrix}3/2&0\\0&6/7\end{pmatrix} +\begin{pmatrix}33/14&-9/14\\-9/14&33/14\end{pmatrix}\\ &=\begin{pmatrix}51/7&-9/14\\-9/14&33/7\end{pmatrix}. \end{aligned}

Dividing by $np=6$ gives

V_1=\begin{pmatrix}17/14&-3/28\\-3/28&11/14\end{pmatrix}.

To check that the next row update uses the inverse of this matrix, its determinant and inverse are

|V_1|=\frac{739}{784},\qquad V_1^{-1}=\frac1{739}\begin{pmatrix}616&84\\84&952\end{pmatrix}.

For a symmetric matrix $W=\begin{pmatrix}a&b\\b&c\end{pmatrix}$, direct multiplication of these three data matrices gives

\sum_{i=1}^3X_iWX_i^T =\begin{pmatrix}5a+2b+2c&a+3b-c\\a+3b-c&2a-2b+2c\end{pmatrix}.

Substituting $W=V_1^{-1}$ and dividing by six therefore yields

U_2=\begin{pmatrix} 2576/2217&-14/739\\ -14/739&1484/2217 \end{pmatrix}.

For the following column update, the inverse is

U_2^{-1}=\frac{2217}{3821020}\begin{pmatrix}1484&42\\42&2576\end{pmatrix}.

If $W=U_2^{-1}=\begin{pmatrix}a&b\\b&c\end{pmatrix}$, then direct multiplication gives $\sum_iX_i^TWX_i=\begin{pmatrix}5a+2b+2c&a+3b-c\\a+3b-c&2a-2b+2c\end{pmatrix}$. Substitution and division by six give $V_2$, approximately

V_2\approx\begin{pmatrix}1.2239&-0.0934\\-0.0934&0.7771\end{pmatrix}.

The matrices are positive definite: $U_1$ has positive diagonal entries, $|V_1|=739/784>0$ with positive diagonal entries, and $U_2$ has positive first diagonal entry and determinant $3821020/2217^2>0$. The column scatter that produces $V_2$ is positive definite as well: for nonzero $z$, its quadratic form is $\sum_i (X_i z)^TU_2^{-1}(X_i z)$, which is positive because the matrices $X_i$ have no common nonzero vector in their null spaces. Thus these displayed updates are valid. A finite number of sweeps produces iterates; it does not by itself establish convergence or that the likelihood has reached its maximum.

Scale Normalization Without Changing the Model

The factors are not separately identified by their Kronecker product. For any $a>0$, $(aU)\otimes(V/a)=U\otimes V$. A convenient convention is to set the determinant of the column factor to one. Applying that convention after a sweep changes the factor matrices, but not the covariance represented by them.

Scale-normalization result: For $U$ of dimension $p\times p$ and $V$ of dimension $q\times q$, let $c=|V|^{-1/q}$, $U^*=U/c$, and $V^*=cV$. If $U,V$ are positive definite, then $|V^*|=1$ and $V^*\otimes U^*=V\otimes U$.

Indeed, $|V^*|=|cV|=c^q|V|=1$, while the scalar factors cancel in the Kronecker product: $(cV)\otimes(U/c)=V\otimes U$. This proves the result. $\square$

Worked Example: Normalize the First Column Factor

For the factors in the first sweep, $|V_1|=739/784$, so the determinant-one multiplier is $c=|V_1|^{-1/2}=28/\sqrt{739}$. The normalized factors are therefore

V_1^*=\frac{28}{\sqrt{739}}V_1,\qquad U_1^*=\frac{\sqrt{739}}{28}U_1.

The determinant check is $|V_1^*|=c^2|V_1|=(784/739)(739/784)=1$. The represented covariance is unchanged because $V_1^*\otimes U_1^*=(cV_1)\otimes(U_1/c)=V_1\otimes U_1$. This normalization fixes a scale convention; it does not make the two factors contain new information or improve the fit.

Evaluating the Likelihood at an Iterate

A likelihood calculation should identify exactly which pair of factors is being evaluated. There is a useful check when $V$ has just been updated with $U$ held fixed. By the column update, $\sum_iX_i^TU^{-1}X_i=npV$. Taking a trace after multiplying by $V^{-1}$ gives

\sum_{i=1}^n\operatorname{tr}(U^{-1}X_iV^{-1}X_i^T) =\operatorname{tr}\left(V^{-1}\sum_{i=1}^nX_i^TU^{-1}X_i\right) =np\,\operatorname{tr}(I_q)=npq.

This identity applies to the pair consisting of the fixed row factor and its just-updated column factor. It is not a claim that the pair is a joint maximum likelihood estimate.

Worked Example: Likelihood After the First Column Update

For $U_1,V_1$, the update identity gives the sum of quadratic terms as $npq=3\cdot2\cdot2=12$. Also, $|U_1|=(7/6)(2/3)=7/9$ and $|V_1|=739/784$. Since the determinant exponents in the joint likelihood are both $-n$ for $p=q=2$, substitution into the matrix-normal density gives

L(U_1,V_1)=(2\pi)^{-6} \left(\frac79\right)^{-3} \left(\frac{739}{784}\right)^{-3} e^{-12/2} =(2\pi)^{-6}\left(\frac{1008}{739}\right)^3e^{-6}.

Equivalently, $\ell(U_1,V_1)=-6\log(2\pi)-3\log(739/1008)-6$, approximately $-16.10$. This evaluates the likelihood at the first full sweep, not at a claimed converged estimate. Because normalization leaves $V\otimes U$ unchanged, using $U_1^*,V_1^*$ gives the same likelihood.

Common Pitfalls and Why the Checks Matter

Calling one sweep the MLE. Alternating conditional maximizations produce successive iterates. A calculated pair is not automatically a converged joint maximum; convergence and the likelihood behavior must be checked separately.
Normalizing only one factor. Replacing $V$ by $cV$ while leaving $U$ fixed changes $V\otimes U$. To preserve the model, replace $U$ by $U/c$ at the same time.
Using the wrong dimensions in the divisors. The row update divides by $nq$, while the column update divides by $np$. These factors count observations and the dimension of the opposite axis.
Treating an indefinite scatter as a valid covariance update. The update must be positive definite to be an admissible factor. In a numerical example, check symmetry and positive definiteness rather than assuming them.

These checks matter in applications because factor scale conventions affect reported parameters, while the likelihood depends on the represented covariance. Reporting the normalization and the iteration status makes estimates interpretable and prevents a finite computational step from being mistaken for a completed optimization.

Check Your Understanding

Use the definitions and calculations above to answer the following.

For a $p\times q$ matrix observation, why does the row update divide by $nq$, rather than $np$?
What pair of factor transformations preserves $V\otimes U$ when $V$ is multiplied by a positive constant $c$?
After the first sweep in the example, what is the sum of the three quadratic terms in the log likelihood, and how can it be checked from the column update?
Why does a positive determinant alone not establish that a symmetric covariance factor is positive definite?
What additional conclusion is not justified merely by calculating $U_1$ and $V_1$?