Bayesian Level Set Sampling Methods

• Bayesian level set sampling is a framework that infers interfaces and partitions by thresholding latent functions combined with Bayesian updating.
• It integrates Gaussian priors and surrogate models with likelihood-based MCMC and sequential design to robustly handle inverse and sampling problems.
• The approach has been successfully applied to reconstruct sharp boundaries in geometric inverse problems and to sample nested convex slices in high-dimensional settings.

Bayesian level set sampling denotes a family of Bayesian methodologies in which the object of interest is a set defined by thresholding a function. In geometric inverse problems, the unknown interface between regions is represented implicitly as the level set of a latent function, and Bayesian inference is performed on that latent function; in black-box level-set estimation, the target is a super-level set or a reliable level set of an expensive response; and in distributional sampling, nested level sets of a target density are used to construct samplers over convex slices (Iglesias et al., 2015, Huang et al., 2021, Iwazaki et al., 2019, Cheng et al., 2024, Foster et al., 2012). The resulting methods combine priors on latent fields or surrogates, likelihood-based posterior updating, and either Markov chain Monte Carlo or sequential evaluation policies to quantify uncertainty about interfaces, partitions, and threshold-defined regions.

1. Core formulations

A canonical geometric formulation introduces a real-valued level-set function $u:D\to\mathbb R$ or $\phi:B_R\to\mathbb R$, together with thresholds

$$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$

The physical domain is partitioned by the sublevel bands

$$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$

and a piecewise-constant coefficient is recovered through a level-set map

$$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$

This representation is used for discontinuous coefficients, refractive indices, and related geometric inverse problems (Iglesias et al., 2015, Huang et al., 2021).

A second formulation treats the level set itself as the primary inferential target for an expensive black-box function $f$. For a threshold $c\in\mathbb R$, the super-level set is

$L=\{x\in\mathcal X:f(x)\ge c\},$

while under input uncertainty one instead studies the reliable level set

$$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$

For scalar equality constraints, level sets may also be written as

$A_c=\{\theta\in\Theta:f(\theta)=c\},$

with multivariate generalization to $\phi:B_R\to\mathbb R$0 (Cheng et al., 2024, Iwazaki et al., 2019, Edwards et al., 2024).

A third formulation uses level sets of a density rather than a forward model. For a quasi-concave density $\phi:B_R\to\mathbb R$1, the upper level set

$\phi:B_R\to\mathbb R$2

is convex for every $\phi:B_R\to\mathbb R$3. The sampler operates across a nested sequence $\phi:B_R\to\mathbb R$4, estimates relative slice volumes, and reweights pooled samples to recover the target density (Foster et al., 2012).

Setting	Level-set object	Target of inference or sampling
Geometric inverse problems	$\phi:B_R\to\mathbb R$5 or $\phi:B_R\to\mathbb R$6	Interfaces and piecewise-constant fields
Black-box level-set estimation	$\phi:B_R\to\mathbb R$7 or a surrogate for $\phi:B_R\to\mathbb R$8	Super-level or reliable level sets
Smoothed ABC level sets	$\phi:B_R\to\mathbb R$9 through $$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$0	Posterior mass concentrated near $$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$1
Quasi-concave density sampling	$$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$2	Samples from the target density

This suggests that “Bayesian level set sampling” is not a single algorithmic template. Rather, it is a shared strategy in which threshold-induced geometry is made inferentially explicit.

2. Latent representations and prior models

In the inverse-problem literature, the standard prior is Gaussian on the level-set function. One formulation uses

$$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$3

with $$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$4 for $$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$5 under Neumann-zero-mean conditions, or an integral operator

$$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$6

chosen so that draws lie almost surely in $$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$7 (Iglesias et al., 2015). In acoustic inverse medium scattering, the prior is taken as

$$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$8

