Infinite-dimensional spherical-radial decomposition for probabilistic functions, with application to constrained optimal control and Gaussian process regression

Kewei Wang; Georg Stadler

math.OCMar 20, 2026

Infinite-dimensional spherical-radial decomposition for probabilistic functions, with application to constrained optimal control and Gaussian process regression

Authors: Kewei Wang, Georg Stadler

Organizations: Courant Institute School of Mathematics, Computing and Data Science, New York University, New York, USA

Abstract

The spherical-radial decomposition (SRD) is an efficient method for estimating probabilistic functions and their gradients defined over finite-dimensional elliptical distributions. In this work, we generalize the SRD to infinite stochastic dimensions by combining subspace SRD with standard Monte Carlo methods. The resulting method, which we call hybrid infinite-dimensional SRD (hiSRD) provides an unbiased, low-variance estimator for convex sets arising, for instance, in chance-constrained optimization. We provide a theoretical analysis of the variance of finite-dimensional SRD as the dimension increases, and show that the proposed hybrid method eliminates truncation-induced bias, reduces variance, and allows the computation of derivatives of probabilistic functions. We present comprehensive numerical studies for a risk-neutral stochastic PDE optimal control problem with joint chance state constraints, and for optimizing kernel parameters in Gaussian process regression under the constraint that the posterior process satisfies joint chance constraints.

Figures & tables

Panel of Figure 1 — Figure 1: Comparison of root mean squared error (RMSE; $y$ -axis) versus stochastic truncation dimension ( $x$ -axis) for standard Monte Carlo (gray horizontal line), finite-dimensional SRD (dashed) and for the proposed method (solid). Results are for Gaussian process example detailed in section 5 , but with 1000 repeats to estimate the RMSE.

Panel of Figure 2 — Figure 2: RMSE difference between standard MC and SRD ( $y$ -axis) versus approximation dimension ( $x$ -axis). The dashed line indicates the slope theoretically expected from theorem 4 . Results are for Gaussian process example detailed in section 5 , but with 1000 repeats.

Figure 3: Illustration of standard SRD (left) and hybrid SRD (right) to estimate the probability for the three-dimensional set given by the intersection of half spaces defined by hyperplanes (light blue). The standard SRD integrates the three-dimensional

\chi

-density along 3D rays (black lines in left figure). The hybrid method uses Monte Carlo in the

z

-direction and SRD in the

x,y

-direction, integrating the two-dimensional

\chi

-density along rays (black lines in right figure).

Figure 4: Geometric illustration of SRD rays

\tilde{\xi}_{i}+r\boldsymbol{v}_{i}

(

i=1,2,3

), which start at

\tilde{\xi}_{i}

. The affine constraints defining the feasible set (white area) are shown as half spaces (in blue). The rays demonstrate three different intersection scenarios: a ray entering and exiting the feasible set from the outside, with

0<\tilde{\rho}_{\mathrm{in}}<\tilde{\rho}_{\mathrm{out}}

(

i=1

), a ray originating inside the feasible set, with

0=\tilde{\rho}_{\mathrm{in}}<\tilde{\rho}_{\mathrm{out}}

(

i=2

), and a ray that completely misses the feasible set, illustrating the case in which

\tilde{\rho}_{\mathrm{out}}<\tilde{\rho}_{\mathrm{in}}

