mmes_control: Modular Numerical Controls for Mixed Models

In sommer: Solving Mixed Model Equations in R

Modular Numerical Controls for Mixed Models

Description

Creates a strictly validated sommer_mmes_control object for the numerical methods used by mmes. Model formulas remain separate. Every module accepts only its documented fields; unknown module or field names are errors. Constructing controls does not fit a model or establish that requested methods are eligible for a particular model.

Usage

mmes_control(formulation = "henderson", task = list(),
             representation = list(), solve = list(), logdet = list(),
             score = list(), curvature = list(), update = list(),
             proposal = list(), convergence = list(), output = list(),
             execution = list(), diagnostics = list())

Arguments

formulation	Equation domain: "henderson" (default) or "direct". Henderson forms coefficient-space mixed-model equations; direct uses observation-space inversion.
task	Named list with mode, objective, and pql. Defaults: mode="fit", objective="reml", and pql=list(maxit=20L, tol=1e-5, estimateTheta=FALSE, warmStart=TRUE). mode is "fit" or "solve"; objective is "reml" or "ml". The nested pql list accepts only maxit, tol, estimateTheta, and warmStart, with the defaults shown. These fields own task selection and outer PQL iteration, not inner covariance iterations.
representation	Named list with compression="off", factorScores="none", fixedRank="error", equations="normal", sparseFactors="auto", and matrixFree="auto" by default, unless corresponding sommer.mme.* options override experimental modes. sparseFactors="auto" enables exact sparse innovation kernels for modified-Cholesky factors of dimension at least 64 and maximum innovation lag three. It also enables sparse SAR/CAR kernels when the factor dimension is at least 500 and at most one percent of adjacency weights are nonzero. Other factors retain their existing path. Use "off" to disable or "force" to enable eligible sparse factor kernels regardless of size. If sommer.mme.saddlePoint resolves to "force", the default equations becomes "saddle-point". compression is "off", "groups", or "classes"; factorScores is "none", "fixed-shape", or "profile"; equations is "normal", "saddle-point", or "conditional" (experimental latent FA equations with independently selected AI/EM updates). fixedRank="error" retains the default error for rank-deficient fixed designs. These errors show the opt-in syntax before potentially long term and column labels so it remains visible when error output is truncated. Opt in with fixedRank="drop" to select a full-column-rank subset after observation filtering, before either the Henderson or direct engine runs. A shared C++ sparse QR removes zero columns and detects linear dependencies on unit-norm columns, using pivot threshold max(nrow(X), ncol(X)) * .Machine$double.eps. COLAMD ordering and numerical pivoting choose the basis; retained columns are then restored to their original order. This need not select the same constraints as ASReml. No ridge penalty or generalized inverse is substituted for the fitted model. Omitted coefficients are returned as zero identification constraints, with unavailable standard errors and tests in summaries. fit$fixedRank records retained/omitted indices, original column names, rank, tolerance, and a sparse relation satisfying X_original = X_retained %*% relation to numerical precision. Direct-engine VarBeta remains in retained-column coordinates to avoid allocating a dense covariance over all original columns. Supported workflows are Gaussian identity-link fits, with normal or saddle-point Henderson equations or direct inversion, without PQL, solve mode, conditional FA, factor scores, observation rotation, or compression. Henderson predictions accept only contrasts estimable under the original fixed model; non-estimable rows produce an explicit error. Wald term tests, likelihood-ratio comparisons, Satterthwaite/Kenward-Roger inference, and emmeans currently require an explicitly full-rank refit. Existing direct-engine prediction/VarU limitations remain unchanged. REML likelihoods from different choices of fixed-effect basis can differ by an additive constant; compare covariance estimates and fitted values, not raw REML likelihoods across different bases. For rank-reduced fits, character classify predictions stop before allocating a Cartesian grid. Use predict(fit, D="locale:cultivar", observed=TRUE) for observation-weighted fixed predictions at combinations present in the retained data, or supply explicit estimable contrasts. The prediction maxMemoryMB safeguard checks dense solve/covariance workspace before allocation; see predict.mmes for weighting and limitations. The remaining fields use experimental mode values "auto", "off", or "force" (logical TRUE/FALSE are accepted). This module owns equation, compression, covariance-factor, and factor-score representations.
solve	Named list. method="auto" selects among "auto", "ldlt", "cholmod", and "pcg". PCG fields default to tol=1e-8, maxIters=0L, preconditioner="diagonal", nystromRank=32L, warmStart=FALSE, recycleRank=0L, recycleMethod="coarse", recycleState=NULL, exportRecycle=FALSE, rhsMethod="independent", blockSize=16L, and residualRefresh=0L. Also accepts the experimental solve requests blockSchur="auto", blockBanded="auto", factorSchurMultiFactor="auto", and dFirstML="auto"; global sommer.mme.* options may override these mode defaults. The method owns checked linear solves, preconditioning, reuse, and eligible elimination. PCG settings apply only when an eligible PCG path uses them.
logdet	Named list with scaling="none" and lanczosSteps=20L by default. Scaling is "none" or "diagonal". This module owns determinant approximation and its correction, not linear-solve accuracy.
score	Named list with traceProbes=8L, estimator="hutchinson", rank=16L, covarianceChecks="auto", and traceBatch="off" by default. The estimator is "hutchinson" or "control-variate". covarianceChecks and traceBatch accept experimental mode values. This module owns covariance derivatives, trace estimation, and covariance checks.
curvature	Named list with method="exact-ai" and reduced="off" by default. The method is "exact-ai" or "hybrid-ai"; reduced accepts experimental mode values. This module owns the AI metric and reduced-direction assembly.
update	Named list with method="ai", recovery="em", blendWeight=NULL, and commonScale="off" by default. recovery is "em" or "pxem"; commonScale accepts experimental mode values. This module owns covariance-coordinate updates, EM/PXEM recovery, blending, and common-scale profiling. method="ai" remains the default optimizer. With representation$equations="conditional", it uses exact REML scores, average-information curvature, safeguarded steps and EM recovery in latent FA coordinates. method="em" selects conditional-moment EM and requires that representation. The earlier method="conditional-fa" selector continues to select conditional EM, not AI. Conditional equations support one FA relationship term with diagonal or FA residuals, including missing and unbalanced observations, using sparse LDLT posterior solves and selected inverse panels. They require Gaussian identity Henderson REML, solve$method="auto" or "ldlt", output=list(inverse=0, covariance="none", coefficientMatrix=FALSE), and no weights, rotation, vcc, other representation requests, or acceleration. Each source must be entirely free or entirely fixed. Conditional AI bounds its dense derivative-RHS and curvature workspace by execution$denseMemoryMB; output$panelSize also bounds projection solve panels. A rank-revealing eigensolve removes numerically null curvature directions. Accepted AI steps satisfy an observed-REML ascent check; failed proposals use an audited EM recovery step, never a silent estimator change. AI convergence requires a small diagonally scaled score using convergence$parameterTol; at numerical bounds it uses the projected (KKT) score. Conditional AI bounds free physical specific variances below by max(1e-8, execution$inverseTol) * var(y), and shared genetic allocations to [max(1e-8, execution$inverseTol), 1-max(1e-8, execution$inverseTol)]. This prevents nearly deterministic latent equations from producing unreliable likelihoods. Starting values must respect these bounds; execution$inverseTol must be less than 0.5 on this path. The floors and active coordinate indices are reported in conditionalFA$ai. EM retains its existing machine-precision floor and likelihood/parameter stabilization; it is not silently changed to this constrained AI optimizer. See the conditional FA section of vignette("sommer.mmes.options").
proposal	Named list with stepWeight=NULL and extrapolation="none" by default. Extrapolation is "none" or "aitken". This module owns proposal scaling and extrapolation. Aitken extrapolation cannot be combined with curvature$method="hybrid-ai".
convergence	Named list with maxIters=30L, likelihoodTol=1e-4, parameterTol=1e-4, accuracy="fixed", maxTraceProbes=NULL, maxLanczosSteps=NULL, logDetTol=1e-3, and scoreTol=1e-3. Accuracy is "fixed" or "adaptive"; NULL stochastic limits resolve to the initial traceProbes and lanczosSteps. This module owns inner iteration convergence and stochastic accuracy requests.
output	Named list with inverse=0L, covariance="full", coefficientMatrix=TRUE, setup=FALSE, selectedInverse="auto", and panelSize=64L by default. inverse is 0 (no final coefficient inverse or fit-time PEV computation), 1 (selected inverse and random-effect PEVs), or 2 (full inverse and random-effect PEVs). Modes 1 and 2 compute PEVs during fitting; mode 1 does not form the complete coefficient-matrix inverse. With mode 0, use postPEV after fitting to compute PEVs or the full inverse, provided coefficientMatrix=TRUE; covariance is "full" or "none"; selectedInverse accepts experimental mode values. This module owns coefficient uncertainty, covariance-parameter uncertainty, retained coefficient-matrix output, and setup output.
execution	Named list with denseMemoryMB=1000, memoryPolicy="speed", gpu=FALSE, gpuMinSize=256L, blasThreads=NULL, inverseTol=1e-6, and seed=0L by default. A non-NULL sommer.mme.denseMemoryMB, sommer.mme.memoryPolicy, sommer.mme.gpu, or sommer.mme.gpuMinSize option supplies its corresponding default. This module owns resource policy, hardware, threading, numerical inverse tolerance, and stochastic seed.
diagnostics	Named list with verbose=TRUE, dateWarning=TRUE, and profile="off" by default. profile accepts experimental mode values. This module owns ordinary fit messages, package-age warnings, and phase profiling.