with Neumann boundary conditions, where $$-\infty=c_0<c_1<\cdots<c_n=\infty \quad\text{or}\quad -\infty=c_0<c_1<\cdots<c_L=\infty.$$9 controls smoothness, $$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$0 is an inverse length-scale parameter, and $$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$1 is a marginal variance (Huang et al., 2021). The same Whittle–Matérn prior may be generated by the SPDE

$$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$2

Several extensions retain the level-set geometry but alter the latent state. For multiple level sets, one writes

$$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$3

and places independent Gaussian priors on both the level-set functions $$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$4 and the region magnitudes $$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$5 (Reese et al., 2021). In point-source identification for the heat equation, the positive set of $$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$6 determines the support of a discrete source term,

$$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$7

with $$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$8 on $$D_i=\{x\in D:c_{i-1}\le u(x)<c_i\}, \qquad B_i=\{x\in B_R:c_{i-1}\le \phi(x)<c_i\},$$9 (Deng et al., 4 Sep 2025). In level-set Cox processes, a latent Gaussian process $$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$0 with thresholds $$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$1 induces regions

$$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$2

and the intensity becomes

$$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$3

(Gonçalves et al., 2020).

In sequential level-set estimation, the prior is typically placed on the unknown response rather than on an interface field. The standard choice is a Gaussian process,

$$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$4

with noisy observations $$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$5, $$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$6 (Cheng et al., 2024, Iwazaki et al., 2019). For high-dimensional settings, a Bayesian neural network surrogate replaces the GP: a prior $$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$7 is placed on the network weights, the posterior is approximated by MC-dropout variational inference $$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$8, and predictive mean and variance are computed by Monte Carlo over dropout samples (Ha et al., 2020).

3. Likelihoods, posteriors, and well-posedness

For geometric inverse problems, the observational model has the form

$$F(u)(x)=\sum_{i=1}^{n}\kappa_i\,\mathbf1_{D_i}(x), \qquad (F\phi)(x)=\sum_{i=1}^{L}b_i\,\mathbf1_{c_{i-1}\le \phi(x)<c_i}(x).$$9

or, more generally,

$f$0

The negative log-likelihood is written as

$f$1

and Bayes’ theorem yields a posterior absolutely continuous with respect to the Gaussian prior: $f$2 The principal well-posedness result is that the level-set map is discontinuous only on functions for which a threshold value is attained on a set of positive measure; under the Gaussian prior these events have zero probability, so the forward map is $f$3-almost surely continuous, the posterior is well-defined, and the data-to-posterior map is locally Lipschitz in the Hellinger distance (Iglesias et al., 2015, Huang et al., 2021).

Approximate Bayesian formulations modify the likelihood while retaining the level-set target. In Smoothed Approximate Bayesian Computation, one begins from the formal posterior

$f$4

and replaces the Dirac constraint by a Gaussian kernel $f$5. The approximate posterior becomes

$f$6

with $f$7 replaced in practice by the GP predictive density. As $f$8, the marginal $f$9 converges to $c\in\mathbb R$0 (Edwards et al., 2024).

Other formulations emphasize exactness rather than approximation. For level-set Cox processes with piecewise constant intensity, the obstacle is the intractable factor

$c\in\mathbb R$1

because the areas $c\in\mathbb R$2 are induced by a random level-set partition. The proposed solution is an almost-surely positive unbiased Poisson estimator $c\in\mathbb R$3, embedded in a pseudo-marginal infinite-dimensional MCMC algorithm with retrospective sampling. The resulting inference is “exact” in the sense that no space discretization approximation is used and MCMC error is the only source of inaccuracy (Gonçalves et al., 2020).

A distinct approximation route is the Gauss–Newton Laplace approximation for piecewise constant reconstructions. With

$c\in\mathbb R$4

the MAP point $c\in\mathbb R$5 is found by Gauss–Newton, and the Hessian

$c\in\mathbb R$6

defines the Gaussian approximation $c\in\mathbb R$7 (Reese et al., 2021).

