EnvelopeBuild: GPU-Accelerated Envelope Construction for Posterior...

In glmbayes: Bayesian Generalized Linear Models (IID Samples)

GPU-Accelerated Envelope Construction for Posterior Simulation

Description

GPU-Accelerated Envelope Construction for Posterior Simulation

Usage

EnvelopeBuild(bStar,A,y,x,mu,P,alpha,wt,family = "binomial",link = "logit", 
Gridtype = 2L,n = 1L,n_envopt=NULL,sortgrid = FALSE,use_opencl = FALSE,verbose = FALSE)

EnvelopeSetGrid(GridIndex, cbars, Lint)

EnvelopeSetLogP(logP, NegLL, cbars, G3)

Arguments

bStar	Point at which envelope should be centered (typically posterior mode).
A	Diagonal precision matrix for the log-likelihood in standard form.
y	A vector of observations of length m.
x	A design matrix of dimension m * p.
mu	A vector giving the prior means of the variables.
P	Prior precision matrix of the variables (positive-definite).
alpha	Offset vector.
wt	A vector of weights.
family	Family for the envelope: binomial, quasibinomial, poisson, quasipoisson, or Gamma.
link	Link function ("logit", "probit", "cloglog" for binomial; "log" for Poisson/Gamma).
Gridtype	Method to determine the number of subgradient densities in the grid.
n	Number of draws from the posterior (used for grid sizing).
n_envopt	Effective sample size passed to EnvelopeOpt for grid construction. Defaults to match n. Larger values encourage tighter envelopes.
sortgrid	Logical; if TRUE, sort the envelope descending by component probability.
use_opencl	Logical; if TRUE, use OpenCL for gradient evaluations.
verbose	Logical; if TRUE, print progress messages.
GridIndex	A matrix indicating, for each grid component, whether the component lies in the left tail, center, or right tail of the density. Rows correspond to grid components; columns correspond to standardized variables.
cbars	A matrix containing the subgradient of the (adjusted) negative log-likelihood at each grid component.
Lint	A matrix storing the lower and upper bounds for each grid component, depending on whether sampling is from the left, center, or right.
logP	A matrix (typically two columns) with information for each grid component. The first column usually holds the output from EnvelopeSetGrid(), corresponding to the restricted normal density.
NegLL	A vector of negative log-likelihood evaluations at each grid component.
G3	A matrix of tangency points used in the grid.

Details

Constructs an enveloping function for posterior simulation using a grid of tangency points. The envelope is used in accept-reject sampling to guarantee iid draws from the posterior distribution. The implementation follows \insertCiteNygren2006glmbayes, with extensions for GPU acceleration (via OpenCL), dynamic grid optimization, and parallelized evaluation.

The envelope is typically built around the posterior mode \theta^\star for a model in standard form (which in this context means a model with a diagonal posterior precision matrix and prior identity precision matrix - glmb_Standardize_Model). It uses dimension-specific width parameters \omega_i derived from the precision matrix. Tangency points are selected per dimension, and the full grid is formed via Cartesian expansion. Negative log-likelihood and gradient values are computed at each grid point, either on CPU or GPU depending on the use_opencl flag. These values are used to construct a piecewise envelope function that dominates the posterior density.

Value

EnvelopeBuild()

A list of envelope components used for accept-reject sampling:

GridIndex

Integer matrix encoding sampling type (tail, center, line) per dimension and region.

thetabars

Matrix of tangency points \bar{\theta}_j for each grid region.

cbars

Matrix of subgradients c(\bar{\theta}_j) of the negative log-likelihood at tangency.

loglt

Matrix of log left-tail probabilities per dimension and region.

logrt

Matrix of log right-tail probabilities per dimension and region.

logU

Matrix of selected per-dimension log-density contributions (tail/center) for each region.

logP

Matrix of total log-probabilities per region (first column); used to derive mixture weights.

PLSD

Vector of normalized mixture weights over grid regions used to draw region indices.

