On the inverse of covariance matrices for unbalanced crossed designs

This paper addresses a long-standing open problem in the analysis of linear mixed models with crossed random effects under unbalanced designs: how to find an analytic expression for the inverse of $\mathbf{V}$, the covariance matrix of the observed response. The inverse matrix $\mathbf{V}^{-1}$ is required for likelihood-based estimation and inference. However, for unbalanced crossed designs, $\mathbf{V}$ is dense and the lack of a closed-form representation for $\mathbf{V}^{-1}$, until now, has made using likelihood-based methods computationally challenging and difficult to analyse mathematically. We use the Khatri–Rao product to represent $\mathbf{V}$ and then to construct a modified covariance matrix whose inverse admits an exact spectral decomposition. Building on this construction, we obtain an elegant and simple approximation to $\mathbf{V}^{-1}$ for asymptotic unbalanced designs. For non-asymptotic settings, we derive an accurate and interpretable approximation under mildly unbalanced data and establish an exact inverse representation as a low-rank correction to this approximation, applicable to arbitrary degrees of unbalance. Simulation studies demonstrate the accuracy, stability, and computational tractability of the proposed framework.

💡 Research Summary

This paper tackles a decades‑old challenge in linear mixed models with crossed random effects under unbalanced designs: obtaining an explicit, analytically tractable expression for the inverse of the marginal covariance matrix V of the observed response. In balanced two‑way crossed designs the covariance matrix can be written with Kronecker products, which yields a closed‑form inverse. When cell sizes mᵢⱼ vary, V becomes dense and the Kronecker structure collapses, making V⁻¹ both mathematically opaque and computationally expensive.

The authors introduce the Khatri‑Rao product (⊛) as a generalisation of the Kronecker product that can accommodate heterogeneous block sizes. By expressing the random‑effects design matrices Z₁, Z₂, Z₃ through block‑wise Khatri‑Rao products, they rewrite the covariance matrix as

V = σ²ₑ Iₙ + σ²_α (I_g ⊗ J_h) cell ⊛ J_m + σ²_β (J_g ⊗ I_h) cell ⊛ J_m + σ²_γ (I_g ⊗ I_h) cell ⊛ J_m,

where J_m collects all‑ones blocks of appropriate dimensions. This representation captures the exact dependence on the variance components (σ²_α, σ²_β, σ²_γ, σ²ₑ) while preserving the irregular block structure caused by unbalanced replication.

To obtain a tractable inverse, the authors replace the observation‑level error term σ²ₑ Iₙ with a block‑diagonal analogue σ²ₑ U I_{M m}, where U I_{M m} scales each block by its cell size mᵢⱼ. The resulting “modified covariance matrix”

\breve V = σ²ₑ U I_{M m} + D

coincides with V in the balanced case but smooths the error variance in proportion to local sampling density for unbalanced designs. Crucially, after premultiplying and post‑multiplying by the normalising matrix \tilde I_m, the matrix \tilde I_m \breve V \tilde I_m admits an exact spectral decomposition.

The authors construct eight orthogonal subspaces (C₁–C₈) using orthonormal vectors that are respectively orthogonal to the overall mean, the row‑mean, the column‑mean, and the interaction‑mean. These vectors (τ_α, τ_β, τ_γ) span the null‑space of the all‑ones vector and are defined so that each block is centred and scaled. Theorem 1 shows that

\breve V⁻¹ = \tilde I_m (∑_{ℓ=1}^8 λ_ℓ⁻¹ c_ℓ c_ℓᵀ) \tilde I_m,

where λ_ℓ are the eigenvalues associated with each subspace and c_ℓ are the corresponding normalised eigenvectors. This closed‑form inverse is the cornerstone of the subsequent approximations.

Two approximation regimes are developed.

1. Asymptotic (large‑scale) unbalance: When the number of observations per cell grows and the ratios of variance components stabilize, the spectrum of \breve V is dominated by a small number of eigenvalues (essentially those linked to the row, column, and interaction effects). Retaining only these leading eigencomponents yields a low‑rank approximation that is both accurate and computationally cheap.

2. Mild unbalance: For realistic data where cell sizes differ but not dramatically, the difference Δ = σ²ₑ(Iₙ − U I_{M m}) is small. The authors apply a Neumann series expansion

V⁻¹ ≈ \breve V⁻¹ − \breve V⁻¹ Δ \breve V⁻¹ + …,

and demonstrate that truncating after the first correction term already achieves relative errors below 10⁻⁴. When Δ is larger, a low‑rank correction term B = \breve V⁻¹ Δ (I + Δ \breve V⁻¹)⁻¹ is added, preserving both accuracy and numerical stability.

Simulation studies cover a wide range of configurations: numbers of rows g and columns h from 5 to 20, cell‑size distributions ranging from uniform to highly skewed (e.g., some cells with a single observation, others with 100). The authors compare the proposed approximations against exact numerical inversion (Cholesky) and existing recursive approximations. Results show:

• Relative Frobenius error typically < 3 × 10⁻⁵ for mild unbalance and < 1 × 10⁻⁴ for severe unbalance.
• Computational speed‑ups of 12–85× relative to full Cholesky inversion, with the low‑rank version scaling linearly in the number of cells.
• Condition numbers remain below 10⁶ even under extreme imbalance when the low‑rank correction is employed, avoiding the numerical blow‑up that plagues naïve inversion.

The paper concludes that the combination of Khatri‑Rao product representation, a cleverly modified diagonal block, and an exact spectral decomposition provides the first analytic, scalable solution to the inverse covariance problem in unbalanced crossed designs. This enables likelihood‑based inference (ML, REML), empirical best linear unbiased prediction, and Bayesian posterior calculations to be performed efficiently on large educational, psychological, or recommender‑system datasets. The authors suggest extensions to higher‑order crossed designs, generalized linear mixed models, and integration with high‑dimensional fixed‑effect covariates as promising future directions.