REML through packed ssSNPBLUP equations

The experimental standalone YAML interface can estimate unstructured genetic and residual covariance matrices using the same conditional packed PCG equations as prediction. Natural covariance matrices in YAML are starting values. Numerical computation belongs to greml; the Python frontend validates inputs, prepares files and labels outputs. Existing R interfaces are unchanged.

Model and computation

The initial route supports an ordinary pedigree, one genetic channel per trait, and a shared unstructured residual covariance across independent observation blocks. Missing traits are integrated using the observed submatrix of the current residual covariance. Fixed effects and additional identity random effects are retained; covariances of those additional effects remain fixed. Marker centering, normalization, genomic scaling, ridge/blend and any pedigree adjustment are held fixed throughout REML. This is conditional on the supplied genomic construction, rather than simultaneous estimation of its base parameters.

The method is EM-REML. Conditional genetic coordinates have independent prior blocks with row precisions Qnn, the genomic residual precision and marker identity, all sharing the genetic covariance K. Each update uses posterior second moments of those coordinates, including coefficient covariance from the fixed-effect equations. Dividing their contracted moments by the total number of genetic and marker levels gives the genetic EM update. Residual updates use posterior residual moments, with missing components integrated conditionally. Fixed covariance blocks retain their supplied natural values exactly.

The exact pedigree factor, ancestral reduction and packed genotype provider are reused across covariance states. The current implementation copies the sparse topology/maps when making a new numeric state; it does not refactor the pedigree. Record precisions, right-hand sides and the PCG preconditioner are refreshed. Qnn reuse avoids another pedigree factor for each new preconditioner; the shifted-pedigree option does build its covariance-dependent preconditioner factor. The existing final prediction solve independently checks the fitted model.

The fit does not assemble a global equation matrix, full inverse, dense genomic relationship or marker cross-product matrix. The optional exact likelihood reference materializes bounded small coefficient matrices. Source records and trait-pattern metadata, sparse maps and factors, coefficient vectors and the selected record cache still consume memory. Old/new numeric states can coexist during refresh; this is not a bound on process RSS.

YAML controls

Use the complete packed REML example (models/multitrait-packed-reml.yml). The main estimation controls are:

analysis:
  task: reml
  estimate: [genetic, residual]
  optimizer:
    maximum_iterations: 200
    score_tolerance: 1.0e-6
    step_fraction: 1.0
  traces:
    method: exact
    maximum_solves: 10000

estimate selects whole covariance blocks; omit a block to fix it. Packed EM does not yet support fixing individual Cholesky coordinates. step_fraction is the EM update fraction in (0,1], not a likelihood line-search control.

Exact traces solve every conditional coefficient unit column at each covariance state. They are intended for small reference problems, not industrial data. For a bounded stochastic calculation, replace the trace controls explicitly:

  traces:
    method: stochastic
    probes: 64
    seed: 4401
    maximum_solves: 10000
    maximum_score_mcse: 0.001

The stochastic option uses a common deterministic Rademacher probe bank across covariance states and streams one query at a time. The example numbers above are syntax examples, not a calibrated industrial recommendation. Score error is in natural variance/covariance coordinates and depends on data and scaling; the probe count alone does not establish adequate precision.

Each state uses one response solve plus probes stochastic solves, or one plus the conditional coefficient dimension for exact traces. maximum_solves caps their total across all states. It excludes the initial setup and subsequent independent final prediction solve, which have separate timing/iteration metrics. Every REML solve must pass both transformed and original augmented equation residual checks. The internal solve uses a numerical margin; the requested acceptance tolerances are unchanged.

Convergence requires the maximum absolute free score plus twice the maximum free-score Monte Carlo standard error to be no larger than score_tolerance. The MCSE must also meet maximum_score_mcse. Exact traces have zero sampling MCSE, subject to PCG numerical error. A stochastic moment matrix that is not SPD is rejected without clipping, jitter, silent probe increases or an exact fallback.

A stochastic candidate must pass one final score/precision check using an independent probe bank. Its additional response and probe solves count against the same total cap, and history marks the validation state. Failed independent validation rejects the fit; its probes never guide another update or select another validation bank. This check can add a qualification state after the last permitted covariance-update state. Score convergence does not establish uniqueness or statistical identifiability of an arbitrary covariance model.

Results and qualification

Results include variance_components.csv, parameters.csv, reml_history.csv, likelihood.json and convergence.json, plus native natural-score and score-MCSE vectors. History records scores, MCSE and cumulative solve counts. Converged fits then produce the ordinary prediction tables. Failed fits retain diagnostic parameter states and history but publish no qualified prediction tables.