4. Sampling mechanisms

The most common function-space sampler in Bayesian level-set inversion is preconditioned Crank–Nicolson. Given the current state $c\in\mathbb R$8 or $c\in\mathbb R$9, one draws $L=\{x\in\mathcal X:f(x)\ge c\},$0 and proposes

$L=\{x\in\mathcal X:f(x)\ge c\},$1

with acceptance probability

$L=\{x\in\mathcal X:f(x)\ge c\},$2

Because the proposal preserves the Gaussian prior, the method remains well-behaved under mesh refinement, and no explicit velocity law or reinitialization of the level set is required (Iglesias et al., 2015, Huang et al., 2021). In the acoustic scattering experiments, $L=\{x\in\mathcal X:f(x)\ge c\},$3 was chosen $L=\{x\in\mathcal X:f(x)\ge c\},$4, the mesh size was $L=\{x\in\mathcal X:f(x)\ge c\},$5, and $L=\{x\in\mathcal X:f(x)\ge c\},$6 samples were generated, with the last $L=\{x\in\mathcal X:f(x)\ge c\},$7 used to compute conditional-mean estimates (Huang et al., 2021).

For quasi-concave densities, the sampling mechanism is different. The level-set hit-and-run sampler constructs nested convex slices $L=\{x\in\mathcal X:f(x)\ge c\},$8, runs hit-and-run within each slice, enforces a warm-start condition through the volume-ratio rule

$L=\{x\in\mathcal X:f(x)\ge c\},$9

estimates slice masses, and finally reweights the pooled samples. For exponentially-tilted quasi-concave posteriors $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$0, the method augments the state with an auxiliary variable $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$1 and samples along one-dimensional chords with density proportional to $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$2 (Foster et al., 2012). Theoretical mixing guarantees cited there yield $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$3 hit-and-run mixing within a warm convex slice and overall cost $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$4 when the number of slices is $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$5 (Foster et al., 2012).

Smoothed ABC level-set MCMC uses a Metropolis–Hastings chain on $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$6. At each step one proposes $$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$7, draws

$$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$8

and accepts with

$$H^*=\{x\in D:p^*(x)\ge \alpha\}, \qquad p^*(x)=P_{s\sim g(\cdot|x)}[f(s)\le h].$$9

thereby targeting the smoothed posterior $A_c=\{\theta\in\Theta:f(\theta)=c\},$0 (Edwards et al., 2024).

When a Laplace approximation is adopted, posterior sampling is often performed without a Cholesky factorization. The construction in (Reese et al., 2021) applies a preconditioned Lanczos method to the preconditioned Hessian $A_c=\{\theta\in\Theta:f(\theta)=c\},$1, approximates $A_c=\{\theta\in\Theta:f(\theta)=c\},$2 through a low-dimensional tridiagonal projection, and returns

$A_c=\{\theta\in\Theta:f(\theta)=c\},$3

All steps require only matrix-vector products with $A_c=\{\theta\in\Theta:f(\theta)=c\},$4, $A_c=\{\theta\in\Theta:f(\theta)=c\},$5, and applications of $A_c=\{\theta\in\Theta:f(\theta)=c\},$6 and $A_c=\{\theta\in\Theta:f(\theta)=c\},$7 (Reese et al., 2021).

5. Sequential design for discovering level sets

In Bayesian experimental design, “sampling” often refers to adaptive data acquisition rather than posterior MCMC. Under input uncertainty, the Input-Uncertain Reliable Level-Set Estimation framework models the actual query point by $A_c=\{\theta\in\Theta:f(\theta)=c\},$8 and defines the reliability

$A_c=\{\theta\in\Theta:f(\theta)=c\},$9

A pointwise credible interval

$\phi:B_R\to\mathbb R$00