Details

Module ownership describes mathematical responsibility, not a promise that every cross-module combination is supported. For example, formulation and solver choices are model-dependent; PCG has REML and output restrictions; compression, saddle-point equations, factor-score augmentation, common-scale profiling, and reduced AI each have additional prerequisites. Consult vignette("sommer.mmes.options", package="sommer") and the internal .mmes_capabilities() catalog for the current conceptual decision graph and capability inventory. The catalog is descriptive; fitting validation remains authoritative. The modular API does not imply a kernel rewrite or arbitrary method composition.

task$mode="solve" selects a known-covariance path rather than ordinary variance-component fitting; it requires an eligible Gaussian identity model and has narrower PCG/output options. task$pql$estimateTheta=TRUE requires a negative-binomial family. representation$compression needs eligible matching designs and independent residuals. Saddle-point equations, factor-score representation, sparse factors, and matrix-free PCG each have separate model and backend gates. PCG fitting is a Henderson REML path; selected-inverse output mode 1 is not supported by PCG. The vignette and capability catalog describe these prerequisites, but the fit-time validator decides eligibility.

Within PCG, solve$preconditioner="nystrom" uses solve$nystromRank; solve$recycleState and solve$exportRecycle require solve$recycleMethod="krylov" and a positive solve$recycleRank. logdet$scaling and Lanczos depth apply to stochastic determinant evaluation. score$rank is used by the control-variate estimator. Adaptive trace and Lanczos budgets apply to stochastic PCG accuracy. output$panelSize applies to full inverse panels and conditional FA posterior-moment/projection panels, not ordinary selected-inverse mode.