LLconst

Vector of acceptance-test constants per region used in the inequality for rejection sampling.

EnvelopeSetGrid()

A list of matrices computed for grid-based log-density evaluation:

Down

Lower bounds for truncated-normal evaluation per dimension and region.

Up

Upper bounds for truncated-normal evaluation per dimension and region.

lglt

Log left-tail probabilities (from (-\infty, \mathrm{Up}]) per dimension and region.

lgrt

Log right-tail probabilities (from [\mathrm{Down}, \infty)) per dimension and region.

lgct

Log central-interval probabilities (from [\mathrm{Down}, \mathrm{Up}]) per dimension and region.

logU

Selected log-probability per grid cell based on GridIndex (tail or center).

logP

Matrix with row-wise sums of logU (first column) used to form mixture weights.

EnvelopeSetLogP()

A list with updated mixture-weight and acceptance constants:

logP

Input logP with its second column populated by the log of unnormalized visit probabilities per region (mixture denominators).

LLconst

Vector of acceptance constants -\log f(y \mid \bar{\theta}_j) - c(\bar{\theta}_j)^{T}\bar{\theta}_j used in the accept-reject test.

Models in standard form

The standard-form restriction and its closed-form truncated-normal integrals follow \insertCiteNygren2006glmbayes. See \insertCiteglmbayesChapterA08glmbayes for the full theoretical details (standard form restriction, closed-form truncated-normal integrals, and the resulting log-scale tractability).

In the implementation, these standard-form quantities determine the grid-based tangency shifts and the precomputed region constants (e.g., the log-CDF pieces) that drive the mixture weights used by the envelope sampler.

Construction of restricted subgradient densities

For the full restricted density construction and the resulting envelope constants, see \insertCiteglmbayesChapterA08glmbayes.

In the implementation, these theory objects become the precomputed region log-constants (via closed-form CDF pieces) that are used to normalize the envelope mixture weights for the accept-reject sampler.

Mixture construction and tractable probabilities

The mixture construction and its tractable region probabilities are derived in \insertCiteglmbayesChapterA08glmbayes. In the implementation, these theory objects become the precomputed region log-constants and the mixture weights (PLSD) used by the envelope-based accept-reject sampler.

Log-scale properties of the envelope function

The log-scale form of the envelope factor and the subgradient inequality that imply envelope dominance are given in \insertCiteglmbayesChapterA08glmbayes. In the implementation, these properties allow pointwise evaluation in the log-domain and provide the theoretical basis for the rejection test inside the sampler.

Use of the envelope during sampling

The standardized sampler called through .rNormalGLM_std_cpp() uses the envelope to generate posterior samples via rejection sampling. Although not exported, this routine is called internally by .rNormalGLM_cpp(), which in turn is invoked by the user-facing function rNormal_reg(). Together, these routines implement envelope-based sampling for generalized linear models with log-concave likelihood functions and multivariate normal priors.

The envelope provides a mixture of restricted likelihood-subgradient densities, each defined over a region A_i, with associated mixture weights \tilde{p}_i stored in PLSD. The sampling proceeds as follows:

1. A region index J(i) is drawn from the discrete distribution defined by PLSD.

2. A candidate \theta_i is drawn from the restricted density q^{\bar{\theta}_{J(i)}}_{A_{J(i)}}, using the normal CDF bounds loglt and logrt, and subgradient vector cbars. Simulation for each dimension uses the internal C++ function ctrnorm_cpp(), which explicitly uses these inputs.

3. The log-likelihood \log f(y \mid \theta_i) is computed and stored in testll[0] using the appropriate likelihood function f2.

The acceptance test is performed using the inequality

\log(U_2) \le \mathrm{LLconst}[J(i)] + \mathrm{cbars}[J(i), ]^{T} \theta_i + \log f(y \mid \theta_i),

which is equivalent to