(

i=3

Figure 5: Left: Root mean squared error (RMSE, $y$ -axis) of probability estimation ( 6 ) for different number $N$ of MC samples ( $x$ -axis). Compared are the finite-dimensional SRD for various KL-mode truncations $K$ , and the proposed hybrid method. The “exact” probability, computed with $10^{8}$ standard Monte Carlo samples, is $p\approx 0.648926$ . Right: Gradient at the nominal control, which is largest close to $\partial{\mathcal{D}}_{2}$ , where the stochastic Neumann boundary force is applied.

Panel of Figure 6 — Figure 6: Gaussian process regression: Shown on the left is a comparison of the root mean squared error (RMSE; $y$ -axis) of probability estimation with $N=10^{4}$ samples versus KL dimension $K$ ( $x$ -axis) for finite-dimensional SRD (dashed) and for the proposed method (solid), with both MC and QMC sampling. Shown for comparison is the result with standard Monte Carlo (gray horizontal line). Shown on the right are, for $N=500$ samples, the individual components discussed in section 3.2 contributing to the variance of hiSRD -MC sampling. The larger variance compared to the left figure is due to the smaller $N$ .

Panel of Figure 7 — Figure 7: Gaussian process regression: Companion plot to fig. 6 , but for estimation of the gradient. Probability estimation via standard MC sampling is non-differentiable and therefore not included, since its RMSE is unbounded.

$\underline{\xi}$	$l$	$\sigma$	$\sigma_{n}$	$\mathcal{J}(u)$
$-\infty$	$1.16\text{\times}{10}^{-01}$	$4.53\text{\times}{10}^{-01}$	$5.55\text{\times}{10}^{-05}$	$3.19\text{\times}{10}^{00}$
$-0.2$	$8.54\text{\times}{10}^{-02}$	$1.18\text{\times}{10}^{-01}$	$5.57\text{\times}{10}^{-05}$	$3.60\text{\times}{10}^{01}$
$-0.05$	$6.94\text{\times}{10}^{-02}$	$3.07\text{\times}{10}^{-02}$	$5.71\text{\times}{10}^{-05}$	$6.49\text{\times}{10}^{02}$
$-0.01$	$6.71\text{\times}{10}^{-02}$	$7.10\text{\times}{10}^{-03}$	$1.23\text{\times}{10}^{-05}$	$1.24\text{\times}{10}^{04}$

Table 1: Optimal kernel parameters and objective function ( 7 ) values for different lower bounds

\underline{\xi}

. Results use hiSRD with

N=10^{4}

samples.

\underline{\xi}=-\infty

corresponds to the unconstrained Gaussian process regression.

Figure 8: Posterior Gaussian processes with optimal kernel parameters, with different lower bounds

\underline{\xi}

for the chance constraint. Each figure shows the posterior mean (blue), the 90% credible interval of the Gaussian process (dashed), and 10 samples from the Gaussian process (gray). Left:

\underline{\xi}\equiv-0.2

. Right:

\underline{\xi}\equiv-0.05

Explore similar work

Riemannian Stochastic Optimization for Sufficient Dimension Reduction
May 29, 2026Thibault Pautrel, François PortierStochastic Optimization Riemannian Optimization
Laboratoire des Signaux et Systèmes (L2S), CentraleSupélec, Université Paris-Saclay, Gif-sur-Yvette, France · ENSAI, CREST, Bruz, France
Stochastic Gradient Descent over P2
Sep 11, 2026Maria Oprea, Qin Li, Yunan YangStochastic Optimization Stochastic Gradient Descent
Center for Applied Mathematics, Cornell University · Department of Mathematics, University of Wisconsin–Madison, USA · Department of Mathematics, Cornell University, USA
Statistical Inference for Stochastic Gradient Descent Beyond Finite Variance
May 25, 2026Jose Blanchet, Peter Glynn, Wenhao YangStochastic Optimization Confidence Region Estimation
Management Science and Engineering, Stanford University

May 29, 2026J/K move · Enter open · S save

Riemannian Stochastic Optimization for Sufficient Dimension Reduction

Thibault Pautrel, François Portier

Laboratoire des Signaux et Systèmes (L2S), CentraleSupélec, Université Paris-Saclay, Gif-sur-Yvette, France · ENSAI, CREST, Bruz, France

Stochastic Optimization Riemannian Optimization Inverse Dimensionality Reduction Dimensionality Reduction

Open paper AskPDF arXiv

Sufficient dimension reduction (SDR) makes high-dimensional regression tractable by projecting the covariates onto a low-dimensional subspace that preserves the conditional mean of the response. Existing gradient-based estimators either operate in the ambient space and suffer from the curse of dimensionality, or localize in the reduced space at a per-outer-iteration cost at least quadratic in the sample size. We show that minimizers of the population Minimum Average Variance Estimation (MAVE) risk approximate the same Grassmannian target as the Outer Product of Gradients (OPG), and recast the empirical criterion as a smooth maximization on the Stiefel manifold with closed-form Riemannian gradient. The resulting algorithm, SMAVE, combines sparse projected-space nearest-neighbor localization with Riemannian stochastic gradient ascent. A simplified version comes with almost-sure convergence and a non-asymptotic rate matching the standard non-convex stochastic first-order scaling. Empirically, SMAVE matches or improves on RMAVE's synthetic subspace recovery at moderate-to-high ambient dimension, and on four real datasets it uniformly improves over OPG and is competitive with or outperforms RMAVE at orders of magnitude lower runtime.

Abstract

Explore similar work

Riemannian Stochastic Optimization for Sufficient Dimension Reduction

Stochastic Gradient Descent over P2

Statistical Inference for Stochastic Gradient Descent Beyond Finite Variance