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

1. Problem Statement

The paper addresses stochastic optimization problems under uncertainty, specifically those involving probabilistic (chance) constraints where the underlying random variables are infinite-dimensional.

• Mathematical Formulation: The core problem is minimizing an objective function $J(u)$ subject to a probabilistic constraint:
  $$\min_{u \in U_{ad}} J(u) \quad \text{s.t.} \quad \phi(u) := P(g(u, \xi) \le 0) \ge p$$
  Here, $\xi$ is an elliptical random variable (e.g., a Gaussian random field) taking values in an infinite-dimensional Hilbert space $\mathcal{H}$. The function $g$ is typically convex but non-smooth with respect to $\xi$.
• The Challenge:
  1. Infinite Dimensions: The random variable $\xi$ is often represented via an infinite series (e.g., Karhunen-Loève expansion), making direct integration impossible.
  2. Truncation Bias: Standard approaches truncate the series to a finite dimension $K$ and apply Spherical-Radial Decomposition (SRD). However, this introduces a truncation-induced bias.
  3. Variance Degeneration: While SRD is known to reduce variance compared to standard Monte Carlo (MC) in finite dimensions, the paper demonstrates that as the truncation dimension $K$ increases, the variance reduction benefit of SRD deteriorates, eventually converging to the variance of standard MC.
  4. Differentiability: Efficient optimization requires computing gradients of the probability function $\phi(u)$, which is difficult with standard MC and complex with truncated SRD when the mean lies outside the feasible set.

2. Methodology: Hybrid Infinite-dimensional SRD (hiSRD)

The authors propose a novel method called hybrid infinite-dimensional SRD (hiSRD) to eliminate bias while preserving the variance reduction and differentiability properties of SRD.

Core Concept

Instead of simply truncating the random variable, hiSRD decomposes the infinite-dimensional space $\mathcal{H}$ into two orthogonal subspaces:

1. Subspace $\mathcal{H}_K$: Spanned by the first $K$ basis functions (the "active" subspace).
2. Remainder Space $\mathcal{H}_R$: Spanned by the remaining infinite basis functions.

The random variable is decomposed as $\xi = \bar{\xi} + \xi_K + \xi_R$, where $\xi_K \in \mathcal{H}_K$ and $\xi_R \in \mathcal{H}_R$.

The Algorithm

The probability estimation is formulated as a nested integration:
$$\phi(u) = \int_{\mathcal{H}_R} \left[ \int_{S^{K-1}} \mu_{\chi_K}(\{r \ge 0 : g(u, \bar{\xi} + z + r L_K v) \le 0\}) \, d\mu_U(v) \right] d\mu_{\xi_R}(z)$$

1. Outer Loop (Monte Carlo): Samples $z$ are drawn from the distribution of the remainder $\xi_R$. This accounts for the infinite-dimensional tail, ensuring the estimator is unbiased.
2. Inner Loop (SRD): For each fixed remainder realization $z$, the problem is projected onto the $K$-dimensional subspace. SRD is applied to integrate over the unit sphere $S^{K-1}$ along radial directions.
   • Unlike standard finite-dimensional SRD, the "center" of the ray is shifted by $z$ ($\bar{\xi} + z$).
   • The method calculates the entry distance ($\rho_{in}$) and exit distance ($\rho_{out}$) of the ray relative to the feasible set defined by $g$.
   • The probability contribution for a specific ray is $F_{\chi_K}(\rho_{out}) - F_{\chi_K}(\rho_{in})$.

Differentiability

The paper extends existing differentiability results for finite-dimensional SRD to this hybrid setting.

• It establishes conditions (Rank-2 Constraint Qualification) under which the estimator $\tilde{\phi}_N(u)$ is differentiable.
• An explicit gradient formula is derived (Eq. 3.7), which involves the probability density function of the Chi-distribution and gradients of the constraint function $g$. This allows the use of gradient-based optimization algorithms (e.g., SQP).

3. Key Contributions

1. Generalization to Infinite Dimensions: The first rigorous extension of SRD to infinite stochastic dimensions that eliminates truncation bias.
2. Theoretical Variance Analysis: A proof showing that the variance reduction advantage of standard SRD decays as $O(K^{-1/2})$ as the dimension $K$ increases. This explains why simply increasing $K$ in standard SRD is inefficient.
3. Unbiased, Low-Variance Estimator: The hiSRD method provides an unbiased estimator with lower variance than standard Monte Carlo, regardless of the truncation dimension $K$.
4. Gradient Computation: The derivation of explicit gradient formulas for the infinite-dimensional case, enabling efficient gradient-based optimization.
5. Numerical Validation: Comprehensive studies on two distinct problem types:
   • PDE-Constrained Optimal Control: Controlling an elliptic PDE with uncertain Neumann boundary conditions under joint chance constraints.
   • Gaussian Process (GP) Regression: Optimizing kernel hyperparameters such that the posterior GP satisfies joint pointwise chance constraints.

4. Numerical Results

The authors conducted extensive numerical experiments comparing Standard MC, Finite-dimensional SRD (with truncation), and hiSRD.

• Bias Elimination: In the PDE control example, finite-dimensional SRD showed significant bias that dominated the error for large sample sizes. hiSRD remained unbiased.
• Variance Reduction:
  • In the GP regression example, hiSRD consistently achieved lower Root Mean Squared Error (RMSE) than standard MC.
  • The variance of hiSRD was decomposed into $V_{SRD}$ (conditional variance) and $V_{rem}$ (remainder variance). Results showed that $V_{rem}$ decays rapidly as $K$ increases, while $V_{SRD}$ remains the dominant but controlled component.
• Gradient Accuracy: hiSRD provided accurate gradient estimates for the kernel parameters, whereas standard MC is non-differentiable.
• Optimization Performance:
  • In the GP regression task, the algorithm successfully optimized kernel parameters to satisfy strict chance constraints (e.g., $p=0.95$).
  • As the constraint became stricter (lower bound increased), the optimizer naturally reduced the posterior variance ($\sigma$) and length scale ($l$) to ensure the constraint was met, demonstrating the method's utility in risk-averse design.

5. Significance and Limitations

Significance:

• Bridging Theory and Practice: The method resolves the trade-off between bias (from truncation) and variance (from high dimensions) in stochastic optimization.
• Applicability: It enables the solution of complex, infinite-dimensional chance-constrained problems in fields like engineering (PDE control) and machine learning (robust GP regression) that were previously intractable or required heuristic approximations.
• Efficiency: By combining SRD with MC, it leverages the fast convergence of SRD on low-dimensional subspaces while handling the infinite tail via MC, offering a computationally efficient path to high-accuracy solutions.

Limitations:

• Distribution Type: The method is currently limited to elliptical distributions (e.g., Gaussian, Student's t, Laplace).
• Convexity: The feasible sets defined by the constraints must be convex (typically intersections of half-spaces).
• Convergence Rate: While variance is reduced, the method still relies on (quasi) Monte Carlo sampling for the remainder, meaning the convergence rate is limited to the standard Monte Carlo rate ($O(N^{-1/2})$), though with a smaller constant.

In conclusion, this work provides a robust mathematical framework and practical algorithm for handling uncertainty in infinite-dimensional optimization, significantly advancing the state-of-the-art in chance-constrained optimization and probabilistic machine learning.