Regularized Linear Randomize-then-Optimize (RLRTO)

• RLRTO is an optimization-based sampling technique that encodes implicit priors by embedding regularization into a randomized least-squares framework for linear-Gaussian inverse problems.
• It enables the incorporation of hard constraints like nonnegativity, total variation, and sparsity directly into the optimization, enhancing computational efficiency.
• The method circumvents the need for explicit density formulations while offering a trade-off between speed and theoretical challenges in uniqueness and data dependence.

Searching arXiv for the cited papers to ground the article in the provided literature. Regularized Linear Randomize-then-Optimize (RLRTO) is an optimization-based sampling methodology for Bayesian inverse problems in which prior information is encoded through a randomized regularized optimization problem rather than through an explicit prior density. In the framework developed for implicit priors in Bayesian inverse problems, RLRTO is presented as “enhancing a Gaussian distribution with the deterministic effects of regularization by adding explicit regularization to an optimization-based characterization of the posterior distribution.” It is targeted at the linear-Gaussian setting, where posterior samples are obtained by perturbing data and prior terms and then solving a regularized least-squares problem. The resulting prior information is therefore implicit: it is encoded through Gaussian mean and covariance together with regularization or constraints, and the corresponding density is generally computationally intractable to write down, often improper, and in some settings data-dependent (Everink et al., 15 Sep 2025).

1. Conceptual position within Bayesian inverse problems

RLRTO arises in settings where prior information is hard or impossible to express as a tractable density, where one wishes to impose constraints such as nonnegativity or monotonicity, where regularity or sparsity should be promoted efficiently, or where optimization-based sampling is computationally preferable to generic Markov chain Monte Carlo. In this sense, it belongs to the broader class of implicit priors: the prior enters through a computational mechanism rather than a closed-form density (Everink et al., 15 Sep 2025).

For an explicit prior, the density is directly specified, for example

$\pi(x) \propto \exp(-R(x)).$

Bayes’ rule then yields an explicit posterior,

$$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$

RLRTO does not proceed this way. Instead, prior information is encoded in a randomized optimization problem whose optimizer is treated as a sample. The corresponding prior may exist, but it is implicit and typically not available in tractable density form (Everink et al., 15 Sep 2025).

The method also differs from Plug-and-Play (PnP) priors. PnP priors modify iterative sampling or optimization algorithms by replacing a proximal operator or prior-gradient step with a denoiser or restorator; in the cited framework, PnP-ULA is built from a restoration operator $D_\alpha$, and may or may not correspond to an actual underlying prior. RLRTO does not replace a proximal operator in a Langevin step. Instead, it inserts regularization directly into the randomized optimization characterization of a Gaussian posterior, so the implicit information enters through the objective function rather than through a denoiser inside a Markov chain step (Everink et al., 15 Sep 2025).

Standard MCMC methods such as MH, Langevin, and HMC require the posterior density, or at least its log-density and gradients. RLRTO avoids this by drawing randomized perturbations, solving a regularized optimization problem, and treating the optimizer output as the sample. This suggests a different computational tradeoff: when structured optimization is easier than running long Markov chains, RLRTO can be attractive, especially in high dimensions (Everink et al., 15 Sep 2025).

Framework	Prior encoding	Sampling mechanism
Explicit prior	Density specified directly	Bayes’ rule plus posterior sampling
PnP prior	Denoiser or restorator in an iterative scheme	Langevin or optimization with replaced proximal/prior step
RLRTO	Randomized regularized optimization problem	Draw perturbations, solve optimization, use optimizer as sample

2. Linear-Gaussian formulation and randomized optimization

The basic inverse problem considered for RLRTO has unknown $x \in \mathbb{R}^n$, data $y \in \mathbb{R}^m$, and observation model

$y = Ax + e, \qquad e \sim \mathcal{N}(0,\Sigma_e).$

The likelihood is

$$L_{y_{\mathrm{obs}}}(x)=\pi(y_{\mathrm{obs}}\mid x)\propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2\right).$$

With an explicit Gaussian prior

$$\pi(x)\propto \exp\!\left(-\frac12\|x-\mu\|_{\Sigma_p^{-1}}^2\right),$$

the posterior is Gaussian: $$\pi(x\mid y_{\mathrm{obs}})\propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2-\frac12\|x-\mu\|_{\Sigma_p^{-1}}^2\right).$$ RLRTO starts from this linear-Gaussian setting and alters the posterior representation by defining samples as optimizers of randomized regularized problems (Everink et al., 15 Sep 2025).

