NExt-LMM: Scalable GWAS Mixed Model

• NExt-LMM is a linear mixed model framework for GWAS that leverages a shared heritability ratio and HODLR approximation to efficiently handle large-scale genetic similarity matrices.
• It minimizes computational bottlenecks by avoiding per-SNP variance component re-estimation, reducing the cost of repeated matrix inversions while retaining low approximation error.
• The method guarantees near-exact inference through convergence of REML estimators and vanishing Kullback–Leibler divergence as the approximation tolerance is tightened.

Searching arXiv for the specified NExt-LMM paper and closely related LMM scalability work. Near-Exact Linear Mixed Model (NExt-LMM) is a linear mixed model framework for genome-wide association studies (GWAS) that targets the computational bottlenecks of variance-component estimation and repeated large-matrix operations while retaining near-exact statistical behavior. It is introduced in “A near-exact linear mixed model for genome-wide association studies” (Pu et al., 7 Aug 2025), where the core design combines two elements: a shared heritability ratio across SNP tests and a Hierarchical Off-Diagonal Low-Rank (HODLR) approximation of the genetic similarity matrix (GSM) and its induced covariance. In the paper’s formulation, these choices reduce the cost of GWAS-scale inference while preserving low approximation error, with theoretical guarantees stated in terms of convergence of REML estimators and vanishing Kullback–Leibler divergence as the HODLR approximation tolerance tends to zero (Pu et al., 7 Aug 2025).

1. Model formulation and statistical setting

NExt-LMM is built on the standard single-SNP linear mixed model for a quantitative trait $y$ across $n$ individuals,

$$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$

where $X$ includes an intercept, covariates, and the tested SNP, and $K$ is the $n\times n$ genetic similarity matrix. The marginal covariance is parameterized as

$$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$

with narrow-sense heritability

$$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$

This is the parameterization adopted by NExt-LMM for continuous-trait GWAS (Pu et al., 7 Aug 2025).

For fixed covariance $\Sigma$, the fixed effects are estimated by generalized least squares,

$$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$

and the Gaussian log-likelihood is

$n$0

NExt-LMM reparameterizes $n$1 using a shared heritability-derived ratio $n$2, writing $n$3, and then uses

$n$4

The article’s central claim is that the combination of a shared $n$5 and fast HODLR inversion makes these estimators computationally tractable at GWAS scale without materially altering inferential output (Pu et al., 7 Aug 2025).

A GSM in this setting is typically constructed from a standardized genotype matrix $n$6 as

$n$7

while the paper also writes $n$8 when $n$9 denotes the standardized genotype matrix. The stated results are intended to apply under either normalizing convention (Pu et al., 7 Aug 2025).

2. Shared heritability ratio and maximum-likelihood acceleration

A defining feature of NExt-LMM is the use of a single population-level heritability estimate for a phenotype, shared across all SNP tests. The paper estimates heritability using PCGC and converts it to a variance-ratio parameter through

$$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$0

This removes the need to re-estimate variance components for each SNP, which is identified as one of the dominant costs in conventional LMM-based GWAS (Pu et al., 7 Aug 2025).

For standardized $$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$1 under the additive model for continuous traits, the paper states

$$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$2

PCGC estimation is then defined by least squares over off-diagonal entries,

$$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$3

with $$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$4 complexity. The paper emphasizes that this estimation uses only off-diagonal elements of $$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$5 (Pu et al., 7 Aug 2025).

Once $$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$6 is fixed, each SNP test proceeds by forming the corresponding design matrix, computing $$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$7, estimating $$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$8, and then constructing a Wald statistic for the SNP coefficient. The variance of $$y \;=\; X\beta \;+\; u \;+\; \varepsilon,\quad u \sim \mathcal{N}(0,\sigma_g^2 K),\quad \varepsilon \sim \mathcal{N}(0,\sigma_e^2 I),$$9 is

$X$0

and the Wald statistic for SNP $X$1 is given in the paper as

$X$2

with a p-value obtained via the standard normal, or a $X$3-distribution with $X$4 degrees of freedom. Score and likelihood-ratio tests are described as available in the same framework, but the implementation is stated to center on fast GLS/Wald testing enabled by HODLR (Pu et al., 7 Aug 2025).

This design directly contrasts with approaches that repeatedly optimize variance components during the scan. A plausible implication is that NExt-LMM trades per-SNP variance-component flexibility for a shared, phenotype-level estimate whose adequacy depends on the empirical stability of heritability across tested markers.

3. HODLR approximation and numerical linear algebra

The second core component of NExt-LMM is the use of the HODLR format to represent matrices whose off-diagonal blocks are well approximated by low-rank factors. The motivation is empirical: GSMs are reported to exhibit strong low-rankness in off-diagonal blocks, with population stratification, relatedness, and linkage disequilibrium producing correlation patterns that vary slowly across individuals when ordered by ancestry or principal components (Pu et al., 7 Aug 2025).

