Eigenvalue Decomposition for Exact LMM Scans

Why Rotate the Relationship Matrix?

In a genome-wide association scan, the same phenotype and relatedness structure are used to test many genetic variants. A dense relationship matrix makes direct covariance calculations expensive: factoring or inverting a new dense covariance matrix for every variant would repeat much of the same work. Eigenvalue decomposition exposes a reusable structure. It changes the coordinate system so that the covariance is diagonal, while preserving the exact model and its likelihood-based and association calculations.

The key is to apply the same rotation to every subject-level quantity. Rotating the relationship matrix alone does not transform the analysis; the phenotype, fixed covariates, and each tested genotype must be expressed in that rotated coordinate system too.

Model and Assumptions

Let \(y\) be an \(n\)-vector of phenotypes, \(X\) an \(n \times p\) full-column-rank matrix of fixed covariates, and \(z\) the \(n\)-vector of dosages for one tested variant. Under the null model, use the Gaussian LMM

Here \(K\) is a known, symmetric positive-semidefinite relationship matrix; \(\sigma_g^2\) and \(\sigma_e^2\) are nonnegative variance components; and \(g\) and \(\varepsilon\) are independent. Thus \(y\) has covariance \(V=\sigma_g^2K+\sigma_e^2I\). For the association model, the mean becomes \(X\beta+z\gamma\), where \(\gamma\) is the variant effect per dosage unit. The rotation below assumes that \(K\), the sample, and the residual covariance form used in this model remain fixed. In particular, the stated diagonalization of \(V\) uses an identity residual covariance.

Because \(K\) is positive semidefinite, \(\lambda_i\geq 0\). If \(\sigma_e^2>0\), then \(d_i>0\), so \(D\) is invertible. Orthogonal rotation preserves the information in the data: it changes coordinates, not the observations or the model. The covariance of \(y^*\) is \(D\), and its components are independent under the Gaussian model. The transformed covariates and genotype are not generally independent across components; their purpose is to represent the same design in the new coordinates.

How the Diagonal Weights Arise

Since \(U\) is orthogonal, \(V=UDU^\mathsf{T}\), so \(V^{-1}=UD^{-1}U^\mathsf{T}\). Consequently, generalized least-squares products can be computed using diagonal weights:

The determinant also simplifies: \(|V|=|D|=\prod_i d_i\), because an orthogonal change of coordinates has determinant magnitude one. Thus, for any candidate variance components, both the covariance inverse and log determinant can be evaluated from \(n\) diagonal entries rather than by a new dense factorization.

For a fixed variance-component choice, define \(W=D^{-1}\), \(A=X^{*\mathsf{T}}WX^*\), and \(b=X^{*\mathsf{T}}Wy^*\). The GLS estimate of the fixed covariate coefficients is \(A^{-1}b\), provided \(X\) has full column rank. The generalized residual quadratic is

These formulas are exact identities under the model, not approximations to the dense-matrix calculations. In restricted maximum likelihood (REML), the variance-dependent terms in the negative twice log likelihood are \(\log|V|+\log|X^\mathsf{T}V^{-1}X|+y^\mathsf{T}Py\), with an additional \((n-p)\log(2\pi)\) constant. The rotation gives \(\log|V|=\sum_i\log d_i\) and \(X^\mathsf{T}V^{-1}X=X^{*\mathsf{T}}WX^*\). Variance components may therefore be optimized using diagonal weights; the eigendecomposition is reused as the candidate variances change.

Association Calculations for Each Variant

With the null covariance specified or fitted, the rotated coordinates provide an efficient conditional association calculation. Let \(A=X^{*\mathsf{T}}WX^*\), \(c=X^{*\mathsf{T}}Wz^*\), \(b=X^{*\mathsf{T}}Wy^*\), \(q=z^{*\mathsf{T}}Wz^*\), and \(r=z^{*\mathsf{T}}Wy^*\). After accounting for the fixed covariates, the genotype information and phenotype-genotype product are

If \(\widetilde q>0\), the GLS estimate of the variant effect is \(\widehat\gamma=\widetilde r/\widetilde q\), with standard error \(SE(\widehat\gamma)=1/\sqrt{\widetilde q}\) when the covariance \(V\) is treated as known. If variance components were estimated, the usual plug-in standard error uses their fitted values. The score statistic for testing \(H_0:\gamma=0\) is \(U^2/I=\widetilde r^2/\widetilde q\), where \(U=\widetilde r\) and \(I=\widetilde q\). Under standard large-sample conditions, it is compared with a chi-squared distribution with one degree of freedom. The rotation makes these calculations efficient; it does not make a large-sample reference distribution exact in finite samples.

Why the Decomposition Can Be Reused

The eigendecomposition costs approximately \(O(n^3)\) operations for a dense \(n \times n\) matrix, but it is performed once for a fixed \(K\). For a given variance-component choice, evaluating \(d_i\) and the weights costs \(O(n)\). Products for covariates and a rotated genotype then use diagonal weighted sums. Rotating each new genotype directly still costs \(O(n^2)\); if many genotypes are processed together, matrix multiplication can perform these rotations efficiently. Importantly, the method does not make all per-variant work constant-time. Its benefit is avoiding a fresh dense covariance factorization for each scan test.

A common alternative is to factor or solve with the dense covariance directly. That approach can be numerically convenient for a single analysis and does not require explicitly storing eigenvectors, but repeating the factorization for every variant is wasteful when the covariance is shared. The eigenvalue route is algebraically equivalent under the stated model. Score, Wald, and likelihood-ratio tests remain distinct inferential procedures; diagonalization accelerates their required covariance calculations but does not make them the same test.

Diagnostics, Failure Modes, and Limits

Several checks help distinguish a valid reuse of the decomposition from a coding or modeling error. Numerically, verify that \(U^\mathsf{T}U\) is close to \(I\), that \(U\Lambda U^\mathsf{T}\) is close to \(K\), and that all \(d_i\) are positive before taking inverses or logarithms. Recompute a small number of covariance products directly and compare them with the rotated weighted sums. Material disagreement suggests an eigenvector-ordering, sample-ordering, or transformation error. Tiny negative eigenvalues caused by floating-point error may occur for a theoretically positive-semidefinite matrix; their treatment should be based on scale and numerical precision, not silently interpreted as substantial negative relatedness variance.

The decomposition ceases to apply unchanged if the analysis changes the subjects, the relationship matrix, or the residual covariance structure. For example, a non-identity residual covariance \(R\) generally gives \(\sigma_g^2K+\sigma_e^2R\), which is not necessarily diagonalized by the eigenvectors of \(K\). Likewise, a leave-one-chromosome-out matrix requires its own decomposition if it differs from the matrix used previously. These are model and data alignment issues, not merely computational details.

A further pitfall is to call the resulting scan “exact” and infer that every p-value is an exact finite-sample probability. The eigenvalue transformation is exact linear algebra, conditional on the specified model and inputs. Variance-component estimates, plug-in standard errors, and asymptotic test references have their own inferential properties. Boundary estimates, such as a fitted variance component of zero, also require care in interpretation and testing. Finally, a dense eigendecomposition and storage of \(U\) can be impractical at very large \(n\); the exact dense method is most useful when its one-time cost and memory requirements are manageable.