Stochastic Galerkin Method and Hierarchical Preconditioning for PDE-constrained Optimization

1. Problem Statement

The paper addresses optimal control problems governed by Partial Differential Equations (PDEs) with uncertain coefficients.

• Context: In many applications (fluid mechanics, structural optimization), PDE coefficients (e.g., diffusion $k(x, \xi)$) are not deterministic but modeled as random variables or stochastic processes due to measurement errors or environmental variations.
• Objective: Minimize a cost functional $J(y, u)$ that typically includes a tracking term (difference between state $y$ and target $y_d$), a control cost ($u$), and a variance penalty term ($\sigma^2(y)$) to ensure robustness against uncertainty.
• Challenge: The resulting optimization problems lead to large-scale, ill-conditioned linear systems (Karush-Kuhn-Tucker or KKT systems) when discretized. Solving these systems efficiently is difficult due to the "curse of dimensionality" introduced by the stochastic discretization and the ill-conditioning caused by the coupling of state, control, and adjoint variables.

2. Methodology

The authors adopt a "discretize-then-optimize" framework, integrating three main components:

A. Discretization

1. Spatial: Finite Element Method (FEM).
2. Stochastic: Stochastic Galerkin (SG) Method. The solution is expanded using Generalized Polynomial Chaos (gPC) basis functions (Hermite polynomials for Gaussian inputs). This transforms the stochastic PDE into a coupled deterministic system.
3. Temporal (for time-dependent problems): An all-at-once discretization scheme (coupling all time steps simultaneously) using implicit Euler or trapezoidal rules.

B. The Linear System

The discretization yields a large, symmetric, indefinite KKT system of the form:
$$\begin{bmatrix} M_\gamma & 0 & -A^T \\ 0 & \beta M & M^T \\ -A & M & 0 \end{bmatrix} \begin{bmatrix} y \\ u \\ \lambda \end{bmatrix} = \begin{bmatrix} M y_d \\ 0 \\ g \end{bmatrix}$$
where $M$ is the mass matrix, $A$ is the stiffness matrix (dependent on the random field), and $\lambda$ is the Lagrange multiplier. The system is solved using Krylov subspace methods (specifically GMRES).

C. Hierarchical Preconditioning (The Core Innovation)

To accelerate GMRES convergence, the authors propose a Hierarchical Gauss-Seidel (hGS) preconditioner.

• Schur Complement Approximation: The ideal preconditioner involves the exact Schur complement ($S_{exact}$), which is computationally prohibitive to invert. The authors approximate $S_{exact}$ using a matrix $Z$ derived from the stochastic expansion.
• Truncation Strategy: Instead of using the full stochastic expansion (which is expensive), they utilize a truncated hierarchical expansion. The matrix $Z$ is decomposed into blocks based on the polynomial order of the gPC basis.
• Algorithm: The preconditioner applies a block Gauss-Seidel iteration. It solves smaller, deterministic sub-problems (involving the mean stiffness matrix) sequentially while updating the right-hand side using matrix-vector products with truncated off-diagonal blocks of the stochastic matrices.
  • Mean-based ($r=1$): Uses only the mean coefficient.
  • Hierarchical ($r=m_\xi+1$): Includes first-order stochastic terms.
  • Full ($r=n_A$): Includes all terms.
• Time-Dependent Extension: For time-dependent problems, the method is extended to a Parallel-in-Time (PINT) framework. Because the all-at-once system has a block-diagonal structure with respect to time steps, the steady-state hGS preconditioner is applied simultaneously to all time steps, enabling massive parallelization.

3. Key Contributions

1. Novel Preconditioner for Stochastic Control: Development of a hierarchical preconditioner specifically tailored for the KKT systems of stochastic optimal control problems, extending previous work on forward stochastic PDEs.
2. Spectral Analysis: Rigorous proof that the proposed preconditioner is spectrally equivalent to the exact Schur complement. The eigenvalues of the preconditioned system are bounded within a fixed interval independent of the spatial mesh size ($h$) and stochastic discretization parameters ($m_\xi, p$), provided the expansion coefficients decay sufficiently.
3. Time-Dependent Extension: Successful adaptation of the method to time-dependent problems using an all-at-once formulation, creating a scalable parallel-in-time solver.
4. Efficiency vs. Accuracy Trade-off: Demonstration that a truncated expansion (specifically retaining terms up to order $m_\xi + 1$) offers performance nearly identical to the full expansion but with significantly reduced computational cost.

4. Numerical Results

The authors conducted extensive numerical experiments on both steady-state and time-dependent problems using the IFISS toolbox and GMRES.

• Mesh Independence: Iteration counts remained stable (sub-linear growth) as the spatial mesh was refined, confirming the theoretical spectral bounds.
• Parameter Robustness: The preconditioner performed robustly across varying regularization parameters ($\beta$) and variance levels ($\sigma_k$). Unlike mean-based preconditioners, which degraded significantly as uncertainty increased, the hierarchical method maintained stable convergence.
• Truncation Efficiency:
  • The $n_\tau = m_\xi + 1$ truncation strategy consistently achieved iteration counts comparable to the full expansion ($n_\tau = n_A$).
  • However, the CPU time for the truncated version was significantly lower because it avoided the expensive assembly and application of high-order stochastic coupling terms.
• Scalability: The parallel-in-time approach for time-dependent problems showed linear scaling with respect to the number of degrees of freedom (DoF), handling systems with over 12 million DoF efficiently.

5. Significance

This work provides a robust and scalable solution for a class of problems that are notoriously difficult to solve: large-scale stochastic PDE-constrained optimization.

• Practical Impact: It enables the use of high-fidelity stochastic models in optimization without being bottlenecked by the computational cost of solving the resulting linear systems.
• Theoretical Value: The spectral analysis provides a mathematical guarantee of mesh and parameter independence, which is crucial for designing reliable solvers.
• Methodological Advance: By bridging the gap between forward stochastic PDE solvers and optimization KKT systems, and extending this to time-dependent parallel-in-time frameworks, the paper sets a new standard for solving uncertainty quantification problems in control theory.

In summary, the paper proves that hierarchical truncation combined with Gauss-Seidel smoothing creates an optimal balance between computational cost and solver convergence, making large-scale stochastic optimal control problems tractable.