\log(U_2) \le \log f(y \mid \theta_i) - \left( \log f(y \mid \bar{\theta}_{J(i)}) - c(\bar{\theta}_{J(i)})^{T}(\theta_i - \bar{\theta}_{J(i)}) \right),

where:

• LLconst[J(i)] stores the precomputed quantity -\log f(y \mid \bar{\theta}_{J(i)}) - c(\bar{\theta}_{J(i)})^{T} \bar{\theta}_{J(i)}, computed during envelope construction via EnvelopeSet_LogP_C2().

• cbars[J(i), ] is the precomputed subgradient vector c(\bar{\theta}_{J(i)}), extracted via cbars(J(i), _). It defines the exponential tilt direction used to evaluate the envelope.

• testll[0] is the log-likelihood at the candidate draw \theta_i, evaluated using the model specified by family and link.

• -\log(U_2) is the threshold from a uniform draw U_2 \sim \mathrm{Unif}(0,1).

The right-hand side of this inequality is always non-positive, and equals zero when \theta_i = \bar{\theta}_{J(i)}. This reflects the fact that the envelope is tangent to the log-likelihood at each \bar{\theta}_j, and lies above it elsewhere.

This procedure guarantees that accepted samples are drawn from the posterior \pi(\theta \mid y). The envelope ensures bounded rejection probability, and the mixture structure allows efficient sampling across regions. The output out contains accepted draws, and draws records the number of attempts per sample.

The components returned by EnvelopeBuild() are used in specific steps of the sampling procedure as follows:

• PLSD is used to randomly select a region index J(i) from the envelope mixture.

• loglt and logrt define the truncated normal bounds for each dimension, used together with cbars to generate candidate values \theta_i.

• cbars provides the subgradient vectors c(\bar{\theta}_j) used both for candidate generation and for computing the acceptance test.

• LLconst stores precomputed constants used in the acceptance inequality, avoiding recomputation of posterior terms at tangency points.

• logU stores the per-dimension log-density contributions for each region, computed during envelope setup. These values are summed to produce logP, which determines the mixture weights PLSD.

• logP contains the total log-probabilities for each grid component, which are normalized to form the mixture weights PLSD.

• thetabars stores the tangency points \bar{\theta}_j used to define subgradients and region-specific densities.

• GridIndex encodes the sampling type (tail, center, line) used for each dimension and region, guiding how each coordinate is simulated.

Algorithmic steps (linked to theory)

The implementation of EnvelopeBuild follows the envelope construction in \insertCiteNygren2006glmbayes for models in standard form (see Section 3–3.3 there). Each computational step corresponds to a theoretical guarantee:

1. Compute width parameters \omega_i from the diagonal precision matrix. In particular, let \theta^{\ast} denote the unique posterior mode. For each dimension i, define

   \omega_{i} := \frac{\sqrt{2} - \exp\!\big(-1.20491 - 0.7321\,\sqrt{0.5 - \partial^{2}\log f(\theta^{\ast}\mid y)/\partial\theta_{i}^{2}}\big)} {\sqrt{1 - \partial^{2}\log f(\theta^{\ast}\mid y)/\partial\theta_{i}^{2}}}.

   As seen from the above, the widths \omega_i are derived from the local curvature of the log-likelihood at the posterior mode. This ensures that the three-interval construction per dimension below yields an envelope whose efficiency does not deteriorate with sample size.

2. Use the width parameters to construct intervals around the posterior mode \theta^\star. Specifically, we set

   \ell_{i,1} = \theta^{\ast}_{i} - 0.5\,\omega_{i}, \quad \ell_{i,2} = \theta^{\ast}_{i} + 0.5\,\omega_{i},

   and construct three intervals per dimension:

   A_{i,1} = (-\infty,\ell_{i,1}), \quad A_{i,2} = [\ell_{i,1},\ell_{i,2}], \quad A_{i,3} = (\ell_{i,2},\infty).

For each dimension i, let J_{i} = \{1,2,3\} and define J = \prod_{i=1}^{p} J_{i}, which has 3^{p} elements. Each j \in J is a vector (j_{1},\ldots,j_{p}), and we define

