Canonical LMM architecture

This section defines the LMM architecture. Current dependency versions are recorded in package metadata; the original architecture baseline used greml 0.66.0. Capability status belongs in the model catalog, reusable machinery in the primitive registry, evidence in model qualification, and established equivalences in model reductions. Historical development notes do not supersede the current capability catalog.

Current authority and status

Current public headers and R exports define what is implemented. Focused tests and committed qualification evidence determine whether a specific model/backend/inference combination is qualified. The catalog uses exactly:

  • qualified: implemented, publicly reachable where claimed, and covered by focused numerical or contract evidence for that combination;
  • implemented_unqualified: code exists, but the complete claimed combination lacks focused qualification;
  • unsupported: the combination rejects explicitly or has no implementation;
  • not_applicable: the concept does not apply to the route.

Code presence alone never establishes qualified. A limitation on one route must not be generalized to the suite.

Canonical scientific specification

The backend-independent specification contains caller-ordered observed records and stable record IDs; response values; compact fixed-coefficient references; random component, level, channel/trait, and basis references; deterministic external domains and coefficient order; relationship covariance or precision contracts; typed covariance layouts; residual blocks and groups; prediction targets; uncertainty requests; and separate resource/execution policy.

Caller record order may change only record-shaped output order. Scientific domains, parameter positions, likelihoods, active sets, coefficients, and uncertainty results are deterministic after alignment by external identity. Production stores short record-level coefficient-reference lists and never materializes a global fixed- or random-effect incidence matrix, dense or sparse. Pedigree routes reuse the committed ordering, inbreeding, Henderson precision, relationship diagonal, and determinant.

Ordering contract

  • Component LMM coefficients are component-major, then level-major, then channel-major.
  • Pedigree random-regression coefficients are animal-major, then basis axis.
  • Pedigree multitrait coefficients retain declared trait order for every external animal.
  • Fixed effects retain validated caller column order.
  • Structured residual values use canonical variance-group order followed by rho or phi; they are not packed SPD matrices.
  • SPD covariance values retain declared channel/trait order and use explicit natural and log-Cholesky transformed layouts.

Every public result returns identities or names; callers do not infer them from undocumented offsets.

Execution and retained-state boundary

Scientific syntax is independent of PCG, CHOLMOD, PARDISO, batching, threading, and memory budgets. Gsolve owns numerical sessions and selected inversion; greml owns model assembly, likelihoods, derivatives, optimization, BLUP, and uncertainty; gsuite owns R validation, lowering, typed serialization, and external-order restoration.

Backend and trace method are distinct. CHOLMOD and PARDISO consume the same sparse MME, but only CHOLMOD supplies selected inversion. Coordinate traces are an exact bounded route. Retained-factor selected_solve target-block uncertainty is another exact bounded route and is not selected inversion. PCG has no retained direct factor. Stochastic Hutchinson/SLQ routes exist only where the catalog identifies them; they are never silent fallbacks for exact sparse-direct fits.

Incomplete-multitrait pedigree stochastic REML constructs deterministic counter-Rademacher trace and determinant banks once per fit and retains common probe identities across accepted and trial covariance states. Greml inference owns phase-aware precision: finite raw trace and score standard errors are diagnostics; nonstationary step resolution is assessed on curvature/Newton scales; convergence requires score compatibility with zero and adequate dimensionless log-Cholesky-scale uncertainty. Natural and transformed AI are validated and canonically symmetrized at the greml inference boundary before one identical matrix is used for rank diagnostics and the gsolve moment solve.

Stochastic covariance convergence, approximate likelihood availability, and BLUP finalization are distinct states. Final SLQ imprecision does not revoke an accepted score/AI optimum or prevent retained-response BLUP, but it makes the restricted likelihood unavailable for reporting or model comparison. Solver diagnostics likewise separate logical PCG systems from attempts, recovered warm-start retries, and unrecovered failures.

Successful direct optimization retains the accepted likelihood, factor, response, score, AI state, and relevant selected-inverse values. BLUP reuses the retained response. Exact coefficient PEV is the random-effect block of the inverse of the complete MME, conditional on fitted covariance parameters. Existing CHOLMOD selected inversion, pedigree trace contractions, longitudinal pedigree uncertainty, and component-LMM uncertainty are current capabilities. Exact PARDISO selected_solve PEV is qualified for the three sparse pedigree routes; PCG exact PEV and PARDISO selected inversion remain unsupported.

Boundaries and structured residuals

The component LMM’s opt-in active-set policy supports an exact zero scalar random component and rank-one two-channel covariance with KKT checks and reactivation. Larger covariance blocks and all residual parameters remain positive-interior. AR(1) and continuous exponential residual blocks share one typed residual-block abstraction; this does not imply ARMA, heterogeneous temporal scales, spatial covariance, or mixtures of residual structures.