is formed around $\phi:B_R\to\mathbb R$01 using $\phi:B_R\to\mathbb R$02, and points are classified as reliable, unreliable, or uncertain according to whether $\phi:B_R\to\mathbb R$03, $\phi:B_R\to\mathbb R$04, or neither condition holds (Iwazaki et al., 2019). The acquisition function seeks the expected increase in the size of the reliably classified set, and the approximation $\phi:B_R\to\mathbb R$05 reduces the cost from $\phi:B_R\to\mathbb R$06 to $\phi:B_R\to\mathbb R$07. The same work proves an accuracy guarantee via Chebyshev’s inequality and finite-time termination with probability $\phi:B_R\to\mathbb R$08 under mild regularity and a vanishing randomization schedule with $\phi:B_R\to\mathbb R$09 (Iwazaki et al., 2019).

Posterior-sampling-based design adopts a simpler policy. In PS-BAX, one samples a function realization $\phi:B_R\to\mathbb R$10, forms the sampled super-level set

$\phi:B_R\to\mathbb R$11

and then chooses

$\phi:B_R\to\mathbb R$12

The final estimate is $\phi:B_R\to\mathbb R$13, and on a finite domain the method is asymptotically convergent under the complement-independence condition (Cheng et al., 2024).

High-dimensional level-set estimation replaces GP surrogates with Bayesian neural networks. For explicit LSE with threshold $\phi:B_R\to\mathbb R$14, the acquisition function is the mutual information

$\phi:B_R\to\mathbb R$15

while for implicit LSE with threshold $\phi:B_R\to\mathbb R$16, where $\phi:B_R\to\mathbb R$17, the acquisition is

$\phi:B_R\to\mathbb R$18

Their stated evaluation costs over $\phi:B_R\to\mathbb R$19 candidates are $\phi:B_R\to\mathbb R$20 and $\phi:B_R\to\mathbb R$21, respectively (Ha et al., 2020).

A common source of confusion is the relation between reliable and standard level-set estimation. The IU-rLSE formulation shows that when there is no input uncertainty, $\phi:B_R\to\mathbb R$22 and $\phi:B_R\to\mathbb R$23, so the method reduces to the usual GP-based LSE (Iwazaki et al., 2019).

6. Applications, empirical behavior, and limitations

The acoustic inverse medium-scattering study provides a representative function-space Bayesian level-set workflow. For a love-shaped obstacle with $\phi:B_R\to\mathbb R$24 and $\phi:B_R\to\mathbb R$25, a cross-shaped domain with $\phi:B_R\to\mathbb R$26 and $\phi:B_R\to\mathbb R$27, and two disjoint circles with $\phi:B_R\to\mathbb R$28 and $\phi:B_R\to\mathbb R$29, both regular Bayesian and level-set Bayesian methods recover the shape, but the level-set posterior mean shows a much sharper boundary or yields crisper interfaces (Huang et al., 2021). Trace-plots of $\phi:B_R\to\mathbb R$30 and selected coefficients are stationary after $\phi:B_R\to\mathbb R$31 steps, autocorrelation functions decay rapidly under pCN, and Gaussian noise with $\phi:B_R\to\mathbb R$32 produces mild widening of uncertainty belts while MAP and conditional-mean estimates remain accurate (Huang et al., 2021).

In the earlier geometric inverse-problem study, the same level-set methodology was demonstrated on an inverse potential problem and on discontinuous permeability in Darcy flow (Iglesias et al., 2015). The implementation used an $\phi:B_R\to\mathbb R$33 finite-difference grid, pCN chains of $\phi:B_R\to\mathbb R$34–$\phi:B_R\to\mathbb R$35 steps, and tuning of $\phi:B_R\to\mathbb R$36 for $\phi:B_R\to\mathbb R$37 acceptance. A recurring observation was that the prior correlation length $\phi:B_R\to\mathbb R$38 or smoothness $\phi:B_R\to\mathbb R$39 markedly influences the posterior mean and variance, with a “critical” length-scale above which the main geometric features are accurately recovered and posterior variance localizes near interfaces (Iglesias et al., 2015). This indicates that prior length-scale selection is not a secondary implementation detail.

