Multivariate sparse-pedigree evaluation and REML

gs_multitrait_pedigree_records() stores one row per observed scalar response. Rows sharing an observation-unit label belong to one animal and may contain any nonempty subset of the declared traits. A trait can occur once per unit; missing traits are absent and never imputed. Units are residual-independent, while traits observed together use the corresponding principal submatrix of the full SPD residual covariance R.

ped <- gs_pedigree_data(
  individual = c("child", "sire", "dam"),
  sire = c("sire", NA, NA), dam = c("dam", NA, NA))
records <- gs_multitrait_pedigree_records(
  ped,
  animal = c("sire", "sire", "dam", "child", "child"),
  trait = c("height", "weight", "height", "height", "weight"),
  observation_unit = c("sire", "sire", "dam", "child", "child"),
  phenotype = c(1.2, 2.0, .9, 1.8, 2.7),
  fixed = Matrix::Matrix(cbind(
    height_intercept = c(1, 0, 1, 1, 0),
    weight_intercept = c(0, 1, 0, 0, 1)), sparse = TRUE),
  trait_levels = c("height", "weight"))
start <- gs_multitrait_pedigree_start(records)
control <- gs_multitrait_pedigree_reml_control(
  linear_solver = "sparse_direct", direct_backend = "cholmod",
  exact_trace_method = "selected_inverse",
  maximum_selected_inverse_workspace_bytes = 2 * 1024^3,
  trace_mode = "complete_deterministic",
  maximum_direct_coefficient_dimension = 256,
  max_exact_trace_systems = 256, trace_batch_size = 32,
  requested_threads = 8)
fit <- gs_multitrait_pedigree_fit(
  records, start = start, covariance = "optimized", control = control)

The deterministic start fits the active fixed columns separately by trait and sets both starting covariance matrices to one half of the diagonal residual phenotypic variances. It does not estimate correlations and is not RHE. Explicit SPD G and R override it.

Three capability levels are intentional:

  • Large-scale fixed SPD G,R: sparse-PCG MME/BLUP, fixed effects, all-animal breeding values, record fitted/residual/Py, and stochastic fixed-parameter REML diagnostics.
  • Exact-work-limited sparse-direct estimation: one lower-CSC pattern and symbolic PARDISO or CHOLMOD analysis, deterministic coordinate traces in bounded retained-factor batches or one CHOLMOD supernodal selected-inverse derivative-support request per full trace evaluation, exact factor log determinants, and likelihood-monotone log-Cholesky optimization. Preflight reports coefficient dimension, input/factor nonzeros, trace systems, RHS workspace, fill, and factor memory, and enforces the configured limits.
  • Scalable stochastic estimation: PCG, fixed caller-controlled counter-Rademacher trace/determinant probes and seeds, common probes across accepted and trial states, standard batched AI, and noise-aware log-Cholesky optimization. Raw trace/score SEs are reported diagnostics; step admission and convergence use scale-aware curvature/Newton and standardized-score rules. An unresolved state returns insufficient_stochastic_precision without direct fallback.

The direct accepted state already contains the exact final likelihood. Rejected trials perform no trace or AI work, and final fixed effects and BLUPs reuse the accepted response coefficients without another factorization or solve. The result identifies fixed versus optimized covariance, trace mode, optimizer qualification, exact-work preflight, numerical and stochastic status, starting/final covariances, scores and AI, likelihood and SE, all-animal BLUPs, optional record values, retained-state reuse, observed patterns, sparse storage, PCG work, or direct symbolic/numerical factorization, solve, fill, memory, exact-trace, and thread-domain diagnostics.

Stochastic optimization returns separate optimizer, likelihood, and BLUP statuses. Score/AI convergence promotes the accepted covariance state and allows retained-response all-animal BLUP independently of final SLQ. If the requested log-determinant SE is not met, likelihood_available is false and the typed likelihood status is likelihood_precision_insufficient; the actual estimate, SE, probe/Lanczos configuration, and work remain diagnostic only. Likelihood-dependent model comparisons must require likelihood_available = TRUE. Exact sparse-direct REML is preferred for small or moderate problems within its resource limits, especially when an exact likelihood is required.

The qualified stochastic route retains standard batched AI. The separate augmented-AI experiment, component-LMM stochastic routes, and factor-analytic qualification boundaries are recorded in model qualification.

exact_trace_method = "coordinate" is the unchanged default and works with PARDISO or CHOLMOD. "selected_inverse" requires a resolved CHOLMOD backend and a sufficient maximum_selected_inverse_workspace_bytes limit. Derivative support is compiled once from a deterministic positive-definite planning state; one trace plan, derivative mapping, and selected-inverse workspace are reused across covariance states. It computes the exact inverse values on the coalesced covariance-derivative support, reports request/factor/workspace/peak-memory accounting, and performs zero coordinate trace systems. It does not remove the separately counted AI derivative-response systems. Explicit PARDISO, or automatic resolving to PARDISO, returns not_supported without fallback.

For direct fits, automatic prefers PARDISO. PARDISO retains qualified internal threading; CHOLMOD orchestration and multitrait support work remain 1/1 while selected-inverse dense kernels use their independently requested oneMKL thread count. PCG supports fixed-covariance solves and the stochastic REML route described above; both report zero direct work. An incompatible PARDISO selected-inverse request performs one symbolic capability resolution and then returns not_supported before numerical factorization.

The accepted factor, response, determinant, likelihood, and selected values are retained so BLUP finalization performs zero additional response solves. The shared typed serializer only translates native fields and performs no factorization, solve, trace, likelihood, score, AI, selected inversion, or BLUP work. In direct diagnostics, sparse_direct$requested is the logical direct-route flag; requested_backend and resolved_backend are backend strings.

The combined Windows process loads official Intel libiomp5md and no GNU libgomp; the CHOLMOD-only process loads no OpenMP runtime.

No route constructs dense pedigree, observation covariance, incidence, Kronecker, projection, or MME matrices. Exact covariance boundaries, PARDISO selected inversion, MUMPS, fixed-effect rank repair, general multiple relationship components, longitudinal direct REML through this API, RHE, marker effects, genomic or ssGBLUP models, and serial dependence between observation units are excluded.

Exact conditional PEV/reliability is specified separately in the prediction-error contract, including CHOLMOD selected inversion and bounded PARDISO retained-factor target solves.