Implementation Details: glmb_Standardize_Model()

The function glmb_Standardize_Model() implements this transformation in practice. It performs two key eigenvalue decompositions:

1. The first decomposition diagonalizes the Hessian of the log-posterior at the posterior mode. This yields a transformation matrix L2 that is used to normalize the posterior precision. The transformed design matrix and prior precision are updated accordingly.

2. The second decomposition identifies a diagonal matrix \(\varepsilon\) such that the difference \(P - \varepsilon\) remains positive definite. This ensures that the transformed model has a diagonal prior precision and can be reparameterized into standard form. A scaling loop is used to find the largest possible \(\varepsilon\) that satisfies this condition.

The function then applies a second transformation L3 to finalize the standardization. This results in a model where both the prior and the data precision matrices are diagonal and that the prior precision matrix is the identity matrix, satisfying the theoretical conditions for standard form.

The output includes:

• Transformed posterior mode (bstar2)
  
• Transformed design matrix (x2)
  
• Prior contribution to the modified log-likelihood (P2): the part of the prior precision that is shifted to the log-likelihood
  
• Transformed prior mean (mu2, typically zero)
  
• Precision matrix for standardized data (A)
  
• Inverse transformation matrices (L2Inv, L3Inv) for undoing the standardization

This implementation ensures that the model satisfies the theoretical conditions for standard form, enabling efficient envelope construction and posterior sampling. It also provides a modular and transparent framework for debugging and extending the standardization process in future versions of the package.

1. Gridtype 1: Static Threshold Test

• If \(\sqrt{1 + a_i} \leq \frac{2}{\sqrt{\pi}} \approx 1.128379\), then the expected number of candidate draws per accepted sample under a single tangent envelope is less than the limiting upper bound provided by the three-tangent envelope.
  Therefore, building more points does not meaningfully reduce rejections.

  • Use a single-point envelope at the mode:
    \(G2[i] = \{\theta_i^{\ast}\}\), \(GIndex1[i] = \{4\}\)

• Else: build a three-point envelope at
  \(\{\theta_i^{\ast} - \omega_i,\ \theta_i^{\ast},\ \theta_i^{\ast} + \omega_i\}\),
  \(GIndex1[i] = \{1,2,3\}\)

2. Gridtype 2: Dynamic Envelope via EnvelopeOpt(a, n)

• Rather than a fixed-bound test, this approach solves (per dimension) the cost tradeoff between:

  • \(T_{\text{build}}(g_i) \propto g_i\) (linear grid-construction cost)
    
  • \(T_{\text{sample}}(n, \text{acc}_i(g_i)) \propto n / \text{acc}_i(g_i)\), where \(\text{acc}_i(g_i)\) is the approximate acceptance rate given \(g_i\) grid points

• EnvelopeOpt returns \(\text{gridindex}[i] \in \{1, 3\}\) by minimizing \(T_{\text{build}} + T_{\text{sample}}\) under approximations from subgradient theory.

• This allows more envelope points up front when \(n\) is large, because reduced rejections during sampling more than pay for the extra build cost.

3. Gridtype 3: Always Three-Point Grid

We now consider the sampling performance using the foregoing approach when the data are multivariate normal with a diagonal precision matrix given by NP for a positive definite diagonal matrix P and positive scalar N.

Theorem 3: Consider a k-dimensional multivariate normal data model in standard form with precision as above, and let \(\tilde{a}^*(N)\) denote the value of \(\tilde{a}\) at N; then
\[ \lim_{N \to \infty} \tilde{a}^*(N) = \left( \frac{2}{\sqrt{\pi}} \right)^k \]

Proof: This follows immediately from Theorem 2 in Nygren & Nygren (2006).

4. Gridtype 4: Always Single-Point Grid (Mode Only)

Following Example 1 in Nygren & Nygren (2006), the expected number of candidate draws per accepted draw in the multivariate normal case is:
\[ \prod_{i=1}^{k} \sqrt{1 + a_i} \]

where k is the number of dimensions. More generally, if the data are “normal-like,” then we expect similar acceptance rates even for other generalized linear models.

Deep Dive: Gridtype 1 — Static Threshold Envelope

Gridtype 1 uses a simple threshold test to decide whether a single-point or three-point envelope should be constructed for each dimension. The key quantity is the posterior precision \(1 + a_i\), and the threshold \(\frac{2}{\sqrt{\pi}} \approx 1.128379\) arises from bounding the expected number of candidate draws per accepted sample.

If \(\sqrt{1 + a_i} \leq \frac{2}{\sqrt{\pi}}\), then the expected number of candidate draws per accepted sample under a single tangent envelope is less than the limiting upper bound provided by the three-tangent envelope. This implies that adding more grid points does not meaningfully reduce rejection rates, and a single tangent at the mode suffices.

This logic is especially useful for low-variance dimensions, where the posterior is sharply peaked and additional tangents offer little marginal benefit. The resulting envelope is minimal and fast to construct:

