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
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
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.
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 χ -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 χ -density along rays (black lines in right figure).
Figure 4: Geometric illustration of SRD rays ξ~i+rvi ( i=1,2,3 ), which start at ξ~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<ρ~in<ρ~out ( i=1 ), a ray originating inside the feasible set, with 0=ρ~in<ρ~out ( i=2 ), and a ray that completely misses the feasible set, illustrating the case in which ρ~out<ρ~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 108 standard Monte Carlo samples, is p≈0.648926 . Right: Gradient at the nominal control, which is largest close to ∂D2 , where the stochastic Neumann boundary force is applied.
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=104 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 .
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.
ξ
l
σ
σn
J(u)
−∞
1.16\times10−01
4.53\times10−01
5.55\times10−05
3.19\times1000
−0.2
8.54\times10−02
1.18\times10−01
5.57\times10−05
3.60\times1001
−0.05
6.94\times10−02
3.07\times10−02
5.71\times10−05
6.49\times1002
−0.01
6.71\times10−02
7.10\times10−03
1.23\times10−05
1.24\times1004
Table 1: Optimal kernel parameters and objective function ( 7 ) values for different lower bounds ξ . Results use hiSRD with N=104 samples. ξ=−∞ corresponds to the unconstrained Gaussian process regression.
Figure 8: Posterior Gaussian processes with optimal kernel parameters, with different lower bounds ξ 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: ξ≡−0.2 . Right: ξ≡−0.05 .
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.
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 gradient descent (SGD) admits diffusion approximations that replace the complicated randomness of stochastic gradients by Gaussian noise, providing a powerful tool for understanding its dynamics and long-time behavior. We investigate whether an analogous approximation principle holds for optimization over probability measures, where the objective is a functional defined on the Wasserstein space P2. The nonlinear geometry and infinite-dimensional nature of P2 prevent a direct extension of the classical Euclidean theory. Using Lions differentiability, we lift the problem to a linear Hilbert space, where higher-order differential calculus becomes available. We then construct a Gaussian random-field approximation whose velocity field matches the mean and covariance of the original stochastic gradient. By exploiting this moment matching through higher-order Taylor expansions, we show that the Gaussian approximation captures the SGD dynamics with second-order weak accuracy. Our result provides a rigorous foundation for replacing sample-driven randomness by analytically tractable Gaussian fluctuations in stochastic optimization over probability measures.
Maria Oprea, Qin Li, Yunan Yang
Center for Applied Mathematics, Cornell University · Department of Mathematics, University of Wisconsin–Madison, USA · Department of Mathematics, Cornell University, USA
Stochastic gradient descent (SGD) is a foundational algorithm for large-scale statistical learning and stochastic optimization. However, statistical inference based on SGD iterates remains challenging when stochastic gradients have infinite variance, as the relevant limiting distributions depend on unknown nuisance parameters. In this paper, we develop an efficient, model-agnostic methodology for constructing confidence regions from SGD trajectories that applies in both finite- and infinite-variance regimes. The procedure is based on a joint weak convergence result for the Polyak-Ruppert averaged estimator and an empirical second-moment normalizer constructed from stochastic gradients along the SGD trajectory. This joint limit yields a self-normalized statistic in which the leading tail-dependent scaling terms cancel. We then use a subsampling calibration scheme to estimate the relevant critical values, avoiding explicit estimation of tail indices, slowly varying functions, or stable-law parameters. The resulting confidence regions are straightforward to implement and are asymptotically valid under both the finite- and infinite-second-moment regimes. Simulation studies show reliable coverage in various settings, supporting the proposed method as a practical tool for uncertainty quantification in stochastic optimization.
Jose Blanchet, Peter Glynn, Wenhao Yang
Management Science and Engineering, Stanford University