Stable LD Matrix Computation for Summary-Statistic Methods

Why LD Matrix Stability Matters

Many summary-statistic methods estimate variant weights by relating marginal association statistics to joint effects through linkage disequilibrium (LD). In a common standardized formulation, \(m\) is the vector of marginal effect estimates and \(R\) is the LD correlation matrix for the included variants. A basic linear relation is \(m \approx Rw\), where \(w\) is a vector of joint effects or weights. Solving this relation directly requires an invertible, numerically stable matrix. If two variants are nearly redundant, small perturbations in the entries of \(R\) or \(m\) can produce large changes in the resulting weights.

The issue is not only whether software returns an answer. An ill-conditioned matrix can yield weights that depend strongly on rounding, small changes in the reference sample, or details of the numerical solver. Regularization can make computation more stable, but it generally changes the estimated weights as well. It is therefore important to identify the source of instability and state what the chosen remedy does.

Where Instability Comes From

An LD matrix is usually computed as a correlation matrix of genotype columns after allele alignment and consistent variant filtering. In exact arithmetic, a correlation matrix is symmetric and positive semidefinite: for every vector \(a\), \(a^\mathsf{T}Ra\geq 0\). It is positive definite only when the included genotype columns are linearly independent in the relevant sample. If the reference panel has \(n_{\mathrm{ref}}\) individuals, centered genotype columns have rank at most \(n_{\mathrm{ref}}-1\). Thus, when the number of included variants exceeds that rank, the matrix must be singular. Even when there are fewer variants, duplicates, near-duplicates, or strong LD can make the smallest eigenvalue very small.

These facts assume that the matrix is built from the same aligned variants and the same reference individuals, with valid genotype measurements and consistent scaling. They do not guarantee that the reference panel represents the target population; that is a separate question. Nor does a matrix that looks symmetric and positive semidefinite guarantee good conditioning. A matrix can satisfy both properties and still have eigenvalues spanning many orders of magnitude.

Before solving a system, inspect the dimensions and variant order, check that \(R\) is approximately symmetric and has the expected diagonal, and examine its eigenvalues or singular values. Tiny negative eigenvalues can arise from rounding, but substantial negative eigenvalues suggest a construction, alignment, or data problem rather than ordinary near-singularity. For a proposed solution \(\hat w\), also check the residual \(r=m-R\hat w\). A small residual is useful but not conclusive: a sensitive system can have a small residual and still yield weights that change substantially under small input perturbations.

Stabilizing the System: Filtering and Regularization

One response is to remove variants that are redundant or poorly measured, then recompute the matrix and solve the reduced system. This can improve rank and conditioning, but changes the set of modeled variants. Variant removal should follow defensible rules, such as excluding duplicate records or choosing one representative from a highly correlated group, rather than being driven solely by a desire to obtain convenient weights. Clumping is one familiar way to reduce correlated variants for particular score constructions; it is not identical to stabilizing a dense LD matrix for a shrinkage method.

A second response is ridge regularization: replace the system \(Rw=m\) with \((R+\lambda I)w=m\), where \(I\) is the identity matrix and \(\lambda>0\). If \(R\) is symmetric positive semidefinite, its eigenvalues are nonnegative, and those of \(R+\lambda I\) are \(d_j+\lambda\). The regularized matrix is positive definite, with condition number \((d_{\max}+\lambda)/(d_{\min}+\lambda)\). This is no larger than \(d_{\max}/d_{\min}\) when \(d_{\min}>0\), and is strictly smaller when \(d_{\max}>d_{\min}\). If the original matrix is singular, adding positive \(\lambda\) makes the smallest eigenvalue positive. Regularization reduces sensitivity but shrinks the solution along low-eigenvalue directions; \(\lambda\) is therefore a modeling choice as well as a numerical one.

For the positive-semidefinite case, the reason is direct: adding \(\lambda\) increases every eigenvalue by the same amount, while the ratio of largest to smallest becomes less extreme. If \(0\leq d_{\min}\leq d_{\max}\), then \((d_{\max}+\lambda)/(d_{\min}+\lambda)\leq d_{\max}/d_{\min}\) when \(d_{\min}>0\), because cross-multiplication reduces the inequality to \(d_{\min}\lambda\leq d_{\max}\lambda\). Equality holds when \(d_{\min}=d_{\max}\); if \(d_{\min}=0\), the unregularized matrix is singular and has no finite condition number.

A Small Derivation of Weight Sensitivity

