Eigenvalue Decomposition for Efficient LMM Likelihoods

Why diagonalize the covariance?

In a single-kernel linear mixed model, genetic similarity between individuals makes their phenotypes correlated. The covariance matrix is usually dense, so evaluating likelihoods by repeatedly manipulating that matrix can be costly. An eigendecomposition offers a useful alternative: when the covariance consists of a genetic relationship matrix multiplied by a variance component, the kernel’s eigenvectors also diagonalize the full covariance. In those coordinates, the likelihood is a sum of contributions from scalar variances rather than a calculation requiring a general matrix inverse at every parameter value.

This change of coordinates does not alter the statistical model or discard observations. It separates the likelihood into directions of genetic variation, each with its own variance. The same representation also helps explain how fixed effects are fitted and how a candidate SNP is tested while accounting for polygenic covariance.

The model and its eigendecomposition

After integrating out \(g\), the marginal model is \(y\sim N(X\beta,V)\), with \(V=\sigma_g^2K+\sigma_e^2I\). The assumption \(\sigma_e^2>0\) ensures that \(V\) is positive definite even when \(K\) has zero eigenvalues. Because \(K\) is symmetric, it has an orthogonal eigendecomposition \(K=U\Lambda U^\mathsf{T}\), where \(U^\mathsf{T}U=I\) and \(\Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n)\). Positive semidefiniteness means each \(\lambda_i\geq 0\).

Rotate the observations and design matrix into the eigenvector coordinates: \(z=U^\mathsf{T}y\) and \(W=U^\mathsf{T}X\). Orthogonal rotation preserves the information in the data, but changes how that information is expressed. In particular, \(z\sim N(W\beta,D)\), so the coordinates \(z_i\) are independent conditional on the fixed effects, with means \(w_i^\mathsf{T}\beta\) and variances \(d_i\). Here \(w_i^\mathsf{T}\) is row \(i\) of \(W\).

Likelihood as a sum of weighted terms

For any fixed \(\beta,\sigma_g^2,\sigma_e^2\), the Gaussian log likelihood is \(-\tfrac12[n\log(2\pi)+\log|V|+(y-X\beta)^\mathsf{T}V^{-1}(y-X\beta)]\). The decomposition gives \(|V|=|D|=\prod_i d_i\) and \(V^{-1}=UD^{-1}U^\mathsf{T}\). Substituting \(z=U^\mathsf{T}y\) and \(W=U^\mathsf{T}X\) yields:

Each direction contributes a log-variance term and a squared residual weighted by the reciprocal variance. Directions with larger \(d_i\) contribute less to the quadratic term for the same residual size. The formula is exact, not an approximation. Its direct use requires the covariance to have the stated single-kernel form and the kernel eigendecomposition to be valid; if a separate covariance component is added, the same eigenvectors need not diagonalize the resulting covariance.

Fitting fixed effects and restricted likelihood

For fixed variance components, maximizing the likelihood over \(\beta\) is weighted least squares in the transformed coordinates. Define \(D=\operatorname{diag}(d_1,\ldots,d_n)\); then the generalized least-squares estimate is:

This requires \(W\) to have full column rank, which follows from full column rank of \(X\) because \(U\) is invertible. The expression weights each transformed residual according to its variance. Substituting \(\hat{\beta}\) into the log likelihood gives the profile likelihood for the variance components; optimizing it over \(\sigma_g^2\geq0\) and \(\sigma_e^2>0\) is one route to maximum-likelihood estimation.

Restricted maximum likelihood (REML) accounts for estimating fixed effects by using error contrasts: linear combinations of the data that remove \(X\beta\). If \(X\) has rank \(p\), the restricted log likelihood, up to constants independent of the variance parameters, is:

where \(P=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}\). In eigenvector coordinates, the matrix terms are obtained from \(D\) and \(W\): \(|V|=|D|\), \(X^\mathsf{T}V^{-1}X=W^\mathsf{T}D^{-1}W\), and \(y^\mathsf{T}Py=z^\mathsf{T}[D^{-1}-D^{-1}W(W^\mathsf{T}D^{-1}W)^{-1}W^\mathsf{T}D^{-1}]z\). Thus REML also uses the scalar variances \(d_i\), while still accounting for the fixed-effect design.

Using the representation for SNP association

A candidate SNP can be included as an additional fixed effect while the single-kernel covariance accounts for polygenic correlation. Let \(X\) contain the baseline fixed effects under the null hypothesis, and let \(x\) be the candidate genotype dosage vector. Under the alternative, the mean is \(X\beta+x\gamma\), where \(\gamma\) is the phenotype change per dosage unit, conditional on the model. For fixed covariance parameters, define \(P\) as above. The generalized least-squares estimate and its model-based standard error are:

These expressions require \(x^\mathsf{T}Px>0\): after accounting for the baseline design, the genotype must contain information not explained by \(X\). With the variance components treated as known, the Wald statistic \(\hat{\gamma}/SE(\hat{\gamma})\) has a standard normal distribution under \(H_0:\gamma=0\), and its square has a chi-square distribution with one degree of freedom. In practice, variance components are often estimated under the null and substituted into these calculations; the resulting test uses an approximation whose calibration depends on the model assumptions.

In transformed coordinates, replace \(x,y,X\) by \(U^\mathsf{T}x,z,W\), and use \(D^{-1}\) for the covariance weighting. For example, the null-model score numerator is \(x^\mathsf{T}Py\), equivalently the weighted transformed association after removing the baseline fixed effects. The eigendecomposition therefore supports association testing without treating genetically correlated individuals as independent.

Why the transformation matters—and a common pitfall

Once \(U\) and the kernel eigenvalues are available, likelihood evaluation at each proposed pair of variance components reduces to updating \(d_i=\sigma_e^2+\sigma_g^2\lambda_i\), summing \(\log d_i\), and forming weighted sums. The same diagonal weights feed the fixed-effect and association calculations. This is especially useful when the same kernel is used across many likelihood evaluations or candidate-variant tests. The transformation does not make the initial eigendecomposition free; its practical value is that the decomposition can be reused as the variance components change.

A key distinction is between the transformed observations and residuals after fitting fixed effects. The transformed observations \(z_i\) are independent under the Gaussian model conditional on their means. But estimating \(\beta\) imposes shared constraints across coordinates, so the fitted residuals \(z-W\hat{\beta}\) are generally correlated. It is incorrect to treat those fitted residuals as independent observations with no adjustment for fixed effects.

• Wrong reasoning: “The covariance is diagonal after rotation, so fitting the intercept can be done separately for each coordinate.”
• Correction: The likelihood has independent coordinate contributions before fitting, but the intercept is a shared parameter. Estimate it jointly by weighted least squares using \(W^\mathsf{T}D^{-1}W\).
• Wrong reasoning: “A SNP test in the rotated data is just an unweighted correlation of \(U^\mathsf{T}x\) and \(z\).”
• Correction: The test must use the weights \(D^{-1}\) and remove the baseline fixed effects. The relevant quantities are \(x^\mathsf{T}Px\) and \(x^\mathsf{T}Py\).