The immediate precursor is linear Randomize-then-Optimize (LRTO). For Gaussian likelihood and Gaussian prior, LRTO characterizes a posterior sample as the optimizer of

$$\arg\min_{x \in \mathbb{R}^n} \left\{ \frac12\|Ax-\hat y\|_{\Sigma_e^{-1}}^2 + \frac12\|x-\hat\mu\|_{\Sigma_p^{-1}}^2 \right\},$$

with randomized quantities

$$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$0

This yields an efficient sampler for Gaussian posteriors (Everink et al., 15 Sep 2025).

RLRTO modifies LRTO by adding a regularizer $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$1: $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$2 Within this construction, the randomization $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$3 introduces posterior variability, the optimization computes the sample, and the regularization encodes structure such as constraints or sparsity. A single sample is obtained by drawing $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$4 and $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$5, solving the regularized problem, and repeating this process for many independent draws (Everink et al., 15 Sep 2025).

A central interpretive point is that the posterior is no longer represented through an explicit prior density. Instead, it is represented through a randomized optimization procedure. The cited framework therefore treats RLRTO as an implicit prior mechanism rather than merely as a numerical solver (Everink et al., 15 Sep 2025).

3. Regularization, constraints, and prior interpretation

The regularization term $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$6 is the main vehicle through which RLRTO imposes structure. The implementation described for Bayesian inverse problems supports nonnegativity, box constraints, monotonicity-related constraints in one dimension, $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$7 sparsity, total variation (TV), and user-defined proximal or projection operators. The constraints specifically listed are "nonnegativity", "box", "increasing", "decreasing", "convex", and "concave", while the regularization options include "l1" and "TV" in 1D and 2D (Everink et al., 15 Sep 2025).

Several special cases clarify how the induced prior should be understood. When

$$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$8

with $$\pi(x \mid y_{\mathrm{obs}}) \propto \exp\!\left(-\frac12\|Ax-y_{\mathrm{obs}}\|_{\Sigma_e^{-1}}^2 - R(x)\right).$$9 an affine subspace and $D_\alpha$0 the indicator of $D_\alpha$1, the prior is a Gaussian restricted to $D_\alpha$2. When

$D_\alpha$3

the posterior gives positive probability to coordinates being exactly zero, and more generally to many lower-dimensional subsets of $D_\alpha$4. The cited implementation also presents the combined nonnegativity-plus-TV example

$D_\alpha$5

In hierarchical examples, the Gaussian prior term

$D_\alpha$6

may be removed by setting $D_\alpha$7, interpreted as an unbounded uniform distribution (Everink et al., 15 Sep 2025).

These constructions motivate the description of RLRTO priors as implicit, and in some cases improper or data-dependent. The method can represent structures that are awkward to encode in conventional density form, but this flexibility comes with theoretical subtleties. The cited framework explicitly notes that the resulting prior may be improper or depend on the data, and that uniqueness of the optimizer can be subtle when $D_\alpha$8 has a null space and $D_\alpha$9 is not positive definite (Everink et al., 15 Sep 2025).

A common misconception is to equate the regularizer $x \in \mathbb{R}^n$0 with an explicit prior density term. In the RLRTO setting this identification is generally not valid. The regularizer participates in the optimization that defines samples, but the induced prior is not normally written as a tractable $x \in \mathbb{R}^n$1. A plausible implication is that posterior interpretation must be tied to the sampling construction rather than to a closed-form prior density.

4. Algorithmic realization and CUQIpy implementation

The implementation described for Computational Uncertainty Quantification in Inverse Problems makes a design distinction between implicit priors used via Langevin-based samplers and implicit priors used via RLRTO, where the sampler itself is part of the prior representation. For RLRTO, the base implicit prior class is a regularized Gaussian prior that stores the Gaussian mean, Gaussian covariance or precision, and the regularization or constraints. Once that prior is included in a posterior, it is intended to be sampled using an RLRTO sampler (Everink et al., 15 Sep 2025).

The main class names are RegularizedGaussian for the implicit prior and RegularizedLinearRTO for the sampler. The implementation also includes the special cases ConstrainedGaussian, NonnegativeGaussian, RegularizedGMRF, ConstrainedGMRF, NonnegativeGMRF, and RegularizedUnboundedUniform for cases with no Gaussian prior term (Everink et al., 15 Sep 2025).