A^{\ast}_{j} = \prod_{i=1}^{p} A_{i,j_{i}}.

The collection A^{\ast} = \{A^{\ast}_{j} : j \in J\} forms a partition of \Theta.

1. For each member of the partition, select tangency points \theta^\star \pm \omega_i.

For each j \in J, define index sets

C_{j1} = \{i : j_{i} = 1\}, \quad C_{j2} = \{i : j_{i} = 2\}, \quad C_{j3} = \{i : j_{i} = 3\}.

The tangency points \bar{\theta}_{j} are then defined componentwise by

\bar{\theta}_{j,i} = \begin{cases} \theta^{\ast}_{i} - \omega_{i}, & i \in C_{j1}, \\ \theta^{\ast}_{i}, & i \in C_{j2}, \\ \theta^{\ast}_{i} + \omega_{i}, & i \in C_{j3}. \end{cases}

The tangency points are hence chosen so that the envelope touches the log-likelihood at representative points in each interval, guaranteeing dominance and tightness.

1. Build the full grid of tangency points (Cartesian product across dimensions).

   The Cartesian product of per-dimension partitions yields the 3^p restricted densities described in the paper, ensuring coverage of the full parameter space.

2. Evaluate negative log-likelihood and gradients at each grid point to construct the likelihood subgradient densities and to facilitate accept rejection sampling

   The subgradients c(\bar{\theta}) enter the likelihood-subgradient density construction (\insertCiteNygren2006glmbayes; see also \insertCiteglmbayesChapterA08glmbayes), and both subgradients and negative log-likelihoods (through h_{\bar{\theta}}(\cdot)) are used in the accept-reject procedure. CPU and GPU routines compute these values efficiently.

   • On CPU: via f2_f3_non_opencl.

   • On GPU: via f2_f3_opencl, which computes these in parallel across faces

3. Call EnvelopeSet_Grid_C2_pointwise to evaluate restricted multivariate normal log-densities. Each restricted density corresponds to a subset of the partition, normalized as in Remark 5 of \insertCiteNygren2006glmbayes.

4. Call EnvelopeSet_LogP_C2 to compute component log-probabilities and constants. The constants \tilde{a} and mixture weights \tilde{p}_i are computed explicitly as in Remark 6 of the paper, ensuring that the mixture envelope is properly normalized.

5. Normalize probabilities (PLSD) and optionally sort grid components. Normalization implements Claim 2 of the paper so the mixture forms a valid dominating density for the posterior. Sorting is an implementation detail to improve sampling efficiency.

Theory reference (JASA paper and vignette)

Definitions, claims, theorems, remarks, and examples through Remark 16 (including standard form, the 3^p partition, and sampling remarks) are in \insertCiteNygren2006glmbayes. An expanded narrative is in vignette("Chapter-A08", package = "glmbayes").

Subgradient density formulation

Each grid component corresponds to a tilted multivariate normal density, normalized using the moment-generating function (MGF). In the single-point case, centered at the posterior mode \theta^\star, the density is:

f(\theta) = \frac{1}{(2\pi)^{p/2} |A|^{-1/2} \cdot \text{MGF}_A(c)} \exp\left( -\frac{1}{2} (\theta - \mu)^T A (\theta - \mu) + c^T (\theta - \theta^\star) \right)

where:

• A is the precision matrix,

• \mu is the prior mean vector,

• c is the gradient of the log-likelihood at \theta^\star,

• \text{MGF}_A(c) is the moment-generating function:

  \text{MGF}_A(c) = \exp\left( \frac{1}{2} c^T A^{-1} c \right)

This closed-form density dominates the posterior locally and is used when Gridtype = 1. For richer envelopes, multiple such components are constructed at tangency points \theta_j, each with its own gradient c_j, and combined into a mixture:

f_{\text{env}}(\theta) = \sum_{j=1}^{K} p_j f_j(\theta)

