Why Computation Becomes a Modeling Choice
A genomic linear mixed model can be statistically straightforward and computationally demanding at the same time. With many individuals, even storing the covariance matrix of the observations can be burdensome; factoring it may require substantially more work. The practical question is not merely which formula is correct, but how to evaluate it without changing the intended model—or, if an approximation is used, how to understand the change.
This tutorial compares three approaches: dense matrix factorization, low-rank representations or approximations, and iterative linear solvers. They are computational strategies for quantities already defined by the model, not alternative inferential targets. The appropriate choice depends on the number of individuals, the number and structure of markers, the effective rank, and whether the calculation needs exact likelihoods or only selected covariance solves.
For a standardized \(n\times m\) marker matrix \(W\), a common genomic relationship matrix is \(K=WW'/m\). The matrix \(V\) is then positive definite because, for any nonzero vector \(a\), \(a'Va=\sigma_g^2a'Ka+\sigma_e^2a'a>0\). Dense likelihood and REML calculations commonly require solving systems involving \(V\) and evaluating \(\log|V|\). Generalized least squares and REML also require operations involving \(X\), so a fast covariance solve is not necessarily a complete likelihood algorithm.
Dense Factorization: Direct and Predictable
A direct approach forms \(V\) as an \(n\times n\) matrix and computes its Cholesky factorization \(V=LL'\), where \(L\) is lower triangular. Since \(V\) is positive definite under the assumptions above, this factorization exists. It supports stable solves by forward and backward substitution, and \(\log|V|=2\sum_{i=1}^n\log(L_{ii})\). The factor can be reused for multiple right-hand sides at the same covariance parameters.
For a general dense matrix, Cholesky factorization takes order \(n^3\) arithmetic operations and storage is order \(n^2\). These costs are predictable and the calculation is exact up to floating-point rounding. If variance parameters change during likelihood optimization, however, the factorization generally must be recomputed for each new \(V\). A dense approach is attractive for moderate \(n\), particularly when exact likelihood quantities and reliable numerical behavior matter.
Worked Example: Estimating Dense Storage and Factorization Work
Suppose a hypothetical analysis has \(n=6{,}000\) individuals. Using double-precision entries, estimate the storage needed for one dense \(V\), and use the leading Cholesky cost of about \(n^3/3\) floating-point operations.
The matrix has \(6{,}000^2=36{,}000{,}000\) entries. At 8 bytes per entry, its storage is \(36{,}000{,}000\times8=288{,}000{,}000\) bytes, or about 288 MB (about 274.7 MiB), before accounting for the factor, data, and other objects.
The leading factorization work is \(6{,}000^3/3=216{,}000{,}000{,}000/3=72{,}000{,}000{,}000\), or about 72 billion floating-point operations. The order-of-magnitude calculation shows why a mathematically simple dense procedure can become expensive. It is not a precise prediction of run time, which also depends on hardware and implementation.
Low-Rank Structure: Exact Representation or Approximation
The marker representation can avoid explicitly forming \(K\). Set \(a=\sigma_e^2\) and \(\tau=\sigma_g^2/m\), so \(V=aI_n+\tau WW'\). If \(m<n\), the marker-space matrix \(W'W\) is only \(m\times m\). Matrix identities then give an exact alternative to factoring \(V\), provided the smaller computations are feasible. This is an exact re-expression of the model, not a rank approximation.
The inverse identity is the Woodbury matrix identity; the determinant identity is the matrix determinant lemma. Both preserve the full covariance \(V\). If the smaller matrix is factored directly, the cost includes forming \(W'W\), which takes order \(nm^2\) operations, and factoring an \(m\times m\) matrix, which takes order \(m^3\). This can be appealing when \(m\) is much smaller than \(n\). If \(m\) is large, forming and factoring the marker-space matrix may be no better than working in individual space. The dimensions alone do not guarantee a speed-up.
A Short Derivation of the Marker-Space Formulas
For the inverse, apply the Woodbury identity to \(A+UC U'\), with \(A=aI_n\), \(U=W\), and \(C=\tau I_m\). Since \(A^{-1}=a^{-1}I_n\) and \(C^{-1}=\tau^{-1}I_m\), the identity gives the displayed expression for \(V^{-1}\). For the determinant, the matrix determinant lemma states that \(|A+UC U'|=|C^{-1}+U'A^{-1}U|\,|C|\,|A|\). Substituting these matrices and collecting factors yields
Taking logarithms gives the stated log determinant. The identities require no approximation. If \(\sigma_g^2=0\), use \(V=aI_n\) directly rather than the inverse formula containing \(\tau^{-1}\).
A different use of low rank is to approximate \(K\) by retaining only its leading \(r\) eigencomponents: \(K\approx U_r\Lambda_rU_r'\), where the columns of \(U_r\) are orthonormal eigenvectors and \(\Lambda_r\) contains the retained eigenvalues. This changes the covariance model unless the omitted eigenvalues are zero. Retaining the residual term \(\sigma_e^2I_n\) keeps the approximated covariance positive definite, but it does not restore the omitted genetic covariance. The rank \(r\) is therefore both a computational choice and an approximation choice.
Worked Example: Quantifying a Truncated-Rank Approximation
Consider a hypothetical \(K\) with eigenvalues \(4,\ 1,\ 0.1\), and let \(\sigma_g^2=0.5\) and \(\sigma_e^2=1\). The eigenvalues of \(V\) are \(\sigma_e^2+\sigma_g^2\lambda_i\), namely \(3,\ 1.5,\ 1.05\). Thus
If the smallest eigenvalue of \(K\) is discarded, the approximated covariance eigenvalues are \(3,\ 1.5,\ 1\). Its determinant is \(4.5\), and its log determinant is \(\log(4.5)\approx1.5041\). The log-determinant difference is about \(1.5529-1.5041=0.0488\). The approximation has reduced the determinant by \(0.225\), or about \(4.76\%\) relative to \(4.725\). These calculations show the effect on this covariance quantity; they do not establish that the approximation is acceptable for every phenotype or inferential target.
Iterative Solvers: Compute the Solve You Need
When explicitly forming or factoring \(V\) is too costly, an iterative method can solve \(Vx=b\) using repeated matrix-vector products. Conjugate gradients (CG) applies when \(V\) is symmetric positive definite, as it is under the model assumptions above. In exact arithmetic it reaches the solution in at most \(n\) iterations, but large problems are normally stopped earlier when the residual is sufficiently small. Practical convergence depends on the spectrum and conditioning of \(V\), as well as the stopping tolerance and any preconditioner.
The marker representation makes a product with \(V\) possible without forming \(K\) or \(V\):
For dense \(W\), one such product takes order \(nm\) operations. CG then costs roughly the number of iterations times the matrix-vector cost, in addition to preconditioning and data-handling costs. The method can be attractive when the required solves converge in far fewer than \(n\) iterations. A preconditioner aims to improve convergence by making the transformed system better conditioned, but its construction has a cost and must not change the target solution.
An iterative solve is not by itself a complete likelihood calculation. In particular, CG does not directly return \(\log|V|\). Likelihood or REML procedures that need a log determinant require an additional method, such as a stochastic approximation, or a suitable exact factorization. Those additional methods have their own computational cost and, when stochastic, approximation error. Similarly, GLS calculations involving several right-hand sides require multiple solves or a block strategy.
Worked Example: Two Conjugate-Gradient Iterations
For a small hypothetical covariance matrix, let \(V=\begin{pmatrix}2&1\\1&2\end{pmatrix}\) and solve \(Vx=b\) for \(b=(3,0)'\). This matrix is symmetric positive definite: its eigenvalues are 1 and 3. Start CG with \(x_0=(0,0)'\), residual \(r_0=b-Vx_0=b\), and search direction \(p_0=r_0\).
The first step length is \(\alpha_0=(r_0'r_0)/(p_0'Vp_0)=9/18=0.5\), since \(Vp_0=(6,3)'\). Thus \(x_1=x_0+\alpha_0p_0=(1.5,0)'\), and \(r_1=r_0-\alpha_0Vp_0=(0,-1.5)'\). The direction update uses \(\beta_0=(r_1'r_1)/(r_0'r_0)=2.25/9=0.25\), giving \(p_1=r_1+\beta_0p_0=(0.75,-1.5)'\).
Now \(Vp_1=(0,-2.25)'\), so \(\alpha_1=(r_1'r_1)/(p_1'Vp_1)=2.25/3.375=2/3\). Therefore \(x_2=x_1+\alpha_1p_1=(2,-1)'\). Direct substitution verifies the answer: \(V(2,-1)'=(3,0)'=b\). The example illustrates a covariance solve; it does not imply that two iterations will suffice for a large genomic model.
Choosing Among the Approaches
| Approach | Computational feature | Useful setting | Main caution |
|---|---|---|---|
| Dense factorization | Order \(n^3\) factorization work and order \(n^2\) storage | Moderate \(n\); exact solves and log determinants are needed | Repeated covariance updates may require repeated factorizations |
| Exact marker-space calculation | Uses \(W'W\) and an \(m\times m\) factorization | \(m\) is substantially smaller than \(n\) | Forming \(W'W\) costs order \(nm^2\) |
| Truncated low-rank approximation | Work can fall with retained rank \(r\) | A compact representation is adequate for the intended analysis | Omitted components alter the covariance model |
| Iterative solve | Repeated matrix-vector products, often without forming \(V\) | Large problems where solves converge quickly | Convergence depends on conditioning; log determinants need extra work |
Determine whether the analysis needs covariance solves, a log determinant, many right-hand sides, or all of these.
Assess \(n\), \(m\), and any plausible retained rank \(r\), along with memory and the cost of forming intermediate matrices.
For a truncated representation, assess the consequences of omitted covariance components; for an iterative solve, monitor residuals and convergence under the chosen tolerance.
Common Pitfall and Practical Interpretation
A common mistake is to call every low-rank method an approximation, or to assume that an iterative solver makes the entire likelihood calculation fast. The correction is to distinguish the operations carefully: marker-space identities can be exact; eigenvalue truncation changes the covariance unless discarded components are zero; and CG computes a linear solve, not a log determinant. A computational shortcut is statistically harmless only when it preserves the intended calculation or when its approximation error is understood and judged acceptable.
These distinctions matter in genetic analyses because computational limits can influence which cohorts, covariance components, and likelihood evaluations are feasible. Reported estimates do not become more reliable merely because an algorithm is faster. The calculation must still correspond to the specified mixed model, and approximation or numerical tolerances must be sufficiently controlled for the intended use.
Check Your Understanding
Try these questions before reviewing the main points.
- Under what model conditions is \(V=\sigma_g^2K+\sigma_e^2I_n\) positive definite?
- Why is use of the marker-space determinant identity not necessarily a low-rank approximation?
- If the number of markers \(m\) is close to or greater than \(n\), why might forming and factoring \(W'W\) fail to offer a computational advantage?
- What changes statistically when positive eigenvalues of \(K\) are omitted from a truncated approximation?
- Why does a converged CG solve not, on its own, provide the log determinant required by a likelihood calculation?