Solver choice is structure-dependent. The cited implementation uses FISTA for nonnegativity constraints, ADMM for TV regularization, and CGLS for the basic linear RTO case; other solvers such as L-BFGS-B may be used depending on the structure. The paper emphasizes that the optimization algorithm matters substantially for speed, so solver selection is part of the practical identity of RLRTO rather than a peripheral implementation detail (Everink et al., 15 Sep 2025).

The implementation is described as providing a convenient high-level API, easy switching between constraints and regularizers, support for hierarchical models, and efficiency for structured problems. It can also avoid the scalability issues of trans-dimensional samplers like RJMCMC for sparsity-like behavior. At the same time, the cited limitations remain essential: the implicit prior is generally not available in explicit density form; the induced prior may be improper or data-dependent; optimizer uniqueness can be subtle; and user-defined constraints or regularization have limitations in hierarchical formulations unless custom conjugacy rules are provided (Everink et al., 15 Sep 2025).

5. Empirical behavior in inverse problems and constrained regression

The inverse-problem examples used to illustrate RLRTO cover toy linear systems, hierarchical formulations, deconvolution, and PDE parameter inference. In a simple linear problem,

$x \in \mathbb{R}^n$2

with $x \in \mathbb{R}^n$3, the implementation uses RegularizedGaussian(np.zeros(2), 10, constraint="nonnegativity") together with RegularizedLinearRTO(posterior).sample(1000). The reported behavior is that nonnegativity-constrained samples are restricted to the positive orthant and TV regularization yields samples reflecting component similarity. The paper emphasizes that RLRTO can easily enforce hard constraints and produce interpretable samples (Everink et al., 15 Sep 2025).

Hierarchical RLRTO examples infer hyperparameters jointly with the unknown field, for example

$x \in \mathbb{R}^n$4

within a TV-regularized constrained inverse problem. The implementation supports hierarchical formulations with RLRTO for all implemented constraint and regularization options, and the regularization and hyperparameters can be sampled jointly using HybridGibbs (Everink et al., 15 Sep 2025).

The deconvolution example with a staircase-like signal is particularly informative about constraint misspecification. A Gaussian prior with a monotonicity constraint is compared with an unconstrained Gaussian and latent-variable step expansion priors. The constrained samples look much more like the ground truth, but an incorrect constraint structure can overconfidently misrepresent uncertainty. The reported monotonicity-constrained credible intervals peak near jumps, reflecting the imposed structure (Everink et al., 15 Sep 2025).

In PDE settings, RLRTO is used for 1D Poisson source term inference with a nonnegativity constraint and for 2D Poisson boundary value inference with GMRF plus TV regularization. In the source-term problem, RLRTO strictly enforces nonnegativity for each sample and the posterior variance is significantly reduced relative to an unconstrained GMRF reference, especially where the true source is exactly zero. In the boundary-value problem, TV-regularized RLRTO recovers the main trend and jump locations well, produces a clear staircase effect, and substantially reduces oscillations relative to a plain GMRF prior. By contrast, for the Poisson conductivity problem the mapping is nonlinear, so RLRTO is not applicable and MYULA is used instead. This delineates the intended scope of the method: linear-Gaussian problems (Everink et al., 15 Sep 2025).

A distinct application appears in monotonicity-constrained Gaussian process regression, where virtual derivative observations produce a constrained posterior over derivative values. There, RLRTO samples from the constrained posterior by solving

$x \in \mathbb{R}^n$5

The objective is strictly convex and the feasible set is convex, so the solution exists and is unique. The optimizer may land on the boundary $x \in \mathbb{R}^n$6, which permits flat monotone regions. In synthetic functions and differential-equation surrogate models, the reported results show that constrained methods improve over the unconstrained GP, that RLRTO is usually the best or among the best in MSE, and that it has the highest ESS per second among the compared constrained samplers (Zhang et al., 9 Jul 2025).

6. Extensions via prior transformation and relation to non-Gaussian priors

Although the core RLRTO framework is formulated for linear-Gaussian problems with regularization added directly in the optimization objective, related work extends the randomize-then-optimize idea to non-Gaussian priors through prior transformations. For $x \in \mathbb{R}^n$7-type priors of the form

$x \in \mathbb{R}^n$8

an exact transformation maps a standard Gaussian variable $x \in \mathbb{R}^n$9 to the physical variable $y \in \mathbb{R}^m$0 through

$y \in \mathbb{R}^m$1

After transformation, the posterior in $y \in \mathbb{R}^m$2 takes the RTO form

$y \in \mathbb{R}^m$3

and proposals are generated by solving a randomized optimization problem followed by a Metropolis-Hastings correction (Wang et al., 2016).