• Grid: \(\{\theta_i^{\ast}\}\)
  
• Index: \(\{4\}\) (single tangent at the mode)

If the threshold is exceeded, a symmetric three-point grid is used:

• Grid: \(\{\theta_i^{\ast} - \omega_i,\ \theta_i^{\ast},\ \theta_i^{\ast} + \omega_i\}\)
  
• Index: \(\{1, 2, 3\}\) (left, center, right tangents)

The width parameter \(\omega_i\) is derived from the local curvature of the log-likelihood at the posterior mode. Specifically:

\[ \omega_i := \frac{ \sqrt{2} - \exp\left( -1.20491 - 0.7321 \sqrt{0.5 + a_i} \right) }{ \sqrt{1 + a_i} } \]

Since \(\sqrt{1 + a_i}\) is the square root of the precision for dimension i, the inverse is essentially the standard deviation. When \(a_i\) is large, the second term approaches zero, so:

\[ \omega_i \approx \frac{ \sqrt{2} }{ \sqrt{1 + a_i} } \]

and hence is essentially proportional to the standard deviation.

Deep Dive: Gridtype 2 — Adaptive Envelope via EnvelopeOpt()

Gridtype 2 uses a dynamic strategy to select the number of grid points per dimension based on a cost-benefit analysis. The goal is to minimize total runtime by balancing envelope construction cost against sampling efficiency.

The tradeoff is formalized as:

• \(T_{\text{build}}(g_i) \propto g_i\): cost of building \(g_i\) tangents in dimension \(i\)

• \(T_{\text{sample}}(n, \text{acc}_i(g_i)) \propto n / \text{acc}_i(g_i)\): cost of drawing \(n\) samples given acceptance rate \(\text{acc}_i(g_i)\)

The function EnvelopeOpt(a_i, n) solves this tradeoff using approximations from subgradient theory. Dimensions are ranked by decreasing posterior standard deviation — equivalently, by increasing values of \(1 / (1 + a_i)\). Dimensions with larger a_i (i.e., sharper curvature) are more likely to benefit from additional tangents.

The algorithm proceeds as follows:

1. Rank dimensions by \(1 / (1 + a_i)\)
2. For each possible number of three-point dimensions (from 0 to l1):
   • Assign three-point grids to the top-ranked dimensions
   • Estimate total build cost and sampling cost
   • Compute total evaluation cost: evalest = T_build + T_sample
3. Select the configuration with minimal evalest

This results in a vector dimcount where each entry is either 1 (single-point) or 3 (three-point), depending on the optimal tradeoff.

This adaptive logic allows the simulation routine to scale efficiently across models of varying size and complexity, and reflects the theoretical insights from Nygren & Nygren (2006) on envelope efficiency.

In the implementation of EnvelopeOpt(), the line

intest[i] <- 3^(i-1)

represents the grid construction cost at iteration i. This term estimates the number of log‑likelihood and gradient evaluations required to build the envelope as more dimensions are promoted from single‑point to three‑point tangents:

• At the start (\(i = 1\)), all dimensions use single‑point tangents, so the build cost is minimal.
  
• Each time a new dimension is switched to a three‑point tangent, the number of required evaluations grows multiplicatively.
  
• The factor \(3^{\,i-1}\) captures this exponential growth: with each additional three‑point dimension, the grid size triples, and so does the number of evaluations needed to construct the envelope.

In addition to the build cost term intest[i], the function also tracks the expected sampling cost through the quantity

slopeest[i] <- prod(scaleest[i,])

This term represents the expected number of likelihood/gradient (slope) evaluations required to obtain a valid draw under the current grid configuration:

• At initialization, scaleest[1,] <- sqrt(1 + a1), so
  \[ \text{slopeest}[1] \;=\; \prod_j \sqrt{1 + a_j}, \]
  which corresponds to the baseline rejection burden when all dimensions use single‑point tangents.

• Each time a dimension is promoted to a three‑point tangent, its entry in scaleest is replaced by the constant
  \[ \tfrac{2}{\sqrt{\pi}}, \]
  reflecting the improved acceptance bound for that dimension.

• Consequently, as more dimensions are switched to three‑point tangents, the product
  \[ \prod_j \text{scaleest}[i,j] \]
  decreases, lowering the expected number of slope evaluations per accepted draw.

Thus, slopeest[i] captures the sampling efficiency of the envelope:

• Large slopeest values → more rejections, higher per‑draw cost.
  
• Smaller slopeest values → fewer rejections, more efficient sampling.

Together with the build cost term intest[i], the total evaluation cost is

\[ \text{evalest}[i] \;=\; \underbrace{3^{\,i-1}}_{\text{build cost}} \;+\; \underbrace{n \cdot \text{slopeest}[i]}_{\text{sampling cost}}. \]

This balance ensures that the algorithm adaptively chooses the grid size that minimizes the combined cost of envelope construction and sampling.

