This vignette is a technical tour of ai_mme_sp2(), the C++ engine used by
mmes() for Henderson-based mixed-model fitting. It explains the equations
the engine solves, how covariance parameters are represented and estimated,
why there are three sparse-system solvers, and where the major computational
savings come from. Small examples connect the equations to actual R output.
The discussion is specific to the Henderson engine (mmes(..., formulation="henderson"), the default). formulation="direct" selects a different,
direct-inversion engine and is outside the solver comparisons below. Some
implementation details evolve; the source comments in src/MNR.cpp are the
authoritative record of the current implementation.
Reading map
ai_mme_sp2() solvesFor one response, write the model as
$$ y = X\beta + \sum_{i=1}^m Z_i u_i + e, \qquad u_i\sim N(0,G_i),\qquad e\sim N(0,R). $$
Here \(X\) and \(Z_i\) describe how observations connect to fixed and random
effects. The unknown covariance matrices are not arbitrary dense matrices.
They are generated from a small parameter vector: for example, a random
intercept has one variance, while an unstructured \(q\)-trait covariance has
\(q(q+1)/2\) covariance coordinates. vsm() and its covariance constructors
describe those structures to the engine.
The marginal covariance of the observations is
$$ V = \sum_i Z_iG_iZ_i’ + R. $$
One possible fitting strategy forms and factors \(V\), whose dimension is the
number of observations. ai_mme_sp2() instead uses Henderson’s equations,
whose dimension is the total number of fixed- and random-effect
coefficients. In grouped-data problems that coefficient system can be
substantially smaller and sparser than \(V\).
The regression coefficients and random-effect predictions are solved as linear-system quantities at each covariance-parameter iterate. The iterative optimization estimates the covariance parameters: the scale parameters, variance ratios, correlations, loadings, and any other free coordinates in the supplied covariance descriptors. Fixed covariance parameters are held at their specified values and omitted from the free update.
For REML, the target is the restricted likelihood, which accounts for the
fixed-effect design while estimating covariance. For ML, it is the ordinary
likelihood. task$objective="reml" is the default. task$objective="ml" is supported by the
ldlt and cholmod solvers; the current pcg path is REML-only.
Stack all design matrices as \(W=[X\ Z]\), and stack all coefficients as \(b=(\beta',u')'\). Let \(G\) denote the block-diagonal covariance of the random effects. Henderson’s mixed-model equations are
$$ C b = r, \qquad C = W’R^{-1}W + \begin{pmatrix}0&0\0&G^{-1}\end{pmatrix}, \qquad r=W’R^{-1}y. $$
The solution gives the generalized least-squares estimate \(\hat\beta\) and
the best linear unbiased predictions (BLUPs) \(\hat u\). In code, the random
prior contribution is assembled as structured precision blocks; if a
relationship inverse \(A_i\) is supplied, the corresponding contribution has
the form \(K_i^{-1}\otimes A_i\), with the exact factor order determined by the
vsm() descriptor.
The fitted-data quadratic can be computed without constructing the dense projection matrix \(P\):
$$ y’Py = y’R^{-1}y - r’C^{-1}r. $$
This identity is central to the Henderson route. The engine solves for \(b\) and obtains the quadratic from a scalar correction, rather than forming the observation-sized \(P=V^{-1}-V^{-1}X(X'V^{-1}X)^{-1}X'V^{-1}\).
Ignoring constants that do not depend on covariance parameters, the restricted log likelihood has the schematic form
$$ \ell_R = -\frac12\left[ \log|C| + \log|R| + y’Py
The multipliers \(a_i,b_i\) depend on the dimensions of the structured random term and its relationship precision. They account for the covariance determinants represented by the descriptors. The expression emphasizes the three numerical ingredients that matter on each iteration:
The code evaluates \(y'Py\) using the Henderson identity above and computes log determinants from factorizations or structure-specific formulas. It does not generally build \(P\) or an inverse of the full observation covariance.
For ML, the fixed effects are not removed by a restricted-likelihood determinant. The objective replaces \(\log|C|\) by \(\log|D|\), where \(D\) is the random-effects-only block of \(C\):
$$ D=Z’R^{-1}Z+G^{-1}. $$
The quadratic form is unchanged. This is why the ML and REML paths share most of the work but differ in the determinant term and corresponding score traces.
If \(V_j=\partial V/\partial\phi_j\) is the derivative of the marginal covariance with respect to a covariance parameter \(\phi_j\), the familiar REML score can be written
$$ U_j = \frac12\left[y’PV_jPy-\operatorname{tr}(PV_j)\right]. $$
The average-information (AI) matrix is commonly written
$$ \mathcal I_{ij}^{AI} =\frac12 y’PV_iPV_jPy. $$
These formulas are useful for understanding the statistics, but they are
not a recipe the C++ implementation follows literally. ai_mme_sp2() uses
equivalent Henderson-system identities, structured precision derivatives,
selected inverse entries, and linear solves. Avoiding explicit \(P\) and dense
\(V_j\) matrices is a major part of the speed advantage.
An optimizer that takes arbitrary real steps can propose a negative variance, a correlation outside its valid interval, or a non-positive definite covariance. The covariance descriptors solve this by optimizing working coordinates \(\eta\) and mapping them to legal natural parameters \(\theta=g(\eta)\).
Examples include:
| Quantity | Working coordinate | Natural value |
|---|---|---|
| Positive scale or variance ratio | \(\eta\in\mathbb R\) | \(\exp(\eta)>0\) |
| AR(1) correlation | \(\eta\in\mathbb R\) | \(\rho=\tanh(\eta)\) |
| Uniform correlation for \(q\) levels | \(\eta\in\mathbb R\) | logistic map into \((-1/(q-1),1)\) |
| Unstructured covariance | Cholesky-factor coordinates | \(K=LL'\) |
The product scale for a vsm() term is kept separate from covariance
shape. Schematically,
$$ G = \sigma^2(K_1\otimes\cdots\otimes K_s), $$
with one overall scale and dimensionless factor shapes. This avoids confounding multiple free scales inside one product. Internally, the overall scale is represented on a log scale and scaled relative to the response variance; results are transformed back to the data scale for reporting.
Each compiled covariance descriptor supplies the factor evaluator and derivative information needed by the optimizer. For built-in structures, these derivatives are computed natively. A user-defined structure can also provide a derivative callback; if it does not, the small factor matrix can be centrally finite-differenced.
The statistical optimizer therefore does not need a separate branch for every named covariance model. It requests a covariance, precision, log determinant, and derivatives from the descriptor interface. This is both a software design choice and a computational one: derivatives can be calculated at the compact factor level before any full coefficient-space work is attempted.
An explicitly fixed parameter is never included as an unknown in the information-system solve. Positive scales have lower bounds, and structured factors use their own valid working ranges. The engine also checks that proposed covariance matrices remain positive definite where required.
The coordinate transforms make many invalid values unreachable, but they do not guarantee that every numerical step is useful. A very large working step can produce an ill-conditioned covariance or reduce the likelihood. That is handled separately by trust caps, positive-definiteness repairs, and likelihood backtracking (Chapter 4).
At the current covariance parameter vector, the engine computes a score vector and an information matrix. The score measures the local direction in which the likelihood changes. The information matrix scales that direction according to the curvature of the likelihood surface. An information system is solved for a parameter correction; conceptually,
$$ \text{information}\times\text{correction}\approx\text{score}. $$
In ai_mme_sp2(), the information calculation uses the Henderson equations
and derivatives of the covariance precision. The AI terms require
sensitivity solves of the form
$$ C,\frac{\partial b}{\partial\phi_j} =-\frac{\partial C}{\partial\phi_j}b -\frac{\partial r}{\partial\phi_j}. $$
Those solutions feed the AI curvature without numerically differentiating the complete likelihood. Cross-derivatives between distinct covariance structures are zero where the structures are parameter-independent; within a structure, second derivative corrections are included as needed.
AI steps are usually much faster than expectation-maximization (EM) near the
optimum, but can be aggressive when a parameter is near a boundary or when
the local curvature is poorly behaved. EM updates provide a stable alternative.
update$blendWeight now combines covariance-update directions, not information matrices:
$$ \Delta_{blend} = (1-w)\Delta_{AI}+w\Delta_{EM}. $$
update$blendWeight is \(w\): 1 selects the EM direction, 0 selects AI, and intermediate
values interpolate proposals. The default update$blendWeight=NULL instead starts
with AI and computes no EM/PX moments while AI makes useful progress.
At exact explicit EM
endpoints, iterations skip AI construction except for final inference. Generic
covariance structures use increasing-Q generalized EM updates.
A single whole-number update$blendWeight greater than 1 sets the decay length.
For example, update$blendWeight=5 decays exponentially from 1 to 0.03 over five
iterations, then holds at 0.03. The decay length is capped at convergence$maxIters.
Fractional counts and vectors containing values greater than 1 are errors.
proposal$stepWeight damps the resulting proposal and defaults to 1 when blending is
active or automatic recovery is selected, preserving exact standalone EM
updates and allowing full AI proposals subject to trust caps and backtracking.
Explicit pure-AI fits retain their
previous damping schedule. Set update$blendWeight=0 to request the existing pure-AI path.
Positive blend weights support Gaussian Henderson or direct REML/ML
fitting; Henderson PCG, PQL, acceleration, factor-score augmentation,
and saddle-point fitting do not support them. These modes default to pure AI
when weights are omitted; explicitly positive weights are rejected.
Diagnostics, including numerical M-step convergence, appear in
engineDiagnostics$blend.
The same AI-first policy applies to task$objective="ml". At each accepted covariance,
the engine profiles the fixed effects by GLS. REML E-steps use the joint
coefficient inverse, including fixed-effect uncertainty. ML instead conditions
on that GLS estimate and uses random-effect posterior covariance
\(S=(Z'R^{-1}Z+G^{-1})^{-1}=D^{-1}\).
The random moment is \(mm'+S\); the residual moment is
\(rr'+ZSZ'\), where \(r=y-X\hat\beta-Zm\).
Weights apply the existing residual transformation to errors and designs.
Fixed effects stay deterministic during covariance Q and PX optimization,
including reduced-rank loading probes. Every likelihood trial then profiles
fixed effects again at its proposed covariance. Acceptance uses the observed
profile ML likelihood, not the auxiliary Q objective. The scalar residual
update divides by \(n\), not \(n-p\); its uncertainty trace includes only random
coefficients for ML. Final coefficient uncertainty still uses the joint
inverse. Diagnostics identify the likelihood mode as "ML" or "REML".
With verbose=TRUE, all mmes() optimizer paths share a compact table:
iter LogLik deltaLL update policy trust line notPD time(s)
update identifies AI, scaled regularized AI (AI-reg), EM, PXEM, or a
weighted mixture such as AI+EM:0.40.
policy shows auto:em or auto:pxem while AI-first monitoring is active.
Recovery rows show backtracking, low-gain, line-search, or info-solve;
the update column identifies the selected recovery method. scheduled means
an explicit positive blend weight; off means automatic recovery is disabled.
Trigger labels denote attempts, not successful recovery or convergence.
Detailed acceptance and stop reasons remain in engineDiagnostics$recovery.
deltaLL compares successive accepted likelihoods. time(s) is cumulative
elapsed time, including rejected trials; redundant clock and CPU columns are
removed. Trust, line-search and failed-structure diagnostics remain visible.
- denotes unavailable/unattempted scales or no applicable local notPD repair.
Q gains and step norms remain available in stored diagnostics.
See ?mmes for the timing of safeguard scales on each optimizer path.
Automatic recovery uses update$recovery="em" by default; update$recovery="pxem" selects the
scalar expanded alternative. Recovery requires a non-small projected score
and one of these events:
convergence$likelihoodTol.For an unreliable information solve, ordinary AI first retries with a
scaled, ridge-regularized direction in the feasible parameter tangent.
Trust-cap scaling and an increasing ridge bound each working-coordinate step,
preventing a nearly flat coordinate from shrinking all other updates.
An otherwise usable ordinary AI direction also takes this retry when its
undamped trust multiplier would be below 0.1; user-requested step damping
does not trigger it.
The AI-reg / info-solve row identifies this retry. It must remain feasible
and improve observed likelihood by more than 1e-10 through backtracking.
Only if that direction or its line search fails are EM/PXEM moments computed.
Other recovery triggers still use EM/PXEM directly.
An EM/PXEM recovery step also must improve observed likelihood by more than
1e-10. The next iteration always retries ordinary AI. This sequence may
repeat up to convergence$maxIters; repeated recoveries have no separate iteration cap.
If EM/PXEM recovery also stalls, the fit returns the
last accepted parameters, warns, and reports convergence=FALSE. An approaching
variance floor alone does not trigger recovery. These are initial fixed
heuristics, not guarantees that recovery will improve speed on every model.
Variance floors apply to free product-level log variances and to log variance,
specific-variance and innovation-variance ratios, not to log loadings,
Cholesky factors or spatial range parameters. The floor is
max(1e-8, execution$inverseTol) in each parameter’s natural units (absolute
variance for product scales, relative variance for ratios).
Adaptive AI uses a boundary-constrained solve: a joint step that crosses a
floor, supported by an outward marginal score/curvature step, fixes that
coordinate at the floor and re-solves the remaining directions. Correlated
overshoots without this marginal support retain the trust safeguard instead.
Boundary moves do not shrink unrelated coordinates through the global trust
multiplier. Fixed parameters and equality/scaling constraints remain exact;
a floor-bound parameter can re-enter the solve when its score points inward.
All such proposals still require observed-likelihood acceptance.
Experimental common-scale profiling rejects proposals at any variance floor;
refit with common-scale profiling off when a boundary is needed.
Convergence requires both a small projected score and either a small likelihood
change or relative parameter step. The score is projected through the free/tied
parameter tangent and removes outward directions at variance floors. Each
component is scaled by the square root of its projected AI diagonal, floored
at execution$inverseTol. Inward-pointing directions at a variance floor use
the unfloored curvature (bounded away from numerical zero), so a vanishing
log-variance derivative cannot falsely certify convergence. The maximum must not exceed
sqrt(max(convergence$likelihoodTol, convergence$parameterTol)).
engineDiagnostics$recovery records the method, attempted and accepted step
counts, source iterations,
trigger reasons, stop reason, and final score certificate. Those step counts
refer to EM/PXEM. regularizedAttempts counts stabilized-direction solves,
including failures and active-set retries; regularizedAcceptedSteps counts
accepted regularized AI updates. Paired regularizedIterations and
regularizedRidges record source iterations and ridges for usable stabilized
directions, whether or not the subsequent proposal is accepted. The existing blend
counters report actual AI, moment, and M-step work. Explicit numeric
update$blendWeight schedules disable automatic recovery and retain their semantics.
AI, EM, and PXEM remain local methods. A feasible covariance and a small score do not establish a global maximum. Compare converged fits from alternative starting values when competing optima are a concern.
Set update$recovery="pxem" to replace the EM/GEM direction with a scalar
parameter-expanded generalized EM (PX-GEM) direction:
$$ \Delta_{blend}=(1-w)\Delta_{AI}+w\Delta_{PX}. $$
Each random term uses \(u_j=a_jv_j\), \(a_j>0\), and \(G_{*j}=G_j(\phi)/a_j^2\). The current E-step uses \(a_j=1\). An ordinary EM/GEM update seeds joint numerical optimization over admissible original covariance parameters and log expansion factors. Reduction recovers \(G_j=a_j^2G_{*j}\) before blending; expansion resets to identity each iteration. Keeping constraints on the original parameters preserves fixed values and equality/ratio ties without clipping a reduced proposal.
# Standalone scalar PX-GEM
fit <- mmes(y ~ 1,
random = ~id,
data = dat,
control = mmes_control(update = list(recovery = "pxem", blendWeight = 1)))
# PX-GEM/AI with a five-iteration exponential taper
fit <- mmes(y ~ 1,
random = ~id,
data = dat,
control = mmes_control(update = list(recovery = "pxem", blendWeight = 5)))
This is increasing-Q numerical PX-GEM, not full matrix expansion. All existing covariance families supported by EM/GEM can use it, including latent reduced-rank terms. Residual-only models have no expansion and retain their ordinary update. Solver and fitting-mode restrictions remain unchanged. Full weight skips iterative AI; zero weight skips all EM and PX work. Non-reduced-rank AI directions use native exact AI kernels even within a blended schedule. Reduced-rank AI uses the same latent kernels as standalone RR AI. Accepted trial fits can supply the next zero-weight iteration and final exact inference without repeated coefficient factorization. Existing blend constraints and likelihood guards remain active.
The expanded residual objective includes cross-term posterior uncertainty.
Non-reduced-rank models cache quadratic expansion moments, with storage
quadratic in the number of random terms. Reduced-rank models cache exact
quadratic moments in a normalized loading/design basis once per E-step;
objective evaluations require no additional coefficient solves.
Cache storage is quadratic in the basis size. Extra workspace obeys
denseMemoryMB. PX proposals always pass observed-likelihood backtracking;
rejected proposals can fall back to the ordinary EM/GEM seed.
Both EM and PX reuse exact Schur inverse blocks or selected-inverse contractions for random moments. Diagonal covariance updates omit unused off-diagonal moments. Scalar residual updates use a trace identity, including fixed-effect uncertainty. Scalar non-reduced-rank PX expansion moments also reuse exact coefficient contractions when their sparse products fit the memory budget. Generic bounded panel solves remain the exact fallback.
Diagnostics report pxQGains, pxLogExpansion, numerical PX convergence and
evaluation counts, and pxEMFallbacks. These describe proposals before
backtracking. Fewer iterations need not mean less wall time: compare PX-AI
against EM-AI at the same final accuracy before choosing a method.
The default remains update$recovery="em".
The direct engine (formulation="direct") accepts the same blend and update$blendWeight
controls. Its default is AI-first recovery, explicit zero selects AI only, and
positive scalar/vector schedules blend genuine posterior-moment EM/GEM or
scalar PX-GEM directions with AI. These updates do not call Henderson.
For a random term with covariance \(\Sigma\otimes K\) and \(m\) relationship levels, cache \(B_{ab}=Z_aKZ_b'\). At an accepted fit, let \(r=Py\) and \(A=P\) for REML or \(A=V^{-1}\) for ML. The sufficient statistic is
$$ C_{ab}=r’B_{ab}r-\operatorname{tr}(AB_{ab}),\qquad S=m\Sigma+\Sigma C\Sigma. $$
Eligible scalar, diagonal, and unstructured random M-steps use these moments directly. Other families use numerical increasing-Q GEM, including custom covariances, constrained parameters, and structured or incomplete residual blocks. Zero-specific RR moves loadings into a latent design with nonsingular prior covariance. Candidate loading/PX changes retain exact cross-term posterior uncertainty in residual Q without a coefficient-sized posterior matrix. ML profiles beta inside the residual Q objective; REML integrates fixed-effect uncertainty.
Full EM proposals bypass AI trust caps and still pass covariance and observed likelihood checks. The direct driver retains projected-score convergence certification and clean final AI standard errors. Healthy AI iterations prepare no EM moments; failed or persistently difficult AI steps can request one EM/PX recovery, then resume AI. Direct information reliability uses diagonal-normalized curvature, avoiding recovery solely because a variance approaches its numerical floor.
engineDiagnostics$direct reports AI update solves and accepted-state
likelihood work. engineDiagnostics$blend reports moment preparations, Q gains,
numerical M-step status, actual weights, expansions, and fallbacks.
engineDiagnostics$recovery records automatic triggers and accepted recoveries.
Exact EM workspace obeys denseMemoryMB, separately from existing direct
caches. Structured/PX numerical updates may cost more than AI.
The candidate update is not accepted merely because the information system was solvable. The engine applies several layers of protection:
This distinction matters: a valid covariance is not necessarily a likelihood-improving update, and an improving update is not necessarily a well-conditioned one. The protections address both issues.
Convergence is assessed using likelihood and parameter-change criteria. A signed likelihood decrease is not interpreted as convergence; it triggers backtracking. Iteration counts describe accepted optimization iterations, not every internal line-search retry.
All three choices solve the same statistical model. The difference is how they handle the coefficient matrix \(C\) and the traces/log determinants needed for covariance estimation.
solve$method="ldlt": sparse direct factorizationThe default for the usual sparse MME is Eigen’s simplicial sparse LDLT factorization:
$$ C = LDL’, \qquad \log|C|=\sum_k\log(D_{kk}). $$
The factor supports accurate solves for the BLUP, fixed effects, and parameter sensitivities. If the sparsity pattern is unchanged between iterations, the symbolic analysis (ordering and fill structure) is reused; the numerical factorization is redone because covariance values change.
For score traces, the code can use a selected subset of \(C^{-1}\) instead of
forming all of \(C^{-1}\). In the LDLT case, a Takahashi recursion computes the
inverse entries associated with the factor’s nonzero pattern. If a needed
entry is outside the selected subset, the code falls back to batched solves.
output$inverse=1 uses this sparse selected-inverse route for requested
uncertainty calculations.
Strengths: deterministic direct solves and log determinants; useful selected-inverse structure; generally effective when \(C\) remains sparse.
Cost: a direct factorization can become expensive when a dense relationship precision causes extensive fill-in. Sparse input does not guarantee a sparse factor.
solve$method="cholmod": supernodal direct factorizationCHOLMOD is accessed through the Matrix package’s registered CHOLMOD
callables. Its supernodal sparse Cholesky factorization groups dense parts
of the factor into supernodes and uses matrix-matrix kernels. In broad terms,
this spends more arithmetic on dense blocks to use cache-friendly,
BLAS-3 operations effectively. That can be a better fit than simplicial
LDLT when the MME has substantial fill-in, particularly with dense
relationship matrices.
The fitted model and likelihood remain exact up to floating-point
round-off. For the score traces, CHOLMOD’s factor representation is not
converted to the LDLT/Takahashi selected-inverse representation in this
implementation. It uses cached inverse blocks when available and otherwise
batched factor solves. If output$inverse=1 is requested, a final LDLT
factorization is used to obtain the selected inverse subset.
There is also an exact block-Schur engine for eligible REML patterns. When random-effect blocks are uncoupled except through a relatively small border, each block is factored separately and the border is eliminated via a Schur complement. For disjoint blocks \(D_g\) and border Schur complement \(S\), the determinant identity is
$$ \log|C|=\sum_g\log|D_g|+\log|S|. $$
This changes the size of the expensive dense factorizations without
approximating the answer. Eligibility is structural; it is not a guarantee
that every model using cholmod activates this engine.
Strengths: often strong when dense fill makes a supernodal factorization more efficient; may activate the exact block-Schur path.
Cost: still a direct factorization in the general case. Dense factors need memory, and trace calculations may require solves when the inverse cache does not cover a requested entry.
solve$method="pcg": iterative solves and stochastic tracesPreconditioned conjugate gradients (PCG) solves \(Cx=b\) using products with \(C\) rather than a sparse direct factorization. The current preconditioner is diagonal. This avoids storing a large fill-heavy factor and is attractive when a matrix-vector product is cheaper than factorization.
REML also needs \(\log|C|\) and inverse traces, not just a BLUP solve. PCG approximates them as follows:
The probes are deterministic Rademacher vectors in this implementation. Reusing them at successive parameter values gives common random numbers: likelihood changes are less confounded by fresh Monte Carlo variation at every iteration. The estimates are still approximations, so probe count, Lanczos steps, and solve tolerance influence the speed-accuracy tradeoff.
Current defaults are solve$tol=1e-8, score$traceProbes=8, and
logdet$lanczosSteps=20; solve$maxIters=0 chooses an automatic iteration limit.
The defaults are practical settings, not a universal error bound.
solve$method="pcg" currently requires task$objective="reml". output$inverse=1 is not
available because it requests the LDLT/Takahashi selected inverse. Use
output$inverse=0 for the factorization-free fit, or request
output$inverse=2 only when an explicit full inverse is genuinely needed and
the system is small enough.
Strengths: can avoid direct factorization and its fill/memory cost; often benefits from parallel SLQ probes.
Cost: log determinants and traces are stochastic approximations, and the iteration may need many matrix-vector products or PCG solves. It is not automatically faster for small systems, poorly conditioned systems, or matrices for which the preconditioner is weak.
| Feature | ldlt |
cholmod |
pcg |
|---|---|---|---|
| Main system method | Sparse simplicial LDLT | Supernodal sparse Cholesky | Iterative PCG |
| `B0(\log | C | )B0` | Factor pivots |
| Trace strategy | Selected inverse, then solve fallback | Cached inverse blocks, then solves | Hutchinson probes and solves |
| Exact likelihood calculations | Yes, up to round-off | Yes, up to round-off | Approximate determinant/traces |
ML (task$objective="ml") |
Yes | Yes | No |
Selected inverse (output$inverse=1) |
Native Takahashi subset | Final LDLT subset | Not supported |
| Typical fit | Sparse pedigree-like precision | Dense/fill-heavy precision | Very large systems where solves dominate |
The default solve$method="auto" in mmes() selects cholmod if any supplied
random-effect relationship inverse has density greater than 0.2, and
otherwise selects ldlt. This is a heuristic for separating dense genomic
relationship matrices from sparse pedigree-style matrices. An explicit
solver choice always overrides it. pcg must be requested explicitly.
There is no single trick that makes every model fast. ai_mme_sp2() first
recognizes algebraic structure, then tries to perform the same statistical
calculation with smaller factors, fewer entries, or fewer repeated passes.
For a structured covariance
$$ G=s(K_1\otimes K_2\otimes\cdots\otimes K_d), $$
the inverse and determinant obey
$$ G^{-1}=s^{-1}(K_1^{-1}\otimes\cdots\otimes K_d^{-1}), $$
and
$$ \log|G|=q\log s+ \sum_{j=1}^d\frac{q}{q_j}\log|K_j|, \qquad q=\prod_jq_j. $$
Computing and retaining each small factor inverse can be much cheaper than assembling and factoring the full product. The same factor-wise strategy is used for first derivatives and log-determinant derivatives. Selected second derivative contractions are also evaluated factor by factor.
This is exact algebra, not an approximation. Generic factors that do not have a specialized formula use a small dense factorization. Some structures have additional native formulas; for example, AR(1), compound symmetry, and antedependence can avoid a generic factorization for their precision or determinant.
If residual covariance is diagonal, applying \(R^{-1}\) is elementwise and its log determinant is the sum of log diagonal entries. No sparse residual factorization is needed.
If observations form complete repeated blocks with a common Kronecker residual covariance, the engine can factor the small block once and reuse its inverse and log determinant across blocks. Trace contributions are accumulated block by block. Where the derivative and inverse blocks are symmetric, a trace can be computed as an elementwise inner product, \(\operatorname{tr}(AB)=\sum_{ij}A_{ij}B_{ij}\), rather than forming \(AB\).
Missing cells or irregular observation blocks can break the repeated Kronecker pattern. The code then uses sparse residual factorization or a more general path; it does not pretend that an incomplete block has the complete-block inverse.
Most score traces need only entries of \(C^{-1}\) corresponding to nonzeros of a sparse derivative matrix. For LDLT, the Takahashi recursion computes the inverse entries associated with the factor sparsity. The trace is then an accumulation over derivative nonzeros. This avoids the \(O(n^2)\) storage of a dense inverse when only a sparse subset is needed.
If a requested inverse entry is unavailable from the selected subset, the
implementation solves for batches of right-hand sides rather than
immediately materializing all of \(C^{-1}\). A full inverse is reserved for
the explicit output$inverse=2 request.
Covariance values change at every accepted parameter update, but many structural objects do not. When the sparse pattern is unchanged, LDLT and CHOLMOD reuse symbolic analysis and repeat only the numeric factorization. Other examples include cached residual block mappings, reusable derivative topology, and cached inverse blocks for repeated trace access. This is a common sparse-optimization pattern: distinguish the graph/pattern from the numeric values stored on that graph.
For a subset of REML models with compatible factor-wise random covariance
precision and output$inverse=0, the PCG iteration can apply the random prior
through Kronecker/tensor operations. It does not need to assemble the full
random precision contribution into \(C\) on each iteration. The full \(C\) is
materialized once at the end when needed for the fitted-object interface.
This path is intentionally narrower than PCG itself. It does not currently
cover ML or output$inverse>0; direct LDLT and CHOLMOD use an explicit coefficient
matrix. Compatibility depends on the model’s covariance structure, not only
on its size.
A dense relationship precision can make every random-effect block dense.
That raises both assembly work and factor fill. cholmod often helps because
its supernodal factorization can exploit dense blocks efficiently; eligible
models may also use the exact block-Schur engine. PCG can avoid factor
storage, but then each iteration’s repeated matrix-vector products and
stochastic trace calculations matter.
The best solver depends on the whole system: number of effects, sparsity, fill pattern, conditioning, covariance structure, requested uncertainty outputs, and available numerical-library threading. Solver labels alone do not predict elapsed time.
The speed comes from both the algebra above and the numerical kernels that execute it. The main components are:
| Component | Role in the engine | Performance relevance |
|---|---|---|
Eigen (RcppEigen) |
Sparse LDLT, PCG, sparse matrix operations, dense factor kernels | The main sparse-system machinery for ldlt and iterative solves; supports reusable symbolic analysis and ordering choices. |
| CHOLMOD, through Matrix | Supernodal sparse Cholesky and solves for solve$method="cholmod" |
Dense supernodes use efficient matrix-matrix kernels; runtime access is through Matrix’s registered C callables rather than a separate package link. |
Armadillo (RcppArmadillo) |
Dense covariance-factor arithmetic and much of the descriptor/derivative work | Compact factor matrices are small enough for dense linear algebra; interfaces with R data and Eigen-backed sparse work. |
| BLAS/LAPACK | Dense matrix and factorization kernels selected by the R build | Can have a major impact on supernodal/dense blocks and dense covariance factors. Actual implementation and threading depend on the R installation. |
| METIS (optional) | Nested-dissection ordering for some sparse factorizations | Can reduce fill for suitable large sparse graphs. If not detected at build time, Eigen AMD ordering is used instead. |
| OpenMP (when enabled in R) | Parallel independent probes and selected block/trace loops | Helps tasks with independent blocks or probes; availability depends on how R/package was built. |
DESCRIPTION declares Rcpp, RcppArmadillo, RcppEigen, and Matrix;
Matrix is also a package dependency. CHOLMOD is provided by the Matrix
installation, so a separate direct SuiteSparse link is not required by this
package. METIS is optional and detected at configure time. OpenMP and BLAS
threading are build/runtime properties; having a parallel-capable library
does not mean every loop is parallel or that increasing thread counts always
helps.
overhead can dominate and solver timing is not informative.
This section follows one covariance-parameter iterate of the model
$$ {yield ~ env},\qquad {random = vsm(usm(env), ism(id), Gu = Ai)}, $$
with Henderson’s equations and the CHOLMOD direct solver. The purpose is to
show the numerical mechanics of one likelihood evaluation and one MME solve,
not to claim that the parameter values chosen below are the converged REML
estimates. During fitting, ai_mme_sp2() repeats these calculations at
successive parameter values and uses the likelihood score and information
to update them.
The usual rcov = ~ units residual model is assumed. There is one record in
each environment-by-entry cell, so the residual covariance at the worked
iterate is \(R=I_6\). formulation="henderson" is the mmes() default; the solver choice
is explicitly cholmod.
Use the observation order
$$ (A,X),(A,Y),(A,Z),(B,X),(B,Y),(B,Z). $$
The response vector is
$$ y=\begin{pmatrix}10&11&12&20&19&21\end{pmatrix}’. $$
With treatment contrasts and A as the reference environment, the fixed effects are an intercept and the B-versus-A contrast:
$$ X=\begin{pmatrix} 1&0\1&0\1&0\1&1\1&1\1&1 \end{pmatrix}, \qquad \beta=\begin{pmatrix}\beta_0\\beta_B\end{pmatrix}. $$
Each environment-entry pair occurs once. In the stated random-effect ordering
$$ u=(u_{A,X},u_{A,Y},u_{A,Z},u_{B,X},u_{B,Y},u_{B,Z})’, $$
the random design is consequently \(Z=I_6\). This makes the example especially transparent: every observation has its own random coefficient, but those coefficients are correlated according to the environment and entry relationship structures.
AiLet the entry relationship covariance be
$$ A=\begin{pmatrix} 1&0.5&0.25\ 0.5&1&0.25\ 0.25&0.25&1 \end{pmatrix}, $$
so X and Y have relationship 0.5, while X and Z and Y and Z each have relationship 0.25. This matrix is positive definite. Its inverse is
$$ Ai=A^{-1}=\frac1{11} \begin{pmatrix} 15&-7&-2\ -7&15&-2\ -2&-2&12 \end{pmatrix} =\begin{pmatrix} 1.363636&-0.636364&-0.181818\ -0.636364&1.363636&-0.181818\ -0.181818&-0.181818&1.090909 \end{pmatrix}. $$
The distinction is important in this API: Gu in vsm(..., Gu=Ai) is the
relationship inverse/precision for the final main-effect levels, not
the covariance \(A\) itself. The model description may naturally start from
\(A\), but Henderson’s equations add its inverse to the random-effect
precision block.
To make one iteration numerically concrete, set the current product scale to \(\sigma_g^2=1\) and the current environment covariance shape to
$$ K_{env}=\begin{pmatrix}1&0.5\0.5&1\end{pmatrix}. $$
This is a valid two-level usm(env) shape: its first diagonal is fixed to
one, and its off-diagonal is 0.5. With ism(id) and the supplied entry
relationship, the random-effect covariance is
$$ G=\sigma_g^2(K_{env}\otimes A)=K_{env}\otimes A. $$
The first Kronecker factor is the slow index and the second is the fast index, so the coefficient order is environment first, then entry. Expanded,
$$ G=\begin{pmatrix} 1&0.5&0.25&0.5&0.25&0.125\ 0.5&1&0.25&0.25&0.5&0.125\ 0.25&0.25&1&0.125&0.125&0.5\ 0.5&0.25&0.125&1&0.5&0.25\ 0.25&0.5&0.125&0.5&1&0.25\ 0.125&0.125&0.5&0.25&0.25&1 \end{pmatrix}. $$
The first \(3\times3\) block is the entry relationship within environment A; the second diagonal block is the same covariance within B. Cross-environment covariance blocks are \(0.5A\). For example, \(\operatorname{cov}(u_{A,X},u_{B,X})=0.5\), while \(\operatorname{cov}(u_{A,X},u_{B,Z})=0.125\).
At this iterate, the random precision is
$$ G^{-1}=K_{env}^{-1}\otimes A^{-1},\qquad K_{env}^{-1}=\begin{pmatrix}4/3&-⅔\-⅔&4/3\end{pmatrix}, $$
or numerically
$$ G^{-1}=\begin{pmatrix} 1.818182&-0.848485&-0.242424&-0.909091&0.424242&0.121212\ -0.848485&1.818182&-0.242424&0.424242&-0.909091&0.121212\ -0.242424&-0.242424&1.454545&0.121212&0.121212&-0.727273\ -0.909091&0.424242&0.121212&1.818182&-0.848485&-0.242424\ 0.424242&-0.909091&0.121212&-0.848485&1.818182&-0.242424\ 0.121212&0.121212&-0.727273&-0.242424&-0.242424&1.454545 \end{pmatrix}. $$
Notice the negative off-diagonal entries in the precision. They do not mean negative relationships: the covariance is \(G\), whose corresponding cross-environment and cross-entry covariances above are positive. Precision encodes conditional coupling, and its signs need not match covariance signs.
The usm working coordinates at this exact shape are also easy to see.
Its lower Cholesky factor is
$$ L=\begin{pmatrix}1&0\0.5&\sqrt{0.75}\end{pmatrix}, \qquad K_{env}=LL’. $$
Thus the two free shape coordinates are the unconstrained lower entry \(L_{21}=0.5\) and the log diagonal coordinate \(\log(L_{22})=\log(\sqrt{0.75})\approx-0.143841\). Together with \(\log(\sigma_g^2)\) and the residual log variance, this example has four free covariance working coordinates at an ordinary unconstrained iterate. The optimizer updates these coordinates, not the entries of \(G\) directly.
At the selected iterate \(R=I_6\), so \(R^{-1}=I_6\). Since \(Z=I_6\), the Henderson matrix and right-hand side simplify to
$$ C=\begin{pmatrix} X’X&X’\ X&I_6+G^{-1} \end{pmatrix}, \qquad r=\begin{pmatrix}X’y\y\end{pmatrix}. $$
The fixed-effect crossproducts are
$$ X’X=\begin{pmatrix}6&3\3&3\end{pmatrix},\qquad X’y=\begin{pmatrix}93\60\end{pmatrix}. $$
Consequently,
$$ r=\begin{pmatrix}93&60&10&11&12&20&19&21\end{pmatrix}’. $$
Substituting the precision matrix from the preceding section gives this 8-by-8 symmetric positive-definite coefficient matrix (rounded to six decimal places):
$$ C=\begin{pmatrix} 6&3&1&1&1&1&1&1\ 3&3&0&0&0&1&1&1\ 1&0&2.818182&-0.848485&-0.242424&-0.909091&0.424242&0.121212\ 1&0&-0.848485&2.818182&-0.242424&0.424242&-0.909091&0.121212\ 1&0&-0.242424&-0.242424&2.454545&0.121212&0.121212&-0.727273\ 1&1&-0.909091&0.424242&0.121212&2.818182&-0.848485&-0.242424\ 1&1&0.424242&-0.909091&0.121212&-0.848485&2.818182&-0.242424\ 1&1&0.121212&0.121212&-0.727273&-0.242424&-0.242424&2.454545 \end{pmatrix}. $$
The upper-left \(2\times2\) block comes entirely from the fixed-effect design. The upper-right block records which observations belong to each random coefficient. The lower-right block combines the observation information \(Z'R^{-1}Z=I_6\) with the prior precision \(G^{-1}\). This is the key regularization: without \(G^{-1}\), the six cell effects would simply reproduce the six observations and would not be shrunk toward their covariance model.
CHOLMOD first applies a fill-reducing permutation \(P\) and computes a sparse Cholesky factorization of the permuted system,
$$ PCP’=LL’. $$
It then solves, in order, the triangular systems
$$ Lq=Pr,\qquad L’t=q, $$
and unpermutes \(t\) to obtain \(\hat b=(\hat\beta',\hat u')'\). On this small example the matrix is dense enough that sparse storage has little advantage; the same procedure is shown because it is the production solver path requested here. CHOLMOD obtains the log determinant from its factor, without constructing \(C^{-1}\):
$$ \log|C|=2\sum_j\log(L_{jj}) $$
for the permuted Cholesky factor (permutation does not change the determinant). Numerically, for this iterate, \(\log|C|\approx5.751654\).
Solving the displayed system gives
$$ \hat\beta_0=11.055556,\qquad \hat\beta_B=9.000000, $$
and
$$ \hat u=\begin{pmatrix} -0.433333&-0.233333&0.500000&-0.233333&-0.433333&0.500000 \end{pmatrix}’. $$
Thus the fitted fixed means are \(11.055556\) for A and \(20.055556\) for B. The six fitted values after adding random effects are approximately
$$ \hat y=\begin{pmatrix} 10.622222&10.822222&11.555556&19.822222&19.622222&20.555556 \end{pmatrix}’, $$
with residuals
$$ y-\hat y=\begin{pmatrix} -0.622222&0.177778&0.444444&0.177778&-0.622222&0.444444 \end{pmatrix}’. $$
The fixed effects absorb the broad environment difference. The random effects then account for entry-specific deviations, while their covariance structure borrows information both among entries and across environments. For instance, the estimated effect for entry X in B is informed not only by the B-X observation but also by the A-X effect through the 0.5 environment covariance and by Y/Z effects through \(A\).
The engine does not need to construct the full projection matrix \(P\). It uses
$$ y’Py=y’R^{-1}y-r’C^{-1}r. $$
Here \(y'R^{-1}y=y'y=1436\), and the solved correction is \(r'C^{-1}r\approx1433.866667\). Therefore
$$ y’Py\approx2.133333. $$
The random covariance determinant is computed factor-wise:
$$ \log|G| =3\log|K_{env}|+2\log|A| =3\log(0.75)+2\log(0.6875) \approx-1.612433. $$
The multipliers are the opposite factor dimensions: each environment covariance determinant is repeated for the three entries, and each entry relationship determinant is repeated for the two environments. Since \(R=I_6\), \(\log|R|=0\). Omitting the Gaussian constant that is independent of covariance parameters, the restricted log likelihood contribution at this iterate is
$$ \ell_R =-\frac12\left(\log|C|+\log|G|+\log|R|+y’Py\right) \approx-3.136277. $$
This is one evaluated point on the likelihood surface, not the answer to
the optimization. During fitting, the score asks how this objective changes
as the overall random scale, the two usm coordinates, and the residual
scale move.
For the chosen usm shape, write
$$ L=\begin{pmatrix}1&0\a&e^d\end{pmatrix},\qquad K=LL’=\begin{pmatrix}1&a\a&a^2+e^{2d}\end{pmatrix}. $$
At the current point, \(a=0.5\) and \(d=\log(\sqrt{0.75})\). The shape derivative matrices are
$$ \frac{\partial K}{\partial a} =\begin{pmatrix}0&1\1&2a\end{pmatrix} =\begin{pmatrix}0&1\1&1\end{pmatrix}, \qquad \frac{\partial K}{\partial d} =\begin{pmatrix}0&0\0&2e^{2d}\end{pmatrix} =\begin{pmatrix}0&0\0&1.5\end{pmatrix}. $$
Consequently, the covariance derivatives for the full random effect are
$$ \frac{\partial G}{\partial a} =\sigma_g^2\left(\frac{\partial K}{\partial a}\otimes A\right), \qquad \frac{\partial G}{\partial d} =\sigma_g^2\left(\frac{\partial K}{\partial d}\otimes A\right), \qquad \frac{\partial G}{\partial\log\sigma_g^2}=G. $$
The precision derivatives follow from differentiating an inverse:
$$ \frac{\partial G^{-1}}{\partial\phi} =-G^{-1}\frac{\partial G}{\partial\phi}G^{-1}. $$
The determinant derivatives can also be computed without finite-differencing the full six-dimensional covariance:
$$ \frac{\partial\log|G|}{\partial\phi} =\operatorname{tr}\left(G^{-1}\frac{\partial G}{\partial\phi}\right). $$
The descriptor works on the 2-by-2 environment factor and then lifts the result through the Kronecker product. The residual log-scale derivative is similarly simple because \(R=\sigma_e^2I_6\). These derivatives feed the likelihood score and the AI/EM information calculation described earlier.
The optimizer then solves its small covariance-parameter information system, applies working-coordinate step caps and covariance-validity checks, and evaluates the trial likelihood. A likelihood-decreasing proposal is shortened. The next accepted parameter point rebuilds or rescales the necessary numeric matrix values and repeats the solve. The fixed design, incidence pattern, and symbolic factorization structure can be reused when their sparsity pattern remains unchanged.
For comparison with the deliberately chosen iterate above, fitting the
specified six observations with solve$method="cholmod" and
formulation="henderson" (using the usual independent units residual) gives the
following native-scale covariance estimates on this run:
$$ \widehat{\operatorname{var}}(u_{A,\cdot})=0.664395,\qquad \widehat{\operatorname{cov}}(u_{A,\cdot},u_{B,\cdot})=0.683457, \qquad \widehat{\operatorname{var}}(u_{B,\cdot})=0.706016, \qquad \widehat{\sigma}_e^2=0.482605. $$
Here the three random-effect quantities describe the environment covariance factor; the entry relationship \(A\) still multiplies it by the Kronecker construction. The fitted environment covariance is approximately
$$ \widehat{\Sigma}_{env}= \begin{pmatrix} 0.664395&0.683457\ 0.683457&0.706016 \end{pmatrix}. $$
Its correlation is about \(0.998\) and its determinant is only about \(0.00196\). This is a near-boundary covariance estimate from a deliberately tiny dataset, not a recommendation to infer a stable two-environment covariance from six records. It also illustrates why positive-definiteness checks, transformed coordinates, AI/EM blending, and likelihood step acceptance matter: covariance estimates can be weakly identified even when the linear MME solve itself is straightforward.
These converged values are not the values used to assemble the 8-by-8 matrix earlier. That matrix used the simple illustrative point \(\sigma_g^2=1\), \(K_{env,AB}=0.5\), and \(\sigma_e^2=1\) so each intermediate calculation could be inspected by hand. A real fit evaluates many such systems while moving from its initialized covariance parameters toward the accepted likelihood optimum.
Even at this tiny size, the fitting work has distinct layers:
The example is deliberately small enough that every matrix can be written down. In a large model, the same statistical steps remain, but sparse ordering, fill, selected/cached inverse information, repeated-block structure, Kronecker identities, and numerical-library kernels determine whether those steps are affordable.
solve$method="auto" or solve$method="ldlt" for sparse pedigree-like
relationship precision. Try solve$method="cholmod" when supplied relationship
precision is dense or the LDLT factor has substantial fill.solve$method="pcg" when factor memory or direct factorization is the
bottleneck and approximate log determinants/traces are acceptable. Check
convergence and sensitivity to solve$tol, score$traceProbes, and
logdet$lanczosSteps; do not assume solver agreement from BLUP convergence
alone.output$inverse=0 when inverse-based prediction-error output is not
required. output$inverse=1 requests a sparse selected inverse and is a natural
choice with LDLT. output$inverse=2 requests the full inverse and can be very
expensive in both time and memory.The central design principle is simple: ai_mme_sp2() tries to preserve the
statistical calculation while changing the linear algebra representation.
It uses Henderson’s equations to avoid observation-space dense matrices,
descriptor-level formulas to exploit covariance structure, and a solver
chosen for the sparsity and scale of the resulting coefficient system.