A key result for this transformed setting is that if the original forward model is linear, $y \in \mathbb{R}^m$4, and the $y \in \mathbb{R}^m$5-prior is transformed via $y \in \mathbb{R}^m$6, then the transformed problem satisfies the conditions required by RTO. In that case, the method targets the correct posterior after MH correction. This establishes a precise connection between linear RTO and regularized non-Gaussian inference: the regularization is not inserted directly as $y \in \mathbb{R}^m$7 in the original variable, but is absorbed into a prior transformation that restores a Gaussian-prior-like structure in transformed coordinates (Wang et al., 2016).

A closely related construction appears for Besov priors. There, the discrete posterior is

$y \in \mathbb{R}^m$8

which is not directly compatible with classical RTO when $y \in \mathbb{R}^m$9. The remedy is a componentwise transform from a standard Gaussian variable $y = Ax + e, \qquad e \sim \mathcal{N}(0,\Sigma_e).$0 to generalized Gaussian wavelet coefficients, yielding

$y = Ax + e, \qquad e \sim \mathcal{N}(0,\Sigma_e).$1

The transformed posterior becomes

$y = Ax + e, \qquad e \sim \mathcal{N}(0,\Sigma_e).$2

so the optimization is again a regularized stochastic least-squares problem, but now in the transformed variable and with a nonlinear transformed forward map. Because of that nonlinearity, Metropolis-Hastings correction is required (Horst et al., 20 Jun 2025).

These transformed-prior approaches are not identical to the direct RLRTO construction in linear Gaussian coordinates. Nevertheless, they show that the randomize-then-optimize paradigm extends beyond explicit Gaussian priors. This suggests a broader methodological family in which optimization-based sampling is preserved while prior complexity is shifted either into a regularizer in the objective or into a transformation of variables.

7. Advantages, limitations, and interpretive cautions

The principal strengths attributed to RLRTO are computationally efficient sampling with hard structural constraints and regularization, particularly in linear problems, and the ability to exploit fast deterministic optimization solvers. In constrained GP regression, the absence of a Markov chain means no burn-in, independent samples, and natural parallelizability; in the virtual-point monotonicity setting, the integrated autocorrelation time is reported to remain close to 1. In inverse problems, the method supports constraints, sparsity-promoting regularizers, and hierarchical models within a unified implementation (Zhang et al., 9 Jul 2025, Everink et al., 15 Sep 2025).

The main drawbacks are theoretical rather than merely practical. The induced prior or posterior is often implicit, possibly improper, and harder to analyze theoretically. It may depend on the data, and its density is generally unavailable in explicit form. In addition, optimizer uniqueness can become subtle when the forward operator has a null space or when the Gaussian precision is not positive definite. In transformed non-Gaussian settings such as Besov or $y = Ax + e, \qquad e \sim \mathcal{N}(0,\Sigma_e).$3-type priors, exactness generally requires MH correction, and the determinant or Jacobian terms in those corrections can create scalability limitations, with $y = Ax + e, \qquad e \sim \mathcal{N}(0,\Sigma_e).$4 determinant costs explicitly noted in the transformed-RTO literature (Everink et al., 15 Sep 2025, Horst et al., 20 Jun 2025, Wang et al., 2016).

A further caution concerns uncertainty quantification under strong structural assumptions. In the monotonicity-constrained deconvolution example, incorrect constraint structure can overconfidently misrepresent uncertainty. In the GP setting, the truncated-prior formulation forces strict positivity almost surely, whereas the RLRTO formulation can place samples on the boundary and therefore represent flat monotone regions with zero gradient. These contrasts show that the choice of constrained formulation affects not only computation but also the qualitative geometry of posterior uncertainty (Everink et al., 15 Sep 2025, Zhang et al., 9 Jul 2025).

Taken together, the cited work places RLRTO between explicit Bayesian modeling and algorithmic regularization. Unlike explicit priors, it does not require a tractable prior density; unlike PnP priors, it does not insert a denoiser into a Langevin step; unlike generic MCMC, it can exploit randomized optimization with structure-aware constraints. Its natural domain is the linear-Gaussian inverse problem, but related transformed constructions show how the same general logic can be adapted to Besov and $y = Ax + e, \qquad e \sim \mathcal{N}(0,\Sigma_e).$5-type priors, and to constrained Gaussian process regression, provided that the additional nonlinearities are accounted for appropriately (Everink et al., 15 Sep 2025, Horst et al., 20 Jun 2025, Wang et al., 2016, Zhang et al., 9 Jul 2025).