• Build cost grows rapidly as more dimensions are promoted.
  
• Sampling cost decreases because acceptance rates improve with additional tangents.

The optimizer selects the configuration that minimizes this total. Intuitively:

• For small sample sizes (\(n\) small), the build cost dominates, so the algorithm favors fewer three‑point tangents.
  
• For large sample sizes (\(n\) large), the sampling cost dominates, so the algorithm invests more in upfront grid construction to reduce rejections later.

When GPU acceleration is available through OpenCL, the Grid construction can be performed in parallel, reducing the amount of time that components takes. To that end, we modify the call to the function to reflect the lower Grid contruction cost by passing a scaled up sample size (which yields the same optimization result).

  int core_CNT=get_opencl_core_count();
  
  if (verbose) {
    Rcpp::Rcout << "[INFO] OpenCL core count = " << core_CNT << "\n";
  }
  
  // Second row in G1b here is the posterior mode
  
  
  
  NumericVector gridindex(l1);
  
  if(Gridtype==2){
    // Temporarily scale n to encourage richer envelopes when GPU parallelism is available
    int scaled_n = std::max(1, n * core_CNT);
    
    if (verbose) {
      Rcpp::Rcout << "[INFO] Scaling n from " << n << " to " << scaled_n
                  << " for envelope optimization.\n";
    }
    
    gridindex = EnvelopeOpt(a_2, scaled_n);
    
    //gridindex=EnvelopeOpt(a_2,n);
  }

1. Compute width parameters \(\omega_i\) from the diagonal precision matrix

Let \(\theta^{\ast}\) denote the unique posterior mode. For each dimension \(i\), define:

\[ \omega_{i} := \frac{\sqrt{2} - \exp\left(-1.20491 - 0.7321 \sqrt{0.5 - \frac{\partial^{2} \log f(\theta^{\ast} \mid y)}{\partial \theta_{i}^{2}}} \right)}{\sqrt{1 - \frac{\partial^{2} \log f(\theta^{\ast} \mid y)}{\partial \theta_{i}^{2}}}} \]

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 yields an envelope whose efficiency does not deteriorate with sample size.

2. Construct intervals around the posterior mode \(\theta^\star\)

Set:

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

Then define 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) \]

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

3. 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\} \]

Tangency points \(\bar{\theta}_{j}\) are 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} \]

These tangency points ensure the envelope touches the log-likelihood at representative points in each interval, guaranteeing dominance and tightness.

4. Build the full grid of tangency points

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

5. Evaluate negative log-likelihood and gradients at each grid point

This step constructs the likelihood subgradient densities and facilitates accept–reject sampling. The subgradients \(c(\bar{\theta})\) are part of the definitions of likelihood subgradient densities (see Definition 2). Both subgradients and negative log-likelihoods (via \(h_{\bar{\theta}}(.)\)) are used during the accept–reject procedure.

• On CPU: via f2_* and f3_* routines
  
• On GPU: via f2_f3_opencl, which computes both together

6. Call Set_Grid_C2_pointwise to evaluate restricted multivariate normal log-densities

(Claim 2; Remark 5)
Each restricted density corresponds to a subset of the partition, normalized as in Remark 5.

7. Call setlogP_C2 to compute component log-probabilities and constants

(Remark 6)
The constants \(\tilde{a}\) and mixture weights \(\tilde{p}_i\) are computed explicitly as in Remark 6, ensuring proper normalization of the mixture envelope.

8. Normalize probabilities (PLSD) and optionally sort grid components

(Claim 2)
Normalization ensures the mixture of restricted densities forms a valid dominating density for the posterior. Sorting improves sampling efficiency.

8.1 Reverse-transformation pipeline

• Undo L3 transform: left-multiply the matrix of standardized draws by L3Inv.
  
• Undo L2 transform: left-multiply the result by L2Inv.
  
• Restore prior mean: add the original prior mean vector mu to every draw.
  
• Optional diagnostics: recompute the user-supplied log-posterior at the transformed draws when requested.

8.2 Returned object (fields you can rely on)

• coefficients: n × k matrix of posterior draws on the original parameter scale (returned as trans(out2) in the C++ routing function).
  
• coef.mode: k-vector containing the posterior mode on the original scale (optimizer output + mu).
  
• dispersion: final dispersion used (1 for Poisson/Binomial; user value otherwise).
  
• Prior: list with components mean and Precision (original mu and P).
  
• offset: offset vector passed through.
  
• prior.weights: original weights (before dispersion scaling).
  
• y, x: data echoed for convenience.
  
• fit: the full optim() return list used to obtain the posterior mode and Hessian.
  
• iters: integer vector of attempts per accepted draw (sampler output).
  
• Envelope: envelope artifact returned by EnvelopeBuild_c (mixture weights, tangency points, LLconst, cbars, loglt, logrt, GridIndex, etc.).