Future-development rule

Proposed extensions must begin with an unsupported or implemented_unqualified catalog cell, identify whether the gap is representation, exposure, backend contract, or qualification only, and reuse the registry and reductions. They must not duplicate existing solvers, selected inversion, PEV calculations, component composition, residual blocks, or boundary optimization.

Three layers

Scientific model recipe

The recipe defines the response, fixed terms, random components and incidence, coefficient domains, covariance and residual structures, missingness, requested predictions, identifiability conditions, and exclusions. It must not encode a backend, trace route, batch width, or memory policy.

Canonical compiled representation

The retained compiled route is selected with control = gs_fit_control(...); formula and formulas are aliases and do not select different routes. It requires explicit gvc() components with complete starts and no separate pedigree argument. The general data-frame adapter is described in the pedigree tutorial. For the compiled route, base R model.frame() and model.matrix() construct trait-specific fixed designs. Supported random terms are lowered recordwise as component/member, entity, trait, residual-unit, and unit-weight metadata; no R random-effect incidence matrix is constructed. Scalar lowering now composes at most one pedigree and one independent identity random-intercept domain. Formula order fixes component blocks and parameter order, while data-induced cross-domain MME entries are retained. Canonical scalar records, caller mappings, covariance-component metadata, positive-interior starts, resource limits, and execution controls cross the R/C++ boundary (gsuite/R/zz-formula-reml.R, gsuite/src/compiled_formula.cpp).

Greml’s installed CompiledMmeSpecification, immutable CompiledMmePlan, handle-based CompiledMmeParameterState, and move-only CompiledMmeSession own coefficient ordering, sparse structure, evaluation, optimization, and retained results (greml/include/greml/compiled_mme.hpp).

Execution

Gsolve owns PCG and sparse-direct numerical machinery (gsolve/include/gsolve/pcg.hpp, gsolve/include/gsolve/sparse_direct.hpp). The same compiled representation feeds each qualified backend. Execution owns symbolic analysis, numerical factorization, response/information solves, coordinate traces or selected inversion, determinants, likelihood-monotone optimization, retained state, and BLUP finalization. Gsuite maps owned results to one stable R object and focused S3 methods.

The component PCG route is a second execution of the existing compiled component representation, not a PCG-only model. It applies \(W^{\mathsf T}R^{-1}W\) by gathering and scattering the canonical short record references, adds compiled component precision, and uses a deterministic fixed/per-level channel block-Jacobi preconditioner derived from that same numerical state. The RHS, parameter and coefficient order, member projections, residual units, observed channel patterns, and caller-record restoration are shared with exact CHOLMOD. No global W, Z, incidence matrix, dense relationship, or full MME is constructed for operator application.

Fixed-budget stochastic inference constructs fixed counter-Rademacher trace and determinant plans once per session. Coefficient-space stochastic traces, the existing explicit derivative-response AI, natural-to-log-Cholesky transformation, symmetry/rank/ conditioning checks, phase-aware convergence, optional SLQ, and retained BLUP all consume the canonical plan. Parameter-step Monte Carlo covariance is \(\mathrm{AI}^{-1}\Sigma_g\mathrm{AI}^{-1}\). Under-resolved nonstationary updates use a scale-invariant diagonal Levenberg term, while convergence remains assessed against the original canonical AI. Optimizer convergence, likelihood availability, BLUP finalization, and solver status remain distinct. Missing-channel inference uses this same path for arbitrary observed channel subsets: missing channels are absent scalar records, and channel-specific fixed references and output mappings remain plan-owned. No inference layer reconstructs those layouts.

An optional adaptive execution policy is available at the component-inference layer while leaving the fixed-budget path as the default. Exploratory, accepted-state, and certification tiers choose deterministic nested counter-Rademacher prefixes and progressively tighter PCG tolerances. Extending a prefix at one covariance state solves only newly added probes and updates retained trace moments and score covariance; matching probe indices preserve common random numbers between accepted and trial states. Certification always uses the predetermined fixed final prefix, never the adaptively stopped intermediate prefix.

Step resolution and convergence remain distinct. Nonstationary directions use the parameter-step covariance \(\mathrm{AI}^{-1}\Sigma_g\mathrm{AI}^{-1}\) and may extend, damp, or fail explicitly. A stationary candidate is certified through deterministic curvature/Newton tolerances, simultaneous standardized scores, the numerical- rank-aware statistic \(g^{\mathsf T}\Sigma_g^{+}g\), and an absolute transformed-parameter uncertainty envelope. It is not rejected merely because uncertainty divided by a vanishing proposed update is large. Trials omit SLQ and result finalization; a potentially acceptable trial is compared at a matched prefix with paired common-probe merit uncertainty. Optimizer, optional likelihood, BLUP, and solver statuses remain separate, and no adaptive route falls back to direct inference.

