Why the Computational Representation Matters
A linear mixed model (LMM) can account for correlation among individuals through a relationship matrix, but computing with that matrix can become a major cost in a large cohort. The same statistical model may be implemented using a dense matrix, a sparse matrix, or a low-rank approximation. These choices determine storage, runtime, and—when the representation is approximate—the extent to which computed association results differ from those under the intended model.
The central decision is not simply “which method is fastest?” It is whether the chosen representation preserves the covariance structure needed for the analysis while fitting within available memory and computation time. Dense calculations can provide a useful exact reference. Sparse calculations are attractive when sparsity is genuine and the required factorization remains sparse. Low-rank methods can reduce cost when the relationship structure is well represented by a small number of components, but introduce approximation error that must be assessed.
Model and Computational Choices
Let \(y\) be the \(n\)-vector of phenotypes, \(X\) the fixed-effect design matrix, and \(K\) a positive semidefinite relationship matrix for the \(n\) individuals. A standard single-random-effect LMM is
Here \(u\) and \(e\) are independent, and \(\sigma_e^2>0\). Thus the marginal covariance matrix is \(V=\sigma_g^2K+\sigma_e^2I\). The choice of computational strategy concerns how \(K\), \(V\), or operations involving them are represented and calculated. It does not change the scientific model unless an approximation replaces the intended \(K\).
A full dense \(n\times n\) double-precision matrix uses about \(8n^2\) bytes. This estimate excludes workspace, factorization storage, and other objects used by the software. Exploiting symmetry can reduce storage, but does not remove the steep growth in computation: a general dense factorization has cost on the order of \(n^3\), while storing the matrix remains on the order of \(n^2\).
For a sparse matrix, storage depends on the number of nonzero entries, but runtime also depends on how those entries are arranged and on fill-in. Fill-in occurs when a factorization creates nonzero entries where the original matrix had zeros. Consequently, the fraction of nonzero entries alone does not establish that a sparse strategy will be efficient.
For a low-rank representation \(K\approx U_r\Lambda_rU_r'\), with \(U_r\) having \(r\) columns, storage for the factors is on the order of \(nr\), rather than \(n^2\). A common alternative is \(K\approx BB'\), where \(B\) is \(n\times r\). When \(r\ll n\), matrix identities can reduce computations involving \(V\) to smaller \(r\times r\) systems. The representation is useful only if its approximation is adequate for the task.
Assumptions and a Practical Decision Rule
Before choosing a strategy, check that the relationship matrix corresponds to the intended set and ordering of individuals, is symmetric, and is positive semidefinite up to small numerical error. The covariance model also assumes a valid variance-component specification; in particular, \(V\) is positive definite when \(\sigma_e^2>0\). The fixed effects and phenotype must be aligned with the rows of \(K\). Violations of these conditions are data or model problems, not reasons to switch matrix formats.
The following process separates computational feasibility from statistical adequacy.
Calculate \(8n^2\) bytes for one full matrix, then allow additional memory for factorization and working objects. If this exceeds available memory, a dense implementation is not feasible without a different computational resource or approach.
If \(K\) is exactly sparse, determine whether the software can exploit that pattern and whether expected factorization fill-in is manageable. If sparsity comes from dropping small entries, treat the resulting matrix as an approximation.
Estimate a feasible rank \(r\), its storage cost, and the approximation error relevant to the planned association analysis. A small rank is computationally appealing but is not, by itself, evidence of adequate accuracy.
Where possible, compare it with a more accurate calculation on a manageable subset or selected variants. Examine covariance approximation and changes in association estimates or test statistics, using tolerances specified before final analysis.
Worked Examples
Worked Example: Is a Dense Calculation Feasible?
Suppose a hypothetical cohort has \(n=12{,}000\) individuals, and the available memory for the analysis is 64 GB. The relationship matrix is not known to have useful exact sparsity. Estimate the storage for one dense double-precision matrix and decide whether a dense strategy is plausible.
The matrix has \(12{,}000^2=144{,}000{,}000\) entries. At 8 bytes per entry, its storage is
That figure is for one matrix only. The analysis may need additional arrays and factorization workspace, and a software implementation may temporarily hold several large objects. Still, a 1.15 GB dense matrix is plausible within 64 GB, provided those additional costs fit. The cohort size alone does not rule out dense computation. A practical choice is to benchmark the intended software and retain a dense calculation as an exact reference if its full memory requirements are acceptable.
Worked Example: When Exact Sparsity Supports a Sparse LMM
Consider a hypothetical cohort of \(n=60{,}000\) individuals in 300 independent groups of 200. Suppose the specified \(K\) is exactly block diagonal, with a full \(200\times200\) block for each group and zeros between groups. Estimate storage for its nonzero entries using 8-byte values and 4-byte column indices, plus 4-byte row pointers.
There are \(300\times200^2=12{,}000{,}000\) nonzero entries. The values require \(12{,}000{,}000\times8=96{,}000{,}000\) bytes, and the column indices require \(12{,}000{,}000\times4=48{,}000{,}000\) bytes. The \(n+1=60{,}001\) row pointers use \(60{,}001\times4=240{,}004\) bytes. The total is approximately \(144.2\) MB, before workspace.
For comparison, one full dense matrix would require \(8\times60{,}000^2=28{,}800{,}000{,}000\) bytes, or 28.8 GB. Here the block structure is exact, not created by discarding small entries. If the implementation can process the independent blocks and its factorization does not create cross-block fill, a sparse strategy is an especially natural choice. The model still needs checking: the group structure must genuinely match the specified covariance, and the individuals must be ordered consistently.
Worked Example: Evaluating a Low-Rank Approximation
Suppose a hypothetical cohort has \(n=50{,}000\) and a proposed relationship representation uses \(r=300\) columns in \(B\), so \(K\approx BB'\). Storing \(B\) in double precision requires
A full dense \(K\) would require \(8\times50{,}000^2=20{,}000{,}000{,}000\) bytes, or 20 GB. Thus the factor representation can substantially reduce storage. Suppose, for illustration, that a covariance approximation \(\widetilde V\) has been checked to satisfy \(\|V-\widetilde V\|_2\leq0.02\), and that \(\sigma_g^2=0.6\), \(\sigma_e^2=0.4\). Because both relationship matrices are positive semidefinite, the smallest eigenvalue of each covariance is at least \(0.4\). The inverse-matrix bound derived below then gives
This bound is a worst-case guarantee on inverse-covariance error, not a guarantee that every association statistic changes by a small amount. The rank-300 representation is therefore a promising computational candidate, but should be evaluated by comparing association quantities with a more accurate calculation on relevant test cases. If it fails a pre-specified tolerance, increase the rank or use another strategy.
Why Approximation Error Needs Context
A useful justification for checking covariance error is the identity
Taking the matrix 2-norm and using the fact that the norm of a product is no greater than the product of the norms gives
If both \(K\) and its approximation are positive semidefinite and \(\sigma_e^2>0\), the smallest eigenvalue of each covariance is at least \(\sigma_e^2\). Hence \(\|V^{-1}\|_2\leq1/\sigma_e^2\) and \(\|\widetilde V^{-1}\|_2\leq1/\sigma_e^2\), which yields
This bound explains why a covariance approximation can affect calculations involving the inverse covariance. Association procedures under a fitted LMM use such covariance-dependent quantities, so an approximation should be assessed for its effect on the actual outputs of interest, not only by the percentage of variance captured by retained eigenvectors. A small eigenvalue-truncation error can still matter for a particular phenotype or variant; conversely, a conservative matrix bound may be too loose to describe typical changes.
Diagnostics, Alternatives, and Limitations
| Strategy | Most suitable when | Main checks and limitations |
|---|---|---|
| Dense | The cohort and available memory permit it; exact calculations or a reference result are valuable. | Estimate full memory and workspace. Computation and storage grow rapidly with \(n\). |
| Sparse | The specified covariance structure is genuinely sparse and the software exploits its structure. | Check nonzero pattern, factorization fill-in, and runtime. Thresholding a dense matrix changes the model unless the approximation is justified. |
| Low-rank | A small number of components adequately represents the relevant covariance structure, or dense computation is infeasible. | Check approximation error and its effect on association results. A small \(r\) is not proof of adequate accuracy. |
The main alternative to an approximate sparse or low-rank calculation is an exact dense calculation, when feasible. It can serve as a benchmark, but may be too costly at large \(n\). Sparse and low-rank approaches are not interchangeable: sparsity describes where entries are zero, whereas low rank describes dependence among columns or directions. A matrix can be sparse without being low rank, or low rank without being sparse.
Useful warning signs include unexpectedly high memory use, sparse factorization fill-in that approaches dense storage, numerical failures in covariance operations, and material changes in association estimates or test statistics as the approximation is refined. For a thresholded sparse approximation, also check whether dropped entries collectively alter important relationship patterns. For a low-rank approximation, compare more than a single global summary: evaluate representative variants and individuals, including cases where the retained components appear least adequate.
A frequent mistake is to choose the method with the smallest nominal matrix representation and assume it is therefore both faster and statistically acceptable. A sparse matrix can produce a costly factorization, and a low-rank approximation can omit covariance structure relevant to association. The correction is to separate three questions: can the representation be stored, can the required computations exploit it, and is its approximation error acceptable for the analysis? A favorable answer to only the first question is not enough.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- For \(n=10{,}000\), how many bytes does one full dense double-precision relationship matrix require? Give the result in GB using \(10^9\) bytes per GB.
- A matrix has few nonzero entries, but sparse factorization creates extensive fill-in. Why might a sparse implementation still be inefficient?
- Does retaining 99% of the eigenvalue sum in a low-rank approximation guarantee accurate association statistics? Explain.
- State two checks that should precede selecting a sparse or low-rank strategy.
Solutions. (1) \(8\times10{,}000^2=800{,}000{,}000\) bytes, or \(0.800\) GB. (2) Fill-in increases the storage and computation required by factorization, potentially eliminating the advantage suggested by the original nonzero count. (3) No. A global eigenvalue summary does not guarantee that approximation error is small in directions relevant to a phenotype or variant; compare association outputs with a more accurate calculation. (4) Examples include confirming that \(K\) is correctly aligned with the individuals and is positive semidefinite, checking whether sparsity is exact or induced by thresholding, estimating factorization fill-in, and assessing approximation error against a pre-specified tolerance.