Laplace Approximations with adlaplaceExample
Long-form lab notebook. This vignette is intentionally detailed: it records end-to-end development checks (simulation, shard construction, derivative validation, inner Laplace solves, and outer optimization) rather than a minimal getting-started tutorial. Skim the section headings if you only need a quick overview; work through the chunks in order when reproducing or extending the example backend.
This vignette walks through a hierarchical skew-normal model using
the adlaplaceExample backend: simulate data, build an
automatic-differentiation (AD) representation as
shards, validate derivatives, solve the inner Laplace
problem for random effects, and fit outer parameters with a profile
likelihood.
Model
Observations follow a skew-normal with location \(\xi = \eta\), scale \(\omega\), and shape \(\alpha\). The log density for one observation is
\[ \ell(y;\xi,\omega,\alpha) = -\log\omega - \frac12\log(2\pi) -\frac12\left(\frac{y-\xi}{\omega}\right)^2 + \log\left[\operatorname{erfc}\left(-\frac{\alpha(y-\xi)}{\sqrt{2}\,\omega}\right)\right]. \]
The linear predictor is \(\eta = X\beta +
A\gamma\), with two independent random-effect groups in \(A\), log-standard-deviation parameters for
each group, plus skew-normal scale \(\omega\) (log scale when
transform_theta = TRUE) and shape \(\alpha\). Densities are implemented in
src/objectiveFunction.cpp; tapes for custom shards are
built via adlaplace::ad_fun_ptr() using the
package slot on each ad_data shard
(adlaplaceExample for obs/extra, adlaplace for
built-in random effects).
Simulated data
We generate 5000 observations, fixed effects, two random-effect
groups, and skew-normal scale and shape. Group SDs and \(\omega\) are stored on the log scale
(transform_theta = TRUE); \(\alpha\) is untransformed.
##
## Attaching package: 'adlaplace'## The following objects are masked from 'package:stats':
##
## binomial, gaussianlibrary(adlaplaceExample)
Nobs <- 5000
Nrandom1 <- 10
Nrandom2 <- 25
set.seed(0)
X <- Matrix::Matrix(cbind(1, rbinom(Nobs, 1, prob = 0.5)))
Adf <- data.frame(
r1 = sample(Nrandom1, Nobs, replace = TRUE),
r2 = sample(Nrandom2, Nobs, replace = TRUE)
)
AmatList <- list(
r1 = Matrix::sparseMatrix(
i = seq_len(Nobs),
j = Adf$r1
),
r2 = Matrix::sparseMatrix(
i = seq_len(Nobs),
j = Adf$r2
)
)
Amat <- do.call(cbind, AmatList)
beta <- rep_len(c(1, 0.5), ncol(X))
thetaOrig <- c(sd1 = 0.1, sd2 = 0.1, omega = 0.25, alpha = 0.8)
gamma <- rnorm(
ncol(Amat),
sd = rep(thetaOrig[1:2], c(Nrandom1, Nrandom2))
)
eta_true <- as.vector(X %*% beta + Amat %*% gamma)
if (requireNamespace("sn", quietly = TRUE)) {
y <- sn::rsn(Nobs, xi = eta_true, omega = thetaOrig["omega"], alpha = thetaOrig["alpha"])
} else {
y <- stats::rnorm(Nobs, mean = eta_true, sd = thetaOrig["omega"])
}
config <- list(
beta = beta,
theta = c(log(thetaOrig[setdiff(names(thetaOrig), "alpha")]), thetaOrig["alpha"]),
transform_theta = TRUE,
gamma = rep(0, ncol(Amat)),
shards = adlaplace::ad_shards(Amat, num_shards = 100),
num_threads = 4L,
verbose = FALSE
)
n_beta <- length(config$beta)
n_gamma <- length(config$gamma)
n_theta <- length(config$theta)Building AD shards
Each model piece is compiled separately with
ad_fun_ptr(). The full model has four shards: observation
likelihood, hyperparameter (extra) terms, and one shard per
random-effect group. Shards are merged later with
ad_fun().
Observations and hyperparameters
The observation shard uses \(\omega\) and \(\alpha\) (last two \(\theta\) components). The extra shard contributes the \(-\sum_i \log\omega\) normalization terms.
data_obs <- adlaplace::ad_data(
y = y,
A = Amat,
X = X,
beta_map = Matrix::Diagonal(n_beta),
gamma_map = Matrix::Diagonal(n_gamma),
theta_map = list(c(
which(names(thetaOrig) == "omega"),
which(names(thetaOrig) == "alpha")
), n_theta),
ad_kind = "observations",
ad_fun = "skewnormal_obs",
package = "adlaplaceExample"
)
data_extra <- adlaplace::ad_data(
y = data_obs@y,
beta_map = data_obs@beta_map,
gamma_map = data_obs@gamma_map,
theta_map = data_obs@theta_map,
ad_kind = "parameters",
ad_fun = "skewnormal_extra",
package = "adlaplaceExample"
)
ad_fun_obs <- adlaplace::ad_fun_ptr(data = data_obs, config = config)
ad_fun_extra <- adlaplace::ad_fun_ptr(data = data_extra, config = config)Random effects
Each group gets a diagonal normal prior with its own \(\theta\) index.
nr1 <- Nrandom1
nr2 <- Nrandom2
model_r1 <- adlaplace::ad_data(
beta_map = n_beta,
gamma_map = Matrix::sparseMatrix(
i = seq_len(nr1),
j = seq_len(nr1),
dims = c(n_gamma, nr1)
),
theta_map = c(1L, n_theta),
ad_kind = "random",
ad_fun = "random_diagonal",
package = "adlaplace",
precision = rep(1, nr1)
)
ad_fun_r1 <- adlaplace::ad_fun_ptr(data = model_r1, config = config)
model_r2 <- adlaplace::ad_data(
beta_map = n_beta,
gamma_map = Matrix::sparseMatrix(
i = seq.int(nr1 + 1L, length.out = nr2),
j = seq_len(nr2),
dims = c(n_gamma, nr2)
),
theta_map = c(2L, n_theta),
ad_kind = "random",
ad_fun = "random_diagonal",
package = "adlaplace",
precision = rep(1, nr2)
)
ad_fun_r2 <- adlaplace::ad_fun_ptr(data = model_r2, config = config)Observation likelihood check
Scan \(\log(\omega)\) while holding
\((\beta, \gamma, \theta_{1:2},
\alpha)\) fixed at the data-generating values. The AD joint log
density (observation + extra shards) should match sn::dsn
on the full sample when the sn package is
available.
if (!requireNamespace("sn", quietly = TRUE)) {
stop("Install package 'sn' to run the observation likelihood check.")
}
xx <- c(beta, gamma, config$theta)
log_omega_true <- config$theta["omega"]
alpha_true <- config$theta["alpha"]
omega_idx <- n_beta + n_gamma + which(names(thetaOrig) == "omega")
eta_true <- as.vector(X %*% beta + Amat %*% gamma)
sn_log_dens <- function(log_omega) {
omega <- exp(log_omega)
sum(sn::dsn(y, xi = eta_true, omega = omega, alpha = alpha_true, log = TRUE))
}
ad_sn_log_dens <- function(log_omega) {
xx_scan <- xx
xx_scan[omega_idx] <- log_omega
adlaplace::joint_log_dens(ad_fun_obs, xx_scan, negative = FALSE)
}
ad_sn_log_dens_extra <- function(log_omega) {
xx_scan <- xx
xx_scan[omega_idx] <- log_omega
adlaplace::joint_log_dens(ad_fun_extra, xx_scan, negative = FALSE)
}
# skewnormal_obs shard (logDensObs): sum_i -z_i^2 + log(erfc(-t_i)), z = (y - eta)/(omega sqrt(2)),
# t = alpha * z; log(erfc(-t)) = log(2) + log Phi(alpha * (y - eta) / omega)
sn_obs_log_dens <- function(log_omega) {
omega <- exp(log_omega)
z <- (y - eta_true) / (omega * sqrt(2))
tt <- alpha_true * z
sum(
-z^2 +
log(pracma::erfc(-tt))
)
}
# skewnormal_extra shard (logDensExtra): N * (-log omega - 1/2 log(2 pi))
sn_extra_log_dens <- function(log_omega) {
Nobs * (-log_omega - 0.5 * log(2 * pi))
}
seq_log_omega <- log_omega_true + seq(-0.5, 0.5, length.out = 7)
cmp <- data.frame(
omega = exp(seq_log_omega),
ad_obs = vapply(seq_log_omega, ad_sn_log_dens, numeric(1)),
obs_ref = vapply(seq_log_omega, sn_obs_log_dens, numeric(1)),
ad_extra = vapply(seq_log_omega, ad_sn_log_dens_extra, numeric(1)),
extra_ref = vapply(seq_log_omega, sn_extra_log_dens, numeric(1)),
sn = vapply(seq_log_omega, sn_log_dens, numeric(1))
)
cmp$ad <- cmp$ad_obs + cmp$ad_extra
cmp$sum_ref <- cmp$obs_ref + cmp$extra_ref
# cmp$diff <- cmp$ad - cmp$sn
cmp## omega ad_obs obs_ref ad_extra extra_ref sn ad
## 1 0.1516327 -6116.7116 -6116.7116 4836.7791 4836.7791 -1279.93244 -1279.93244
## 2 0.1791328 -4125.7560 -4125.7560 4003.4458 4003.4458 -122.31024 -122.31024
## 3 0.2116204 -2724.1456 -2724.1456 3170.1125 3170.1125 445.96688 445.96688
## 4 0.2500000 -1745.6223 -1745.6223 2336.7791 2336.7791 591.15685 591.15685
## 5 0.2953401 -1069.7282 -1069.7282 1503.4458 1503.4458 433.71761 433.71761
## 6 0.3489031 -609.1811 -609.1811 670.1125 670.1125 60.93138 60.93138
## 7 0.4121803 -300.8549 -300.8549 -163.2209 -163.2209 -464.07574 -464.07574
## sum_ref
## 1 -1279.93244
## 2 -122.31024
## 3 445.96688
## 4 591.15685
## 5 433.71761
## 6 60.93138
## 7 -464.07574## [1] -3.410605e-13 1.989520e-12Combining shards
Shard handles are combined with c(). This uses move
semantics: after ad_fun() is called, the individual
ad_fun_ptr objects are cleared and must not be reused.
ad_fun_plain <- c(ad_fun_obs, ad_fun_r1, ad_fun_r2, ad_fun_extra)
ad_pack <- adlaplace::ad_fun(ad_fun_plain, num_threads = config$num_threads)The per-shard log densities sum to the joint log density at a fixed parameter vector.
x_full <- c(config$beta, config$gamma, config$theta)
shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_pack@ptr))
by_shard <- vapply(
shards,
function(s) {
adlaplace::joint_log_dens(
ad_pack, x_full,
shards = s, negative = FALSE
)
},
numeric(1)
)
names(by_shard) <- c("obs", "extra", "r1", "r2")
by_shard## obs extra r1 r2 <NA> <NA>
## -37.372876 -36.185846 -44.788895 -35.690614 -49.604208 -44.050728
## <NA> <NA> <NA> <NA> <NA> <NA>
## -19.155657 -13.529529 -18.266421 -25.635001 -46.733058 -22.311467
## <NA> <NA> <NA> <NA> <NA> <NA>
## -42.500017 -30.568011 -20.087328 -22.665720 -7.600922 -32.571815
## <NA> <NA> <NA> <NA> <NA> <NA>
## -16.922736 -23.608143 -34.183124 -40.922141 -55.964326 -18.846272
## <NA> <NA> <NA> <NA> <NA> <NA>
## -22.931841 -20.324322 -52.319801 -43.722951 -14.799815 -24.201110
## <NA> <NA> <NA> <NA> <NA> <NA>
## -26.534741 -7.357971 -30.646356 -27.029511 -19.670840 -69.997814
## <NA> <NA> <NA> <NA> <NA> <NA>
## -73.601382 -18.912684 -42.515551 -12.472141 -28.442488 -17.317243
## <NA> <NA> <NA> <NA> <NA> <NA>
## -6.548248 -27.146111 -43.147575 -18.198781 -36.894560 -68.080997
## <NA> <NA> <NA> <NA> <NA> <NA>
## -13.765659 -38.301258 -30.409036 -29.901397 -22.223958 -26.956642
## <NA> <NA> <NA> <NA> <NA> <NA>
## -33.202936 -24.635023 -15.689118 -39.750814 -37.737883 -20.130089
## <NA> <NA> <NA> <NA> <NA> <NA>
## -28.201987 -22.942585 -38.170597 -44.441230 -18.121639 -14.531615
## <NA> <NA> <NA> <NA> <NA> <NA>
## -43.544988 -40.761923 -60.267746 -26.332572 -30.518310 -21.706637
## <NA> <NA> <NA> <NA> <NA> <NA>
## -22.256423 -30.495266 -22.109757 -46.143945 -23.180770 -13.996295
## <NA> <NA> <NA> <NA> <NA> <NA>
## -42.672070 -32.741418 -13.376656 -63.108027 -8.139921 -12.404384
## <NA> <NA> <NA> <NA> <NA> <NA>
## -17.849283 -14.927238 -16.427269 -38.140472 -41.776411 -25.658926
## <NA> <NA> <NA> <NA> <NA> <NA>
## -19.476649 -14.045354 -16.640235 -15.979835 -12.251755 -16.259670
## <NA> <NA> <NA> <NA> <NA> <NA>
## -36.835316 -38.960456 -30.642088 -9.294910 13.836466 34.591164
## <NA>
## 2336.779140## [1] -525.4094## [1] -525.4094Inner optimization
Given outer parameters \((\beta,
\theta)\), inner_opt() maximizes over \(\gamma\) (equivalently minimizes \(-\log p(y, \gamma \mid \beta, \theta)\)).
The returned solution should track the simulated random effects.
x_outer <- c(config$beta, config$theta)
inner_res <- adlaplace::inner_opt(
parameters = x_outer,
gamma = config$gamma,
ad_fun = ad_pack,
control = list(
maxit = 100,
report.level = 0,
report.freq = 0
),
deriv = TRUE,
verbose = FALSE
)
plot(
inner_res$opt$solution, gamma,
xlab = "estimated", ylab = "true",
main = "Inner mode vs simulated gamma"
)
abline(0, 1, col = "red", lty = 2)
Profile likelihood
log_lik_laplace() runs the inner optimizer and returns
the Laplace approximation to the marginal log likelihood. With
deriv = TRUE, it also returns the profile gradient and
third-order trace terms used in the determinant correction.
res_noderiv <- adlaplace::log_lik_laplace(
x = x_outer,
ad_fun = ad_pack,
config = config,
deriv = FALSE
)## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -642.693699 242.051935 Continuing 1.000000 8 Reached tolerance
## 2 -645.211652 0.560731 Continuing 1.000000 7 Reached tolerance
## 3 -645.211667 0.000028 Continuing 1.000000 7 Reached tolerance
## 4 -645.211667 0.000000 Continuing 1.000000 7 Reached tolerance
##
## Iteration has terminated
## 4 -645.211667 0.000000 Successres_deriv <- adlaplace::log_lik_laplace(
x = x_outer,
ad_fun = ad_pack,
config = config,
deriv = TRUE
)## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -642.693699 242.051935 Continuing 1.000000 8 Reached tolerance
## 2 -645.211652 0.560731 Continuing 1.000000 7 Reached tolerance
## 3 -645.211667 0.000028 Continuing 1.000000 7 Reached tolerance
## 4 -645.211667 0.000028 Continuing - TR contract 0.500000 7 Reached tolerance
## 5 -645.211667 0.000028 Continuing - TR contract 0.250000 7 Reached tolerance
## 6 -645.211667 0.000028 Continuing - TR contract 0.125000 7 Reached tolerance
## 7 -645.211667 0.000028 Continuing - TR contract 0.062500 7 Reached tolerance
## 8 -645.211667 0.000028 Continuing - TR contract 0.031250 7 Reached tolerance
## 9 -645.211667 0.000028 Continuing - TR contract 0.015625 7 Reached tolerance
## 10 -645.211667 0.000028 Continuing - TR contract 0.007812 7 Reached tolerance
##
## iter f nrm_gr status radCG iter CG result
## 11 -645.211667 0.000028 Continuing - TR contract 0.003906 7 Reached tolerance
## 12 -645.211667 0.000028 Continuing - TR contract 0.001953 7 Reached tolerance
## 13 -645.211667 0.000028 Continuing - TR contract 0.000977 7 Reached tolerance
## 14 -645.211667 0.000028 Continuing - TR contract 0.000488 7 Reached tolerance
## 15 -645.211667 0.000028 Continuing - TR contract 0.000244 7 Reached tolerance
## 16 -645.211667 0.000028 Continuing - TR contract 0.000122 7 Reached tolerance
## 17 -645.211667 0.000028 Continuing - TR contract 0.000061 7 Reached tolerance
## 18 -645.211667 0.000028 Continuing - TR contract 0.000031 7 Reached tolerance
## 19 -645.211667 0.000028 Continuing - TR contract 0.000015 7 Reached tolerance
## 20 -645.211667 0.000028 Continuing - TR contract 0.000008 7 Reached tolerance
##
## iter f nrm_gr status radCG iter CG result
## 21 -645.211667 0.000028 Continuing - TR contract 0.000004 7 Reached tolerance
## 22 -645.211667 0.000028 Continuing - TR contract 0.000002 7 Reached tolerance
## 23 -645.211667 0.000028 Continuing - TR contract 0.000001 7 Reached tolerance
## 24 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 25 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 26 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 27 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 28 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 29 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 30 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
##
## Iteration has terminated
##
## iter f nrm_gr status
## 30 -645.211667 0.000028Radius of trust region is less than stop.trust.radiusdata.frame(
log_lik = c(res_noderiv$log_lik, res_deriv$log_lik),
neg_log_lik = c(res_noderiv$neg_log_lik, res_deriv$neg_log_lik),
inner_fval = c(res_deriv$opt$fval, res_deriv$opt$fval),
row.names = c("deriv = FALSE", "deriv = TRUE")
)## log_lik neg_log_lik inner_fval
## deriv = FALSE 527.9328 -527.9328 -645.2117
## deriv = TRUE 527.9328 -527.9328 -645.2117## [1] -44.737371 53.750083 -3.191366 4.015087 81.246438 -31.442109Outer optimization
outer_fn() and outer_gr() wrap
log_lik_laplace() for use with stats::optim().
A small cache environment stores the current inner \(\gamma\) start.
cache <- new.env(parent = emptyenv())
cache$gamma <- rep(0, n_gamma)
control_inner <- list(
maxit = 100,
report.level = 0,
report.freq = 0
)
optim_control <- list(maxit = 2, trace = 0)
bounds <- list(lower = c(-2, -2, -5, -5, -2, -2), upper = c(2, 2, 0, 0, 3, 2))
common_optim_args <- c(
list(
par = x_outer,
fn = adlaplace::outer_fn,
method = "L-BFGS-B",
control = optim_control,
config = config,
ad_fun = ad_pack,
cache = cache,
control_inner = control_inner
),
bounds
)
outer_fit_no_g <- do.call(stats::optim, common_optim_args)
outer_fit <- do.call(stats::optim, c(common_optim_args, list(gr = adlaplace::outer_gr)))
true_outer <- c(beta, log(thetaOrig[1:3]), thetaOrig[4])
fit_cmp <- rbind(
optim = outer_fit$par,
optim_no_grad = outer_fit_no_g$par,
true = true_outer
)
colnames(fit_cmp) <- c(
paste0("beta[", seq_len(n_beta), "]"),
names(thetaOrig)
)
fit_cmp## beta[1] beta[2] sd1 sd2 omega alpha
## optim 1.015456 0.4971901 -2.290575 -2.316767 -1.395689 0.8099855
## optim_no_grad 1.015456 0.4971907 -2.290575 -2.316767 -1.395689 0.8099854
## true 1.000000 0.5000000 -2.302585 -2.302585 -1.386294 0.8000000Profile derivative checks
Finite-difference checks along one outer coordinate validate the profile gradient, the determinant derivative, and the chain rule through the inner mode \(u(\beta, \theta)\).
if (!requireNamespace("abind", quietly = TRUE)) {
stop("Install suggested package 'abind' to run derivative checks.")
}
par(mfrow = c(3, 2), mar = c(2, 2, 2, 0), mgp = c(1, 0.5, 0))
x1 <- outer_fit$par
Dpar <- 3 # length(x1) - 1
Ngrid <- 7L
par_grid <- matrix(x1, nrow = Ngrid, ncol = length(x1), byrow = TRUE)
Sx <- x1[Dpar] + seq(-0.5, 0.5, length.out = Ngrid)
SxD <- Sx[-1] - diff(Sx) / 2
par_grid[, Dpar] <- Sx
res_scan <- lapply(
split(par_grid, row(par_grid)),
adlaplace::log_lik_laplace,
ad_fun = ad_pack,
config = config,
deriv = TRUE,
gamma = cache$gamma
)## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -642.103294 0.028812 Continuing 1.000000 10 Reached tolerance
## 2 -642.103294 0.000003 Continuing 1.000000 9 Reached tolerance
##
## Iteration has terminated
## 2 -642.103294 0.000003 Success
##
## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -644.665312 0.009224 Continuing 1.000000 10 Reached tolerance
## 2 -644.665312 0.000001 Continuing 1.000000 10 Reached tolerance
##
## Iteration has terminated
## 2 -644.665312 0.000001 Success
##
## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -646.122073 0.001802 Continuing 1.000000 10 Reached tolerance
## 2 -646.122073 0.000000 Continuing 1.000000 11 Reached tolerance
##
## Iteration has terminated
## 2 -646.122073 0.000000 Success
##
## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -646.756455 0.000000 Continuing - TR contract 0.500000 11 Reached tolerance
##
## Iteration has terminated
## 1 -646.756455 0.000000 Success
##
## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -646.779757 0.001014 Continuing 1.000000 10 Reached tolerance
## 2 -646.779757 0.000000 Continuing 1.000000 11 Reached tolerance
##
## Iteration has terminated
## 2 -646.779757 0.000000 Success
##
## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -646.349904 0.002760 Continuing 1.000000 10 Reached tolerance
## 2 -646.349904 0.000000 Continuing 1.000000 11 Reached tolerance
##
## Iteration has terminated
## 2 -646.349904 0.000000 Success
##
## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -645.585102 0.004623 Continuing 1.000000 10 Reached tolerance
## 2 -645.585102 0.000000 Continuing 1.000000 11 Reached tolerance
##
## Iteration has terminated
## 2 -645.585102 0.000000 SuccessSnegLik <- vapply(res_scan, `[[`, numeric(1), "neg_log_lik")
Sdet <- vapply(res_scan, function(r) r$extra$hessian$half_log_det, numeric(1))
grad_mat <- do.call(rbind, lapply(res_scan, `[[`, "grad"))
extra_df <- do.call(abind::abind, c(lapply(res_scan, `[[`, "deriv"), along = 3))
dU <- do.call(
abind::abind,
c(lapply(res_scan, function(r) as.matrix(r$extra$dU)), along = 3)
)
u_hat <- do.call(rbind, lapply(res_scan, function(r) r$opt$solution))
res_mid <- res_scan[[(Ngrid + 1L) %/% 2L]]
plot(Sx, SnegLik, type = "l", xlab = expression(alpha), ylab = "-log lik")
plot(Sx, grad_mat[, Dpar], type = "l", xlab = expression(alpha), ylab = "grad")
points(SxD, diff(SnegLik) / diff(Sx), pch = 16)
legend("bottomright", legend = c("AD", "finite diff"), lty = c(1, NA), pch = c(NA, 1), bty = "n")
Du <- 2L
plot(Sx, u_hat[, Du], type = "l", xlab = expression(alpha), ylab = expression(u[Du]))
diff_here <- diff(u_hat[, Du]) / diff(Sx)
plot(
Sx, dU[Du, Dpar, ],
type = "l",
ylim = range(c(dU[Du, Dpar, ], diff_here), na.rm = TRUE),
xlab = expression(alpha),
ylab = expression(d * u[Du] / d * alpha)
)
points(SxD, diff_here, pch = 16)
plot(Sx, Sdet, type = "l", xlab = expression(alpha), ylab = "half log det")
diff_here <- diff(Sdet) / diff(Sx)
plot(
Sx, extra_df[Dpar, "d_det", ],
type = "l",
ylim = range(c(extra_df[Dpar, "d_det", ], diff_here), na.rm = TRUE),
xlab = expression(alpha),
ylab = expression(d / d * alpha ~ log ~ det)
)
points(SxD, diff_here, pch = 16)
hessian_outer <- Matrix::forceSymmetric(res_mid$extra$hessian$outer)
hessian_plain <- adlaplace::hessian(
ad_fun_plain,
res_mid$full_parameters,
inner = FALSE,
negative = TRUE
)
max_abs_hess_diff <- max(abs((hessian_outer - hessian_plain)@x), na.rm = TRUE)
max_abs_hess_diff## [1] 1.455192e-11Joint density derivative checks
The same finite-difference idea applies to the full joint log density \(-\log p(\beta, \gamma, \theta \mid y)\) and its AD gradient and Hessian.
par(mfrow = c(2, 3), mar = c(2, 2, 2, 0), mgp = c(1, 0.5, 0))
x_here <- res_mid$full_parameters
Dpar_dens <- length(x_here) - 1
Ngrid <- 11L
shards <- seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_fun_plain))
par_grid <- matrix(x_here, nrow = Ngrid, ncol = length(x_here), byrow = TRUE)
Sx <- x_here[Dpar_dens] + seq(-0.1, 0.1, length.out = Ngrid)
SxD <- Sx[-1] - diff(Sx) / 2
par_grid[, Dpar_dens] <- Sx
x_list <- split(par_grid, row(par_grid))
dens <- vapply(
x_list,
adlaplace::joint_log_dens,
numeric(1),
ad_fun = ad_fun_plain,
shards = shards,
negative = FALSE
)
grad <- do.call(
cbind,
lapply(
x_list,
adlaplace::grad,
ad_fun = ad_fun_plain,
inner = FALSE,
shards = shards,
negative = FALSE
)
)
plot(Sx, grad[Dpar_dens, ], type = "l", ylab = "AD gradient")
points(SxD, diff(dens) / diff(Sx), pch = 16)
legend("topright", legend = c("AD", "finite diff"), lty = c(1, NA), pch = c(NA, 1), bty = "n")
hes <- array(
dim = c(length(x_here), length(x_here), Ngrid),
dimnames = list(NULL, NULL, NULL)
)
for (i in seq_len(Ngrid)) {
hes[, , i] <- as.matrix(adlaplace::hessian(
ad_fun_plain,
x_list[[i]],
inner = FALSE,
shards = shards,
negative = FALSE
))
}
grad_slope <- apply(grad, 1, diff) / mean(diff(Sx))
for (Dpar2 in c(Dpar_dens, 1L, 5, 6, length(x_here))) {
plot(
Sx, hes[Dpar_dens, Dpar2, ],
type = "l",
ylab = paste0("H[", Dpar_dens, ",", Dpar2, "]"),
ylim = range(c(hes[Dpar_dens, Dpar2, ], grad_slope[, Dpar2]), na.rm = TRUE)
)
points(SxD, grad_slope[, Dpar2], pch = 16)
}
Nicer interface
Code verifying the omnibus interface gives the same results as the barebones implementation.
dat <- cbind(
data.frame(y = y),
Adf,
x = X[, 2]
)
md <- adlaplace::model_data(
data = dat,
formula =
skewnormal(y) ~
x + adlaplace::iid(r1, init = 0.1) +
adlaplace::iid(r2, init = 0.1)
)
md$data$info$theta[, c("label", "transform", "id")]## label transform id
## skewnormal(y).1 y_skewnormal_omega TRUE 0
## skewnormal(y).2 y_skewnormal_alpha FALSE 1
## r1 r1_iid TRUE 2
## r2 r2_iid TRUE 3## 4 x 2 sparse Matrix of class "ngCMatrix"
##
## [1,] | .
## [2,] . |
## [3,] . .
## [4,] . .## 4 x 2 sparse Matrix of class "ngCMatrix"
##
## [1,] . .
## [2,] . .
## [3,] | .
## [4,] . |## [1] "sd1" "sd2" "omega" "alpha"# order of parameters is different, previously alpha was last
# also intercept is now second
param_remap <- c(2, 1, 5, 6, 3, 4)
xx_remap <- res_deriv$parameters[param_remap]
config2 <- list(
transform_theta = TRUE,
shards = adlaplace::ad_shards(
md$data$A,
num_shards = 100
),
beta = xx_remap[seq_len(n_beta)],
theta = xx_remap[-seq_len(n_beta)],
gamma = rep(0, n_gamma),
verbose = TRUE
)
xx_all <- c(
res_deriv$parameters[seq_len(n_beta)],
rep(0, n_gamma),
res_deriv$parameters[-seq_len(n_beta)]
)
xx_all2 <- c(
xx_remap[seq_len(n_beta)],
rep(0, n_gamma),
xx_remap[-seq_len(n_beta)]
)
ad_extra <- adlaplace::ad_fun_ptr(data = data_extra, config = config)
ad_extra2 <- adlaplace::ad_fun_ptr(data = md$parameters$y_extra, config = config2)## build_ad_fun_parameters: taping...
## build_ad_fun_parameters: computing sparsity...
## discovering Hessian sparsity pattern
## sparsity: grad 1, grad_inner 0, hes full0, hes upper0, hes_inner 0## ad_fun: building 4 density shard(s) (100 observation group(s) per obs shard)...
## [1/4] y (CppAD tape, 100 groups)...
## build_ad_fun_obs groups 100
## taping observation group 1 / 100
## taping observation group 2 / 100
## taping observation group 3 / 100
## taping observation group 4 / 100
## taping observation group 5 / 100
## taping observation group 6 / 100
## taping observation group 7 / 100
## taping observation group 8 / 100
## taping observation group 9 / 100
## taping observation group 10 / 100
## taping observation group 11 / 100
## taping observation group 12 / 100
## taping observation group 13 / 100
## taping observation group 14 / 100
## taping observation group 15 / 100
## taping observation group 16 / 100
## taping observation group 17 / 100
## taping observation group 18 / 100
## taping observation group 19 / 100
## taping observation group 20 / 100
## taping observation group 21 / 100
## taping observation group 22 / 100
## taping observation group 23 / 100
## taping observation group 24 / 100
## taping observation group 25 / 100
## taping observation group 26 / 100
## taping observation group 27 / 100
## taping observation group 28 / 100
## taping observation group 29 / 100
## taping observation group 30 / 100
## taping observation group 31 / 100
## taping observation group 32 / 100
## taping observation group 33 / 100
## taping observation group 34 / 100
## taping observation group 35 / 100
## taping observation group 36 / 100
## taping observation group 37 / 100
## taping observation group 38 / 100
## taping observation group 39 / 100
## taping observation group 40 / 100
## taping observation group 41 / 100
## taping observation group 42 / 100
## taping observation group 43 / 100
## taping observation group 44 / 100
## taping observation group 45 / 100
## taping observation group 46 / 100
## taping observation group 47 / 100
## taping observation group 48 / 100
## taping observation group 49 / 100
## taping observation group 50 / 100
## taping observation group 51 / 100
## taping observation group 52 / 100
## taping observation group 53 / 100
## taping observation group 54 / 100
## taping observation group 55 / 100
## taping observation group 56 / 100
## taping observation group 57 / 100
## taping observation group 58 / 100
## taping observation group 59 / 100
## taping observation group 60 / 100
## taping observation group 61 / 100
## taping observation group 62 / 100
## taping observation group 63 / 100
## taping observation group 64 / 100
## taping observation group 65 / 100
## taping observation group 66 / 100
## taping observation group 67 / 100
## taping observation group 68 / 100
## taping observation group 69 / 100
## taping observation group 70 / 100
## taping observation group 71 / 100
## taping observation group 72 / 100
## taping observation group 73 / 100
## taping observation group 74 / 100
## taping observation group 75 / 100
## taping observation group 76 / 100
## taping observation group 77 / 100
## taping observation group 78 / 100
## taping observation group 79 / 100
## taping observation group 80 / 100
## taping observation group 81 / 100
## taping observation group 82 / 100
## taping observation group 83 / 100
## taping observation group 84 / 100
## taping observation group 85 / 100
## taping observation group 86 / 100
## taping observation group 87 / 100
## taping observation group 88 / 100
## taping observation group 89 / 100
## taping observation group 90 / 100
## taping observation group 91 / 100
## taping observation group 92 / 100
## taping observation group 93 / 100
## taping observation group 94 / 100
## taping observation group 95 / 100
## taping observation group 96 / 100
## taping observation group 97 / 100
## taping observation group 98 / 100
## taping observation group 99 / 100
## taping observation group 100 / 100
## sparsity observation group 1 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 2 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 3 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 4 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 5 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 6 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 7 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 8 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 9 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 10 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 11 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 12 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 13 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 14 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 15 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 16 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 17 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 18 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 19 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 20 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 21 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 22 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 23 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 24 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 25 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 26 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 27 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 28 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 29 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 30 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 31 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 7
## sparsity observation group 32 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 33 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 34 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 35 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 36 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 37 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 38 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 39 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 40 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 41 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 42 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 43 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full45, hes upper26, hes_inner 5
## sparsity observation group 44 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 45 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 46 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 47 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 48 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 49 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 50 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 51 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 52 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 53 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 54 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 55 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 56 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 57 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 58 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 59 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 60 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 61 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 62 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 63 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 64 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 65 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full63, hes upper36, hes_inner 8
## sparsity observation group 66 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 67 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 68 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 69 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 70 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 71 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 72 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 73 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 74 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 75 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 76 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 10, grad_inner 6, hes full76, hes upper43, hes_inner 9
## sparsity observation group 77 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 7
## sparsity observation group 78 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 79 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 80 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 81 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 82 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 83 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 84 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 85 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 86 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 87 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 88 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 89 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full67, hes upper38, hes_inner 8
## sparsity observation group 90 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 91 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full58, hes upper33, hes_inner 7
## sparsity observation group 92 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 93 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 7, grad_inner 3, hes full47, hes upper27, hes_inner 5
## sparsity observation group 94 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 95 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 96 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 9, grad_inner 5, hes full63, hes upper36, hes_inner 8
## sparsity observation group 97 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 8, grad_inner 4, hes full56, hes upper32, hes_inner 6
## sparsity observation group 98 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 6, grad_inner 2, hes full36, hes upper21, hes_inner 3
## sparsity observation group 99 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 6, grad_inner 2, hes full36, hes upper21, hes_inner 3
## sparsity observation group 100 / 100
## discovering Hessian sparsity pattern
## sparsity: grad 6, grad_inner 2, hes full36, hes upper21, hes_inner 3
## [1/4] done (100 AD group(s)).
## [2/4] r1_iid (CppAD tape)...
## build_ad_fun_random: taping...
## theta index 39 logVariance -2.30259 precision 100
## random_diagonal n_gamma 10 qDet -13.8365 qpart 0
## build_ad_fun_random: computing sparsity...
## using external Hessian sparsity pattern
## sparsity: grad 11, grad_inner 10, hes full31, hes upper21, hes_inner 10
## [2/4] done (1 AD group(s)).
## [3/4] r2_iid (CppAD tape)...
## build_ad_fun_random: taping...
## theta index 40 logVariance -2.30259 precision 100
## random_diagonal n_gamma 25 qDet -34.5912 qpart 0
## build_ad_fun_random: computing sparsity...
## using external Hessian sparsity pattern
## sparsity: grad 26, grad_inner 25, hes full76, hes upper51, hes_inner 25
## [3/4] done (1 AD group(s)).
## [4/4] y_extra (CppAD tape)...
## build_ad_fun_parameters: taping...
## build_ad_fun_parameters: computing sparsity...
## discovering Hessian sparsity pattern
## sparsity: grad 1, grad_inner 0, hes full0, hes upper0, hes_inner 0
## [4/4] done (1 AD group(s)).
## ad_fun: merging density handles...
## ad_fun: attaching Hessian templates (103 AD group(s), 1 thread(s))...
## assigning owner threads...
## collecting sparsity patterns (beta=2, gamma=35, theta=4)...
## sparsity group 1/103
## sparsity group 2/103
## sparsity group 3/103
## sparsity group 4/103
## sparsity group 5/103
## sparsity group 6/103
## sparsity group 7/103
## sparsity group 8/103
## sparsity group 9/103
## sparsity group 10/103
## sparsity group 11/103
## sparsity group 12/103
## sparsity group 13/103
## sparsity group 14/103
## sparsity group 15/103
## sparsity group 16/103
## sparsity group 17/103
## sparsity group 18/103
## sparsity group 19/103
## sparsity group 20/103
## sparsity group 21/103
## sparsity group 22/103
## sparsity group 23/103
## sparsity group 24/103
## sparsity group 25/103
## sparsity group 26/103
## sparsity group 27/103
## sparsity group 28/103
## sparsity group 29/103
## sparsity group 30/103
## sparsity group 31/103
## sparsity group 32/103
## sparsity group 33/103
## sparsity group 34/103
## sparsity group 35/103
## sparsity group 36/103
## sparsity group 37/103
## sparsity group 38/103
## sparsity group 39/103
## sparsity group 40/103
## sparsity group 41/103
## sparsity group 42/103
## sparsity group 43/103
## sparsity group 44/103
## sparsity group 45/103
## sparsity group 46/103
## sparsity group 47/103
## sparsity group 48/103
## sparsity group 49/103
## sparsity group 50/103
## sparsity group 51/103
## sparsity group 52/103
## sparsity group 53/103
## sparsity group 54/103
## sparsity group 55/103
## sparsity group 56/103
## sparsity group 57/103
## sparsity group 58/103
## sparsity group 59/103
## sparsity group 60/103
## sparsity group 61/103
## sparsity group 62/103
## sparsity group 63/103
## sparsity group 64/103
## sparsity group 65/103
## sparsity group 66/103
## sparsity group 67/103
## sparsity group 68/103
## sparsity group 69/103
## sparsity group 70/103
## sparsity group 71/103
## sparsity group 72/103
## sparsity group 73/103
## sparsity group 74/103
## sparsity group 75/103
## sparsity group 76/103
## sparsity group 77/103
## sparsity group 78/103
## sparsity group 79/103
## sparsity group 80/103
## sparsity group 81/103
## sparsity group 82/103
## sparsity group 83/103
## sparsity group 84/103
## sparsity group 85/103
## sparsity group 86/103
## sparsity group 87/103
## sparsity group 88/103
## sparsity group 89/103
## sparsity group 90/103
## sparsity group 91/103
## sparsity group 92/103
## sparsity group 93/103
## sparsity group 94/103
## sparsity group 95/103
## sparsity group 96/103
## sparsity group 97/103
## sparsity group 98/103
## sparsity group 99/103
## sparsity group 100/103
## sparsity group 101/103
## sparsity group 102/103
## sparsity group 103/103
## building Hessian map...
## building trace column map...
## attaching Hessian templates to C++ handle...
## ad_fun: Hessian attach complete.## [1] 2336.779## [1] 2336.779adlaplace::joint_log_dens(ad_pack2, xx_all2,
shards = adlaplace:::n_groups(ad_pack2@ptr) - 1L,
negative = FALSE
)## [1] 2336.779by_shard <- vapply(
seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_pack@ptr)),
function(s) adlaplace::joint_log_dens(ad_pack, xx_all, shards = s, negative = FALSE),
numeric(1)
)
by_shard2 <- vapply(
seq.int(from = 0L, length.out = adlaplace:::n_groups(ad_pack2@ptr)),
function(s) adlaplace::joint_log_dens(ad_pack2, xx_all2, shards = s, negative = FALSE),
numeric(1)
)
rbind(by_shard, by_shard2)[, seq(to = length(by_shard), len = 6)]## [,1] [,2] [,3] [,4] [,5] [,6]
## by_shard -38.96046 -30.64209 -9.29491 13.83647 34.59116 2336.779
## by_shard2 -38.96046 -30.64209 -9.29491 13.83647 34.59116 2336.779## [1] -525.4094 -525.4094## [1] -525.4094log_lik1 <- adlaplace::log_lik_laplace(
x = res_deriv$parameters,
config = config,
ad_fun = ad_pack,
deriv = TRUE
)## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -642.693699 242.051935 Continuing 1.000000 8 Reached tolerance
## 2 -645.211652 0.560731 Continuing 1.000000 7 Reached tolerance
## 3 -645.211667 0.000028 Continuing 1.000000 7 Reached tolerance
## 4 -645.211667 0.000028 Continuing - TR contract 0.500000 7 Reached tolerance
## 5 -645.211667 0.000028 Continuing - TR contract 0.250000 7 Reached tolerance
## 6 -645.211667 0.000028 Continuing - TR contract 0.125000 7 Reached tolerance
## 7 -645.211667 0.000028 Continuing - TR contract 0.062500 7 Reached tolerance
## 8 -645.211667 0.000028 Continuing - TR contract 0.031250 7 Reached tolerance
## 9 -645.211667 0.000028 Continuing - TR contract 0.015625 7 Reached tolerance
## 10 -645.211667 0.000028 Continuing - TR contract 0.007812 7 Reached tolerance
##
## iter f nrm_gr status radCG iter CG result
## 11 -645.211667 0.000028 Continuing - TR contract 0.003906 7 Reached tolerance
## 12 -645.211667 0.000028 Continuing - TR contract 0.001953 7 Reached tolerance
## 13 -645.211667 0.000028 Continuing - TR contract 0.000977 7 Reached tolerance
## 14 -645.211667 0.000028 Continuing - TR contract 0.000488 7 Reached tolerance
## 15 -645.211667 0.000028 Continuing - TR contract 0.000244 7 Reached tolerance
## 16 -645.211667 0.000028 Continuing - TR contract 0.000122 7 Reached tolerance
## 17 -645.211667 0.000028 Continuing - TR contract 0.000061 7 Reached tolerance
## 18 -645.211667 0.000028 Continuing - TR contract 0.000031 7 Reached tolerance
## 19 -645.211667 0.000028 Continuing - TR contract 0.000015 7 Reached tolerance
## 20 -645.211667 0.000028 Continuing - TR contract 0.000008 7 Reached tolerance
##
## iter f nrm_gr status radCG iter CG result
## 21 -645.211667 0.000028 Continuing - TR contract 0.000004 7 Reached tolerance
## 22 -645.211667 0.000028 Continuing - TR contract 0.000002 7 Reached tolerance
## 23 -645.211667 0.000028 Continuing - TR contract 0.000001 7 Reached tolerance
## 24 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 25 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 26 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 27 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 28 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 29 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 30 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
##
## Iteration has terminated
##
## iter f nrm_gr status
## 30 -645.211667 0.000028Radius of trust region is less than stop.trust.radiuslog_lik2 <- adlaplace::log_lik_laplace(
x = xx_remap,
config = config2,
ad_fun = ad_pack2,
deriv = TRUE
)## inner_opt: threads = 1, shards = 103, params = 41 (beta = 2, gamma = 35, theta = 4)
## Beginning optimization
##
## iter f nrm_gr status radCG iter CG result
## 1 -642.693699 242.051935 Continuing 1.000000 8 Reached tolerance
## 2 -645.211652 0.560731 Continuing 1.000000 7 Reached tolerance
## 3 -645.211667 0.000028 Continuing 1.000000 7 Reached tolerance
## 4 -645.211667 0.000028 Continuing - TR contract 0.500000 7 Reached tolerance
## 5 -645.211667 0.000028 Continuing - TR contract 0.250000 7 Reached tolerance
## 6 -645.211667 0.000028 Continuing - TR contract 0.125000 7 Reached tolerance
## 7 -645.211667 0.000028 Continuing - TR contract 0.062500 7 Reached tolerance
## 8 -645.211667 0.000028 Continuing - TR contract 0.031250 7 Reached tolerance
## 9 -645.211667 0.000028 Continuing - TR contract 0.015625 7 Reached tolerance
## 10 -645.211667 0.000028 Continuing - TR contract 0.007812 7 Reached tolerance
##
## iter f nrm_gr status radCG iter CG result
## 11 -645.211667 0.000028 Continuing - TR contract 0.003906 7 Reached tolerance
## 12 -645.211667 0.000028 Continuing - TR contract 0.001953 7 Reached tolerance
## 13 -645.211667 0.000028 Continuing - TR contract 0.000977 7 Reached tolerance
## 14 -645.211667 0.000028 Continuing - TR contract 0.000488 7 Reached tolerance
## 15 -645.211667 0.000028 Continuing - TR contract 0.000244 7 Reached tolerance
## 16 -645.211667 0.000028 Continuing - TR contract 0.000122 7 Reached tolerance
## 17 -645.211667 0.000028 Continuing - TR contract 0.000061 7 Reached tolerance
## 18 -645.211667 0.000028 Continuing - TR contract 0.000031 7 Reached tolerance
## 19 -645.211667 0.000028 Continuing - TR contract 0.000015 7 Reached tolerance
## 20 -645.211667 0.000028 Continuing - TR contract 0.000008 7 Reached tolerance
##
## iter f nrm_gr status radCG iter CG result
## 21 -645.211667 0.000028 Continuing - TR contract 0.000004 7 Reached tolerance
## 22 -645.211667 0.000028 Continuing - TR contract 0.000002 7 Reached tolerance
## 23 -645.211667 0.000028 Continuing - TR contract 0.000001 7 Reached tolerance
## 24 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 25 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 26 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 27 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 28 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 29 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
## 30 -645.211667 0.000028 Continuing - TR contract 0.000000 7 Reached tolerance
##
## Iteration has terminated
##
## iter f nrm_gr status
## 30 -645.211667 0.000028Radius of trust region is less than stop.trust.radius
##
## trace_hinv_t: threads = 1, shards = 103, params = 41
## trace_hinv_t: finished
## inner_opt: finished## [1] 527.9328 527.9328## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 53.75008 -44.73737 81.24644 -31.44211 -3.191366 4.015087
## [2,] 53.75008 -44.73737 81.24644 -31.44211 -3.191366 4.015087