The FEM Matérn prior on the spline weights contributes two terms to
the joint log-density (see vignette("fem_joint_density")):
a quadratic form, random_fem_ssq_*, and a log-determinant,
random_fem_det_*. This vignette is about the first,
\[ f(\boldsymbol\gamma,\boldsymbol\theta) = -\tfrac12\,\boldsymbol\gamma^\top \mathbf{Q}(\boldsymbol\theta)\,\boldsymbol\gamma , \]
where \(\boldsymbol\gamma\in\mathbb{R}^n\) holds the spline weights (the random effects, one per active basis function) and \(\boldsymbol\theta=(\mathrm{range}, \mathrm{sd})\), optionally on the log scale.
Splitting the Gaussian log-density this way matters because the two
halves have opposite cost profiles. The quadratic form is arithmetically
trivial but involves all \(n\) weights,
so its Hessian is an \(n\times n\)
sparse block. The log-determinant needs a sparse factorization but
depends on \(\boldsymbol\theta\) alone,
so its Hessian is \(2\times2\).
adlaplaceFem therefore handles them by different
mechanisms: the determinant is a CppAD atomic operation
(src/fem_logdet_atomic.hpp), while the quadratic form
bypasses the tape entirely with closed-form derivatives
(src/fem_ssq_analytic.hpp and
src/fem_ssq_shard.hpp). The last section explains why the
split falls that way.
Two facts about this term do all the work.
1. \(\mathbf{Q}\) is a linear combination of constant matrices. The SPDE precision assembles as
\[ \mathbf{Q}(\boldsymbol\theta) = \sum_{j=1}^{m} c_j(\boldsymbol\theta)\,\mathbf{M}_j , \qquad \mathbf{M} = \begin{cases} (\mathbf{C},\mathbf{G},\mathbf{G}_2) & \alpha=2,\ m=3,\\ (\mathbf{C},\mathbf{G},\mathbf{G}_2,\mathbf{G}_3) & \alpha=3,\ m=4. \end{cases} \]
The Gram matrices are fixed once the basis is chosen; they carry no \(\boldsymbol\theta\) dependence at all. Every hyperparameter enters through \(m = 3\) or \(4\) scalars.
2. \(f\) is quadratic in \(\boldsymbol\gamma\). So the \(\boldsymbol\gamma\) Hessian block is constant in \(\boldsymbol\gamma\), and third derivatives in \(\boldsymbol\gamma\) alone vanish.
Together these mean every derivative the Laplace approximation needs is a handful of \(O(\mathrm{nnz}(\mathbf{Q}))\) passes over arrays that are already in memory.
With \(\nu = \alpha - 1\) in two dimensions, the public \((\mathrm{range}, \mathrm{sd})\) pair converts to the SPDE parameters as
\[ \kappa = \frac{\sqrt{8\nu}}{\mathrm{range}}, \qquad \tau = \begin{cases} \bigl(\kappa\,\sigma\sqrt{4\pi}\bigr)^{-1} & \alpha=2,\\[2pt] \bigl(\kappa^2\sigma\sqrt{8\pi}\bigr)^{-1} & \alpha=3, \end{cases} \qquad \sigma = \mathrm{sd}, \]
so that sd is the marginal standard deviation of the
field (vignette("fem_matern_methods") derives this). The
coefficients are then
| \(\alpha\) | \(\mathbf{Q}\) | \(\mathbf{c}(\boldsymbol\theta)\) |
|---|---|---|
| 2 | \(\tau^2(\kappa^4\mathbf{C} + 2\kappa^2\mathbf{G} + \mathbf{G}_2)\) | \(\tau^2(\kappa^4,\ 2\kappa^2,\ 1)\) |
| 3 | \(\tau^2(\kappa^6\mathbf{C} + 3\kappa^4\mathbf{G} + 3\kappa^2\mathbf{G}_2 + \mathbf{G}_3)\) | \(\tau^2(\kappa^6,\ 3\kappa^4,\ 3\kappa^2,\ 1)\) |
There is exactly one definition of this map in the package,
fem_Q_coefficients<Alpha>() in
src/ad_pack_random_fem.cpp, and three consumers go through
it: the AD quadratic form, the log-determinant atomic, and the analytic
shard’s coefficient tape. Keeping it single-sourced is deliberate — the
three would otherwise drift apart silently, and a mismatch between the
determinant’s \(\mathbf{Q}\) and the
quadratic form’s \(\mathbf{Q}\)
produces a plausible-looking but wrong likelihood.
Write \(q_j = \boldsymbol\gamma^\top\mathbf{M}_j\boldsymbol\gamma\), so that \(f = -\tfrac12\sum_j c_j q_j\). Everything follows by the product rule, since the \(\mathbf{M}_j\) are constant.
| Block | Closed form |
|---|---|
| \(\partial f/\partial\boldsymbol\gamma\) | \(-\mathbf{Q}\boldsymbol\gamma\) |
| \(\partial f/\partial\theta_k\) | \(-\tfrac12\sum_j \dfrac{\partial c_j}{\partial\theta_k}\,q_j\) |
| \(\partial^2 f/\partial\boldsymbol\gamma\,\partial\boldsymbol\gamma^\top\) | \(-\mathbf{Q}\) |
| \(\partial^2 f/\partial\boldsymbol\gamma\,\partial\theta_k\) | \(-\bigl(\partial\mathbf{Q}/\partial\theta_k\bigr)\boldsymbol\gamma\) |
| \(\partial^2 f/\partial\theta_k\,\partial\theta_l\) | \(-\tfrac12\sum_j \dfrac{\partial^2 c_j}{\partial\theta_k\partial\theta_l}\,q_j\) |
with \(\partial\mathbf{Q}/\partial\theta_k = \sum_j (\partial c_j/\partial\theta_k)\mathbf{M}_j\), another linear combination over the same sparsity pattern.
The third row is the one that pays for the whole design. The \(\boldsymbol\gamma\)-\(\boldsymbol\gamma\) Hessian block is the precision matrix itself, up to sign, and it does not depend on \(\boldsymbol\gamma\). Once \(\mathbf{Q}\) has been assembled for the current \(\boldsymbol\theta\) — which the value and gradient already require — the Hessian block needs no further arithmetic, just a copy with a sign flip. Note also that only \(q_j\) couples to \(\boldsymbol\gamma\) in the two \(\boldsymbol\theta\) rows, so the \(m\) scalars \(q_j\) are the only \(\boldsymbol\gamma\)-dependent quantities needed there.
The Laplace approximation’s derivative needs more than the Hessian.
Correcting \(\tfrac12\log\det\mathbf{H}\) for the
movement of the inner mode requires \(\operatorname{tr}(\mathbf{H}^{-1}\partial\mathbf{H}/\partial\theta_k)\),
which adlaplace obtains through the
trace_hinv_t shard method: for each column \(\mathbf{v}\) of a factor of \(\mathbf{H}^{-1}\), accumulate the third
derivative of the log-density contracted twice with \(\mathbf{v}\).
It is tempting to say the third derivative of a quadratic form is zero, but that is only true of the \(\boldsymbol\gamma\boldsymbol\gamma\boldsymbol\gamma\) block. The mixed block is not zero, and it is the entire answer:
\[ \frac{\partial^3 f}{\partial\theta_k\,\partial\boldsymbol\gamma\,\partial\boldsymbol\gamma^\top} = -\frac{\partial\mathbf{Q}}{\partial\theta_k}. \]
The directions adlaplace supplies are scattered into
\(\boldsymbol\gamma\) slots only, so
\(\mathbf{v}_{\boldsymbol\theta}=\mathbf{0}\).
The order-2 Taylor coefficient of \(f\)
along such a direction is \(-\tfrac12\mathbf{v}^\top\mathbf{Q}(\boldsymbol\theta)\mathbf{v}\),
whose gradient is identically zero on every \(\boldsymbol\gamma\) coordinate, leaving two
numbers per column:
\[ \text{trace contribution to }\theta_k = -\tfrac12\,\mathbf{v}^\top\frac{\partial\mathbf{Q}}{\partial\theta_k}\mathbf{v}. \]
Only first derivatives of the coefficient map are needed for this, and the result is exactly zero on all \(n\) gamma coordinates — a strong invariant, and one the tests assert directly.
This shard is assigned every column of \(\mathrm{half\_H\_inv}\) that meets the spatial field, so the sum over those columns is the Frobenius product of \(\partial\mathbf{Q}/\partial\theta_k\) with the joint inverse Hessian on the pattern of \(\mathbf{Q}\). With the same upper-triangle weights \(w_k\) (1 on the diagonal, 2 off it) the contribution to \(\theta_k\) is
\[ -\frac12\sum_k w_k\Bigl(\frac{\partial\mathbf{Q}}{\partial\theta_k}\Bigr)_k (H^{-1})_{s_k}. \]
The slots \(s_k\) are computed in R
at ad_pack time (hinv_trace_index, stored as
chol_inner_list$hinv_q_index) and the dot is \(O(\mathrm{nnz}(\mathbf{Q}))\). A pair of
\(\mathbf{Q}\) that the symbolic
inverse does not store is a structural zero and is skipped.
inner_opt already holds numeric \(H^{-1}\) when it asks for the trace, so the
production path does not walk columns.
The stamped column loop remains the fallback when that index was not
attached (ad_pack of a bare pointer) or only
half_H_inv was supplied. Because \(\mathbf{v}\) is sparse, the fallback stamps
\(\operatorname{supp}(\mathbf{v})\) and
visits only the pattern entries whose row and column are both stamped,
so the per-column cost is \(O(\sum_{c\,\in\,\mathrm{supp}}\mathrm{nnz}(\mathbf{Q}_{\cdot
c}))\) rather than \(O(n)\).
Every upper-triangle pair \((r,c)\)
with \(r \le c\) lives in CSC column
\(c\), so iterating stamped columns and
testing row stamps visits each pair exactly once.
Taping the quadratic form puts one multiply-add chain on the AD tape per nonzero of \(\mathbf{Q}\), and three consumers then sweep it.
| Quantity | AD tape path | Closed form |
|---|---|---|
| Value | Forward(0) |
one fused pass for all \(q_j\) |
| Gradient | Forward(0) + Reverse(1) |
one extra sparse matvec |
| Hessian | Forward(1)/Reverse(2) per
color; the color count is the distance-2 chromatic number of
\(\mathbf{Q}\)’s graph, in the hundreds
for this stencil |
copy \(\mathbf{Q}\) values, \(O(\mathrm{nnz})\) |
| Trace | Forward(0,1,2) + Reverse(3) per
assigned gamma column, and every gamma column is assigned to
this shard |
dot of \(\partial Q/\partial\theta\) with \(H^{-1}\) at the stored slots, \(O(\mathrm{nnz}(Q))\); column stamp loop if the index is absent |
| Build | record the tape, then four calibration sweeps including the symmetric Hessian coloring | patterns constructed directly, no tape |
The trace row is likely the largest single per-iteration cost in a
taped fit, ahead of even the Hessian coloring. The build row accounts
for most of the long pause with no output before optim()
starts on a large mesh.
Two points are worth separating, because they explain why this is a
shard rather than a CppAD atomic. First, the Hessian and trace losses
happen in CppAD’s drivers — sparse_hes recovers
Hessian values by probing with directional derivatives, and
EvalShard::trace_hinv_t clears and fills a dense direction
vector of length \(n\) once per column.
An atomic sits underneath both and cannot remove either, because
atomic_four exposes forward and reverse products but no
hook for handing over sparse Hessian values. Second, an atomic for this
term would have domain \(n+2\), so
every call marshals dense vectors of that length. Replacing the shard is
the level at which both costs disappear.
All \(m\) Gram matrices share one
sparsity pattern — the upper triangle of \(\mathbf{Q}\) in CSC form, from
fem_precision_payload(). femssq::Data holds
that pattern once (Q_p, Q_i) plus \(m\) aligned value vectors in
M, along with per-entry symmetry weights w (1
on the diagonal, 2 off it) so that a single pass over the stored upper
triangle yields a full symmetric quadratic form.
Three primitives in fem_ssq_analytic.hpp do all the
linear algebra:
| Primitive | Computes | Used for |
|---|---|---|
combine(d, coef, out) |
out[k] \(=\sum_j\)
coef[j] * M[j][k] |
assembling \(\mathbf{Q}\) and each \(\partial\mathbf{Q}/\partial\theta_k\) |
sym_matvec(d, A, x, y) |
\(\mathbf{y}=\mathbf{A}\mathbf{x}\) for symmetric \(\mathbf{A}\) in upper-triangle CSC | \(-\mathbf{Q}\boldsymbol\gamma\) and the cross Hessian blocks |
sym_quad_all(d, x, out) |
out[j] \(=\mathbf{x}^\top\mathbf{M}_j\mathbf{x}\)
for all \(j\) |
the \(m\) scalars \(q_j\), in one pass |
sym_quad_all is fused over \(j\) on purpose: the pattern traversal and
the \(w_k x_r x_c\) product are shared
across all \(m\) Grams, so computing
all \(m\) quadratic forms costs barely
more than computing one.
sym_matvec accumulates the mirrored contribution
separately from the direct one, so a diagonal entry is not counted
twice.
femssq::Workspace holds everything that changes with the
parameters: coefficient values and derivatives (c,
dc, d2c), the combined value vectors
(Q_x, dQ_x[k]), and the derived vectors
(gamma, Qg, dQg[k],
quad). update_theta() takes an
order argument so each entry point pays only for what it
needs — order 0 for the value, 1 to add \(\partial c/\partial\theta\) and \(\partial\mathbf{Q}/\partial\theta_k\), 2 to
add \(\partial^2
c/\partial\theta^2\).
The trace_hinv_t scratch deserves a note. Instead of
clearing an \(n\)-vector per column,
the workspace keeps v_dense alongside a stamp
array and an integer stamp_tag; a row counts only when its
stamp matches the current tag, so stale values are never read and
nothing needs zeroing.
One piece is still automatic differentiation, and correctly so. The
map \(\boldsymbol\theta\mapsto\mathbf{c}\) is
recorded as a CppAD::ADFun<double> with domain 2 and
range \(m\)
(make_fem_coef_fun<Alpha>()), and
update_theta calls Forward,
Jacobian, and \(m\) tiny
Hessian calls on it. This avoids hand-differentiating the
\((\mathrm{range},\mathrm{sd}) \to
(\kappa,\tau)\) composition together with the optional log
transform — four cases, all easy to get subtly wrong — and it costs
nothing next to the \(O(m\cdot\mathrm{nnz})\) combine passes.
This is the same principle that makes the log-determinant a good atomic: AD is cheap and reliable on a two-dimensional domain. It is expensive on an \(n\)-dimensional one. The design applies AD exactly where the domain is small.
Because CppAD sweeps are not thread safe, each shard clone takes its
own copy of this tape in Workspace::init, so threads never
share one.
ad_shard is a virtual interface, and
adlaplaceFem already owns the shard factory through
ADLAPLACE_DEFINE_BACKEND(). FemSsqShard
derives from ad_shard and overrides the evaluation methods.
trace_hinv_from_h is an extra virtual (ABI 5); the numeric
work stays in adlaplaceFem.
| Method | Behaviour |
|---|---|
f |
\(q_j\) pass, then \(-\tfrac12\sum_j c_j q_j\) |
f_grad |
adds one sym_matvec for \(-\mathbf{Q}\boldsymbol\gamma\) |
f_grad_hess |
adds sym_matvec per \(\theta_k\), then fills Hessian nonzeros
from slot maps |
trace_hinv_from_h |
dot \(w\cdot\partial Q/\partial\theta\cdot
H^{-1}\) at hinv_q_index; two scalars on the theta
slots |
trace_hinv_t |
stamped support loop when the index or \(H^{-1}\) is absent |
get_sizes, get_sparse_pattern |
delegate to the free helpers pack_sparsity_sizes /
get_pattern |
assign_memory |
nothing to do: no sweeps means no Taylor workspace |
clone |
copies the shared_ptr to the immutable data |
The bookkeeping that connects closed-form values to the pack’s sparsity patterns is the one genuinely fiddly part. At construction the shard classifies every nonzero of the four patterns it was handed:
GradSlots marks each gradient nonzero as a gamma row, a
theta index, or ignorable.HesSlots marks each Hessian nonzero as
kGammaGamma (carrying a position into Q_i),
kGammaTheta (a gamma row and a theta index),
kThetaTheta, or kIgnore.These are built by reading pack.pattern_grad,
pattern_grad_inner, pattern_hessian, and
pattern_hessian_inner directly, so the nonzero ordering is
correct by construction rather than by assuming how the analytic pattern
happened to be enumerated.
Two conventions must match EvalShard exactly, and they
differ from each other:
out_f and out_grad are
accumulated with +=, since several shards
contribute to the same slots.out_hes[map[D]] is assigned with
=. Structurally present entries whose analytic value is
zero must therefore be written as explicit zeros rather than skipped,
which is what the kIgnore branch does.Because the closed forms are worth trusting only if they are checked, the density keeps an AD reference path behind one config flag.
fem_analytic |
Behaviour |
|---|---|
TRUE (default) |
Closed-form FemSsqShard evaluation; tape never
recorded. |
FALSE |
Standard EvalShard over a CppAD tape that evaluates
\(\gamma^\top M_j\gamma\) with
quadform_atomic (one call per Gram). |
In the default mode get_tape_sizes() reports
domain, size_op, and size_var as
zero. That is a diagnostic only; nothing downstream consumes it, and the
zeros do not indicate a broken pack.
A small model, built both ways from the same FEM ingredients.
ssq_pack <- function(alpha = 2L, analytic = TRUE,
range = sqrt(8) / 1.5, sd = 0.8, nk = 4L) {
degree <- if (alpha == 2L) 2L else 3L
knots_list <- list(
x = seq(0, 1, length.out = nk),
y = seq(0, 1, length.out = nk)
)
sites <- do.call(expand.grid, knots_list)
fem <- fem_bspline(sites, knots_list, degree = degree)
prec <- fem_precision_payload(fem, alpha = alpha)
nr <- nrow(fem$C)
config <- list(
gamma = rep(0, nr),
theta = c(log(range), log(sd)),
transform_theta = TRUE,
fem_analytic = analytic
)
dd <- adlaplace::density_data(
gamma_map = Matrix::Diagonal(nr),
theta_map = list(c(1L, 2L), 2L),
ad_kind = "random",
density = paste0("random_fem_ssq_", alpha),
package = "adlaplaceFem",
precision = prec
)
list(
ptr = adlaplace::ad_pack_ptr(dd, config),
fem = fem,
n_gamma = nr,
theta = c(log(range), log(sd))
)
}
analytic <- ssq_pack(analytic = TRUE)
ad <- ssq_pack(analytic = FALSE)
c(n_gamma = analytic$n_gamma)## n_gamma
## 25The default mode records no tape, while the patterns match the ones the atomic AD tape discovers:
rbind(
analytic = unlist(adlaplace:::get_tape_sizes(analytic$ptr, 0L)),
ad = unlist(adlaplace:::get_tape_sizes(ad$ptr, 0L))
)## domain n_global size_op size_var nnz_grad nnz_grad_inner nnz_hes
## analytic 0 27 0 0 27 25 246
## ad 27 27 131 49 27 25 246
## nnz_hes_inner
## analytic 193
## ad 193domain and size_op collapse to zero on the
analytic pack, but nnz_grad and nnz_hes agree.
The AD tape domain is \(n+2\) with a
small operator count from the Gram atomics rather than an \(O(\mathrm{nnz}(Q))\) product loop.
Both modes agree on the value, gradient, and Hessian:
set.seed(3)
x <- c(rnorm(analytic$n_gamma), analytic$theta)
f_of <- function(p) adlaplace::joint_log_dens(p$ptr, x, negative = FALSE)
g_of <- function(p, inner) {
as.numeric(adlaplace::grad(p$ptr, x, inner = inner, negative = FALSE))
}
H_of <- function(p, inner) {
as.matrix(adlaplace::hessian(p$ptr, x, inner = inner, negative = FALSE))
}
data.frame(
quantity = c("log density", "gradient (outer)", "Hessian (outer)"),
analytic_vs_AD = c(
abs(f_of(analytic) - f_of(ad)),
max(abs(g_of(analytic, FALSE) - g_of(ad, FALSE))),
max(abs(H_of(analytic, FALSE) - H_of(ad, FALSE)))
)
)## quantity analytic_vs_AD
## 1 log density 0.000000e+00
## 2 gradient (outer) 7.105427e-15
## 3 Hessian (outer) 4.263256e-14The sharpest single check is the identity \(\partial^2 f/\partial\boldsymbol\gamma\partial\boldsymbol\gamma^\top = -\mathbf{Q}\). Assembling \(\mathbf{Q}\) independently in R from \(\kappa\) and \(\tau\) must reproduce the gamma block of the inner Hessian exactly, at any \(\boldsymbol\gamma\):
kappa <- 1.5
sd <- 0.8
tau <- 1 / (kappa * sd * sqrt(4 * pi))
nr <- analytic$n_gamma
H_inner <- H_of(analytic, TRUE)
Q <- as.matrix(fem_precision(
kappa = kappa, tau = tau,
analytic$fem$C, analytic$fem$G, analytic$fem$G2,
alpha = 2L
))
c(max_abs_diff = max(abs(H_inner[seq_len(nr), seq_len(nr)] + Q)))## max_abs_diff
## 8.881784e-16Finally the third-order trace. Both modes agree, and the analytic result is exactly zero on every gamma coordinate:
af_analytic <- adlaplace::ad_pack(analytic$ptr, num_threads = 1L)
af_ad <- adlaplace::ad_pack(ad$ptr, num_threads = 1L)
half <- af_analytic@chol_inner_list$half_H_inv
half <- methods::as(methods::as(half, "generalMatrix"), "CsparseMatrix")
set.seed(5)
half@x[] <- rnorm(length(half@x))
cols <- af_analytic@chol_inner_list$trace_columns
t_analytic <- adlaplace::trace_hinv_t(af_analytic, x, half, cols)
t_ad <- adlaplace::trace_hinv_t(af_ad, x, half, cols)
c(
rel_diff_vs_AD = max(abs(t_analytic - t_ad)) / max(abs(t_ad)),
max_abs_on_gamma = max(abs(t_analytic[seq_len(nr)])),
theta_1 = t_analytic[nr + 1L],
theta_2 = t_analytic[nr + 2L]
)## rel_diff_vs_AD max_abs_on_gamma theta_1 theta_2
## 1.669384e-15 0.000000e+00 -1.029183e+03 1.089617e+03The two theta entries are the whole contribution; the gamma entries are structurally zero rather than numerically small.
tests/testthat/test-fem-ssq-analytic.R runs these
comparisons across \(\alpha \in
\{2,3\}\) crossed with transform_theta, adds the
pattern-equality assertion between the hand-built and tape-discovered
patterns, and exercises a clone under OpenMP.
The two FEM terms illustrate the two available extension points, and the choice between them turns on the size of the domain.
random_fem_det_* |
random_fem_ssq_* (analytic) |
random_fem_ssq_* (AD) |
|
|---|---|---|---|
| Mechanism | CppAD::atomic_four |
custom ad_shard |
quadform_atomic on a tape |
| Domain of the hand-written code | \(m = 3\) or \(4\) coefficients | \(n+2\) parameters | fixed Grams; \(n\) on each atomic call |
| Derivative route | trace identity \(\mathrm{d}\log\det\mathbf{Q} = \operatorname{tr}(\mathbf{Q}^{-1}\mathrm{d}\mathbf{Q})\), evaluated with a sparse LDL and a Takahashi selected inverse | direct differentiation of a quadratic form | atomic forward/reverse of \(\gamma^\top M_j\gamma\), then chain through \(c_j(\theta)\) |
| Highest order implemented | reverse order 1 | closed form throughout | Taylor / reverse through order 2 |
| Sits on a tape | yes | no | yes |
The determinant is a natural atomic because the coefficients are the
only inputs: it is called with a 3- or 4-vector, CppAD’s sparsity and
coloring machinery sees a \(2\times2\)
Hessian, and the atomic can decline orders above 1 outright — which it
does. It gets away with that because
trace_columns_from_pattern assigns trace columns by
matching a shard’s grad_inner against the gamma slots, and
the determinant depends on \(\boldsymbol\theta\) only, so it is never
assigned any.
The quadratic form’s production path is the opposite case: it is the
term that gets assigned every gamma column, its Hessian block is an
\(n\times n\) sparse matrix, and the
expensive structure lives at the shard level. The AD reference path
still uses atomics, but only for the fixed Grams \(M_j\) (scalar \(y_j = \gamma^\top M_j\gamma\)), with \(\tau\)/range coefficients kept outside —
the same split random_mult uses for a fixed \(Q\).
One structural precondition makes the analytic shard possible:
random_fem_ssq_* is alone in its pack, so the shard can own
the entire term. A shard replaces a whole tape, so if this term ever had
to share a tape with others, the atomic AD path would become the only
option.