Fast Likelihood Computation for Structured Covariance

Likelihood Computation as Two Linear Algebra Tasks

For a multivariate Gaussian model, evaluating the likelihood does not require forming the inverse covariance matrix. It requires a log-determinant and a quadratic form, and the latter can be obtained by solving a linear system. Structure in the covariance can make both tasks substantially cheaper than dense matrix operations.

Let \(y\in\mathbb{R}^d\) be an observation, with \(y\sim N(\mu,\Sigma)\), where \(\Sigma\) is symmetric positive definite. The log-likelihood for \(\mu\) and \(\Sigma\), evaluated at \(y\), is

To compute the last term, set \(r=y-\mu\) and solve \(\Sigma x=r\). Then \(r^\top\Sigma^{-1}r=r^\top x\). For numerical work, a Cholesky factorization \(\Sigma=LL^\top\) is a standard route: \(\log|\Sigma|=2\sum_i\log L_{ii}\), and two triangular solves give \(x\). Structure is useful when it reduces the cost of the factorization, the solves, or both.

Blocks, Bands, and Their Computational Costs

If \(\Sigma\) is block diagonal, each block can be factored separately. Its determinant is the product of the block determinants, and its solve separates into independent block solves. More generally, block-banded matrices couple only nearby blocks. For a symmetric positive-definite matrix with scalar half-bandwidth \(b\), Cholesky factorization preserves the band. The factorization costs \(O(db^2)\) operations and \(O(db)\) storage; a solve costs \(O(db)\). These costs are advantageous when \(b\) is small relative to \(d\).

For \(N\) blocks of size \(r\) in a block-tridiagonal covariance, block Cholesky uses \(O(Nr^3)\) operations and \(O(Nr^2)\) storage, assuming dense blocks and no additional exploitable structure within them. By comparison, a generic dense factorization of the \(Nr\)-dimensional covariance costs \(O((Nr)^3)\). These are operation-count comparisons; actual run time also depends on implementation and numerical details.

Kronecker Covariances

Suppose \(\Sigma=A\otimes B\), where \(A\) is \(p\times p\) and \(B\) is \(q\times q\), both symmetric positive definite. The full covariance has dimension \(pq\), but its determinant and inverse inherit simple formulas from the factors. For any compatible matrix \(Y\), column-wise vectorization satisfies \(\operatorname{vec}(Y)^\top(A^{-1}\otimes B^{-1})\operatorname{vec}(Y) =\operatorname{tr}(Y^\top B^{-1}Y A^{-1})\). Thus, factor solves on \(A\) and \(B\) can replace a solve with the full \(pq\times pq\) matrix.

Proof. The inverse identity follows from \((A\otimes B)(A^{-1}\otimes B^{-1})=(AA^{-1})\otimes(BB^{-1})=I_p\otimes I_q=I_{pq}\). To obtain the determinant, let \(\lambda_1,\ldots,\lambda_p\) and \(\gamma_1,\ldots,\gamma_q\) be the eigenvalues of \(A\) and \(B\), counting multiplicities. The eigenvalues of \(A\otimes B\) are all products \(\lambda_i\gamma_j\). Therefore \[ |A\otimes B|=\prod_{i=1}^p\prod_{j=1}^q\lambda_i\gamma_j =\left(\prod_{i=1}^p\lambda_i\right)^q \left(\prod_{j=1}^q\gamma_j\right)^p =|A|^q|B|^p. \]

For positive-definite factors, factorizing the two factors costs \(O(p^3+q^3)\), rather than the \(O(p^3q^3)\) cost of dense factorization at dimension \(pq\). This comparison assumes the same factors can be reused as needed and that the calculation can exploit their Kronecker form.

Low-Rank Updates and Matrix Identities

A useful pattern is a covariance that is easy to factor except for a low-rank addition. Suppose \(D\) is \(d\times d\) positive definite, \(U\) is \(d\times r\), and \(C\) is \(r\times r\) positive definite, with \(r\) much smaller than \(d\). Then \(\Sigma=D+UCU^\top\) is positive definite. The determinant and inverse can be reduced to calculations involving an \(r\times r\) matrix.

Proof. Write \(\Sigma=D+UCU^\top=D(I+D^{-1}UC U^\top)\). The identity \(|I+PQ|=|I+QP|\), for rectangular \(P,Q\) of compatible dimensions, gives \[ |\Sigma|=|D|\,|I_r+C U^\top D^{-1}U| =|D|\,|C|\,|C^{-1}+U^\top D^{-1}U|. \] For the inverse, set \(M=C^{-1}+U^\top D^{-1}U\) and \(R=D^{-1}-D^{-1}UM^{-1}U^\top D^{-1}\). Multiplication gives \[ (D+UCU^\top)R =I+U\{C-M^{-1}-C(U^\top D^{-1}U)M^{-1}\}U^\top D^{-1}. \] Since \(M=C^{-1}+U^\top D^{-1}U\), the expression in braces is \(C-(I+C U^\top D^{-1}U)M^{-1}=0\), using \(I+C U^\top D^{-1}U=CM\). Hence \(\Sigma R=I\), so \(R=\Sigma^{-1}\). \(\square\)

This is the Woodbury matrix identity together with the matrix determinant lemma. If \(D\) is diagonal or has a cheap structured factorization, the main additional factorization is only \(r\times r\). The gain depends on being able to apply \(D^{-1}\) efficiently; explicitly forming a dense \(D^{-1}\) would defeat that purpose.

Toeplitz Structure: Exploit a Model, Not Just a Pattern

A Toeplitz matrix has constant entries along each diagonal. Storing its first row and column takes only \(O(d)\) numbers, but that fact alone does not establish that an exact likelihood can be evaluated in \(O(d)\) operations. Specialized algorithms may exploit Toeplitz structure, but their complexity, assumptions, and numerical stability depend on the particular model and algorithm. A transparent example is the covariance of a stationary Gaussian AR(1) process with unit marginal variance: \(\Sigma_{ij}=\rho^{|i-j|}\), where \(|\rho|<1\).

For a vector \(x=(x_1,\ldots,x_d)^\top\) from this model, its density factors into the marginal density of \(x_1\), \(N(0,1)\), and the conditional densities \(x_t\mid x_{t-1}\sim N(\rho x_{t-1},1-\rho^2)\), for \(t=2,\ldots,d\). Consequently,

The formulas give the determinant and quadratic form using a single pass through the observations. They follow from the model's innovation representation, not merely from the fact that its covariance is Toeplitz.

Common Pitfall: Confusing Structure With a Guaranteed Speedup

A common mistake is to see a Toeplitz or block pattern and assume that every likelihood calculation is automatically linear-time. The pattern may reduce storage but still require a costly factorization; alternatively, the right model-specific factorization may make the calculation much faster. Check that the proposed method computes the required log-determinant and solve, and state the assumptions behind its complexity.

Another frequent error is to form \(\Sigma^{-1}\) and then compute its entries. A likelihood needs \(r^\top x\), where \(\Sigma x=r\), not the full inverse. Solving with a Cholesky factor is generally more direct and numerically preferable. For a positive-definite covariance, computing the log-determinant from Cholesky diagonal entries also avoids taking the determinant first and then taking its logarithm.