Experimental mode fields accept "auto", "off", or "force"; logical TRUE and FALSE mean force and off. Package options may supply defaults for these requests. "force" requests a path but does not bypass mathematical, structural, or memory eligibility checks. Nested lists reject unknown names and duplicate or empty names.

Defaults retain the established numerical bundle. Aitken belongs to proposal$extrapolation; hybrid AI belongs to curvature$method; they are mutually exclusive.

Migration from removed numerical arguments

mmes() now takes numerical settings through control. The old formals are not compatibility aliases and must not be passed to mmes(). Map former settings to canonical module paths:

henderson

formulation="henderson" or formulation="direct".

solver, engine

solve$method; pass the complete control object as control=....

nIters, tolParConvLL, tolParConvNorm, tolParInv

convergence$maxIters, convergence$likelihoodTol, convergence$parameterTol, and execution$inverseTol, respectively.

REML, solveOnly, pqlControl

task$objective ("reml" or "ml"), task$mode ("fit" or "solve"), and task$pql.

computeCi, getPEV, returnParam

output$inverse and output$setup. PEV output is controlled by output$inverse; the former getPEV flag has no separate replacement.

blend, blendWeight, stepWeight

update$recovery, update$blendWeight, and proposal$stepWeight.

acceleration="aitken"

proposal$extrapolation="aitken".

