Iterative Likelihood-Based Procedure

• The iterative likelihood-based procedure is a computational method that iteratively refines Gaussian mixture approximations to target densities using likelihood-driven steps.
• It employs local Laplace approximations and residual corrections to effectively capture multimodal, skewed, or heavy-tailed distributions.
• By integrating quasi-Monte Carlo grid evaluation and quadratic programming for weight updates, the approach achieves computational efficiency in high-dimensional settings.

An iterative likelihood-based procedure is a computational technique that constructs, updates, or corrects an approximation to a target probability density (or related quantity) through a sequence of likelihood-driven steps. These procedures are used primarily in Bayesian computation, statistical inference, and applied probability for efficient approximation of multivariate integrals, construction of proposal distributions for Monte Carlo simulation, and adaptive density estimation. A representative example of this paradigm is the iterated Laplace approximation (“iterLap”) for approximating complex multivariate probability densities (1103.3508).

1. Foundations: The Classical Laplace Approximation

The Laplace approximation provides a second-order Taylor expansion to approximate integrals of the form $\int \pi(x) dx$ (with $\pi(x)$ a nonnegative, integrable function, often unnormalized posterior or likelihood). The approximation centers at the mode $\tilde{x}$ of $\log \pi(x)$: $$\log \pi(x) \approx \log \pi(\tilde{x}) + \tfrac{1}{2} (x-\tilde{x})' H(\tilde{x}) (x-\tilde{x}),$$ where $H(\tilde{x})$ is the Hessian (matrix of second derivatives) of $\log \pi(x)$ at $\tilde{x}$. Exponentiating and normalizing, one obtains a multivariate normal density $\varphi(x;\tilde{x},\tilde{\Sigma})$ with mean $\tilde{x}$ and covariance $\tilde{\Sigma} = -H(\tilde{x})^{-1}$, scaled by a constant

$$w = (2\pi)^{p/2} |\tilde{\Sigma}|^{1/2} \pi(\tilde{x}).$$

This basic construction is accurate when $\pi(x)$ is unimodal, symmetric, and locally approximately Gaussian. However, for multimodal, skewed, or heavy-tailed densities, a single Laplace approximation is insufficient.

2. The Iterated Laplace Approximation: Algorithm and Workflow

The iterLap procedure iteratively refines the approximation to $\pi(x)$ by patching the residuals left by the current Gaussian mixture approximation. Each iteration comprises:

1. Initialization (Iteration 0): Detect all modes of $\pi(x)$, compute the Laplace approximation at each mode $\mu_j$ with covariance $\Sigma_j = -H(\mu_j)^{-1}$, form initial mixture

$$\tilde{\pi}_0(x) = \sum_{j=1}^{J^{(0)}} w_j \varphi(x; \mu_j,\Sigma_j),$$

where $w_j$ is as above.