Likelihood reporting is optional, as described below. It does not control EM updates or replace score convergence. Exact EM has its usual likelihood monotonicity, but stochastic moment updates carry no likelihood-ascent certificate. Reported score MCSE describes trace estimation at the evaluated parameter state; it does not include variance-component estimation uncertainty or provide confidence intervals for the fitted parameters.

Focused tests compare natural scores with an independent dense restricted likelihood and finite differences, verify a joint analytic covariance optimum, check missing-trait EM updates and fitted predictions, and exercise stochastic work/precision rejection. This is not a 50K or industrial runtime qualification. EM can require many covariance updates, each with many PCG solves; prediction timings cannot be reused as REML runtime estimates. The native reml stage covers every covariance evaluation, final stochastic check and optional likelihood diagnostic. reml_likelihood reports the determinant diagnostic separately. Existing fine-grained data/operator profile counters are not an aggregate REML profile; use that stage and the explicit solve totals when assessing cost.

FA, serial/grouped residuals, genetic groups, metafounders, coordinate constraints, additional covariance estimation, and joint PEV output are outside this initial packed REML interface. The existing compiled sparse-direct REML implementation retains its broader covariance/constraint support and likelihood reporting.

Optional restricted-likelihood change

Omit analysis.likelihood for zero global determinant/operator-product work. To report the restricted-likelihood change from the supplied starting covariances to the final fitted covariances, add one block under analysis:

  likelihood:
    method: slq
    maximum_operator_products: 6144

Defaults are probes: 32, seed: 4402, lanczos_steps: 32, check_steps: 64, maximum_mcse: 1.0 and maximum_depth_change: 0.5. All can be set explicitly. These illustrate controls, not calibrated production precision. The work cap must cover 2 * probes * (lanczos_steps + check_steps) equivalent vector products: both start/final states at both depths. It is independent of the REML trace-solve cap. Earlier Lanczos breakdown can reduce actual products.

The matrix-free stochastic Lanczos quadrature (SLQ) diagnostic uses common Rademacher probes at the two covariance states and both depths, with a seed stream separate from the score traces. Paired differences give an MCSE for the likelihood change. depth_change is the absolute change in the estimate when increasing Lanczos depth; it is a sensitivity check, not a truncation error bound. A small sampling MCSE alone cannot qualify Lanczos accuracy.

For a small exact reference, use:

  likelihood:
    method: exact
    maximum_dimension: 256

This materializes and factors the start/final conditional coefficient matrices, using 2 * coefficient_dimension products. The explicit dimension cap protects against accidentally selecting a dense reference on a large model. It is a caller-selected work/shape limit, not an application RAM ceiling.

Both methods evaluate -0.5 times the change in log|R_observed| + (n + marker_columns) log|K| + log|C| + y'Py. C is the full conditional coefficient matrix, including fixed and nuisance effects. Observed residual submatrices retain missing-trait covariance. The genetic normalization includes latent marker levels. Other prior and coordinate constants cancel because data, design, genomic construction and fixed additional covariances are unchanged. A stable residual-plus-prior quadratic evaluates y'Py. The initial numeric model is retained for this comparison, so requesting likelihood can increase memory.

likelihood.json and run.json.estimation.likelihood report scale: relative_to_start, value, mcse, depth_change, operator products, seconds, message and a status:

• not_requested: no diagnostic requested.
• not_attempted: covariance fitting did not reach convergence.
• qualified: the diagnostic completed and, for SLQ, met the requested targets.
• unqualified: a finite SLQ estimate missed one or both precision targets.
• failed: a numerical/resource diagnostic failed; no value is published.

Unqualified/failed optional likelihood leaves a separately qualified covariance fit and predictions accepted. No probe increase or exact fallback occurs. This is a start-to-final diagnostic, not an absolute likelihood, per-iteration ascent certificate, or evidence for a global optimum. Do not compare its raw value across models with different starting covariances, fixed designs, data or genomic construction. The compiled sparse-direct route reports its existing absolute restricted likelihood with scale: absolute in the same output block.

The complete exact example (models/multitrait-packed-reml.yml) and SLQ example (models/multitrait-packed-reml-slq.yml) share inputs. Small tests compare exact changes with an independent marginal likelihood (including nonzero markers and missing traits), and paired SLQ estimates/MCSE with an eigenvalue oracle. No industrial likelihood runtime or precision qualification follows from those checks.