For a symmetric matrix $X$5, the HODLR representation recursively partitions $X$6 into $X$7 blocks. At level $X$8,

$X$9

and the diagonal blocks are recursively partitioned further. The rank $K$0 is adapted blockwise to satisfy a user-specified Frobenius-norm tolerance $K$1, and the recursion terminates at a minimum block size $K$2 (Pu et al., 7 Aug 2025).

HODLR supports direct inversion by blockwise Schur complements. With

$K$3

the inverse is written as

$K$4

with off-diagonal terms retained in low-rank form. The paper states that for positive definite matrices such as $K$5, all Schur complements remain well-conditioned (Pu et al., 7 Aug 2025).

The supported computational primitives are central to NExt-LMM’s complexity claims. With off-diagonal rank $K$6,

• inversion and linear solves scale as $K$7;
• matrix–vector products with HODLR matrices scale as $K$8;
• log-determinant evaluation via recursion has complexity comparable to inversion.

NExt-LMM constructs a HODLR approximation of $K$9 and applies a HODLR inverse algorithm directly to $n\times n$0. The paper explicitly notes that it does not rely on Woodbury identities; the speedups are attributed instead to hierarchical block inversions with low-rank off-diagonal approximations (Pu et al., 7 Aug 2025).

Once $n\times n$1 is available in HODLR form, per-SNP quantities such as $n\times n$2, $n\times n$3, residual quadratic forms, and test statistics are computed using $n\times n$4 matrix–vector products. This suggests that the principal numerical advantage is realized after the one-time covariance construction and inversion step.

4. Algorithmic workflow and asymptotic complexity

The paper summarizes NExt-LMM by four high-level steps: compute $n\times n$5 from standardized genotypes; estimate $n\times n$6 via PCGC and set $n\times n$7; form $n\times n$8 and compute $n\times n$9 via HODLR-inverse with tolerance $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$0 and minimum block size $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$1; then, for each SNP, build the corresponding design matrix and compute $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$2, $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$3, and a Wald statistic or p-value (Pu et al., 7 Aug 2025).

The HODLR-inverse subroutine is described recursively. If the matrix size is at most $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$4, it returns a direct inverse. Otherwise it partitions the matrix into four blocks, recursively inverts $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$5, forms the Schur complement $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$6, recursively inverts $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$7, approximates off-diagonal blocks by low-rank factors $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$8 satisfying $$\Sigma \;=\; \sigma_g^2 K \;+\; \sigma_e^2 I \;=\; \sigma_e^2(\lambda K + I), \quad \lambda \equiv \sigma_g^2/\sigma_e^2,$$9, and assembles the inverse from the block formulas (Pu et al., 7 Aug 2025).

The computational profile reported in the paper is summarized below.

Operation	Complexity	Stated note
Build $$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$0	$$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$1 time, $$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$2 memory	Dense storage baseline
HODLR storage	$$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$3	After compression
Invert $$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$4	$$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$5	HODLR inversion
Estimate $$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$6 by PCGC	$$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$7	Off-diagonal entries only
Per-SNP scan	$$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$8	Constant number of HODLR matvecs

Over $$h^2 \;=\; \frac{\sigma_g^2}{\sigma_g^2 + \sigma_e^2} \;=\; \frac{\lambda}{1+\lambda}.$$9 SNPs, the total GWAS cost is stated as approximately

$\Sigma$0

The paper compares this to three established GWAS LMM implementations: FaST-LMM with $\Sigma$1, BOLT-LMM with $\Sigma$2, and Fast-GWA with $\Sigma$3 for sparse GSMs. NExt-LMM is characterized as near-linear in $\Sigma$4 and $\Sigma$5 when $\Sigma$6 remains modest (Pu et al., 7 Aug 2025).

In this respect, NExt-LMM belongs to a broader line of scalable LMM work that replaces exact dense linear algebra with structured approximations. A closely related but methodologically distinct example is “Scalable Algorithms for Learning High-Dimensional Linear Mixed Models” (Tan et al., 2018), which uses dual estimators and the subsampled randomized Hadamard transform (SRHT) to obtain near-exact LMM estimation with explicit approximation guarantees. The two approaches share an emphasis on provable approximation error, but NExt-LMM’s mechanism is HODLR compression of the covariance structure rather than randomized sketching of kernel matrices (Tan et al., 2018).

5. Theoretical guarantees and the meaning of “near-exact”

The term “near-exact” in NExt-LMM is tied to explicit approximation theory rather than to heuristic empirical similarity. The central assumption is that there exists a HODLR approximation $\Sigma$7 such that

$\Sigma$8

with $\Sigma$9 controlled by the user through the low-rank tolerance (Pu et al., 7 Aug 2025).

Under this approximation regime, the paper states that the REML estimators $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$0 computed from the approximate model converge to the exact LMM REML estimators $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$1 as $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$2, where $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$3. Consequently, $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$4 (Pu et al., 7 Aug 2025).