1. Iterative Residual Correction:
   • (a) Residual Construction: Compute $r(x) = \pi(x) - \tilde{\pi}_{t-1}(x)$.
   • (b) Stabilization: Since $r(x)$ is not guaranteed to be positive, stabilize by adding a small constant if needed for logarithmic operations.
   • (c) Mode Localization: Maximize $\log r(x)$, apply Laplace approximation to $r(x)$ at the new mode to define a new Gaussian component $(\mu_{\text{new}},\Sigma_{\text{new}})$.
   • (d) Grid Evaluation: Build a quasi-random grid (e.g., Sobol sequence) in the region encompassing the new component, evaluate $\pi(x)$ and all mixture components on this grid.
   • (e) Weight Update: Solve the quadratic programming problem

   $$\min_w\, (y-F'w)'(y-F'w) \quad \text{subject to} \quad w_j \ge 0\ \forall j,$$

   with $y$ the vector of true density evaluations, $F$ the matrix of component densities. - (f) Update Mixture: Form new

   $$\tilde{\pi}_t(x) = \sum_{j=1}^{J^{(t)}} w_j \varphi(x;\mu_j,\Sigma_j).$$

• (g) Stopping Criteria: End iteration when $\max|y-\tilde{\pi}(y)| < \delta M_t$ for grid maximum $M_t$, or when the normalization constant stabilizes.

This cycle repeats, each time patching regions where the current approximation is deficient and adaptively increasing the complexity and fidelity of the mixture estimate.

3. Theoretical and Algorithmic Advancements over Classical Approaches

The iterative likelihood-based procedure improves upon the standard Laplace approximation in several ways:

• Local Adaptivity: By focusing on the residual $r(x)$, the algorithm automatically identifies regions (e.g., non-Gaussian structure, skewness, multimodality) where the current mixture is insufficient, guiding new components to those domains.
• Structural Flexibility: The mixture-of-Gaussians form

$$\tilde{\pi}(x) = \sum_{j=1}^J w_j \varphi(x;\mu_j,\Sigma_j)$$

can encode skewed, heavy-tailed, and multimodal densities—phenomena that a single Gaussian component cannot represent.

• Computational Efficiency: By employing local Laplace approximations, quasi-Monte Carlo grid evaluation, and least-squares weight updates, the iterative method avoids expensive global optimization or Markov chain Monte Carlo (MCMC) sampling, yielding major computational savings for high-dimensional density approximation tasks.

4. Multivariate Normal Mixtures and Quadratic Programming for Weight Determination

Each mixture component arises from a local, second-order expansion of $\log \pi(x)$ (or $\log r(x)$) at a detected mode, resulting in a multivariate normal. The weight vector $\mathbf{w}$ is determined by constrained least-squares fitting, formally: $$\min_{\mathbf{w}\geq 0}\ \|\mathbf{y} - F'\mathbf{w}\|_2^2$$ where $\mathbf{y}$ is the vector of true density values at grid locations, and $F$ is the matrix of mixture component densities evaluated at these points. This constrained quadratic programming ensures that mixture weights are non-negative and that the global fit minimizes mean squared error between the true and approximated densities on the grid.

This approach resembles nonnegative least squares (NNLS), combining computational tractability with stability and interpretability, and the normalization constant $Z_t = \sum_j w_j$ provides a natural stopping metric.

5. Implementational Details and The iterLap R-package

The iterLap method is implemented in the R-package iterLap, providing:

• Automated mode detection,
• Grid generation via quasi-random sequences (e.g., Sobol),
• Adaptive grid sizing proportional to $50\cdot p^{1.25}$ for dimension $p$,
• Iterative mixture expansion with weight updates through quadratic programming,
• Stopping rules based on maximum grid error and normalization constant stabilization.

In practice, iterLap enables the flexible and efficient construction of proposal distributions for Monte Carlo integration schemes, and is particularly advantageous for:

• Densities with multiple modes or pronounced skewness,
• High-dimensional approximation where analytic integration is intractable,
• Construction of proposals for importance sampling or Metropolis-Hastings algorithms.

6. Empirical Results and Performance Metrics

The effectiveness of the iterative likelihood-based procedure was demonstrated empirically:

• For a univariate skew mixture, initial Laplace fits underestimated tails; adding a residual-based component significantly improved accuracy.
• Bivariate skew and multimodal densities (e.g., skew t, mixture of normals) were fit accurately with a small number of components.
• In a 10-dimensional banana-shaped density, iterLap captured non-linear structure with approximately 19,000 target evaluations, performing comparably or better than adaptive importance sampling via MCMC.

Performance is measured by the terminal maximum grid absolute error and the number of mixture components needed for satisfactory convergence, with resource requirements determined primarily by the dimensionality $p$ and the number of components/grids iterations required for stabilization.

7. Limitations, Trade-Offs, and Practical Deployment

The iterative Laplace-based approach is highly effective in contexts where the target density is smooth and can be locally approximated by Gaussian components. Limitations arise in cases where the density is discontinuous, highly non-smooth, or the Hessian is ill-conditioned near modes (potentially causing negative-definiteness violations).

Trade-offs include:

• Number of Components vs. Approximation Quality: More mixture components increase the approximation's fidelity at the cost of computational overhead.
• Grid Size and Placement: Too coarse a grid impairs accuracy; too fine a grid incurs prohibitive computational cost.
• Optimization for Mode Detection: Local maxima of $r(x)$ must be computationally tractable to find (global maximization is generally hard in high dimensions).

For practical deployment, applications should tailor grid density, stopping thresholds, and the maximum number of components to available computational resources and required precision. The iterLap methodology is particularly suited for Bayesian model evidence evaluation, proposal distribution construction for Monte Carlo methods, and real-time density approximation in high-throughput settings.