acceleration="hybrid-ai"

curvature$method="hybrid-ai". Do not combine with Aitken.

compression, factorScoreAugmentation, covarianceInformation, retainCoefficientMatrix

representation$compression, representation$factorScores, output$covariance, and output$coefficientMatrix.

Former pcg settings

Map solve tolerances, iteration limits, preconditioning, and reuse to solve; determinant scaling and Lanczos depth to logdet; trace probes, estimator, rank, and covariance checks to score; stochastic accuracy budgets to convergence; seed to execution.

Former experimental settings

Map structured representation requests to representation, elimination requests to solve, trace batching to score, reduced AI to curvature, common-scale profiling to update, selected/full inverse requests to output, and profiling to diagnostics.

denseMemoryMB, memoryPolicy, gpu, gpuMinSize, blasThreads

execution fields with the same names.

verbose, dateWarning

diagnostics fields with the same names.

recordDiagnostics, verboseDiagnostics

Removed flags with no direct replacement. Use diagnostics$verbose for fit messages and diagnostics$profile for implemented phase profiling.

With diagnostics$verbose=TRUE, the existing Henderson MME path line includes the coefficient solver, active/inactive latent reduced-rank and Schur variants, and the selected residual engine (including inactive grid and repeated-Kronecker paths). An inactive automatic Schur path is reported as “not selected”, not as proof that the model lacks exploitable structure; disabled controls and unsupported coefficient backends are distinguished.

Value

An object of class "sommer_mmes_control".

Compatibility

Settings are checked against one compatibility table when the control is built and again in mmes. A supplied setting that conflicts with another, or that would have no effect on the selected path, stops the call with a message that names every offending setting and how to fix it. A setting counts as supplied only when it differs from its default, so restating defaults is always accepted. The main rules are:

• formulation="direct" accepts only task, update$recovery, update$blendWeight, proposal$stepWeight, convergence iterations and tolerances, output$setup, output$inverse 0 or 2, memory and threading execution settings, and diagnostics.

• task$mode="solve" accepts only solve$tol, solve$maxIters, solve$warmStart, solve$residualRefresh, a non-Nystrom solve$preconditioner, output$setup, CPU execution settings and diagnostics.

• PCG-only settings (solve tolerances, preconditioner, recycling, logdet, score estimators, adaptive accuracy, representation$matrixFree, execution$seed) require solve$method="pcg"; "auto" never selects PCG.

• solve$method="pcg" rejects EM/PXEM recovery, blending, common-scale profiling, reduced AI, Aitken, hybrid AI, saddle-point equations, factor scores and task$objective="ml".

• Aitken extrapolation, hybrid AI, saddle-point equations, factor scores, blending and output$coefficientMatrix=FALSE each reject the modules they cannot combine with.

Rules that depend on the data, family, vcc or rotation are checked by mmes at fit time; for example a Gaussian identity fit stops when task$pql is changed.

See Also

mmes; vignette("sommer.mmes.options", package="sommer")

Examples

ctrl <- mmes_control(
  formulation = "henderson",
  solve = list(method = "pcg", tol = 1e-8),
  score = list(traceProbes = 16L),
  convergence = list(accuracy = "adaptive")
)
print(ctrl)

sommer documentation built on Oct. 24, 2026, 1:08 a.m.