where the weights p_j are computed using log-CDF differences and constants:

\log p_j = \log \Phi(U_j) - \log \Phi(L_j) - \text{NegLL}_j + \text{LLconst}_j

Gridtype logic

The Gridtype argument controls how many tangency points are used per dimension:

• 1: Threshold rule. If 1 + a_i \le 2/\sqrt{\pi}, use a single-point envelope at the mode; otherwise use three points.

• 2: Dynamic optimization via EnvelopeOpt, which balances grid build cost and expected acceptance rate. Grid size is scaled by n and the number of OpenCL cores when GPU is enabled.

• 3: Always use three points per dimension.

• 4: Always use a single point (mode only).

Supported families and links

The following families and link functions are supported:

• Binomial: logit, probit, cloglog

• Quasibinomial: logit, probit

• Poisson: log

• Quasipoisson: log

• Gamma: log

• Gaussian: identity

GPU acceleration (use_opencl = TRUE) is available for all of the above except Gaussian, which is always evaluated on CPU.

GPU acceleration

When use_opencl = TRUE, likelihood and gradient evaluations are offloaded to the GPU using OpenCL. This can substantially reduce runtime for high-dimensional models or large grids. Results are mathematically equivalent to the CPU version, but small numerical differences may occur due to floating-point arithmetic. If reproducibility across hardware is critical, prefer the CPU path.

If OpenCL support was not detected at compile time, the flag is ignored and the CPU implementation is used. Diagnostic messages are printed when verbose = TRUE.

Verbose output

When verbose = TRUE, the function prints:

• Grid type, number of draws, OpenCL usage, and detected core count.

• Grid size after expansion.

• Time-stamped messages when entering the grid loop, starting likelihood evaluations, starting gradient evaluations, and invoking GPU kernels.

• Messages when setting grid values, computing log-probabilities, and sorting.

Any constants needed by the sampling are added to a list and returned.

References

\insertAllCited

See Also

EnvelopeSize, EnvelopeEval, EnvelopeSort, glmb_Standardize_Model; rNormal_reg, rglmb, glmb. Theory and vignettes: \insertCiteNygren2006glmbayes; \insertCiteglmbayesChapterA08,glmbayesSimmethods,glmbayesChapterA10,glmbayesChapter12glmbayes.

Examples

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

summary(menarche)
plot(Menarche/Total ~ Age, data=menarche)


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))

# 2 standard deviations for prior estimate at age 13 between 0.1 and 0.9
## Specifies uncertainty around the point estimates

V1[1,1]<-((log(0.9/0.1)-log(0.5/0.5))/2)^2 
V1[2,2]=(3*mu[2,1]/2)^2  # Allows slope to be up to 1 times as large as point estimate 

famfunc<-glmbfamfunc(binomial(logit))

f1<-famfunc$f1
f2<-famfunc$f2
f3<-famfunc$f3
f5<-famfunc$f5
f6<-famfunc$f6

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


###### Adjust weight for dispersion

wt2=wt/dispersion2

######################### Shift mean vector to offset so that adjusted model has 0 mean

alpha=x%*%as.vector(mu)+offset2
mu2=0*as.vector(mu)
P2=P
x2=x


#####  Optimization step to find posterior mode and associated Precision

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  ## Posterior mode for adjusted model
bstar
bstar+as.vector(mu)  # mode for actual model
A1=opt_out$hessian # Approximate Precision at mode

## Standardize Model

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
L2Inv=Standard_Mod$L2Inv
L3Inv=Standard_Mod$L3Inv

Env2=EnvelopeBuild(as.vector(bstar2), as.matrix(A),y, as.matrix(x2),
as.matrix(mu2,ncol=1),as.matrix(P2),as.vector(alpha),as.vector(wt2),
family="binomial",link="logit",Gridtype=as.integer(3), n=as.integer(n),
sortgrid=TRUE)

## These now seem to match

Env2

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