The Gaussian-model discrepancy is quantified through the Kullback–Leibler divergence between the exact and approximate distributions with equal mean $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$5,

$$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$6

where

$$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$7

The paper proves that this divergence tends to zero as $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$8. It then interprets the result as implying arbitrarily small discrepancies in the log-determinant and quadratic forms that drive ML, REML, and downstream test statistics (Pu et al., 7 Aug 2025).

A practical interpretation is also stated. With $$\hat{\beta} \;=\; (X^{\!T}\Sigma^{-1}X)^{-1}X^{\!T}\Sigma^{-1}y,$$9 set to a small tolerance, exemplified by $n$00 in experiments, the errors in $n$01, $n$02, $n$03, and $n$04 are described as negligible for GWAS purposes and empirically consistent with exact baselines (Pu et al., 7 Aug 2025).

This notion of near-exactness is narrower than exact algebraic equivalence but stronger than approximate speedup without guarantees. It places NExt-LMM in the same broad category as scalable LMM methods that preserve likelihood geometry up to controlled error, although the specific guarantee here is based on HODLR approximation and KL convergence rather than the SRHT spectral bounds used in (Tan et al., 2018).

6. Empirical performance, practice, and scope

The empirical evaluation in (Pu et al., 7 Aug 2025) includes simulation studies and a real-data GWAS-scale application. In Simulation 1, covariance matrices of the form $n$05 were built for $n$06 and $n$07. Element-wise comparisons between HODLR-based inverses and NumPy-based dense inverses are reported to align on $n$08, with near-identical marginal distributions. Mean absolute errors of inverse entries are reported as below tolerance, specifically $n$09 across all scenarios, with mild growth in MAE as $n$10 increases (Pu et al., 7 Aug 2025).

Runtime in the same simulation is reported to exhibit near-linear scaling. At $n$11, HODLR-inverse is stated to be more than $n$12 faster than a cubic baseline and approximately $n$13 faster than NumPy’s optimized inverse (Pu et al., 7 Aug 2025).

In Simulation 2, NExt-LMM is compared against FaST-LMM, Fast-GWA, and BOLT-LMM for $n$14, $n$15, and $n$16. Time ratios defined as competitor over NExt-LMM exceed 1 across all tested $n$17, and at $n$18, NExt-LMM is reported as $n$19–$n$20 faster. Statistical power, measured by AUC, is reported as near parity, with competitor-to-NExt-LMM AUC ratios approximately 1 across $n$21, which the paper interprets as maintained type I error control and power (Pu et al., 7 Aug 2025).

The real-data demonstration uses SARS-CoV-2 RNA sequences with $n$22 and $n$23. Strong population structure is reported from GSM heatmaps and PCA. NExt-LMM identifies 110 genome-wide significant loci at $n$24, compared with 108 hits for FaST-LMM, and the shared hits with Fast-GWA and BOLT-LMM are given as 100–107, corresponding to 89–97% agreement. Runtime is reported as 4.2 minutes and at least $n$25 faster than the alternatives (Pu et al., 7 Aug 2025).

The paper also provides practical guidance. For HODLR tuning, $n$26 is set through the Frobenius-norm threshold on off-diagonal low-rank approximations, with $n$27 reported to deliver near-exact accuracy in experiments. The off-diagonal rank $n$28 is adapted blockwise until $n$29, and small $n$30, with $n$31, is stated to be typical for GSMs. The minimum block size $n$32 should be chosen so that leaf blocks are small enough for fast direct inversion, for example in the tens to low hundreds depending on hardware (Pu et al., 7 Aug 2025).

Several conditions and limitations are stated explicitly. Because $n$33 is positive definite for $n$34, valid Schur complements are guaranteed. The paper recommends centering and standardizing $n$35 and genotypes and including covariates such as PCs and batch effects to mitigate confounding. It also notes situations where HODLR may be less effective: if off-diagonal blocks are not approximately low-rank, $n$36 may grow enough to erode speedups; extremely dense GSMs without exploitable structure may void the advantage, in which case hybrid sparse–HODLR approximations or re-ordering individuals by ancestry or principal components are suggested (Pu et al., 7 Aug 2025).

The current implementation is limited to a single-component random effect based on a GSM, a single phenotype per scan, and Gaussian quantitative traits. Multi-component models, multi-trait models, binary or survival outcomes, and generalized linear mixed models are described as extensions rather than current capabilities. For human studies, LOCO is recommended to mitigate proximal contamination and linkage disequilibrium effects (Pu et al., 7 Aug 2025).

The authors provide a Python package at the repository specified in the paper, with a typical workflow consisting of loading genotypes and phenotype, building and standardizing $n$37, estimating $n$38 via PCGC, constructing $n$39, computing $n$40 via HODLR, and scanning SNPs with GLS and Wald p-values (Pu et al., 7 Aug 2025).