Low-Rank Approximations to Relationship Matrices

Why Approximate a Relationship Matrix?

A genomic relationship matrix \(K\) summarizes genetic similarity among \(n\) individuals. In a linear mixed model, it contributes to the covariance of the phenotype, so operations involving that covariance can dominate computation when the cohort is large. Storing a dense \(n\times n\) matrix requires \(O(n^2)\) memory, and a generic dense factorization can require \(O(n^3)\) operations. A low-rank approximation replaces \(K\) by a matrix built from a smaller number of directions in the data. This can greatly reduce storage and some repeated computations, but it also changes the fitted covariance and can change association statistics.

The aim is not simply to make a matrix smaller. A useful approximation retains the relationship patterns that matter for the analysis while controlling the effect of discarded patterns on inference. The mixed-model framework and association tests provide the context: the approximation is to the relationship matrix, and its consequences must be assessed in the covariance model and the test statistic.

Spectral Low-Rank Approximation

Assume \(K\) is a symmetric positive semidefinite relationship matrix for the analyzed individuals. Its eigendecomposition is \(K=U\Lambda U^\mathsf{T}\), where the columns of \(U\) are orthonormal eigenvectors and the diagonal entries of \(\Lambda\) are nonnegative eigenvalues. Order them as \(\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_n\geq0\). Retaining the first \(r\) eigencomponents gives a rank-at-most-\(r\) approximation.

The leading eigenvectors represent directions with the largest eigenvalues in \(K\). Truncation preserves those directions and sets the remaining eigenvalues to zero. It does not generally preserve every entry, every diagonal value, or every small-scale relationship in the original matrix. A low-rank factor \(U_r\Lambda_r^{1/2}\) is often stored instead of forming \(K_r\) explicitly.

In a standard linear mixed model, write \(y=X\beta+g+\epsilon\), with \(g\sim N(0,\sigma_g^2K)\), \(\epsilon\sim N(0,\sigma_e^2I)\), and independent random terms. Thus \(y\) has covariance \(V=\sigma_g^2K+\sigma_e^2I\). Under rank truncation, the approximate covariance is \(V_r=\sigma_g^2K_r+\sigma_e^2I\). This description assumes \(K\) is positive semidefinite, the same individuals and scaling of \(K\) are used in both models, and \(\sigma_e^2>0\). The positive residual variance ensures \(V_r\) remains positive definite, even when \(K_r\) is singular.

Computational Savings and the Woodbury Identity

If a factor \(U_r\Lambda_r^{1/2}\) is stored, its storage is about \(nr\) numbers rather than \(n^2\) numbers for a dense matrix. When \(r\ll n\), that is a substantial reduction. Calculations with \(V_r=\sigma_e^2I+\sigma_g^2U_r\Lambda_rU_r^\mathsf{T}\) can use the Woodbury identity, which expresses the inverse through an \(r\times r\) matrix rather than requiring a dense \(n\times n\) inversion:

This form assumes the retained eigenvalues are positive; zero eigenvalues can be omitted from \(U_r\). It is an algebraic identity for positive \(\sigma_e^2\), not an approximation beyond the replacement of \(K\) by \(K_r\). When the factor is fixed and many calculations are needed, the smaller matrix can make repeated solves much cheaper. The time to obtain the factor itself, data movement, and any costs of estimating variance components still matter.

What Truncation Does to the Covariance

In the eigenbasis of \(K\), the exact covariance has eigenvalues \(\sigma_e^2+\sigma_g^2\lambda_j\). The rank-\(r\) covariance has the same values for \(j\leq r\) and \(\sigma_e^2\) for \(j>r\). Therefore \(V-V_r=\sigma_g^2(K-K_r)\) is positive semidefinite: truncation removes covariance in the discarded directions and does not add covariance there.

The largest absolute covariance eigenvalue error is \(\sigma_g^2\lambda_{r+1}\), and the total squared Frobenius error is \(\sigma_g^4\sum_{j=r+1}^{n}\lambda_j^2\). These are useful matrix diagnostics, but neither alone guarantees that an association test is accurate. The phenotype and tested genotype may align especially strongly with a discarded direction, and a small covariance change can matter for a statistic near a decision threshold.

For each discarded direction \(j>r\), the exact inverse-covariance eigenvalue is \(1/(\sigma_e^2+\sigma_g^2\lambda_j)\), while the approximate one is \(1/\sigma_e^2\). Their difference is

This shows why a large residual variance can dampen the effect of a discarded eigenvalue on the inverse covariance, while a relatively large genetic variance or a substantial omitted eigenvalue can make the discrepancy more consequential. It is not, by itself, a bound on a \(p\)-value or on a test statistic with fixed-effect adjustment.

Assessing Association-Test Error

For a single variant with genotype vector \(x\), a score test under the null model uses the covariance \(V\) and fixed-effect design \(X\). Define \(P_V=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}\), assuming \(X\) has full column rank. With the null-model residual vector \(y\), the score statistic is \(Q=(x^\mathsf{T}P_Vy)^2/(x^\mathsf{T}P_Vx)\), provided the denominator is positive. The approximate test substitutes \(V_r\) to form \(P_{V_r}\) and \(Q_r\). Under the null and the model’s regularity conditions, the score statistic is commonly compared with a chi-square distribution with one degree of freedom.

The two statistics need not be close just because \(K_r\) explains a high proportion of the matrix’s total eigenvalue sum. The test depends on the genotype direction, the phenotype residuals, and the fixed-effect adjustment. A direct assessment compares \(Q_r\) with \(Q\) on representative variants for which an exact or higher-rank calculation is feasible. It is useful to compare effect estimates, standard errors, and \(p\)-values as well, especially near the study’s reporting threshold.

Diagnostics, Alternatives, and Limitations

A practical approximation assessment has both matrix-level and test-level checks. At the matrix level, inspect the retained eigenvalue spectrum, the residual eigenvalues, and quantities such as \(\lambda_{r+1}\) and \(\sum_{j>r}\lambda_j^2\). At the test level, compare exact or higher-rank results against the proposed rank on a subset of variants. Include variants with different allele frequencies and association strengths, and examine discrepancies in statistics and standard errors rather than only whether a significance threshold is crossed.

Several failure modes deserve attention. If eigenvalues decay slowly, a small \(r\) may discard meaningful covariance. If the retained subspace is computed for a different set or ordering of individuals, the factor does not represent the intended \(K\). If variance components are estimated using the approximate model, differences between tests reflect both covariance approximation and changed parameter estimates; comparing with fixed common variance components can help isolate the matrix approximation itself. Finally, a high fraction of total eigenvalue mass retained is not a universal guarantee of accurate association inference.

A dense exact approach preserves the specified relationship matrix but may be computationally infeasible. Low-rank approximation is attractive when a relatively small number of directions capture the covariance features needed for the analysis, but it introduces truncation error. Sparse representations instead exploit many zero or negligible entries; that is a different structural approximation, not another name for spectral truncation. These approaches have different error patterns and computational trade-offs, so the suitable method depends on the matrix structure and the required accuracy.

A common mistake is to argue, “The first \(r\) components explain 99% of the eigenvalue sum, so the association \(p\)-values must be accurate.” The eigenvalue sum is a global matrix summary. An individual test depends on how its genotype and phenotype residual align with the discarded directions, so the conclusion does not follow. The correction is to treat eigenvalue summaries as screening diagnostics and directly assess test-level differences on a suitable validation set. A low-rank approximation is a computational model choice, not an automatic guarantee of calibrated inference.