Sparse Relationship Matrices for Large Cohorts

Why Sparsify a Relationship Matrix?

A relationship matrix for \(n\) individuals records genetic similarity across pairs. In a large cohort, storing every entry requires memory proportional to \(n^2\), even when most pairs have only weak estimated relationships. Setting sufficiently small off-diagonal entries to zero can reduce storage and speed some matrix operations. But a weak relationship is not necessarily irrelevant: many small entries can collectively affect covariance calculations.

The goal is therefore not simply to create as many zeros as possible. It is to make a controlled approximation, state what has been discarded, and check whether the approximation is acceptable for the intended mixed-model analysis. This tutorial focuses on that trade-off.

Definition and Mixed-Model Context

Using the absolute value matters when a relationship matrix can contain negative entries: thresholding by \(K_{ij}<\tau\), rather than by \(|K_{ij}|<\tau\), would treat negative and positive relationships asymmetrically. The threshold should be applied symmetrically to the pair \(i,j\) and \(j,i\); otherwise the result may not even be symmetric.

For context, consider the linear mixed model

where \(y\) is the phenotype vector, \(X\) is the fixed-effect design matrix, and \(g\) and \(\varepsilon\) are independent random vectors. Assume \(\sigma_g^2\geq 0\), \(\sigma_e^2>0\), and that \(K\) is positive semidefinite, as a valid covariance relationship matrix should be. Then the phenotype covariance is \(V=\sigma_g^2K+\sigma_e^2 I\). Replacing \(K\) with \(K_\tau\) gives \(V_\tau=\sigma_g^2K_\tau+\sigma_e^2 I\). This replacement approximates the covariance model; it does not merely compress the same matrix without changing its values.

How Much Covariance Information Can Be Lost?

Define the discarded matrix \(\Delta=K-K_\tau\). For symmetric hard thresholding, \(\Delta\) is symmetric, has zero diagonal, and contains the removed off-diagonal entries. A useful summary is its maximum absolute row sum, the matrix infinity norm \(\|\Delta\|_\infty\). It accounts for the combined magnitude of all removed relationships involving each individual, rather than judging each entry in isolation.

The spectral norm \(\|\Delta\|_2\) controls the largest possible change in a covariance quadratic form. For any symmetric matrix, \(\|\Delta\|_2\leq\|\Delta\|_\infty\). Since \(V-V_\tau=\sigma_g^2\Delta\), this gives a bound for all directions in the phenotype space.

Positive Semidefiniteness Is Not Automatic

A relationship matrix used as a covariance matrix must be positive semidefinite: for every vector \(a\), \(a^\mathsf{T}Ka\geq0\). Symmetry and nonnegative diagonal entries alone do not guarantee this property. Hard thresholding does not, in general, preserve positive semidefiniteness. Thus a valid \(K\) can produce a \(K_\tau\) with a negative eigenvalue.

If \(K_\tau\) is positive semidefinite, then \(V_\tau\) is positive definite under the stated model assumptions. For any nonzero \(a\),

If \(K_\tau\) has a negative eigenvalue, the resulting \(V_\tau\) may or may not remain positive definite: the residual variance could outweigh the negative component, but this must be checked rather than assumed. A matrix that is not a valid covariance can cause a factorization to fail or yield an invalid fitted model. In the five-person example, \(K_\tau\) is strictly diagonally dominant with positive diagonal: each row’s off-diagonal absolute sum is less than 1. For a symmetric matrix, that condition here ensures its eigenvalues are positive, so \(K_\tau\) is positive definite. That convenient property of this example is not a general guarantee for thresholded matrices.

Diagnostics and Failure Modes

A threshold should be evaluated against the structure of the matrix and the intended model, not chosen solely to reach a desired storage target. Useful checks include:

• Inspect the retained degree distribution. Record the number of retained neighbors per individual. A small average can conceal individuals whose rows remain dense or others whose relationships have nearly all been removed.
• Summarize discarded row sums. Examine \(\|\Delta\|_\infty\), as well as the distribution of absolute discarded row sums. Many individually small entries can add to a substantial row total.
• Check symmetry and eigenvalues. Confirm that thresholding was symmetric and assess whether the resulting matrix is positive semidefinite. Failed Cholesky factorization can be a warning, though it is not by itself a complete diagnosis of the source of the problem.
• Vary the threshold. Compare the resulting sparsity, covariance summaries, and relevant fitted quantities across plausible thresholds. Large changes under small threshold adjustments indicate that the approximation deserves caution.
• Check computational behavior. Measure actual memory and runtime for the intended operations. Sparse storage alone does not establish that the mixed-model fit will be faster.

A particularly common mistake is to reason that “each deleted relationship is tiny, so their total effect must be tiny.” This ignores accumulation across a row and the direction of the phenotype contrast being considered. The correction is to examine aggregate discarded relationships and, where useful, a bound such as \(\sigma_g^2\|\Delta\|_\infty\), while remembering that it is a worst-direction bound rather than a complete description of every fitted association result.

What Thresholding Does—and Does Not—Replace

A dense relationship matrix retains all estimated pairwise information and avoids this particular thresholding approximation, but its storage scales quadratically with cohort size. A sparse thresholded matrix reduces storage and can accelerate operations that exploit sparsity, at the cost of removing covariance information and potentially creating a non-positive-semidefinite matrix. Low-rank approximation, introduced in the earlier material on low-rank approximations to relationship matrices, compresses the matrix by representing dominant components rather than by deleting selected weak entries. The two approximations discard information in different ways; neither should be assumed to preserve all covariance features.

Thresholding is most defensible when the retained pattern has a meaningful sparse structure, the removed aggregate contributions are acceptably small for the analysis, and the resulting covariance is valid. Its limitations are important: it can distort long-range or diffuse relatedness, is sensitive to the threshold, and does not guarantee faster factorization. The best evidence is a combination of matrix diagnostics and sensitivity analysis, not sparsity alone.