Working with large data

Status — Stable

Mirrors drmTMB's Working with large data. How DRModels.jl stays fast as the number of units grows, and what to reach for when a model is large.

The selling-point model — the q=4 phylogenetic bivariate location–scale fit — is built for sparse tree structure. The marginal likelihood uses a sparse augmented-state Laplace approximation, with an exact sparse outer gradient at a converged inner mode. It never forms the dense p×p phylogenetic covariance; the gradient uses a Takahashi selected inverse rather than differentiating a dense factorisation. The selected-inverse cost depends on Cholesky fill; the measured near-linear scaling below applies to the stated tree shapes and settings.

What the timing evidence shows

Two synthetic timing studies define the current evidence boundary. In a single-thread Julia sweep with four observations per tip, balanced and caterpillar trees at p = 100, 1,000, and 10,000 gave empirical wall-time exponents of 0.90 and 0.91. That is evidence of near-linear growth for those tree shapes and settings, not a runtime guarantee for a new dataset.

A separate paired benchmark used a synthetic near-balanced ultrametric tree, four observations per tip, one discarded warm-up, two timed repetitions, single-threaded linear algebra, and drmTMB 0.6.0 on the Totoro server:

Tips	DRModels.jl median	drmTMB median	drmTMB / DRModels.jl
100	0.529 s	1.752 s	3.31
1,000	6.921 s	5.581 s	0.81
5,000	82.141 s	41.996 s	0.51
10,000	115.178 s	104.869 s	0.91

DRModels.jl was faster at 100 tips in this comparison; drmTMB was comparable or faster at the larger sizes. These figures describe one machine, package version, tree shape, and model. They are not a general speed ratio, and the older extrapolated “N× faster” claim is retired.

Why it scales

• Sparse precision, never dense covariance. The phylogenetic prior precision is sparse (3N − 2 stored non-zeros for a tree with N nodes, about 6p for a binary tree with p tips). The engine factorises that sparse matrix with CHOLMOD; it never materialises the dense Σ.

• Exact sparse outer gradient. At the converged inner mode, the implicit-function gradient reuses the selected inverse entries needed by the sparse Cholesky factor. Its cost depends on factor fill, so the exact gradient does not imply a universal O(p) runtime bound. The timing study above measures near-linear growth for its balanced and caterpillar tree designs.

• A precision sampler for uncertainty. A latent-state draw uses the same sparse precision (Cov(û) ≈ P⁻¹) rather than a dense covariance. Total bootstrap time still depends on the number of refits and their convergence.

Practical tips for large fits

• Stay in ML. ML is the default and is comparable across fixed-effect structures — keep it for model selection on large data. REML is an option, not the default.

• Standardise covariates. Good conditioning matters more as p grows; centre and scale continuous predictors so the optimiser's Hessian stays well-behaved.

• Thread the bootstrap and the profile CIs. Parametric bootstrap replicates are independent refits; profile-likelihood endpoints are independent per coefficient. confint(fit; method = :profile, threads = true) profiles coefficients in parallel when the objective is thread-safe — set JULIA_NUM_THREADS to engage it. On the coupled location–scale random-effect route (NB2 or Gamma with a shared (1 | p | group) term on mu and sigma), each coefficient's lower and upper endpoint searches stay serial, so the gain scales with the number of coefficients profiled; profiling one coefficient there is no faster with threads. On other thread-safe routes a single profiled coefficient runs its two endpoint searches in parallel.

• Parallel backend: threads in one Julia process. threads = true uses Julia threads inside a single process; start Julia with --threads=N or set JULIA_NUM_THREADS. Concurrent multi-threaded BLAS calls can be slower than a serial run, so set BLAS to one thread before threading (using LinearAlgebra; BLAS.set_num_threads(1)); the threaded bootstrap does this for you. DRModels.jl has no Distributed.jl (addprocs / @everywhere / pmap) interface. To spread a simulation grid over many cores or machines, run independent Julia processes (for example one cluster job per cell or seed) and combine their saved results.

• Check the fit cheaply. check_drm reports convergence, the gradient and covariance conditioning without re-fitting — useful before committing to an expensive bootstrap. Its converged flag alone is not proof of a maximum; read max_abs_grad and grad_source with it (see Improving convergence).

• A single fit is single-threaded by design. As in drmTMB, one drm() fit runs on one thread, and drm() has no threads option. The parallel speed-ups are across independent work: bootstrap_ci(...; threads = true) refits replicates in parallel, and confint(...; method = :profile, threads = true) profiles coefficients in parallel where the model allows it (see the profile tip above). Both need Julia started with several threads, for example by setting JULIA_NUM_THREADS.

Beyond the verified engine

The O(p) machinery lives in the verified phylogenetic engine. The non-Gaussian GLMM paths (Poisson/NB2/Beta/Gamma random effects via quadrature) are designed for moderate group counts rather than p = 10,000-scale phylogenies; for very large structured problems, the phylogenetic location–scale engine is the path that has been benchmarked to scale.

Gaussian crossed random effects

Gaussian models with crossed random intercepts, such as y ~ x + (1 | rater) + (1 | item), use a sparse factorisation of the random-effect system, so time grows far more slowly with the number of rows than it did before DRModels.jl#977 was fixed (earlier, 2,000 rows took over five minutes). On a synthetic design shaped like the lme4::InstEval ratings, one machine and one thread, a fit took about 0.3 s at 2,000 rows, 45 s at 20,000 rows and 28 s at 100,000 rows. drmTMB is still faster at 20,000 rows, because the optimiser needs many steps on weakly identified crossed designs, not because each step is slow. For very large crossed data sets, compare timings with drmTMB or glmmTMB before committing.

Related software

Latte.jl and GaussianMarkovRandomFields.jl offer related sparse-precision and Laplace-approximation machinery in Julia. DRModels.jl keeps its own engine and uses them as a comparison and an independent cross-check. On a 50-group Gaussian location-scale test model, Latte.jl's TMB-style Laplace engine, started from suitable values, matched the four DRModels.jl maximum-likelihood estimates to within 0.024 standard errors (the largest absolute difference divided by its standard error). Latte.jl's TMB-style engine reports a posterior mode under hyperpriors rather than a maximum-likelihood fit, so the agreement is a cross-check on one model, not a general equivalence.