A perturbation bound makes the role of conditioning precise. Suppose \(R\) is nonsingular, \(m\ne0\), and \(w=R^{-1}m\). Let the computed or perturbed inputs be \(R+\Delta R\) and \(m+\Delta m\), and let \(\hat w=(R+\Delta R)^{-1}(m+\Delta m)\). Define \(\eta=\|R^{-1}\|_2\|\Delta R\|_2\), and suppose \(\eta<1\). Then \(R+\Delta R\) is nonsingular. Indeed, \(R+\Delta R=R(I+R^{-1}\Delta R)\), and a matrix \(I+A\) is invertible when \(\|A\|_2<1\).

Subtracting the original equation from the perturbed one gives \((R+\Delta R)(\hat w-w)=\Delta m-\Delta Rw\). Also, \(\|(R+\Delta R)^{-1}\|_2\leq\|R^{-1}\|_2/(1-\eta)\), by the geometric-series bound for the inverse of \(I+R^{-1}\Delta R\). Since \(\|m\|_2=\|Rw\|_2\leq\|R\|_2\|w\|_2\), and \(m\ne0\) implies \(w\ne0\), it follows that \[ \frac{\|\hat w-w\|_2}{\|w\|_2} \leq \frac{\kappa_2(R)\,\|\Delta m\|_2/\|m\|_2+\eta}{1-\eta}, \qquad \eta=\kappa_2(R)\frac{\|\Delta R\|_2}{\|R\|_2}<1. \]

Thus a large \(\kappa_2(R)\) can amplify relative perturbations in either input. The relative bound is stated for \(m\ne0\), so its denominators are defined. A version that also covers \(m=0\) is the absolute inequality \(\|\hat w-w\|_2\leq\|R^{-1}\|_2(\|\Delta m\|_2+\|\Delta R\|_2\|w\|_2)/(1-\eta)\), under the same condition \(\eta<1\). These bounds are worst-case guarantees, not predictions of the exact error in a particular data set.

Diagnostics, Alternatives, and Failure Modes

A stable workflow separates checking the matrix from deciding how to modify it. Check that variants appear in the same order in \(R\) and \(m\); that alleles and genotype coding are aligned; that entries are finite; and that the diagonal and symmetry match the construction. Inspect the spectrum and report the condition number or smallest eigenvalues. After solving, evaluate the residual and, where feasible, compare solutions under modest changes in precision or reference sampling. A solver warning, a very large condition number, or weights that change sharply under small perturbations calls for investigation rather than automatic acceptance.

Three common approaches address different problems. Filtering reduces redundancy by changing the variant set. Ridge regularization retains the variants but changes the linear system and shrinks unstable directions. A pseudoinverse sets sufficiently small singular directions to zero, while eigenvalue clipping replaces problematic eigenvalues with a chosen floor; both define altered solutions or matrices and require an explicit threshold. If the matrix is materially indefinite, first investigate its construction. Symmetrizing \(R\) by replacing it with \((R+R^\mathsf{T})/2\) can remove tiny asymmetry from rounding, but cannot repair incorrect variant alignment or a badly constructed reference matrix.

For a singular system with a consistent right-hand side, infinitely many solutions may exist; for an inconsistent right-hand side, no exact solution exists. A pseudoinverse returns a particular minimum-norm least-squares solution, not a uniquely identified biological effect vector. Similarly, adding a ridge penalty produces a unique answer but does not create information about directions unsupported by the data. Regularization parameters, filtering thresholds, and eigenvalue floors can materially affect weights, so they should be documented and checked for sensitivity.

Common Mistakes and Practical Interpretation

• Calling a matrix stable because it is invertible. A smallest eigenvalue close to zero can make an invertible matrix highly sensitive. Report conditioning, not just whether an inverse exists.
• Treating ridge as a numerical fix with no inferential effect. The solution to \((R+\lambda I)w=m\) differs from the unregularized solution and depends on \(\lambda\).
• Interpreting large opposing weights as robust variant-specific evidence. Near-collinear variants can have unstable weights even when their combined contribution is more stable.
• Using a pseudoinverse without stating its threshold. The choice of which small singular values to discard affects the answer.
• Repairing a bad matrix without checking its source. Eigenvalue clipping or forced symmetry cannot correct a variant-order mismatch or allele-coding error.

The practical goal is not to make every weight look reasonable by any available adjustment. It is to identify whether instability comes from redundant variants, limited reference-panel rank, numerical precision, or a matrix construction error; choose a remedy suited to that cause; and describe the resulting estimand honestly.