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(..., henderson=TRUE), the default). henderson=FALSE 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. REML=TRUE is the default. REML=FALSE 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-style information tends to give more conservative, stable movement for simple variance parameters. The implementation combines them:
$$ \mathcal I_{used} = (1-w_{EM})\mathcal I_{AI}+w_{EM}\mathcal I_{EM}. $$
emWeight controls \(w_{EM}\) by iteration. By default it starts relatively
conservatively and tapers toward a small value, so early updates get more
stabilization and later updates rely more heavily on AI. stepWeight
separately scales the score step. They are distinct controls: one blends
information matrices; the other damps the proposed correction.
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.
solver="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.
computeCi=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.
solver="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 computeCi=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.
solver="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 pcgTol=1e-8, pcgTraceProbes=8, and
pcgLanczosSteps=20; pcgMaxIters=0 chooses an automatic iteration limit.
The defaults are practical settings, not a universal error bound.
solver="pcg" currently requires REML=TRUE. computeCi=1 is not
available because it requests the LDLT/Takahashi selected inverse. Use
computeCi=0 for the factorization-free fit, or request
computeCi=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 (REML=FALSE) |
Yes | Yes | No |
Selected inverse (computeCi=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 solver="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 computeCi=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 computeCi=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 computeCi>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 solver="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\). henderson=TRUE 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 solver="cholmod" and
henderson=TRUE (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.
solver="auto" or solver="ldlt" for sparse pedigree-like
relationship precision. Try solver="cholmod" when supplied relationship
precision is dense or the LDLT factor has substantial fill.solver="pcg" when factor memory or direct factorization is the
bottleneck and approximate log determinants/traces are acceptable. Check
convergence and sensitivity to pcgTol, pcgTraceProbes, and
pcgLanczosSteps; do not assume solver agreement from BLUP convergence
alone.computeCi=0 when inverse-based prediction-error output is not
required. computeCi=1 requests a sparse selected inverse and is a natural
choice with LDLT. computeCi=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.