EnvelopeEval: Evaluate Negative Log-Likelihood and Gradients
In glmbayes: Bayesian Generalized Linear Models (IID Samples)

EnvelopeEval

R Documentation

Evaluate Negative Log-Likelihood and Gradients

Description

EnvelopeEval() evaluates the negative log-likelihood and gradients at a grid of parameter values, optionally using OpenCL acceleration.

Usage

EnvelopeEval(G4, y, x, mu, P, alpha, wt,
                    family, link,
                    use_opencl = FALSE, verbose = FALSE)

Arguments

`G4`	Numeric matrix of parameter values (parameters * grid points).
`y`	Numeric response vector.
`x`	Numeric design matrix.
`mu`	Numeric matrix of offsets or prior means.
`P`	Numeric matrix representing the portion of the prior precision shifted into the likelihood.
`alpha`	Numeric offset vector of length `m`
`wt`	Numeric vector of weights.
`family`	Character string; model family (e.g. `"gaussian"`).
`link`	Character string; link function (e.g. `"identity"`).
`use_opencl`	Logical; if `TRUE`, attempt OpenCL acceleration.
`verbose`	Logical; if `TRUE`, print diagnostic output.

Details

The lower-level helpers f2_f3_non_opencl and f2_f3_opencl are internal C++ kernels used by the CPU and OpenCL backends. The internal routine run_opencl_pilot benchmarks OpenCL performance on a pilot subset of the grid to estimate runtime before full evaluation.

These functions implement the grid evaluation logic used in envelope construction for rejection sampling. They make use of the theory described in \insertCiteNygren2006glmbayes and the general implementation outlined in \insertCiteglmbayesSimmethodsglmbayes.

The evaluation workflow has several layers: 1. High-level dispatch (EnvelopeEval)

EnvelopeEval() is the user-facing entry point. It accepts a grid of parameter values (G4) and the data (y, x, mu, P, alpha, wt).
If the grid is large (>= 14 columns), it first calls run_opencl_pilot to benchmark OpenCL performance and optionally report estimated runtime.
It then dispatches to either the CPU or GPU backend:
- If use_opencl = TRUE and the family is not "gaussian", it calls f2_f3_opencl (an internal C++ kernel).
- Otherwise, it calls f2_f3_non_opencl (the CPU kernel).

2. CPU backend (f2_f3_non_opencl)

This function evaluates the negative log-likelihood and gradients using standard CPU routines.
It inspects the family and link arguments and routes to the correct pair of kernels (⁠f2_*⁠ for the likelihood, ⁠f3_*⁠ for the gradient).
For example:
- "binomial" with "logit" calls f2_binomial_logit() and f3_binomial_logit().
- "poisson" calls f2_poisson() and f3_poisson().
- "gaussian" calls f2_gaussian() and f3_gaussian().
These kernels ultimately rely on the same C math routines that R itself uses (from the nmath/rmath libraries), ensuring numerical consistency with base R functions like dnorm, dpois, etc.

3. GPU backend (f2_f3_opencl)

This function mirrors the CPU backend but executes the likelihood and gradient calculations on an OpenCL device (GPU or CPU).
It flattens the input matrices/vectors and allocates output buffers.
It then constructs a full OpenCL program by concatenating:
- a generic OpenCL support header (OPENCL.CL),
- OpenCL ports of R's rmath, nmath, and dpq libraries,
- and the family/link-specific kernel source (e.g. f2_f3_binomial_logit.cl).
The resulting program is compiled and passed to a kernel runner (f2_f3_kernel_runner) which executes the likelihood and gradient calculations in parallel on the device.
This ensures that the GPU backend produces results consistent with the CPU backend, but can scale to much larger grids efficiently.

4. Pilot timing (run_opencl_pilot)

This helper runs a small subset of the grid through the OpenCL backend to estimate runtime.
It is used by EnvelopeEval() to inform users (when verbose = TRUE) whether OpenCL acceleration is likely to be beneficial.

5. Returned values

All backends return a list with:
- NegLL: numeric vector of negative log-likelihood values.
- cbars: numeric matrix of gradients (parameters * grid points).

6. Role of likelihood and gradients in sampling

The outputs of EnvelopeEval() - the negative log-likelihood values (NegLL) and the gradient matrix (cbars) - are not endpoints in themselves. They form the envelope used in the rejection sampler implemented by internal functions such as .rNormalGLM_std_cpp().
This routine is called by .rNormalGLM_cpp(), which underlies the user-facing function rNormal_reg(). Together they implement envelope-based posterior sampling for GLMs with log-concave likelihoods and multivariate normal priors.

7. Simulation execution (accept/reject procedure)

The acceptance test is performed using

\log(U_2) \leq \log f(y \mid \theta_i) - \Big(\log f(y \mid \bar{\theta}_{J(i)}) - c(\bar{\theta}_{J(i)})^T(\theta_i - \bar{\theta}_{J(i)})\Big) \leq 0

Connections between code and notation:

The arguments G4 (in EnvelopeEval) and b (in ⁠f2_f3_*⁠) both represent the grid of tangency points \bar{\theta}_j.
The output NegLL corresponds to -\log f(y \mid \bar{\theta}_{J(i)}), i.e. the negative log-likelihood evaluated at each tangency point.
The output cbars corresponds to the subgradient vectors c(\bar{\theta}_{J(i)}), which define the tangent hyperplanes used in the envelope construction.

Precomputation for efficiency:

Both NegLL and cbars are computed once during envelope construction, prior to the simulation stage.
This means the sampler does not need to recompute likelihoods or gradients at every candidate draw - it simply reuses the stored values (NegLL, cbars, and LLconst) in the acceptance inequality.

This design ensures that the envelope is tangent to the log-likelihood at each \bar{\theta}_j, lies above it elsewhere, and that the accept-reject procedure can run efficiently while still producing samples from the true posterior \pi(\theta \mid y).

Value

EnvelopeEval: List with components NegLL (numeric vector of negative log-likelihood values) and cbars (numeric matrix of gradients).
f2_f3_non_opencl: List with components qf (negative log-likelihood) and grad (gradients) from the CPU kernel.
f2_f3_opencl: List with components qf and grad from the OpenCL kernel.
run_opencl_pilot: Numeric scalar giving estimated runtime (seconds) for OpenCL evaluation on a pilot subset of the grid.

References

\insertAllCited

Examples

############################### Start of EnvelopeEval example ####################

# This example demonstrates EnvelopeEval in isolation. EnvelopeEval evaluates
# the negative log-likelihood and gradients at a grid of parameter values.
# It is called internally by EnvelopeBuild. Here we build the same inputs
# (grid G4, standardized model) using EnvelopeSize and expand.grid, then
# call EnvelopeEval directly. The setup mirrors Ex_EnvelopeBuild through
# the standardization step.

data(menarche, package = "MASS")
Age2 <- menarche$Age - 13

x <- matrix(as.numeric(1.0), nrow = length(Age2), ncol = 2)
x[, 2] <- Age2

y <- menarche$Menarche / menarche$Total
wt <- menarche$Total

mu <- matrix(as.numeric(0.0), nrow = 2, ncol = 1)
mu[2, 1] <- (log(0.9 / 0.1) - log(0.5 / 0.5)) / 3

V1 <- 1 * diag(as.numeric(2.0))
V1[1, 1] <- ((log(0.9 / 0.1) - log(0.5 / 0.5)) / 2)^2
V1[2, 2] <- (3 * mu[2, 1] / 2)^2

famfunc <- glmbfamfunc(binomial(logit))
f2 <- famfunc$f2
f3 <- famfunc$f3

dispersion2 <- as.numeric(1.0)
start <- mu
offset2 <- rep(as.numeric(0.0), length(y))
P <- solve(V1)
n <- 1000

wt2 <- wt / dispersion2
alpha <- x %*% as.vector(mu) + offset2
mu2 <- 0 * as.vector(mu)
P2 <- P
x2 <- x

parin <- start - mu
opt_out <- optim(parin, f2, f3,
  y = as.vector(y), x = as.matrix(x), mu = as.vector(mu2),
  P = as.matrix(P), alpha = as.vector(alpha), wt = as.vector(wt2),
  method = "BFGS", hessian = TRUE
)

bstar <- opt_out$par
A1 <- opt_out$hessian

Standard_Mod <- glmb_Standardize_Model(
  y = as.vector(y), x = as.matrix(x), P = as.matrix(P),
  bstar = as.matrix(bstar, ncol = 1), A1 = as.matrix(A1)
)

bstar2 <- Standard_Mod$bstar2
A <- Standard_Mod$A
x2 <- Standard_Mod$x2
mu2 <- Standard_Mod$mu2
P2 <- Standard_Mod$P2

###############################################################################
# Build grid G4 via EnvelopeSize and expand.grid (as EnvelopeBuild does)
###############################################################################
a <- diag(A)
omega <- (sqrt(2) - exp(-1.20491 - 0.7321 * sqrt(0.5 + a))) / sqrt(1 + a)
b2 <- as.vector(bstar2)
G1 <- rbind(b2 - omega, b2, b2 + omega)

size_info <- EnvelopeSize(a, G1, Gridtype = 3L, n = n)
G2 <- size_info$G2

G3 <- as.matrix(do.call(expand.grid, G2))
G4 <- t(G3)

###############################################################################
# EnvelopeEval: negative log-likelihood and gradients at grid points
###############################################################################
eval_out <- EnvelopeEval(
  G4 = G4,
  y = y,
  x = as.matrix(x2),
  mu = as.matrix(mu2, ncol = 1),
  P = as.matrix(P2),
  alpha = as.vector(alpha),
  wt = as.vector(wt2),
  family = "binomial",
  link = "logit",
  use_opencl = FALSE,
  verbose = FALSE
)

eval_out$NegLL
eval_out$cbars

###############################################################################
# End of EnvelopeEval example
###############################################################################

glmbayes documentation built on Aug. 5, 2026, 1:07 a.m.