For quasi-concave density sampling, empirical comparisons were made against Gibbs sampling on a spike-and-slab mixture and a Cauchy–normal posterior (Foster et al., 2012). In the spike-and-slab example, Gibbs gets “stuck” in one component for $\phi:B_R\to\mathbb R$40, with expected switching time scaling like $\phi:B_R\to\mathbb R$41, whereas the level-set sampler explores both components in $\phi:B_R\to\mathbb R$42 slices. In the multivariate normal example with $\phi:B_R\to\mathbb R$43, Gibbs shows extremely slow ACF decay, whereas the exponentially tilted level-set sampler has ACF near zero after one step (Foster et al., 2012).

For point-source identification in the heat equation, Bayesian level-set sampling is paired with a thinning mechanism that removes statistically unsupported candidate points (Deng et al., 4 Sep 2025). On $\phi:B_R\to\mathbb R$44 with $\phi:B_R\to\mathbb R$45 relative Gaussian noise, the reported reconstruction errors for $\phi:B_R\to\mathbb R$46 are $\phi:B_R\to\mathbb R$47, $\phi:B_R\to\mathbb R$48, $\phi:B_R\to\mathbb R$49, and $\phi:B_R\to\mathbb R$50 in $\phi:B_R\to\mathbb R$51, with source number errors $\phi:B_R\to\mathbb R$52. The chains mix well under pCN with effective sample sizes $\phi:B_R\to\mathbb R$53 for $\phi:B_R\to\mathbb R$54-modes, removal of thinning greatly increases false positives, increasing prior smoothness slightly under-estimates $\phi:B_R\to\mathbb R$55 but maintains location accuracy, and the method tolerates up to $\phi:B_R\to\mathbb R$56 noise with $\phi:B_R\to\mathbb R$57-error increasing linearly (Deng et al., 4 Sep 2025).

In high-dimensional level-set estimation, Bayesian neural networks outperform several GP-based baselines on synthetic and real problems, but the empirical results are explicitly qualified: in low dimensions ($\phi:B_R\to\mathbb R$58), standard GP-based LSE still has the edge, whereas a stated rule of thumb is that BNN-LSE wins when $\phi:B_R\to\mathbb R$59 (Ha et al., 2020). The same study reports that BNN plus hyperparameter tuning takes $\phi:B_R\to\mathbb R$60–$\phi:B_R\to\mathbb R$61 h per experiment on GPUs, while GP methods with batch size $\phi:B_R\to\mathbb R$62 can exceed $\phi:B_R\to\mathbb R$63 h in high dimensions due to cubic scaling in data (Ha et al., 2020). This is a practical limitation rather than a purely statistical one.

The literature also distinguishes clearly between exact and approximate Bayesian constructions. Smoothed ABC targets a tolerance-thickened level set and converges to the true level-set posterior only as $\phi:B_R\to\mathbb R$64 (Edwards et al., 2024); Gauss–Newton Laplace methods approximate the posterior by a Gaussian around the MAP (Reese et al., 2021); exact level-set Cox-process inference avoids spatial discretization error by using retrospective sampling and a pseudo-marginal scheme (Gonçalves et al., 2020). A further misconception addressed directly by the function-space pCN literature is that level-set inference requires evolving a classical level-set PDE with an explicit velocity field. In the Bayesian formulation, the interface is updated implicitly by MCMC on the level-set function, and no explicit velocity field is required (Iglesias et al., 2015).

Taken together, these developments show that Bayesian level set sampling spans well-posed function-space inversion, exact and approximate posterior simulation, sequential design for threshold discovery, and density sampling through nested convex slices. The common structure is the use of threshold-defined geometry as the primary inferential object, with uncertainty represented by a posterior distribution rather than by a single reconstructed interface.