The qualified stochastic component scope is positive-interior identity and Henderson pedigree relationships with independent residuals, for scalar, complete correlated-channel, and arbitrary missing-channel models. The supplied-precision and structured-residual extension reuses the same operator for caller-certified sparse precision relationships and existing patterned, AR(1), and continuous-exponential local residual blocks. That extension is only partially qualified: supplied precision is fixed-state qualified only; AR(1) fixed-state evaluation and stochastic optimization are qualified; continuous exponential coefficient, derivative, trace, AI, and PCG evaluation is fixed-state qualified, but stochastic optimization remains unqualified after the frozen recovery exhausted its score-merit line search. The combined supplied-precision/structured-residual stochastic illustration was not run. Exact CHOLMOD behavior is unchanged. There is no stochastic-to-direct fallback.

A typed fixed-rank factor-analytic covariance layout is available at the model boundary. For declared trait order and rank, it alone owns the identified lower-triangular loading positions, positive pivots, positive specific variances, natural/transformed names, derived covariance/precision/log determinant, and analytic derivatives. The component compiler consumes the derived covariance through the same parameter and coefficient layouts as an unstructured block; neither PCG nor sparse direct contains an FA branch. Structural derivative planning follows the typed incidence pattern and never uses invented loading magnitudes. Numerical covariance-Jacobian rank and conditioning are checked at the actual scientific state.

Canonical FA average information is the conventional first-derivative quantity. Under the FA map it is \(\mathrm{AI}_{\eta} = J^{\mathsf T}\mathrm{AI}_{\theta}J\). Nonlinear covariance and transformed-coordinate second derivatives are retained only as separately named observed-curvature diagnostics; they are not added to AI, required to be positive definite, or repaired with ridge, jitter, or eigenvalue clipping. Exact CHOLMOD FA REML and retained-state BLUP are qualified. Matrix-free PCG execution and stochastic score/AI are qualified at fixed covariance only. Stochastic FA optimization and its post-optimization retained BLUP remain unqualified; there is no fallback to an unstructured layout or direct inference.

Demonstrated invariants

  1. One canonical lowering and one caller-record construction pass occur per fit; later covariance work does not revisit caller data.
  2. Trait, record, observation-unit, pedigree, coefficient, and natural-parameter order are deterministic and externally mapped.
  3. Model semantics do not change with solver choice.
  4. Every supported backend consumes the same compiled representation.
  5. No backend, trace method, batch size, or stochastic fallback is silent.
  6. Resource rejection precedes symbolic analysis and all numerical work.
  7. Failed native fits preserve status and diagnostics; unavailable numerical output is absent rather than invented.
  8. Successful finalization reuses the retained accepted likelihood, factor, and response and performs one BLUP finalization.
  9. Specialized implementations remain reduction oracles until generic replacements are independently qualified.
  10. Public APIs are added only for a concrete model after the abstraction has survived focused validation.

Evidence includes focused formula-interface tests, simulation-based comparisons with specialized model routes, and compiled-MME contract tests. The bounded independent-component composition is qualified against dense observation-space and dense full-MME oracles, the specialized repeatability route, PARDISO and CHOLMOD. See the qualification evidence for the internal test inventory.

Adding one model

  1. Specify the statistical recipe and explicit exclusions.
  2. Match every requirement to the primitive registry and name the smallest missing primitive or composition rule.
  3. Establish an independent dense oracle or qualified specialized reduction.
  4. Implement only that missing general primitive/composition.
  5. Verify canonical lowering and all mappings exactly.
  6. At common fixed parameters, verify covariance/MME, coefficients, projections, score, standard AI, determinants, and likelihood.
  7. Verify positive-interior optimization and retained BLUP behavior.
  8. Verify each backend, resource limit, thread domain, and failure path.
  9. Expose and document R semantics only after focused validation checks pass.

Structural reductions should omit a component when the chosen route does not support an exact zero boundary. The component-LMM scalar/rank-one boundary contract is specified in the model reductions.

Ownership boundaries

  • gsuite: R semantics, base-R fixed design, bundle validation, lowering, result interpretation, and S3 presentation.
  • greml: pedigree/statistical kernels, canonical MME compilation, REML identities and derivatives, optimizer, resources, and retained BLUP state.
  • gsolve: model-independent operators, PCG, direct factorization, solves, selected-inverse capability, and execution resolution.
  • gbits/gmat/gcorr/gsim: retain their committed data, genomic, correlation, and simulation responsibilities; model extensions retain these boundaries.

Possible future interface

Public compile/explain/preflight functions, public R plan/session objects, and an executable component registry are only possibilities. Names such as gs_compile() or gs_explain() are not installed APIs. A future interface requires a concrete model need, ownership/lifetime design, and focused ABI and statistical qualification.