# Stochastic Galerkin Method and Hierarchical Preconditioning for PDE-constrained Optimization ††thanks: Dedicated to the memory of Professor Howard C. Elman, our dear colleague and mentor. We also want to thank Lehigh’s High Performance Computing systems for providing computational resources.
**Authors**: Zhendong Li, Akwum Onwunta, BedĹ™ich SousedĂk
> Department of Industrial and Systems Engineering, Lehigh University, Bethlehem, PA 18015 ( , ).
> Department of Mathematics and Statistics, University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250 ( ).
\headers
A Preconditioner for Stochastic Optimal ControlZ. Li, A. Onwunta, and B. SousedĂk
## Abstract
We develop efficient hierarchical preconditioners for optimal control problems governed by partial differential equations with uncertain coefficients. Adopting a discretize-then-optimize framework that integrates finite element discretization, stochastic Galerkin projection, and advanced time-discretization schemes, the approach addresses challenges of scaling large and ill-conditioned linear systems arising in uncertainty quantification. By exploiting sparsity of linear systems in stochastic Galerkin method, we formulate hierarchical preconditioners based on truncated stochastic expansion that strike an effective balance between computational cost and preconditioning quality. Numerical experiments demonstrate that the proposed preconditioners significantly accelerate the convergence of iterative solvers compared to existing methods, providing robust and efficient solvers for both steady-state and time-dependent optimal control problems under uncertainty.
keywords: Stochastic Galerkin method, preconditioning, iterative solvers, Gauss-Seidel method, hierarchical and multilevel preconditioning
{MSCcodes}
35R60, 65C20, 65F08, 65F10, 60H35, 65N22
## 1 Introduction
Optimal control problems governed by partial differential equations (PDEs) arise in numerous applications, including fluid mechanics, structural optimization, and inverse problems. These problems have been extensively studied over the past decades. For a theoretical overview and computational methods related to deterministic problems, we refer readers to, e.g., [elman_finite_2014, Troltzsch-2010-OCP]. In many practical applications, the PDE coefficients are uncertain. Such uncertainties originate from various sources, including measurement errors, model approximations, and environmental variations, and they are modeled as random variables or stochastic processes. Recently, there has been an increased interest in optimal control problems governed by PDEs with random coefficients, see e.g., [Martinez-Frutos-2018-OCP]. These stochastic problems are inherently more complex than their deterministic counterparts, thus necessitating specialized numerical methods.
Two strategies are commonly used: optimize-then-discretize and discretize-then-optimize. We adopt the latter, which discretizes the objective and PDE first, enabling the use of efficient numerical methods like finite elements [HDO2023, BDOS2016, Benner-2017-SOC, Benner-2016-BDP, elman_finite_2014].
For stochastic discretization, we employ the stochastic Galerkin method [Ghanem-1991-SFE, LeMaitre-2010-SMU, Lee-2013-SGM, Lord-2014-ICS, Rosseel-2012-OCS, Xiu-2010-NMS], which expands the solution using orthonormal polynomials. Although the Monte Carlo method is simpler and the stochastic collocation method [Chen-2014-WRB, Kouri-2013-MSC] decouples the problem, stochastic Galerkin method systematically captures uncertainties and preserves optimal convergence properties, offering computational efficiency for large-scale problems when carefully implemented.
In order to address the PDE-constrained optimization problems, we employ the discretize-then-optimize approach. In practice, it is common to combine temporal discretization (e.g., backward Euler scheme), stochastic expansions (e.g., generalized polynomial chaos expansions), and spatial discretization (e.g., finite element method). This leads to large-scale linear systems obtained via finite element (or possibly finite difference) discretizations, which are then typically solved using Krylov subspace methods, e.g., by the generalized minimal residual method (GMRES). These linear systems are often ill-conditioned, which causes a slow convergence of the iterative method. To address this issue, we construct preconditioners that improve the convergence of iterative solvers.
In this paper, we introduce a hierarchical preconditioning framework specifically tailored for stochastic PDE-constrained optimal control problems. Although the core concepts are inspired by preconditioners for PDE problems [Bespalov-2021-TPS, Sousedik-2014-THP, Sousedik-2014-HSC, Youngnoi-2025-ICS], an application to the Karush-Kuhn-Tucker (KKT) systems arising from optimization problems and an extension to time-dependent problems are nontrivial and constitute a primary contribution of this work. We provide a systematic derivation of the preconditioner for both steady-state and all-at-once formulation of time-dependent problems. The method is supported by a spectral analysis, proving that the proposed preconditioner is spectrally equivalent to the ideal but computationally prohibitive exact Schur complement. The performance is evaluated using a set of numerical experiments.
The paper is organized as follows. In Section 2.1, we introduce the steady-state stochastic optimal control problem and its discretization into a large-scale KKT system. Then, we construct the hierarchical Gauss-Seidel preconditioner for PDE-constrained optimal control problems (hGSoc). In Section 3, we extend the framework to the time-dependent case. We present an all-at-once discretization and develop a corresponding parallel-in-time preconditioner. In Section 4, we provide a spectral analysis of the preconditioners. In Section 5, we demonstrate the efficiency of the methods by a set of numerical experiments. Finally, in Section 6 we conclude and summarize our work.
## 2 Steady-state problem
### 2.1 Problem formulation
Following the standard framework [Ghanem-1991-SFE, Xiu-2010-NMS], let $(\Omega,\mathcal{F},\mathcal{P})$ be a complete probability space with $\xi:\Omega\to\Phi\subset\mathbb{R}^{m_{\xi}}$ a vector of independent random variables. We consider the Hilbert space $L^{2}(\Phi)$ equipped with the inner product $\langle u,v\rangle=\mathbb{E}[uv]$ , where $\mathbb{E}$ denotes the mathematical expectation with respect to the probability measure $\mu$ induced by $\xi$ .
We consider the steady-state optimal control problem given by
$$
\min_{y,u}{J}(y,u):=\,\frac{1}{2}\int_{\Phi}\int_{\mathcal{D}}|y-y_{d}|^{2}dx\,d\mu(\xi)+\frac{\beta}{2}\int_{\Phi}\int_{\mathcal{D}}|u|^{2}\,dx\,d\mu(\xi)+\frac{\gamma}{2}\int_{\mathcal{D}}|\sigma(y)|^{2}\,dx, \tag{1}
$$
subject to
$$
\begin{cases}-\nabla\cdot\left(\Bbbk\left(x,\xi\right)\nabla y\left(x,\xi\right)\right)=u\left(x,\xi\right),&\quad\text{in }\mathcal{D}\times\Phi,\\
y\left(x,\xi\right)=g\left(x\right),&\quad\text{on }\partial\mathcal{D}\times\Phi,\end{cases} \tag{2}
$$
where $y$ is the state, $y_{d}$ is the target state, and $u$ is the (distributed) control. The parameter $\gamma$ penalizes the variance $\sigma^{2}(y)$ of the state $y$ . Observe that both the state $y$ and control $u$ are stochastic. In view of the Doob-Dynkin lemma (see, e.g., [Babuska-2004-GFE]), both $y$ and $u$ admit the same parametric dependence on $\xi$ .
In computations, we work with a finite dimensional subspace $\mathcal{T}_{p}\mathcal{\subset}L^{2}\left(\Phi,\mathcal{B}(\Phi),\mu\right)$ , spanned by a set of generalized polynomial chaos (gPC) functions $\left\{\psi_{\ell}\left(\xi\right)\right\}_{\ell=1}^{n_{A}}$ with
$$
\langle\psi_{\ell},\psi_{k}\rangle=\mathbb{E}\left[\psi_{\ell}\psi_{k}\right]=\int_{\Phi}\psi_{\ell}(\xi)\psi_{k}(\xi)d\mu=\delta_{\ell k}, \tag{3}
$$
where $\left\{\psi_{i}\left(\xi\right)\right\}_{i=1}^{n_{A}}$ is a set of $m_{\xi}-$ dimensional, $p$ -order Hermite polynomials. $\psi_{1}(\xi)=1$ , and $\mathbb{E}[\psi_{k}(\xi)]=0$ for $k>1$ .
It is assumed that $\Bbbk(x,\xi)$ is bounded away from zero and infinity, i.e., $0<\Bbbk_{\min}\leq\Bbbk(x,\xi)\leq\Bbbk_{\max}<\infty$ for some constants $\Bbbk_{\min}$ , $\Bbbk_{\max}$ . We assume the stochastic PDE coefficient $\Bbbk(x,\xi)$ is represented by a gPC expansion
$$
\Bbbk(x,\xi)=\sum_{i=1}^{n_{A}}\kappa_{i}(x)\psi_{i}(\xi), \tag{4}
$$
Similarly, applying the stochastic Galerkin method to (1)–(2), both the state $y$ and control $u$ are expanded as
$$
v=\sum\nolimits_{i=1}^{n_{h}}\sum\nolimits_{k=1}^{n_{\xi}}v_{ik}\phi_{i}(x)\psi_{k}(\xi),\qquad v=y\,\text{or}\,u. \tag{5}
$$
Let $f(x,\xi)=f_{0}(x)+\sum_{i=1}^{n_{\xi}}\sqrt{\theta_{i}}f_{i}(x)\xi_{i}$ be a truncated Karhunen-Loève(KL) expansion of a Gaussian process defined on $\mathcal{D}$ , where $\xi_{i}$ are independent, identically distributed random variables, $f_{0}(x)$ is the mean function, and $(\theta_{i},f_{i}(x))$ is the $i$ -th eigenpair (eigenvalue and eigenfunctions, respectively) of the covariance function $C_{f}(x,y)$ such that
$$
\int_{\mathcal{D}}C_{f}(x,y)\varphi_{i}(y)dy=\theta_{i}\varphi_{i}(x), \tag{6}
$$
where $f_{i}(x)=\sqrt{\theta_{i}}\varphi_{i}(x)$ . In particular, if $\Bbbk(x,\xi)=\exp{[f(x,\xi)]}$ , then following [Ghanem-1999-NGS]
$$
\kappa_{i}{\color[rgb]{0,0,0}(x)}=\frac{\mathbb{E}[\psi_{i}(\xi-f)]}{\mathbb{E}[\psi_{i}^{2}]}\exp{\left[f_{0}(x)+\frac{1}{2}\sum_{j=1}^{m_{\xi}}f^{2}_{j}(x)\right]}.
$$
According to [Matthies-2005-GML], to guarantee a complete representation of the lognormal random field, the variables $y$ and $u$ are defined in the space $\mathcal{T}_{p}$ with dimension $n_{\xi}=\binom{m_{\xi}+p}{p}$ , and the dimension of $\Bbbk(x,\xi)$ is given by $n_{A}=\binom{m_{\xi}+2p}{2p}$ .
Using the stochastic Galerkin finite element framework to discretize problem (1)–(2), we get
$$
\min_{\mathbf{y},\mathbf{u}}J\left(\mathbf{y},\mathbf{u}\right)=\frac{1}{2}\left(\mathbf{y}-\mathbf{y}_{d}\right)^{T}\mathcal{M}\left(\mathbf{y}-\mathbf{y}_{d}\right)+\frac{\gamma}{2}\mathbf{y}^{T}\mathcal{M}_{\sigma}\mathbf{y}+\frac{\beta}{2}\mathbf{u}^{T}\mathcal{M}\mathbf{u},
$$
subject to
$$
-\mathcal{A}\mathbf{y}+\mathcal{M}\mathbf{u}=\mathbf{g}. \tag{7}
$$
Via the Lagrangian formulation, we can derive the first-order optimality conditions
$$
\left[\begin{array}[c]{ccc}\mathcal{M}_{\gamma}&0&-\mathcal{A}^{T}\\
0&\beta\mathcal{M}&\mathcal{M}^{T}\\
-\mathcal{A}&\mathcal{M}&0\end{array}\right]\left[\begin{array}[c]{c}\mathbf{y}\\
\mathbf{u}\\
{\color[rgb]{0,0,0}\boldsymbol{\lambda}}\end{array}\right]=\left[\begin{array}[c]{c}\mathcal{M}\mathbf{y}_{d}\\
0\\
\mathbf{g}\end{array}\right], \tag{8}
$$
where $\boldsymbol{\lambda}$ is a Lagrangian multiplier, $\mathcal{M}_{\gamma}=\mathcal{M+\gamma M}_{\sigma}=H^{\gamma}\otimes M$ , $H^{\gamma}=H_{1}+\gamma H^{\sigma}$ , and similarly, $h^{\gamma}$ is defined as the entry of $H^{\gamma}$ . In particular,
$$
\mathcal{M}=H_{1}\otimes M,\qquad\mathcal{M}_{\sigma}=H^{\sigma}\otimes M, \tag{9}
$$
where the mass matrix, $M$ is given by
$$
M_{a,b}=\int_{D}\phi_{a}(x)\phi_{b}(x)dx,\quad a,b=1,2,\dots,n_{h},
$$
for a suitably chosen finite element basis, $\{\phi_{i}(x)\}_{i=1}^{n_{h}}.$ The global stiffness matrix $\mathcal{A}$ is given by $\mathcal{A}=\sum_{\ell=1}^{n_{A}}H_{\ell}\otimes A_{\ell}$ . Here, the matrices $H_{\ell}$ are defined by their entries $(H_{\ell})_{jk}=h_{\ell jk}$ , with $h_{ijk}=\mathbb{E}[\psi_{i}\psi_{j}\psi_{k}]$ , where we note that all matrices $~H_{\ell}$ are symmetric, and
$$
H^{\sigma}=\operatorname{diag}\left(0,h_{1,jj}\right),\quad j=2,\dots,n_{\xi},
$$
which is obtained from $H_{1}$ by setting $h_{1,11}=0$ . Note also that in our settings $H_{1}=I_{n_{\xi}}$ , i.e., it is an identity matrix of size $n_{\xi}$ . Moreover, all matrices $A_{\ell}$ share the same sparsity pattern and are given by
$$
{\color[rgb]{0,0,0}A_{\ell}=[(A_{\ell})_{ab}],\quad(A_{\ell})_{ab}=\int_{\mathcal{D}}{\color[rgb]{0,0,0}\kappa_{\ell}(x)}\nabla\phi_{a}(x)\cdot\nabla\phi_{b}(x)\,dx,\quad a,b=1,2,\dots,n_{h}.} \tag{10}
$$
We note that the coefficient matrix in (8) is symmetric, indefinite, and in general very large. It is ill-conditioned, and therefore, a good preconditioner is required to solve the system efficiently. Next, we introduce a preconditioner to tackle this problem.
### 2.2 Schur complement-based preconditioner
A block-diagonal preconditioner for (8) (see, e.g. Benner et al. [Benner-2016-BDP]) is given by
$$
\mathcal{P}:=\begin{bmatrix}\mathcal{M}_{\gamma}&0&0\\
0&\beta\mathcal{M}&0\\
0&0&\mathcal{S}_{\text{exact}}\end{bmatrix}, \tag{11}
$$
where $\mathcal{S}_{\text{exact}}$ is the exact Schur complement
$$
\mathcal{S}_{\text{exact}}=\mathcal{A}\mathcal{M}_{\gamma}^{-1}\mathcal{A}^{T}+\frac{1}{\beta}\mathcal{M}. \tag{12}
$$
The first two blocks corresponding to (scaling of) the mass matrix are block diagonal, and the inverses are approximated by Chebyshev semi-iteration. However, forming and applying the inverse of $\mathcal{S}_{\text{exact}}$ is computationally prohibitive. The primary difficulty stems from its additive structure. Therefore, the key to an efficient solution lies in designing an approximation $\mathcal{S}$ that is spectrally equivalent to $\mathcal{S}_{\text{exact}}$ and easy to invert. Following the approach in [Benner-2016-BDP], we employ the approximation
$$
\mathcal{S}=\mathcal{ZM}_{\gamma}^{-1}\mathcal{Z}^{T},\qquad\mathcal{Z}=\mathcal{A}+\sqrt{\frac{1+\gamma}{\beta}}\mathcal{M}=\sum_{\ell=1}^{n_{A}}H_{\ell}\otimes\tilde{A}_{\ell}, \tag{13}
$$
where
$$
\tilde{A}_{1}=A_{1}+\sqrt{\frac{1+\gamma}{\beta}}M,\qquad\tilde{A}_{\ell}=A_{\ell},\quad\ell=2,\dots,n_{A}.
$$
Since $\mathcal{Z}$ is symmetric, we have $\mathcal{S}^{-1}=\mathcal{Z}^{-1}\mathcal{M}_{\gamma}\mathcal{Z}^{-1}$ . In [Benner-2016-BDP], the authors studied a mean-based preconditioner derived from (13) by dropping all the terms of $\mathcal{Z}$ except the first term; that is, $\mathcal{Z}\approx H_{1}\otimes\tilde{A}_{1}$ . However, this mean-based preconditioner is less effective when the standard deviation of the uncertain parameters increases.
In this study, we overcome this shortcoming by hierarchical preconditioning introduced in [Sousedik-2014-THP, Sousedik-2014-HSC] in the context of forward problems. More specifically, unlike [Sousedik-2014-THP, Sousedik-2014-HSC], we extend this strategy to the stochastic optimal control problem, and we also study the spectral properties of the new preconditioner.
To that end, observe first that, by applying the identity
$$
\left(V\otimes W\right)\operatorname{vec}(X)=\operatorname{vec}(WXV^{T}), \tag{14}
$$
to (8) yields the matricized system [Benner-2016-BDP, eq. (53)]:
$$
\displaystyle M\bar{Y}H^{\gamma}-\sum_{\ell=1}^{n_{A}}\tilde{A}_{\ell}^{T}\bar{\Lambda}H_{\ell} \displaystyle=M\bar{Y}_{d}H_{1}, \displaystyle\beta M\bar{U}H_{1}+M^{T}\bar{\Lambda}H_{1} \displaystyle=0, \displaystyle-\sum_{\ell=1}^{n_{A}}\tilde{A}_{\ell}^{T}\bar{Y}H_{\ell}+M\bar{U}H_{1} \displaystyle=\bar{G}. \tag{15}
$$
In addition, the matrix-vector of $\mathcal{Z}\bar{v}$ in (13) reads
$$
\mathcal{Z}\bar{v}=\operatorname{vec}\left(\sum\nolimits_{\ell=1}^{n_{A}}\tilde{A}_{\ell}\bar{V}H_{\ell}\right). \tag{16}
$$
### 2.3 Hierarchical Gauss-Seidel preconditioner
We first recall the preconditioner for the forward PDE problem from [Sousedik-2014-THP]. However, in this work, we present it in the matricized format as Algorithm 1. To set the notation, we note that the matricized format, which utilizes the isomorphism between $\mathbb{R}^{n_{h}n_{\xi}}$ and $\mathbb{R}^{n_{h}\times n_{\xi}}$ , defined via the operators $\operatorname{vec}$ and $\operatorname{mat}$ . Specifically,
$$
\bar{V}=\operatorname{mat}(\bar{v})=\left[v_{1},v_{2},\ldots,v_{n_{\xi}}\right]\in\mathbb{R}^{n_{h}\times n_{\xi}}, \tag{17}
$$
where the column $k$ contains the coefficients associated with the basis function $~\psi_{k}$ , and $\bar{v}=\operatorname{vec}\left(\bar{V}\right)\in\mathbb{R}^{n_{h}n_{\xi}}$ . We will use lowercase letters for the vectorized representation and uppercase letters for the matricized counterpart; so, e.g., $\bar{R}=\operatorname{mat}(\bar{r})$ , etc.
Moreover, we will denote by $~\bar{V}_{\left(i:n\right)}$ a submatrix of $\bar{V}$ containing columns $i,i+1,\dots,n$ , and, in particular, $\bar{V}=\bar{V}_{\left(1:n_{\xi}\right)}$ . There are two components of the preconditioner. The first component consists of block-diagonal solves with blocks of varying sizes. The second component is used in the setup of the right-hand sides for the solves, and consists of matrix-vector products by certain subblocks of the stochastic Galerkin matrix by vectors of corresponding sizes. To this end, we will write $\left[h_{\tau,(\ell)(k)}\right]$ , with $(\ell)$ and $(k)$ denoting a set of (consecutive) rows and columns of matrix $H_{\tau}$ so that, in particular, ${H}_{\tau}=\left[h_{\tau,(1:n_{\xi})(1:n_{\xi})}\right]$ . Let us also denote $\bar{v}_{(\ell)}=\operatorname{vec}(\bar{V}_{(\ell)})$ . Then, the matrix-vector products can be written, cf. (14) and noting the symmetry of $~H_{\tau}$ , as
$$
\bar{v}_{(\ell)}=\sum\nolimits_{\tau\in\mathcal{I}_{\tau}}(\left[h_{\tau,(\ell)(k)}\right]\otimes\tilde{A}_{\tau})\bar{u}_{(k)}\quad\Leftrightarrow\quad\bar{V}_{(\ell)}=\sum\nolimits_{\tau\in\mathcal{I}_{\tau}}\tilde{A}_{\tau}\bar{U}_{(k)}\left[h_{\tau,(k)(\ell)}\right], \tag{18}
$$
where $~\mathcal{I}_{\tau}$ is an index set $~\mathcal{I}_{\tau}\subseteq\left\{1,\dots,n_{\tau}\right\}$ indicating that the matrix-vector products may be truncated. Possible strategies for truncation are discussed in [Sousedik-2014-THP]. In this study, we use $\mathcal{I}_{\tau}=\left\{1,\dots,n_{\tau}\right\}$ with $n_{\tau}=\binom{m_{\xi}+p_{\tau}}{p_{\tau}}$ for some $p_{\tau}\leq p$ . In particular, we set $\tau=\{0,1,2\}$ and with $\mathcal{I}_{\tau}=\emptyset$ both preconditioners in Algorithms 1 and 2 – 3 reduce to mean-based variants. We also note that, since the initial guess is zero in Algorithm 1, the multiplications by $~\mathcal{F}_{1}$ and $\mathcal{F}_{d+1}$ vanish from (19)–(20).
Algorithm 1 [Sousedik-2014-THP, Algorithm 3] Hierarchical Gauss-Seidel preconditioner (hGS)
The preconditioner $\mathcal{Z}_{hGS}:\bar{R}\longmapsto\bar{V}$ is defined as follows.
1: Set the initial solution $~\bar{V}$ to zero and update in the following steps:
2: Solve
$$
\tilde{A}_{1}\bar{V}_{\left(1\right)}=\bar{R}_{\left(1\right)}-\mathcal{F}_{1},\qquad\text{where }\mathcal{F}_{1}=\sum_{\tau\in\mathcal{I}_{\tau}}\tilde{A}_{\tau}\bar{V}_{\left(2:n_{\xi}\right)}\left[h_{\tau,\left(2:n_{\xi}\right)\left(1\right)}\right]. \tag{19}
$$
3: for $d=1,\ldots,p-1$ do
4: Set $\ell=\left(n_{\ell}+1:n_{u}\right),\text{ where }n_{\ell}=\binom{m_{\xi}+d-1}{d-1}\text{ and }n_{u}=\binom{m_{\xi}+d}{d}$ .
5: Solve
$$
\tilde{A}_{1}\bar{V}_{\left(\ell\right)}=\bar{R}_{\left(\ell\right)}-\mathcal{E}_{d+1}-\mathcal{F}_{d+1}, \tag{20}
$$
where
$$
\mathcal{E}_{d+1}=\sum_{\tau\in\mathcal{I}_{\tau}}\tilde{A}_{\tau}\bar{V}_{\left(1:n_{\ell}\right)}\left[h_{\tau,\left(1:n_{\ell}\right)\left(\ell\right)}\right],\quad\mathcal{F}_{d+1}=\sum_{\tau\in\mathcal{I}_{\tau}}\tilde{A}_{\tau}\bar{V}_{\left(n_{u}+1:n_{\xi}\right)}\left[h_{\tau,\left(n_{u}+1:n_{\xi}\right)\left(\ell\right)}\right].
$$
6: end for
7: Set $\ell=\left(n_{u}+1:n_{\xi}\right)$ .
8: Solve
$$
\tilde{A}_{1}\bar{V}_{\left(\ell\right)}=\bar{R}_{\left(\ell\right)}-\mathcal{E}_{p+1},\qquad\text{where }\mathcal{E}_{p+1}=\sum_{t\in\mathcal{I}_{\tau}}\tilde{A}_{\tau}\bar{V}_{\left(1:n_{u}\right)}\left[h_{\tau,\left(1:n_{u}\right)\left(\ell\right)}\right],
$$
9: for $d=p-1,\ldots,1$ do
10: Set $\ell=\left(n_{\ell}+1:n_{u}\right),\text{ where }n_{\ell}=\binom{m_{\xi}+d-1}{d-1}\text{ and }n_{u}=\binom{m_{\xi}+d}{d}$ .
11: Solve (20).
12: end for
13: Solve (19).
Next, we apply the hGS strategy to form a preconditioner for the KKT system. The preconditioner is formulated as Algorithm 2 – 3. It adapts the core idea of Algorithm 1 to handle the coupled variables corresponding to the state $\bar{V}^{Y}$ , control $\bar{V}^{U}$ , and adjoint $\bar{V}^{\Lambda}$ simultaneously at each hierarchical level. A key computational step of Algorithm 2 – 3 entails a solve with $\tilde{P}=\operatorname{blkdiag}\begin{bmatrix}M,\beta M,\tilde{A}_{1}\end{bmatrix}$ , which serves as an auxiliary block-diagonal preconditioner derived from the blocks of the deterministic KKT system. Since this system is relatively small and constant across all hierarchical levels, it can be handled efficiently, for instance, by computing a direct factorization of $\tilde{P}$ once and reusing it for all subsequent solves. The overall performance of the preconditioner is thus determined by the cost of these deterministic solves and the number of truncated off-diagonal matrix-vector products.
Algorithm 2 hGS preconditioner for the optimal control problem (hGSoc)
The preconditioner $\mathcal{P}_{hGSoc}:\left(\bar{R}^{Y},\bar{R}^{U},\bar{R}^{\Lambda}\right)\longmapsto\left(\bar{V}^{Y},\bar{V}^{U},\bar{V}^{\Lambda}\right)$ is defined as follows.
1: Set the initial solution $\left(\bar{V}^{Y},\bar{V}^{U},\bar{V}^{\Lambda}\right)$ to zero and update in the following steps:
2: Solve
$$
\tilde{P}\left[\begin{array}[c]{c}\bar{V}_{\left(1\right)}^{Y}\\
\bar{V}_{\left(1\right)}^{U}\\
\bar{V}_{\left(1\right)}^{\Lambda}\end{array}\right]=\left[\begin{array}[c]{c}\bar{R}_{\left(1\right)}^{Y}-\mathcal{C}_{1}+\mathcal{D}_{1}\\
\bar{R}_{\left(1\right)}^{U}-\mathcal{E}_{1}-\mathcal{F}_{1}\\
\bar{R}_{\left(1\right)}^{\Lambda}+\mathcal{G}_{1}-\mathcal{H}_{1}\end{array}\right], \tag{21}
$$
where
$$
\displaystyle\mathcal{C}_{1} \displaystyle=M\bar{V}_{\left(2:n_{\xi}\right)}^{Y}\left[h_{\left(2:n_{\xi}\right)\left(1\right)}^{\gamma}\right], \displaystyle\mathcal{D}_{1} \displaystyle=\sum_{\tau\in\mathcal{I}_{\tau}}\tilde{A}_{\tau}\bar{V}_{\left(2:n_{\xi}\right)}^{\Lambda}\left[h_{\tau,\left(2:n_{\xi}\right)\left(1\right)}\right], \displaystyle\mathcal{E}_{1} \displaystyle=\beta M\bar{V}_{\left(2:n_{\xi}\right)}^{U}\left[h_{1,\left(2:n_{\xi}\right)\left(1\right)}\right], \displaystyle\mathcal{F}_{1} \displaystyle=M\bar{V}_{\left(2:n_{\xi}\right)}^{\Lambda}\left[h_{1,\left(2:n_{\xi}\right)\left(1\right)}\right], \displaystyle\mathcal{G}_{1} \displaystyle=\sum_{\tau\in\mathcal{I}_{\tau}}\tilde{A}_{\tau}\bar{V}_{\left(2:n_{\xi}\right)}^{Y}\left[h_{\tau,\left(2:n_{\xi}\right)\left(1\right)}\right], \displaystyle\mathcal{H}_{1} \displaystyle=M\bar{V}_{\left(2:n_{\xi}\right)}^{U}\left[h_{1,\left(2:n_{\xi}\right)\left(1\right)}\right]. \tag{1}
$$
3: for $d=1,\ldots,p-1$ do
4: Set $\ell=\left(n_{\ell}+1:n_{u}\right),\text{ where }n_{\ell}=\binom{m_{\xi}+d-1}{d-1}\text{ and }n_{u}=\binom{m_{\xi}+d}{d}$ .
5: Solve
$$
\tilde{P}\left[\begin{array}[c]{c}\bar{V}_{\left(\ell\right)}^{Y}\\
\bar{V}_{\left(\ell\right)}^{U}\\
\bar{V}_{\left(\ell\right)}^{\Lambda}\end{array}\right]=\left[\begin{array}[c]{c}\bar{R}_{\left(\ell\right)}^{Y}-\mathcal{C}_{d+1}+\mathcal{D}_{d+1}\\
\bar{R}_{\left(\ell\right)}^{U}-\mathcal{E}_{d+1}-\mathcal{F}_{d+1}\\
\bar{R}_{\left(\ell\right)}^{\Lambda}+\mathcal{G}_{d+1}-\mathcal{H}_{d+1}\end{array}\right], \tag{22}
$$
where
| | $\displaystyle\mathcal{C}_{d+1}$ | $\displaystyle=M\left(\bar{V}_{\left(1:n_{\ell}\right)}^{Y}\left[h_{\left(1:n_{\ell}\right)\left(\ell\right)}^{\gamma}\right]+\bar{V}_{\left(n_{u}+1:n_{\xi}\right)}^{Y}\left[h_{\left(n_{u}+1:n_{\xi}\right)\left(\ell\right)}^{\gamma}\right]\right),$ | |
| --- | --- | --- | --- |
6: end for
Algorithm 3 hGS preconditioner for the optimal control problem (hGSoc), cont’d
7: Set $\ell=\left(n_{u}+1:n_{\xi}\right)$ .
8: Solve
$$
\tilde{P}\left[\begin{array}[c]{c}\bar{V}_{\left(\ell\right)}^{Y}\\
\bar{V}_{\left(\ell\right)}^{U}\\
\bar{V}_{\left(\ell\right)}^{\Lambda}\end{array}\right]=\left[\begin{array}[c]{c}\bar{R}_{\left(\ell\right)}^{Y}-\mathcal{C}_{p+1}+\mathcal{D}_{p+1}\\
\bar{R}_{\left(\ell\right)}^{U}-\mathcal{E}_{p+1}-\mathcal{F}_{p+1}\\
\bar{R}_{\left(\ell\right)}^{\Lambda}+\mathcal{G}_{p+1}-\mathcal{H}_{p+1}\end{array}\right],
$$
where
| | $\displaystyle\mathcal{C}_{p+1}$ | $\displaystyle=M\bar{V}_{\left(1:n_{u}\right)}^{Y}\left[h_{\left(1:n_{u}\right)\left(\ell\right)}^{\gamma}\right],$ | $\displaystyle\mathcal{D}_{p+1}$ | $\displaystyle=\sum_{\tau\in\mathcal{I}_{\tau}}A_{\tau}\bar{V}_{\left(1:n_{u}\right)}^{\Lambda}\left[h_{\tau,\left(1:n_{u}\right)\left(\ell\right)}\right],$ | |
| --- | --- | --- | --- | --- | --- |
9: for $d=p-1,\ldots,1$ do
10: Set $\ell=\left(n_{\ell}+1:n_{u}\right),\text{ where }n_{\ell}=\binom{m_{\xi}+d-1}{d-1}\text{ and }n_{u}=\binom{m_{\xi}+d}{d}$ .
11: Solve (22).
12: end for
13: Solve (21).
## 3 Time-dependent problem
The time-dependent optimal control problem is given by
$$
\displaystyle\min_{y,u}\mathcal{J}(y,u) \displaystyle= \displaystyle\frac{1}{2}\int_{0}^{T}\int_{\Phi}\int_{\mathcal{D}}|y-y_{d}|^{2}dx\,d\mu(\xi)\,dt+\frac{\beta}{2}\int_{0}^{T}\int_{\Phi}\int_{\mathcal{D}}|u|^{2}dx\,d\mu(\xi)\,dt \displaystyle+\frac{\gamma}{2}\int_{0}^{T}\int_{\Phi}\int_{\mathcal{D}}|\sigma(y)|^{2}dx\,d\mu(\xi)\,dt, \tag{23}
$$
subject to
$$
\begin{cases}\begin{aligned} \frac{\partial y(t,\mathbf{x},\xi)}{\partial t}-\nabla\cdot(\Bbbk(\mathbf{x},\xi)\nabla y(t,\mathbf{x},\xi))&=u(t,\mathbf{x},\xi)\text{ in }(0,T]\times\mathcal{D}\times\Phi,\\
y(t,\mathbf{x},\xi)&=g\text{ on }(0,T]\times\partial\mathcal{D}\times\Phi,\\
y(0,\mathbf{x},\xi)&=y_{0}\text{ in }\mathcal{D}\times\Phi.\end{aligned}\end{cases} \tag{24}
$$
After the application of the stochastic Galerkin finite element discretization to (23), and using the trapezoidal rule for the time discretization, where $n_{t}=T/\tau$ is the number of time steps over the interval $[0,T]$ with time-step size $\tau$ , we obtain
$$
\min_{\mathbf{y},\mathbf{u}}\mathcal{J}\left(\mathbf{y},\mathbf{u}\right)=\frac{\tau}{2}(\mathbf{y}-\mathbf{y}_{d})^{T}\left(D\otimes\mathcal{M}_{\gamma}\right)(\mathbf{y}-\mathbf{y}_{d})+\frac{\tau\beta}{2}\mathbf{u}^{T}\left(D\otimes\mathcal{M}\right)\mathbf{u}, \tag{25}
$$
where $\mathcal{M}$ and $\mathcal{M}_{\gamma}$ are defined in (9), $D$ in (28) below, $\mathbf{y}$ , $\mathbf{y}_{d}$ , and $\mathbf{u}$ are vectors corresponding to the state, desired state, and control, respectively, that contain concatenated vectors $\mathbf{y}_{i},\mathbf{{y}_{d}}_{i},\mathbf{u}_{i}\in\mathbb{R}^{n_{h}n_{\xi}\times 1},i=1,\ldots,n_{t}$ , due to the time-stepping,
$$
\mathbf{y}=\begin{bmatrix}\mathbf{y}_{1}&\dots&\mathbf{y}_{n_{t}}\end{bmatrix}^{T},\,\mathbf{y}_{d}=\begin{bmatrix}\mathbf{{y}_{d}}_{1}&\dots&\mathbf{{y}_{d}}_{n_{t}}\end{bmatrix}^{T},\,\text{and }\mathbf{u}=\begin{bmatrix}\mathbf{u}_{1}&\dots&\mathbf{u}_{n_{t}}\end{bmatrix}^{T}.
$$
After the application of the stochastic Galerkin finite element discretization to (24), and using the implicit Euler method for the time discretization, we obtain
$$
M\mathbf{y}_{k}+\tau A\mathbf{y}_{k}=M\mathbf{y}_{k-1}+\tau M\mathbf{u}_{k}. \tag{26}
$$
Combining all time steps of (26) in all-at-once discretization ([MR3023467, rees2010all]), we can write
$$
\mathcal{A}_{t}\mathbf{y}-\tau\mathcal{N}\mathbf{u}=\begin{bmatrix}\mathcal{M}\mathrm{y}_{0}+\mathbf{g}\,,&\mathbf{g}\,,&\dots,\mathbf{g}\end{bmatrix}^{T}\eqqcolon\mathbf{d},
$$
the vector $\mathbf{g}$ contains the contributions from the Dirichlet boundary data at each time step, and
$$
\mathcal{A}_{t}=\left(I_{n_{t}}\otimes\mathcal{L}\right)-\left(C\otimes\mathcal{M}\right),\,\mathcal{N}=I_{n_{t}}\otimes H_{1}\otimes M,\text{ and }\mathcal{L}=H_{1}\otimes\left(M+\tau A_{1}\right)+\tau\sum_{\ell=2}^{n_{A}}H_{\ell}\otimes A_{\ell}. \tag{27}
$$
where the matrix $C$ and matrix $D$ , used in (25) and also below, are defined as
$$
C=\operatorname{subdiag}(1,1,\dots,1)_{-1},\quad D=\operatorname{diag}\left(\frac{1}{2},1,\dots,1,\frac{1}{2}\right). \tag{28}
$$
Forming the Lagrangean and applying the first-order optimality conditions, we get
$$
\begin{bmatrix}\tau D\otimes\mathcal{M}_{\gamma}&0&-\mathcal{A}_{t}^{T}\\
0&\beta\tau D\otimes\mathcal{M}&\tau\mathcal{N}^{T}\\
-\mathcal{A}_{t}&\tau\mathcal{N}&0\end{bmatrix}\begin{bmatrix}\mathbf{y}\\
\mathbf{u}\\
{\color[rgb]{0,0,0}\boldsymbol{\lambda}}\end{bmatrix}=\begin{bmatrix}\tau\left(D\otimes\mathcal{M}\right)\cdot\left(\mathbf{1}_{n_{t}}\otimes\mathbf{y}_{d}\right)\\
\mathbf{0}\\
\mathbf{d}\end{bmatrix}, \tag{29}
$$
where the boundary-inclusive source vector $\mathbf{d}$ is defined as above, and $\mathbf{1}_{n_{t}}\in\mathbb{R}^{n_{t}\times 1}$ is the column all-ones vector.
### 3.1 PINT-based block-diagonal hierarchical Gauss-Seidel preconditioner
In analogy to (11), we propose a preconditioner for (29) as
$$
\mathcal{P}_{\text{hGSoc-PINT}}\approx\begin{bmatrix}\tau D\otimes\mathcal{M}_{\gamma}&&\\
&\tau\beta D\otimes\mathcal{M}&\\
&&\bar{\mathcal{S}}\end{bmatrix}, \tag{30}
$$
where $\bar{\mathcal{S}}$ is a computationally efficient approximation of the exact Schur complement
$$
\bar{\mathcal{S}}_{\text{exact}}=\frac{1}{\tau}\mathcal{A}_{t}\left(D\otimes\mathcal{M}_{\gamma}\right)^{-1}\mathcal{A}_{t}^{T}+\frac{\beta}{\tau}\mathcal{N}\left(D\otimes\mathcal{M}\right)^{-1}\mathcal{N}^{T}. \tag{31}
$$
The first two blocks that correspond to (scaling of) the mass matrix are approximated by Chebyshev semi-iteration. Since the iteration entails matrix-vector multiplications, we note that using (14) we have
| | $\displaystyle\left(\tau D\otimes H_{1}\otimes M\right)\mathbf{v_{1}}$ | $\displaystyle=\tau\begin{bmatrix}\frac{1}{2}\operatorname{vec}{\left(MV_{11}H_{\gamma}\right)}&\operatorname{vec}{\left(MV_{12}H_{\gamma}\right)}&\cdots&\frac{1}{2}\operatorname{vec}{\left(MV_{1n_{t}}H_{\gamma}\right)}\end{bmatrix}^{T}$ | |
| --- | --- | --- | --- |
where $\mathbf{v_{1}}$ and $\mathbf{v_{2}}$ are the vectors obtained by concatenating $\mathbf{v_{1i}}$ and $\mathbf{v_{2i}}$ , $i=1,\dots,n_{t}$ , respectively, which correspond to the time steps. The matrices $V_{1i},\,V_{2i}\in\mathbb{R}^{n_{h}\times n_{\xi}}$ are then the matricized counterparts of $\mathbf{v_{1i}}$ and $\mathbf{v_{2i}}$ , respectively. Next, since inverting $\bar{\mathcal{S}}_{\text{exact}}$ is computationally prohibitive, we propose an approximation as
$$
\bar{\mathcal{S}}=\frac{1}{\tau}\underbrace{\left(\mathcal{A}_{t}+\tau\sqrt{\frac{1+\gamma}{\beta}}\mathcal{N}\right)}_{\eqqcolon\bar{\mathcal{Z}}}\left(D\otimes\mathcal{M}_{\gamma}\right)^{-1}\left(\mathcal{A}_{t}+\tau\sqrt{\frac{1+\gamma}{\beta}}\mathcal{N}\right)^{T}, \tag{32}
$$
where using (27), we can rewrite $\bar{\mathcal{Z}}$ as
$$
\bar{\mathcal{Z}}=I_{n_{t}}\otimes\left(H_{1}\otimes\left(M+\tau A_{1}\right)+\tau\sum_{\ell=2}^{n_{A}}H_{\ell}\otimes A_{\ell}\right)+\tau\sqrt{\frac{1+\gamma}{\beta}}I_{n_{t}}\otimes H_{1}\otimes M-C\otimes H_{1}\otimes M. \tag{33}
$$
By dropping the last term $C\otimes H_{1}\otimes M$ , we further approximate $\bar{\mathcal{Z}}$ by
$$
\tilde{\mathcal{Z}}=I_{n_{t}}\otimes\left\{H_{1}\otimes\left[(1+\tau\sqrt{\frac{1+\gamma}{\beta}})M+\tau A_{1}\right]+\tau\sum_{\ell=2}^{n_{A}}H_{\ell}\otimes A_{\ell}\right\}. \tag{34}
$$
Observe that $\tilde{\mathcal{Z}}$ is symmetric, and it can be rewritten as
$$
\tilde{\mathcal{Z}}=\tau\sqrt{\frac{1+\gamma}{\beta}}\mathcal{N}+I_{n_{t}}\otimes\mathcal{L}. \tag{35}
$$
The idea is to use
$$
\tilde{\mathcal{S}}\approx\tilde{\mathcal{Z}}\left(D\otimes\mathcal{M}_{\gamma}\right)^{-1}\tilde{\mathcal{Z}}, \tag{36}
$$
and in particular the solves with $\tilde{\mathcal{Z}}$ are approximated by Algorithm 1, similarly to the steady-state case. We remark that by dropping all terms with $\ell>1$ from $\tilde{\mathcal{Z}}$ in (34); that is, considering
$$
\tilde{\mathcal{Z}}_{0}=I_{n_{t}}\otimes H_{1}\otimes\left[(1+\tau\sqrt{\frac{1+\gamma}{\beta}})M+\tau A_{1}\right],
$$
we recover the mean-based preconditioner [Benner-2016-BDP]. Since the application of Algorithm 1 entails matrix-vector multiplications, using (14) we formulate $\tilde{\mathcal{Z}}\mathbf{v_{3}}$ as
$$
\displaystyle\tilde{\mathcal{Z}}\mathbf{v_{3}} \displaystyle=\left(I_{n_{t}}\otimes\left\{H_{1}\otimes\left[(1+\tau\sqrt{\frac{1+\gamma}{\beta}})M+\tau A_{1}\right]+\tau\sum_{\ell=2}^{n_{A}}H_{\ell}\otimes A_{\ell}\right\}\right)\mathbf{v_{3}} \displaystyle=\begin{bmatrix}\sum_{\ell=1}^{n_{A}}\operatorname{vec}{\left(\hat{A}_{\ell}V_{31}H_{\ell}\right)}&\sum_{\ell=1}^{n_{A}}\operatorname{vec}{\left(\hat{A}_{\ell}V_{32}H_{\ell}\right)}&\cdots&\sum_{\ell=1}^{n_{A}}\operatorname{vec}{\left(\hat{A}_{\ell}V_{3n_{t}}H_{\ell}\right)}\end{bmatrix}^{T}, \tag{37}
$$
where
$$
\hat{A}_{\ell}=\begin{cases}\left(1+\tau\sqrt{\frac{1+\gamma}{\beta}}\right)M+\tau A_{1},&\ell=1,\\
\tau A_{\ell},&\ell=2,\dots,n_{A},\end{cases} \tag{38}
$$
and $V_{3i}\in\mathbb{R}^{n_{h}\times n_{\xi}}$ is the matricized form of the $i$ -th block of $\mathbf{v_{3}}$ , with $i=1,\dots,n_{t}$ .
The practical implementation of the preconditioner for the time-dependent system leverages the inherent structure of the all-at-once formulation. As defined in (34), the core operator of the Schur complement preconditioner, $\tilde{\mathcal{Z}}$ , is block-diagonal with respect to the time steps. This structure extends to the entire KKT system, which then makes the preconditioning easily parallelizable. Specifically, an application of the preconditioner $\mathcal{P}_{\text{hGSoc-PINT}}$ entails an application of the steady-state optimal control preconditioner $\mathcal{P}_{\text{hGSoc}}$ from Algorithm 2 – 3 to all time steps simultaneously, and so it represents parallel-in-time (PINT) approach. It is summarized as Algorithm 4.
Algorithm 4 Parallel-in-time hGSoc preconditioner (hGSoc-PINT)
The preconditioner $\mathcal{P}_{\text{hGSoc-PINT}}:\bar{\mathbf{R}}\longmapsto\bar{\mathbf{V}}$ is defined as:
1: for $k=1,\dots,n_{t}$ do
2: Extract $(\bar{R}^{Y}_{k},\bar{R}^{U}_{k},\bar{R}^{\Lambda}_{k})$ . (the subvector $k$ of $\bar{\mathbf{R}}$ )
3: Calculate $(\bar{V}^{Y}_{k},\bar{V}^{U}_{k},\bar{V}^{\Lambda}_{k})=\mathcal{P}_{hGSoc}(\bar{R}^{Y}_{k},\bar{R}^{U}_{k},\bar{R}^{\Lambda}_{k})$ (apply Algorithm 2 – 3)
4: end for
5: Concatenate $\{(\bar{V}^{Y}_{k},\bar{V}^{U}_{k},\bar{V}^{\Lambda}_{k})\}_{k=1}^{n_{t}}$ into $\bar{\mathbf{V}}$ .
## 4 Spectral analysis of the preconditioners
Since the steady-state optimal control problem can be viewed as a special case of the time-dependent formulation (with $n_{t}=1$ ), we focus on analyzing the time-dependent setting; the steady-state results then follow as a direct consequence. The all-at-once discretization, presented in Section 3, couples all time steps simultaneously, yielding a significantly larger KKT system than its steady-state counterpart. Our goal is to prove that the proposed parallel-in-time preconditioner, based on the hGSoc-PINT in Algorithm 4, is spectrally equivalent to the ideal (but computationally more expensive) preconditioner.
**Definition 4.1 (Spectral Equivalence)**
*Two symmetric positive definite (SPD) matrices $A$ and $B$ are said to be spectrally equivalent, denoted $A\sim B$ , if there exist positive constants $a\leq b$ , such that
$$
a\mathbf{v}^{T}B\mathbf{v}\leq\mathbf{v}^{T}A\mathbf{v}\leq b\mathbf{v}^{T}B\mathbf{v}
$$
holds for all non-zero vectors $\mathbf{v}$ . Equivalently, all eigenvalues of the preconditioned matrix $B^{-1}A$ are contained within the fixed interval, which means $\lambda(B^{-1}A)\subset[a,b]$ .*
The proof proceeds by establishing spectral equivalences
$$
\bar{\mathcal{S}}_{\text{exact}}\sim\bar{\mathcal{S}}\sim\tilde{\mathcal{S}}\sim\tilde{\mathcal{S}}_{r}\sim\tilde{\mathcal{S}}_{\text{hGS-PINT}}. \tag{39}
$$
Here, $\bar{\mathcal{S}}_{\text{exact}}$ denotes the exact Schur complement (31). The matrix $\bar{\mathcal{S}}$ , defined in (32), serves as an approximation of $\bar{\mathcal{S}}_{\text{exact}}$ , while $\tilde{\mathcal{S}}$ given in (36) is a block-diagonal approximation of $\bar{\mathcal{S}}$ . The operator $\tilde{\mathcal{S}}_{r}$ represents the truncated hierarchical operator
$$
\tilde{\mathcal{S}}_{r}\;=\;\tilde{\mathcal{Z}}_{r}\,(D\otimes\mathcal{M}_{\gamma})^{-1}\,\tilde{\mathcal{Z}}_{r}^{T},\,\tilde{\mathcal{Z}}_{r}\;=\;I_{n_{t}}\otimes\left\{H_{1}\otimes\left[(1+\tau\sqrt{\tfrac{1+\gamma}{\beta}})M+\tau A_{1}\right]+\tau\sum_{\ell=2}^{r}H_{\ell}\otimes A_{\ell}\right\}, \tag{40}
$$
with $r=1,\dots,n_{A}$ . When $r=1$ , $\tilde{\mathcal{S}}_{r}$ reduces to the mean-based preconditioner employed in [Benner-2016-BDP], and when $r=n_{A}$ , it recovers the full operator $\tilde{\mathcal{S}}$ . Finally, $\tilde{\mathcal{S}}_{\text{hGS-PINT}}=\tilde{\mathcal{Z}}_{\text{hGS-PINT}}(D\otimes\mathcal{M}_{\gamma})^{-1}\tilde{\mathcal{Z}}_{\text{hGS-PINT}}^{T}$ represents the computationally feasible approximation of $\tilde{\mathcal{S}}_{r}$ in which the linear systems $\tilde{\mathcal{Z}}_{r}\mathbf{x}=\mathbf{b}$ are solved approximately via the hierarchical Gauss-Seidel method (Algorithm 1), as implemented in the parallel-in-time framework of Algorithm 4.
To establish the spectral equivalences in this chain, we require several technical results. We state some auxiliary lemmas concerning matrix perturbations and congruence transformations, which will serve as building blocks for the main theorems. Their proofs are provided in Appendix A. In what follows, we denote by $\sigma_{\min}(\cdot)$ and $\sigma_{\max}(\cdot)$ the smallest and largest singular values, respectively.
**Lemma 4.2 (Spectral Equivalence via Error Bound)**
*Let $A$ and $B$ be symmetric positive definite matrices. Define the error matrix $E=A-B$ . If there exists a constant $0<\delta<1$ such that for all non-zero vectors $\mathbf{v}$ ,
$$
|\mathbf{v}^{T}E\mathbf{v}|\leq\delta(\mathbf{v}^{T}B\mathbf{v}),
$$
then $A$ and $B$ are spectrally equivalent with
$$
1-\delta\leq\lambda(B^{-1}A)\leq 1+\delta.
$$*
**Lemma 4.3 (Eigenvalues under Congruence Transformation)**
*Let $C$ and $D$ be symmetric positive definite matrices, and let $Q$ be a nonsingular matrix. The eigenvalues of the pair $(C,D)$ are identical to those of the transformed pair $(Q^{T}CQ,Q^{T}DQ)$ .*
**Lemma 4.4**
*For any nonzero vector $\mathbf{v}$ , we have $\sigma_{\min}(A)\|\mathbf{v}\|_{2}\leq\|A\mathbf{v}\|_{2}.$*
**Lemma 4.5**
*For any matrices $A$ and $B$ of the same dimensions,
$$
\sigma_{\min}(A+B)\geq\sigma_{\min}(B)-\|A\|_{2}.
$$*
With these auxiliary results in place, we now proceed to establish the spectral equivalences in (39). We begin by recalling the following spectral equivalence result from [Benner-2016-BDP].
**Theorem 4.6 ([Benner-2016-BDP], Theorems 4, 6)**
*Let $\bar{\mathcal{S}}_{\text{exact}}$ be the exact Schur complement and $\bar{\mathcal{S}}$ be its approximation as given by (31) and (32), respectively. The eigenvalues of $\bar{\mathcal{S}}^{-1}\bar{\mathcal{S}}_{\text{exact}}$ are given by
$$
\lambda(\bar{\mathcal{S}}^{-1}\bar{\mathcal{S}}_{\text{exact}})\subseteq\left[\frac{1}{2(1+\alpha)},1\right),
$$
where $\alpha$ satisfies $\alpha<\left(\frac{\sqrt{\kappa(\mathcal{A}_{t})}+1}{\sqrt{\kappa(\mathcal{A}_{t})}-1}\right)^{2}-1,$ and $\mathcal{A}_{t}$ is defined in (27).*
Next, we show the relationship between $\bar{\mathcal{S}}$ and $\tilde{\mathcal{S}}$ , as given by (32) and (36), respectively. First, however, we prove the following lemma.
**Lemma 4.7**
*Let $W=(D\otimes\mathcal{M}_{\gamma})^{-1}$ , and assume there exists a constant $\mu>1$ such that
$$
\tau\sqrt{\frac{1+\gamma}{\beta}}\sigma_{\min}(\mathcal{N}W^{\frac{1}{2}})\geq\mu\|(I_{n_{t}}\otimes\mathcal{L})W^{\frac{1}{2}}\|_{2}, \tag{41}
$$
where $\mathcal{L}$ is defined in (27). Then the minimum eigenvalue of $\tilde{\mathcal{S}}$ has the following lower bound
$$
\lambda_{\min}(\tilde{\mathcal{S}})\geq\frac{\tau}{\beta}\left(1-\frac{1}{\mu}\right)^{2}\sigma^{2}_{\min}(M^{\frac{1}{2}}).
$$*
**Proof 4.8**
*Using the definition of $\tilde{\mathcal{S}}$ in (36), and the fact that for any real matrix $A$ , $\lambda(AA^{T})=\sigma^{2}_{\min}(A)$ , we have
$$
\lambda_{\min}(\tilde{\mathcal{S}})=\lambda_{\min}(\tilde{\mathcal{Z}}W\tilde{\mathcal{Z}}^{T})=\lambda_{\min}((\tilde{\mathcal{Z}}W^{\frac{1}{2}})(\tilde{\mathcal{Z}}W^{\frac{1}{2}})^{T})=\sigma^{2}_{\min}(\tilde{\mathcal{Z}}W^{\frac{1}{2}}).
$$
Recall from (35) that $\tilde{\mathcal{Z}}W^{\frac{1}{2}}=\tau\sqrt{\frac{1+\gamma}{\beta}}\mathcal{N}W^{\frac{1}{2}}+(I_{n_{t}}\otimes\mathcal{L})W^{\frac{1}{2}}$ . Now, applying Lemma 4.5 yields
$$
\sigma_{\min}(\tilde{\mathcal{Z}}W^{\frac{1}{2}})\geq\tau\sqrt{\frac{1+\gamma}{\beta}}\sigma_{\min}(\mathcal{N}W^{\frac{1}{2}})-\|(I_{n_{t}}\otimes\mathcal{L})W^{\frac{1}{2}}\|_{2}. \tag{42}
$$
Now using (42) and the assumption (41), we obtain In our experience, this condition is often satisfied numerically when $\tau\gg\sqrt{\beta}$ .
$$
\sigma_{\min}(\tilde{\mathcal{Z}}W^{\frac{1}{2}})\geq\left(1-\frac{1}{\mu}\right)\tau\sqrt{\frac{1+\gamma}{\beta}}\sigma_{\min}(\mathcal{N}W^{\frac{1}{2}}).
$$
Consequently, the lower bound for $\lambda_{\min}(\tilde{\mathcal{S}})$ is
$$
\lambda_{\min}(\tilde{\mathcal{S}})\geq\left(1-\frac{1}{\mu}\right)^{2}\left(\tau\sqrt{\frac{1+\gamma}{\beta}}\right)^{2}\sigma^{2}_{\min}(\mathcal{N}W^{\frac{1}{2}}).
$$
Next, observe that
$$
\sigma^{2}_{\min}(\mathcal{N}W^{\frac{1}{2}})=\sigma^{2}_{\min}\left((I_{n_{t}}\otimes\mathcal{M})(D^{-\frac{1}{2}}\otimes\mathcal{M}_{\gamma}^{-\frac{1}{2}})\right)=\sigma^{2}_{\min}(D^{-\frac{1}{2}})\cdot\sigma^{2}_{\min}(\mathcal{M}\mathcal{M}_{\gamma}^{-\frac{1}{2}}),
$$
where the last equality follows from the property that the singular values of a Kronecker product are the products of the singular values of the factors, i.e., $\sigma(A\otimes B)=\sigma(A)\sigma(B)$ . It is easy to verify that $\sigma_{\min}(D^{-\frac{1}{2}})=1$ and
$$
\sigma^{2}_{\min}(\mathcal{M}\mathcal{M}_{\gamma}^{-\frac{1}{2}})=\sigma^{2}_{\min}((H_{1}\otimes M)(H_{\gamma}^{-\frac{1}{2}}\otimes M^{-\frac{1}{2}}))=\frac{1}{1+\gamma}\sigma^{2}_{\min}(M^{\frac{1}{2}}).
$$
Thus, we have
$$
\sigma^{2}_{\min}(\mathcal{N}W^{\frac{1}{2}})=\frac{1}{1+\gamma}\sigma^{2}_{\min}(M^{\frac{1}{2}}).
$$
Substituting this into the expression for the lower bound of $\lambda_{\min}(\tilde{\mathcal{S}})$ , we finally get
$$
\lambda_{\min}(\tilde{\mathcal{S}})\geq\left(1-\frac{1}{\mu}\right)^{2}\frac{\tau^{2}(1+\gamma)}{\beta}\left(\frac{1}{1+\gamma}\sigma^{2}_{\min}(M^{\frac{1}{2}})\right)=\frac{\tau}{\beta}\left(1-\frac{1}{\mu}\right)^{2}\sigma^{2}_{\min}(M^{\frac{1}{2}}).
$$*
**Theorem 4.9**
*Assume that the conditions of Lemma 4.7 hold. Then, the eigenvalues $\lambda(\bar{\mathcal{S}}^{-1}\tilde{\mathcal{S}})$ satisfy
$$
(1-\theta)^{2}\leq\lambda(\bar{\mathcal{S}}^{-1}\tilde{\mathcal{S}})\leq(1+\theta)^{2},
$$
where $\theta$ is the perturbation parameter, with
$$
\theta\coloneqq\sup_{\mathbf{v}\neq 0}\frac{\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\mathbf{v}\|_{2}}{\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}}\leq\frac{\sqrt{2\beta}}{\sqrt{1+\gamma}\left(1-\frac{1}{\mu}\right)\tau}\kappa(M^{\frac{1}{2}}).
$$*
**Proof 4.10**
*Recall from (33) and (34), $\bar{\mathcal{Z}}=\tilde{\mathcal{Z}}-C\otimes\mathcal{M}$ , which implies for any $\mathbf{v}\neq 0$ ,
$$
W^{\frac{1}{2}}\bar{\mathcal{Z}}^{T}\mathbf{v}=W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}-W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\mathbf{v}.
$$ Using the triangle inequality, we have
$$
\left|\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}-\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\mathbf{v}\|_{2}\right|\leq\|W^{\frac{1}{2}}\bar{\mathcal{Z}}^{T}\mathbf{v}\|_{2}\leq\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}+\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\mathbf{v}\|_{2}.
$$
From the definition of $\theta$ , we have $\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\mathbf{v}\|_{2}\leq\theta\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}$ . Substituting this into the inequality above, we obtain
$$
(1-\theta)\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}\leq\|W^{\frac{1}{2}}\bar{\mathcal{Z}}^{T}\mathbf{v}\|_{2}\leq(1+\theta)\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}.
$$
Squaring these inequalities leads to the spectral bounds
$$
(1-\theta)^{2}\leq\lambda(\bar{\mathcal{S}}^{-1}\tilde{\mathcal{S}})\leq(1+\theta)^{2}.
$$
Next, it remains to analyze the upper bound of $\theta$ . To this end, using Lemma 4.4, we know that $\sigma_{\min}(W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T})\|\mathbf{v}\|_{2}\leq\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}$ .
$$
\theta=\sup_{\mathbf{v}\neq 0}\frac{\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\mathbf{v}\|_{2}}{\|W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T}\mathbf{v}\|_{2}}\leq\frac{\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\|_{2}\|\mathbf{v}\|_{2}}{\sigma_{\min}(W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T})\|\mathbf{v}\|_{2}}=\frac{\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\|_{2}}{\sigma_{\min}(W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T})},
$$
For the numerator, we have
| | $\displaystyle\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\|_{2}$ | $\displaystyle=\|(D^{-\frac{1}{2}}\otimes\mathcal{M}_{\gamma}^{-\frac{1}{2}})(C^{T}\otimes\mathcal{M}^{T})\|_{2}$ | |
| --- | --- | --- | --- |
Furthermore, the spectral norm of the Kronecker product satisfies
$$
\|D^{-\frac{1}{2}}C^{T}\|_{2}\cdot\|\mathcal{M}_{\gamma}^{-\frac{1}{2}}\mathcal{M}^{T}\|_{2}\leq\frac{\sqrt{2}}{\sqrt{1+\gamma}}\|M^{\frac{1}{2}}\|_{2},
$$
which follows from the bound $\|D^{-1/2}C^{T}\|_{2}\leq 1$ and the fact that $\|\mathcal{M}_{\gamma}^{-1/2}\mathcal{M}^{T}\|_{2}$ is bounded by the scaled norm of the square root of the mass matrix. Thus, the numerator satisfies
$$
\|W^{\frac{1}{2}}(C\otimes\mathcal{M})^{T}\|_{2}\leq\frac{\sqrt{2}}{\sqrt{1+\gamma}}\|M^{\frac{1}{2}}\|_{2}.
$$
For the denominator, using lemma 4.7, we have
$$
\sigma_{\min}(W^{\frac{1}{2}}\tilde{\mathcal{Z}}^{T})=\sigma_{\min}(\tilde{\mathcal{Z}}W^{\frac{1}{2}})\geq\left(1-\frac{1}{\mu}\right)\frac{\tau}{\sqrt{\beta}}\sigma_{\min}(M^{\frac{1}{2}}).
$$
Combining the bounds for the numerator and denominator,
| | $\displaystyle\theta$ | $\displaystyle\leq\frac{\frac{\sqrt{2}}{\sqrt{1+\gamma}}\|M^{\frac{1}{2}}\|_{2}}{\left(1-\frac{1}{\mu}\right)\frac{\tau}{\sqrt{\beta}}\sigma_{\min}(M^{\frac{1}{2}})}$ | |
| --- | --- | --- | --- |
thereby completing the proof of the theorem.*
Note that the bound on $\theta$ shows that for sufficiently large $\tau$ (e.g., $\tau\gtrsim\sqrt{\beta}\,\kappa(M)^{1/2}$ ), we have $\theta<1$ ; hence, $\bar{\mathcal{S}}$ and $\tilde{\mathcal{S}}$ are spectrally equivalent. In practice, $\kappa(M)$ can be kept $\mathcal{O}(1)$ via appropriate basis choices and mass lumping, so it suffices to require $\tau\gg\sqrt{\beta}$ .
Next, we establish the spectral equivalence between $\tilde{\mathcal{S}}$ and its truncated form $\tilde{\mathcal{S}}_{r}$ . To this end, we shall rely on the error matrix $\mathcal{E}_{r}=\tilde{\mathcal{S}}-\tilde{\mathcal{S}}_{r}$ , and the result of Lemma 4.2.
**Theorem 4.11 (Truncation Error Bound)**
*Let $\tilde{\mathcal{S}}$ be the full operator and $\tilde{\mathcal{S}}_{r}$ be the truncated hierarchical operator with truncation parameter $r$ defined in (36) and (40), respectively. Assume $\nu\coloneqq\lambda_{\min}(\hat{A}_{1})>0$ . Define the truncation tail sum as
$$
\delta_{r}\coloneqq\frac{1}{\nu}\sum_{\ell=r+1}^{n_{A}}\|H_{\ell}\|_{2}\,\|\hat{A}_{\ell}\|_{2}.
$$
If the decay of the stochastic expansion coefficients is such that $0<\delta_{r}<1$ , then $\tilde{\mathcal{S}}$ and $\tilde{\mathcal{S}}_{r}$ are spectrally equivalent
$$
1-\delta_{r}\leq\lambda(\tilde{\mathcal{S}}_{r}^{-1}\tilde{\mathcal{S}})\leq 1+\delta_{r}.
$$*
**Proof 4.12**
*Consider the error matrix $\mathcal{E}_{r}=\tilde{\mathcal{S}}-\tilde{\mathcal{S}}_{r}$ . Based on the expansion structure of the operator defined in (34), the error term consists of the neglected high-order terms in the stochastic expansion:
$$
\mathcal{E}_{r}=\sum_{\ell=r+1}^{n_{A}}H_{\ell}\otimes\hat{A}_{\ell},
$$
where $\hat{A}_{\ell}$ represents the associated spatial coefficient matrices. We apply Lemma 4.2 by estimating the ratio $|\mathbf{v}^{T}\mathcal{E}_{r}\mathbf{v}|/|\mathbf{v}^{T}\tilde{\mathcal{S}}_{r}\mathbf{v}|$ . First, we bound the numerator using the triangle inequality and the submultiplicativity of the spectral norm:
$$
\left|\mathbf{v}^{T}\mathcal{E}_{r}\mathbf{v}\right|=\left|\sum_{\ell=r+1}^{n_{A}}\mathbf{v}^{T}(H_{\ell}\otimes\hat{A}_{\ell})\mathbf{v}\right|\leq\sum_{\ell=r+1}^{n_{A}}\left|\mathbf{v}^{T}(H_{\ell}\otimes\hat{A}_{\ell})\mathbf{v}\right|\leq\sum_{\ell=r+1}^{n_{A}}\|H_{\ell}\|_{2}\,\|\hat{A}_{\ell}\|_{2}\,\|\mathbf{v}\|^{2}.
$$ Next, we establish a lower bound for the denominator
$$
\mathbf{v}^{T}\tilde{\mathcal{S}}_{r}\mathbf{v}\geq\mathbf{v}^{T}(H_{1}\otimes\hat{A}_{1})\mathbf{v}\geq\nu\|\mathbf{v}\|^{2}.
$$
Combining these inequalities yields
$$
\frac{|\mathbf{v}^{T}\mathcal{E}_{r}\mathbf{v}|}{|\mathbf{v}^{T}\tilde{\mathcal{S}}_{r}\mathbf{v}|}\leq\frac{\sum_{\ell=r+1}^{n_{A}}\|H_{\ell}\|_{2}\,\|\hat{A}_{\ell}\|_{2}\,\|\mathbf{v}\|^{2}}{\nu\|\mathbf{v}\|^{2}}=\frac{1}{\nu}\sum_{\ell=r+1}^{n_{A}}\|H_{\ell}\|_{2}\,\|\hat{A}_{\ell}\|_{2}=\delta_{r}.
$$
Provided that the Karhunen-Loève expansion coefficients decay sufficiently fast (e.g., algebraically or exponentially), the tail sum approaches zero as $r\to n_{A}$ . Thus, for sufficiently large $r$ , we have $\delta_{r}<1$ . The result then follows directly from Lemma 4.2.*
Finally, we proceed to establish the spectral equivalence between the truncated preconditioner $\tilde{\mathcal{S}}_{r}$ and the hierarchical symmetric block Gauss-Seidel approximation $\tilde{\mathcal{S}}_{\text{hGS-PINT}}$ . To this end, following the idea from [Bespalov-2021-TPS], we can rewrite it as $H_{\ell}=L_{\ell}+L_{\ell}^{T},\,\ell=2,3,\dots,n_{A}$ , and matrices $L_{\ell}$ have at most one nonzero entry per row and per column.
Now, define $X_{1}=I_{n_{t}}\otimes H_{1}\otimes\hat{A}_{1}$ , $X_{r}=I_{n_{t}}\otimes\left(\sum_{\ell=2}^{r}L_{\ell}\otimes\hat{A}_{\ell}\right),\,r=2,3,\dots n_{A}$ , so we know $\tilde{\mathcal{Z}}_{\text{hGS-PINT}}=(X_{1}+X_{r})X_{1}^{-1}(X_{1}+X_{r}^{T})=\tilde{\mathcal{Z}}_{r}+X_{r}X_{1}^{-1}X_{r}^{T}$ , the Rayleigh quotient
$$
\frac{\mathbf{v}^{T}\tilde{\mathcal{Z}}_{\text{hGS-PINT}}\mathbf{v}}{\mathbf{v}^{T}\tilde{\mathcal{Z}}_{r}\mathbf{v}}=1+\underbrace{\frac{\mathbf{v}^{T}X_{r}X_{1}^{-1}X_{r}^{T}\mathbf{v}}{\mathbf{v}^{T}(X_{1}+X_{r}+X_{r}^{T})\mathbf{v}}}_{\zeta}.
$$
Introducing the change of variable $\mathbf{u}=X_{1}^{\frac{1}{2}}\mathbf{v}$ and setting $Y=X_{1}^{-\frac{1}{2}}X_{r}X_{1}^{-\frac{1}{2}}$ , we have
$$
\zeta(\mathbf{u})=\frac{\mathbf{u}^{T}YY^{T}\mathbf{u}}{\mathbf{u}^{T}(I+Y+Y^{T})\mathbf{u}},
$$
where $Y=I_{n_{t}}\otimes\left(\sum_{\ell=2}^{r}L_{\ell}\otimes\mathfrak{A}_{\ell}\right)$ , with the scaled coefficient matrices defined as $\mathfrak{A}_{\ell}=\hat{A}_{1}^{-\frac{1}{2}}\hat{A}_{\ell}\hat{A}_{1}^{-\frac{1}{2}}.$ Consequently,
$$
\zeta(\mathbf{u})\leq\max_{\mathbf{u}\neq 0}\frac{\mathbf{u}^{T}YY^{T}\mathbf{u}}{\mathbf{u}^{T}\mathbf{u}}\cdot\max_{\mathbf{u}\neq 0}\frac{\mathbf{u}^{T}\mathbf{u}}{\mathbf{u}^{T}(I+Y+Y^{T})\mathbf{u}}=\frac{\sigma_{\max}^{2}(Y)}{\lambda_{\min}\left(I+Y+Y^{T}\right)}.
$$
The following result holds.
**Lemma 4.13**
*Define
$$
\Delta_{r}\ :=\ \sum_{\ell=2}^{r}\|H_{\ell}\|_{2}\,\rho_{\ell},\qquad\rho_{\ell}\ :=\ \|\mathfrak{A}_{\ell}\|_{2},
$$
and let
$$
Y:=X_{1}^{-1/2}X_{r}X_{1}^{-1/2}=I_{n_{t}}\otimes\Big(\sum_{\ell=2}^{r}L_{\ell}\otimes\mathfrak{A}_{\ell}\Big),\quad E:=Y+Y^{T}=I_{n_{t}}\otimes\Big(\sum_{\ell=2}^{r}H_{\ell}\otimes\mathfrak{A}_{\ell}\Big),
$$
with $H_{\ell}=L_{\ell}+L_{\ell}^{T}$ . Then $\lambda_{\min}\big(I+Y+Y^{T}\big)\ \geq\ 1-\Delta_{r}$ , $\rho_{\ell}\ \leq\ \frac{\|{\color[rgb]{0,0,0}\kappa_{\ell}(x)}\|_{L^{\infty}(\Omega)}}{\Bbbk_{\min}}$ , and
$$
\sigma_{\max}(Y)\ \leq\ \sum_{\ell=2}^{r}\|L_{\ell}\|_{2}\,\rho_{\ell}\ \leq\ \Delta_{r}.
$$*
**Proof 4.14**
*Using
$$
I+Y+Y^{T}=X_{1}^{-1/2}\Big(X_{1}+I_{n_{t}}\otimes\sum_{\ell=2}^{r}H_{\ell}\otimes\hat{A}_{\ell}\Big)X_{1}^{-1/2}=X_{1}^{-1/2}\,\bar{\mathcal{Z}}_{r}\,X_{1}^{-1/2},
$$
the spectrum of $I+Y+Y^{T}$ coincides with the generalized spectrum of the pair $(\bar{\mathcal{Z}}_{r},X_{1})$ . Since $E$ is symmetric, its spectral norm satisfies $\|E\|_{2}=|\lambda_{\max}(E)|\geq|\lambda_{i}(E)|$ for any $i$ . In particular, $\lambda_{\min}(E)\geq-\|E\|_{2}$ , and therefore $\lambda_{\min}(I+E)={\color[rgb]{0,0,0}1+\lambda_{\min}(E)}\geq 1-\|E\|_{2}$ . By the Kronecker product norm rule and the triangle inequality,
$$
\|E\|_{2}=\Big\|I_{n_{t}}\otimes(\sum_{\ell=2}^{r}H_{\ell}\otimes\mathfrak{A}_{\ell})\Big\|_{2}\leq\sum_{\ell=2}^{r}\|H_{\ell}\|_{2}\,\|\mathfrak{A}_{\ell}\|_{2}=\sum_{\ell=2}^{r}\|H_{\ell}\|_{2}\,\rho_{\ell}=\Delta_{r},
$$
which yields $\lambda_{\min}(I+Y+Y^{T})\geq 1-\Delta_{r}$ . Similarly,
$$
\|Y\|_{2}=\Big\|I_{n_{t}}\otimes\sum_{\ell=2}^{r}L_{\ell}\otimes\mathfrak{A}_{\ell}\Big\|_{2}\leq\sum_{\ell=2}^{r}\|L_{\ell}\|_{2}\,\rho_{\ell}.
$$
Since $H_{\ell}=L_{\ell}+L_{\ell}^{T}$ and each $L_{\ell}$ has at most one nonzero per row and per column, we have $\|L_{\ell}\|_{2}\leq\|H_{\ell}\|_{2}$ . Hence $\sigma_{\max}(Y)=\|Y\|_{2}\leq\Delta_{r}$ . For the explicit bound on $\rho_{\ell}$ , recall that $\rho_{\ell}=\|\mathfrak{A}_{\ell}\|_{2}$ is the maximum eigenvalue of the generalized eigenvalue problem $A_{\ell}\mathbf{v}=\lambda\hat{A}_{1}\mathbf{v}$ , where $\hat{A}_{1}$ defined in (38). In terms of the associated finite element function $v_{h}=\sum_{j=1}^{n_{h}}v_{j}\phi_{j}$ , the Rayleigh quotient is given by
$$
\frac{\mathbf{v}^{T}A_{\ell}\mathbf{v}}{\mathbf{v}^{T}\hat{A}_{1}\mathbf{v}}=\frac{\displaystyle\int_{\Omega}\kappa_{\ell}\,|\nabla v_{h}|^{2}\,dx}{\displaystyle\int_{\Omega}\left((1+\tau\sqrt{\tfrac{1+\gamma}{\beta}})|v_{h}|^{2}+\tau\kappa_{0}\,|\nabla v_{h}|^{2}\right)\,dx}.
$$
Since the mass term is non-negative, we can bound this ratio by neglecting the $L^{2}$ -term in the denominator:
$$
\frac{\mathbf{v}^{T}A_{\ell}\mathbf{v}}{\mathbf{v}^{T}\hat{A}_{1}\mathbf{v}}\ \leq\ \frac{\displaystyle\int_{\Omega}\kappa_{\ell}\,|\nabla v_{h}|^{2}\,dx}{\displaystyle\tau\int_{\Omega}\kappa_{0}\,|\nabla v_{h}|^{2}\,dx}\ \leq\ \frac{1}{\tau}\frac{\|{\color[rgb]{0,0,0}\kappa_{\ell}(x)}\|_{L^{\infty}(\Omega)}}{\Bbbk_{\min}},
$$
and taking the supremum gives the stated bound on $\rho_{\ell}$ ; substituting it into $\|Y\|_{2}$ yields the last inequality.*
Observe from above that, with $\Delta_{r}<1$ ,
$$
\zeta(\mathbf{v})\ =\ \frac{\mathbf{v}^{T}YY^{T}\mathbf{v}}{\mathbf{v}^{T}(I+Y+Y^{T})\mathbf{v}}\ \leq\ \frac{\sigma_{\max}^{2}(Y)}{\lambda_{\min}(I+Y+Y^{T})}\ \leq\ \frac{\Delta_{r}^{\,2}}{\,1-\Delta_{r}\,}.
$$
Consequently, we have the bound for the Rayleigh quotient of the factors:
$$
\frac{\mathbf{v}^{T}(\tilde{\mathcal{Z}}_{\text{hGS-PINT}}-\tilde{\mathcal{Z}}_{r})\mathbf{v}}{\mathbf{v}^{T}\tilde{\mathcal{Z}}_{r}\mathbf{v}}=\zeta(\mathbf{v})\leq\frac{\Delta_{r}^{2}}{1-\Delta_{r}}.
$$
To derive the bound for the preconditioner $\tilde{\mathcal{S}}$ , we utilize Lemma 4.3. Let $Q=(D\otimes\mathcal{M})^{1/2}$ . The generalized eigenvalues of $(\tilde{\mathcal{S}}_{\text{hGS-PINT}},\tilde{\mathcal{S}}_{r})$ are identical to those of the pair
$$
\left((Q\tilde{\mathcal{Z}}_{\text{hGS-PINT}}Q)^{2},\quad(Q\tilde{\mathcal{Z}}_{r}Q)^{2}\right).
$$
Since the eigenvalues of the squared operators are simply the squares of the eigenvalues of the base operators (for these symmetric positive definite factors), the spectral bound for the preconditioner is the square of the bound for the factors. Thus, we conclude:
$$
1\leq\lambda(\tilde{\mathcal{S}}_{r}^{-1}\tilde{\mathcal{S}}_{\text{hGS-PINT}})\leq\left(1+\frac{\Delta_{r}^{\,2}}{1-\Delta_{r}}\right)^{2}.
$$
**Remark 4.15 (Summary of Spectral Equivalence for Time-Dependent Case)**
*Under the assumptions of Theorem 4.6, Lemma 4.7, Theorem 4.11, and Lemma 4.13, the hierarchical Gauss-Seidel PINT preconditioner $\tilde{\mathcal{S}}_{\text{hGS-PINT}}$ is spectrally equivalent to the exact Schur complement $\bar{\mathcal{S}}_{\text{exact}}$ . The spectral equivalence is established through the following chain of approximations:
$$
\bar{\mathcal{S}}_{\text{exact}}\sim\bar{\mathcal{S}}\sim\tilde{\mathcal{S}}\sim\tilde{\mathcal{S}}_{r}\sim\tilde{\mathcal{S}}_{\text{hGS-PINT}}.
$$
In particular, provided that the truncation rank $r$ is sufficiently large (to satisfy $\delta_{r}<1$ ) and the mass matrix $M$ is well-conditioned (e.g., via mass lumping), the eigenvalues of the preconditioned system $\tilde{\mathcal{S}}_{\text{hGS-PINT}}^{-1}\bar{\mathcal{S}}_{\text{exact}}$ are contained in a fixed interval independent of the spatial mesh size $h$ and the stochastic discretization parameters $m_{\xi}$ and $p$ (assuming sufficient decay of the expansion coefficients), ensuring mesh robustness in the fixed time step regime.*
**Remark 4.16 (Steady-State Case)**
*Since the steady-state optimal control problem corresponds to the special case $n_{t}=1$ of the time-dependent formulation, all preceding results apply directly with the simplified notation. The spectral equivalence chain for the steady-state Schur complement preconditioner is
$$
\mathcal{S}_{\text{exact}}\sim\mathcal{S}\sim\mathcal{S}_{r}\sim\mathcal{S}_{\text{hGS}},
$$
where $\mathcal{S}_{\text{exact}}$ is defined in (12), $\mathcal{S}$ in (13), and
$$
\mathcal{S}_{r}=\mathcal{Z}_{r}\mathcal{M}_{\gamma}^{-1}\mathcal{Z}_{r}^{T},\qquad\mathcal{Z}_{r}=\sum_{\ell=1}^{r}H_{\ell}\otimes\tilde{A}_{\ell},
$$
with $\tilde{A}_{1}=A_{1}+\sqrt{\frac{1+\gamma}{\beta}}M$ and $\tilde{A}_{\ell}=A_{\ell}$ for $\ell=2,\dots,r$ . When $r=1$ , $\mathcal{S}_{r}$ reduces to the mean-based preconditioner, and when $r=n_{A}$ , it recovers the full operator $\mathcal{S}$ . Finally, $\mathcal{S}_{\text{hGS}}$ represents the computationally feasible approximation of $\mathcal{S}_{r}$ in which the linear systems $\mathcal{Z}_{r}\mathbf{x}=\mathbf{b}$ are solved approximately via the hierarchical Gauss-Seidel method (Algorithm 1), as implemented in Algorithm 2 – 3.*
## 5 Numerical experiments
This section validates the theoretical findings of Sections 2.1 – 3 through comprehensive numerical experiments. We pursue two primary objectives: (i) verifying the mesh-robust for fixed time step and spectral bounds established in the preceding sections, and (ii) demonstrating the computational efficiency of the proposed hierarchical Gauss-Seidel (hGS) preconditioner across varying truncation strategies. Experiments are presented for both steady-state problems (Section 5.1) and time-dependent problems (Section 5.2). The numerical experiments were performed on a system running AlmaLinux-9 with 40GB RAM, and the proposed algorithms were implemented using Matlab 23.2.
The random coefficients $\Bbbk(x,\xi)$ in the problem (2) are constructed as a finite expansion (4), where the spatial modes $\kappa_{i}(x)$ and weights $\theta_{i}$ are eigenpairs of the covariance function
$$
C_{f}(x,y)=\sigma_{\Bbbk}^{2}\exp\left(-\frac{\left|x_{1}-y_{1}\right|}{\ell_{1}}-\frac{\left|x_{2}-y_{2}\right|}{\ell_{2}}\right)\quad\forall(x,y)\in[-1,1]^{2}.
$$
We set the correlation lengths as $\ell_{1}=\ell_{2}=1$ , the mean $\kappa_{1}(x)\equiv 1$ . For the gPC setting in (5), we consider a log-normal random field parameterized by independent Gaussian random variables, for which corresponding Hermite polynomials are employed as the basis. The total number of basis functions is given by $n_{\xi}=\binom{m_{\xi}+p}{m_{\xi}}$ , where $m_{\xi}$ is the number of random variables and $p$ is the polynomial degree. For instance, the case $(m_{\xi},p)=(3,3)$ yields $n_{\xi}=\binom{6}{3}=20$ . This problem has been extensively studied in [Powell-2009-BDP]. Also, we used $\gamma=1$ in both cases, which means we only consider the case with standard deviation. To discretize the spatial domain, we implemented our code based on IFISS 3.7 [ifiss], using $\mathbf{Q}_{1}$ approximation. For temporal discretization, we apply the all-at-once technique proposed in [rees2010all] and set the terminal time as $T=1$ . In all numerical experiments, the spatial mesh size $h$ and the time step $\tau$ are chosen as $2^{-i}$ , with $i=4,5,6,7$ . We solve the linear systems (8) and (29) using the preconditioners given by (11) and (30), respectively, employing the flexible GMRES method (without restarting) [Saad-1993-FIP]. The stopping criterion is defined in terms of the relative residual $\|r_{k}\|/\|b\|$ , with thresholds $10^{-8}$ for the steady-state experiments and $10^{-6}$ or $10^{-4}$ for the time-dependent runs, where $r_{k}$ denotes the residual at iteration $k$ and $b$ is the right-hand side vector. For notational consistency, we set $n_{\tau}\equiv r$ throughout the remainder of the paper. To assess the effectiveness of the hierarchical preconditioning strategy, we systematically compare three truncation settings for the (3,3)-block preconditioner: $r=1$ (mean-based approximation), $r=m_{\xi}+1$ (hGS truncated at the first-order stochastic terms), and $r=n_{A}$ (full expansion retaining all $h_{ijk}$ coefficients). For each configuration, we evaluate the efficiency of the hierarchical Gauss-Seidel method by comparing both the iteration counts and the computational times (in seconds).
We consider homogeneous Dirichlet conditions (i.e., $g(x)=0$ ), corresponding to Example 2 in [elman_finite_2014, Chapter 5]. This example is defined on a square domain $\mathcal{D}=[-1,1]^{2}$ with a discontinuous target function
$$
y_{d}=\begin{cases}1&\text{ in }\Omega_{1}\coloneqq[-1,0]^{2},\\
0&\text{ in }\mathcal{D}\backslash\Omega_{1},\end{cases} \tag{43}
$$
which represents inconsistent Dirichlet boundary data since the target state $y_{d}=1$ differs from the required boundary condition $y=0$ on $\partial\mathcal{D}\cap\partial\Omega_{1}$ .
We subsequently present numerical experiments for both steady-state and time dependent problems to illustrate and verify the efficiency of our proposed hGS method.
Here are some details about the implementation for preconditioners (11), (30). Since time-dependent problems can be seen as a series of steady-state problems, and also because of the diagonal structure of matrix $D$ and matrix $I_{n_{t}}$ , we can just focus on the steady-state preconditioner.
The practical implementation of the preconditioner $\mathcal{P}$ in (11) involves different strategies for its constituent blocks. For the (1,1) and (2,2) blocks, which are based on Kronecker products involving the mass matrix $M$ , applying their inverses requires solving linear systems with $M$ . These solves are handled efficiently by either a direct Cholesky decomposition or the iterative Chebyshev semi-iteration method [golub1961chebyshev, MR2597805].
For the more complex (3,3) block, which represents the approximate Schur complement $\mathcal{S}$ , we employ an outer iterative scheme. Specifically, we use the Preconditioned Richardson method, outlined in Algorithm 5 [MR1386889, Chapter 7], where the core of our proposalâ the hierarchical Gauss-Seidel (hGS) method from Algorithm 1 serves as the preconditioner for each Richardson step.
Algorithm 5 Preconditioned Richardson iteration with hGS preconditioner
1: Given matrix $\mathcal{Z}$ , vector $b$ , and initial guess $x^{(1)}$ .
2: $r^{(1)}=b-\mathcal{Z}x^{(1)}$ . (initial residual)
3: for $k=1,2,\ldots,N$ do
4: Solve $\mathcal{Z}z^{(k)}=r^{(k)}$ . (apply Algorithm 1)
5: $x^{(k+1)}=x^{(k)}+z^{(k)}$ (update solution)
6: $r^{(k+1)}=b-\mathcal{Z}x^{(k+1)}$ . (update residual)
7: end for
8: Return $x^{(N+1)}$
### 5.1 Steady-state case
This subsection focuses on the steady-state optimal control problem (1). We examine the performance of the proposed preconditioner (11) under systematic variations in: (i) the variance $\sigma_{\Bbbk}$ (Tables 1 – 3), (ii) the regularization parameter $\beta$ (Table 4), and (iii) the spatial and stochastic discretization levels. For each configuration, we compare three solvers for the (1,1) and (2,2) blocks—Chebyshev semi-iteration with 5 or 10 steps and direct Cholesky factorization—combined with the three truncation strategies for the (3,3)-block described above. For the (3,3)-block, we apply the Richardson iteration (Algorithm 5) with $N=1$ , i.e., one application of the hGS preconditioner per outer GMRES iteration. All tests use a fixed tolerance of $10^{-8}$ .
Next, by fixing the parameter $\beta$ , we perform further tests summarized in Tables 1 – 3, employing different step settings for the Chebyshev smoother and the Cholesky decomposition for blocks (1,1) and (2,2), as well as various truncation strategies for block (3,3) within the hGS method. The numerical experiments were conducted using multiple mesh sizes and stochastic parameter configurations, with a solver tolerance set to $10^{-8}$ .
Tables 1 – 3 present results for $\beta=10^{-4}$ with $\sigma_{\Bbbk}\in\{0.01,0.1,0.4\}$ , covering a range from near-deterministic to highly stochastic regimes. Tables 1 – 3 demonstrate three key theoretical properties. First, regarding mesh independence, iteration counts grow sub-linearly with spatial refinement for fixed stochastic dimension $n_{\xi}$ , consistent with the spectral bounds established in Section 2.1. Second, concerning truncation efficiency, the $n_{\tau}=m_{\xi}+1$ strategy achieves iteration counts comparable to the full expansion ( $n_{\tau}=n_{A}$ ) while avoiding the computational overhead of summing over all $h_{ijk}$ coefficients, thereby validating the hierarchical approximation framework. Third, regarding smoother comparison, the 5-step Chebyshev semi-iteration balances convergence rate and per-iteration cost more effectively than either the 10-step variant or direct Cholesky factorization. Across all configurations, the $n_{\tau}=m_{\xi}+1$ truncation consistently delivers performance intermediate between the mean-based approximation ( $n_{\tau}=1$ ) and the full expansion, confirming the practical value of the proposed hierarchical preconditioning strategy.
Table 1: Results showing the total number of iterations from low-rank preconditioned GMRES and the total CPU times (in seconds) using preconditioner with $\beta=10^{-4}$ , $\sigma=0.01$ , and selected spatial ( $n_{h}$ ) and stochastic ( $n_{\xi}$ ) degrees of freedom
| $n_{\xi}$ $n_{h}$ $n_{\tau}$ | # iter(t) 289( $h=\frac{1}{2^{4}}$ ) 1 | # iter(t) 1089( $h=\frac{1}{2^{5}}$ ) $m_{\xi}$ +1 | # iter(t) 4225( $h=\frac{1}{2^{6}}$ ) $n_{A}$ | # iter(t) 16641( $h=\frac{1}{2^{7}}$ ) 1 | $m_{\xi}$ +1 | $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\sigma_{a}=0.01$ | | | | | | | | | | | | |
| Chebyshev-5+hGS-1 | | | | | | | | | | | | |
| 20 | 32(5.7) | 32(4.0) | 32(5.0) | 36(9.2) | 35(10.0) | 35(15.6) | 36(23.2) | 27(26.9) | 27(46.5) | 35(173.4) | 33(144.9) | 27(194.8) |
| 70 | 32(129.1) | 32(129.5) | 32(182.9) | 36(243.4) | 35(239.1) | 35(409.2) | 36(752.3) | 34(724.9) | 34(1237.8) | 35(2585.4) | 33(2524.6) | 33(5253.9) |
| 84 | 32(317.5) | 32(318.3) | 32(403.0) | 36(508.9) | 35(498.2) | 35(601.3) | 36(1059.3) | 34(864.8) | 34(1632.5) | 35(3099.0) | 33(3798.1) | 33(5630.2) |
| Chebyshev-10+hGS-1 | | | | | | | | | | | | |
| 20 | 26(2.8) | 26(2.6) | 26(4.4) | 30(6.9) | 29(7.1) | 29(12.0) | 31(20.2) | 31(21.1) | 31(37.6) | 32(145.1) | 31(130.3) | 31(154.5) |
| 70 | 26(111.7) | 26(112.2) | 26(153.7) | 30(191.6) | 29(189.0) | 29(323.8) | 32(550.7) | 31(383.1) | 31(830.2) | 33(1745.9) | 31(2273.0) | 31(3341.3) |
| 84 | 26(263.9) | 26(266.1) | 26(329.2) | 30(457.7) | 29(432.2) | 29(504.1) | 32(959.1) | 31(800.4) | 31(1515.3) | 33(3015.1) | 31(3643.9) | 31(5318.7) |
| Cholesky+hGS-1 | | | | | | | | | | | | |
| 20 | 25(3.5) | 25(4.0) | 25(4.3) | 29(7.4) | 27(6.7) | 27(10.9) | 29(19.1) | 29(26.5) | 29(44.5) | 31(143.6) | 29(162.7) | 29(182.9) |
| 70 | 25(104.2) | 25(107.6) | 25(144.3) | 29(177.3) | 27(168.9) | 27(283.9) | 29(419.0) | 29(480.8) | 29(937.1) | 31(2390.7) | 29(2046.0) | 29(3299.9) |
| 84 | 25(245.4) | 25(247.4) | 25(316.8) | 29(429.9) | 27(393.6) | 27(675.5) | 29(831.9) | 29(1008.3) | 29(1756.8) | 31(3563.3) | 29(3393.6) | 29(6436.3) |
Table 2: Results showing the total number of iterations from low-rank preconditioned GMRES and the total CPU times (in seconds) using preconditioner with $\beta=10^{-4}$ , $\sigma=0.1$ , and selected spatial ( $n_{h}$ ) and stochastic ( $n_{\xi}$ ) degrees of freedom
| $n_{\xi}$ $n_{h}$ $n_{\tau}$ | # iter(t) 289( $h=\frac{1}{2^{4}}$ ) 1 | # iter(t) 1089( $h=\frac{1}{2^{5}}$ ) $m_{\xi}$ +1 | # iter(t) 4225( $h=\frac{1}{2^{6}}$ ) $n_{A}$ | # iter(t) 16641( $h=\frac{1}{2^{7}}$ ) 1 | $m_{\xi}$ +1 | $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\sigma_{a}=0.1$ | | | | | | | | | | | | |
| Chebyshev-5+hGS-1 | | | | | | | | | | | | |
| 20 | 39(3.6) | 33(3.3) | 33(3.7) | 35(11.1) | 35(8.7) | 29(14.2) | 34(29.7) | 36(40.6) | 35(42.3) | 34(162.3) | 28(84.0) | 34(191.2) |
| 70 | 41(199.9) | 34(173.7) | 33(136.4) | 45(218.5) | 35(190.5) | 35(299.8) | 46(841.0) | 36(670.3) | 35(903.4) | 45(2093.8) | 35(2917.1) | 34(4099.5) |
| 84 | 41(460.9) | 34(238.3) | 33(325.3) | 45(475.0) | 35(529.5) | 35(610.9) | 44(1332.7) | 36(1028.1) | 35(1708.7) | 44(4940.4) | 35(4008.7) | 34(5843.5) |
| Chebyshev-10+hGS-1 | | | | | | | | | | | | |
| 20 | 34(2.8) | 26(2.2) | 26(3.1) | 38(7.1) | 30(7.0) | 30(10.5) | 40(25.7) | 31(32.8) | 31(37.8) | 42(108.7) | 31(86.2) | 31(152.0) |
| 70 | 36(191.4) | 26(139.4) | 26(114.9) | 40(186.6) | 30(234.7) | 30(264.6) | 42(488.8) | 31(518.1) | 31(894.2) | 44(1746.1) | 32(1518.3) | 31(3608.8) |
| 84 | 35(411.8) | 26(193.3) | 26(283.1) | 40(436.9) | 30(354.9) | 30(558.3) | 42(1221.2) | 31(970.6) | 31(1529.2) | 43(4988.0) | 32(3801.3) | 31(5417.9) |
| Cholesky+hGS-1 | | | | | | | | | | | | |
| 20 | 33(2.6) | 25(2.2) | 25(3.2) | 37(7.9) | 27(6.9) | 27(11.1) | 39(28.4) | 29(34.0) | 29(34.8) | 39(129.7) | 29(103.2) | 29(150.1) |
| 70 | 35(177.5) | 25(127.6) | 25(177.2) | 39(282.8) | 29(216.6) | 27(320.4) | 41(603.6) | 29(581.4) | 29(818.1) | 41(1928.1) | 29(1582.8) | 29(3003.8) |
| 84 | 35(401.1) | 25(332.4) | 25(407.1) | 39(657.8) | 29(445.5) | 27(665.0) | 41(1156.5) | 29(869.2) | 29(1706.1) | 41(4696.6) | 29(2848.4) | 29(22402.2) |
Table 3: Results showing the total number of iterations from low-rank preconditioned GMRES and the total CPU times (in seconds) using preconditioner with $\beta=10^{-4}$ , $\sigma=0.4$ , and selected spatial ( $n_{h}$ ) and stochastic ( $n_{\xi}$ ) degrees of freedom
| $n_{\xi}$ $n_{h}$ $n_{\tau}$ | # iter(t) 289( $h=\frac{1}{2^{4}}$ ) 1 | # iter(t) 1089( $h=\frac{1}{2^{5}}$ ) $m_{\xi}$ +1 | # iter(t) 4225( $h=\frac{1}{2^{6}}$ ) $n_{A}$ | # iter(t) 16641( $h=\frac{1}{2^{7}}$ ) 1 | $m_{\xi}$ +1 | $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\sigma_{a}=0.4$ | | | | | | | | | | | | |
| Chebyshev-5+hGS-1 | | | | | | | | | | | | |
| 20 | 59(5.7) | 36(3.6) | 27(4.5) | 84(13.3) | 39(10.1) | 29(13.2) | 88(50.9) | 38(37.8) | 36(60.4) | 88(398.2) | 38(163.3) | 28(203.0) |
| 70 | 100(275.7) | 38(198.3) | 34(140.1) | 112(498.4) | 42(299.3) | 37(316.8) | 118(1792.6) | 41(697.2) | 36(929.6) | 117(7712.6) | 41(2655.2) | 36(3742.7) |
| 84 | 86(720.2) | 36(414.0) | 34(463.0) | 96(986.0) | 40(596.4) | 37(624.8) | 102(3051.9) | 40(1452.6) | 37(1838.1) | 102(9946.2) | 41(4576.9) | 36(7289.7) |
| Chebyshev-10+hGS-1 | | | | | | | | | | | | |
| 20 | 69(6.4) | 30(2.8) | 28(3.3) | 77(13.9) | 33(8.2) | 30(9.9) | 81(50.5) | 35(34.4) | 31(38.1) | 83(208.6) | 36(148.7) | 33(162.3) |
| 70 | 92(278.4) | 32(96.5) | 28(125.9) | 105(508.4) | 36(175.3) | 30(269.3) | 110(1318.5) | 37(781.9) | 32(857.9) | 110(5212.5) | 38(2920.1) | 33(4546.2) |
| 84 | 79(931.3) | 30(355.0) | 28(411.9) | 89(1373.9) | 34(529.2) | 30(531.0) | 95(3528.4) | 36(1337.8) | 32(1576.4) | 97(11001.9) | 38(4413.7) | 33(5769.8) |
| Cholesky+hGS-1 | | | | | | | | | | | | |
| 20 | 69(6.8) | 29(3.0) | 27(4.0) | 77(17.5) | 31(8.2) | 29(11.7) | 81(57.5) | 33(35.2) | 29(50.3) | 83(279.8) | 33(180.9) | 31(181.8) |
| 70 | 31(159.1) | 31(156.3) | 29(350.5) | 105(730.2) | 33(238.6) | 29(331.6) | 109(1903.1) | 35(726.1) | 31(953.3) | 109(7600.0) | 35(2737.8) | 31(3350.7) |
| 84 | 33(753.1) | 33(353.2) | 27(425.5) | 89(1325.8) | 33(493.8) | 29(634.8) | 95(3122.6) | 35(1055.9) | 31(1837.7) | 97(10923.5) | 35(4023.8) | 31(6429.3) |
Table 4 examines the sensitivity to the regularization parameter $\beta$ , which balances the tracking term and control cost in the objective functional (1). As $\beta$ decreases from $10^{-2}$ to $10^{-5}$ , the optimization problem becomes increasingly dominated by the tracking term. The iteration counts remain remarkably stable across this range, demonstrating that the hierarchical preconditioner effectively handles varying parameter regimes without requiring problem-specific tuning. The $n_{\tau}=m_{\xi}+1$ truncation consistently performs comparably to the full expansion while maintaining reduced computational cost.
Table 4: Results using the preconditioner with $m_{\xi}$ =3 p=3, $\sigma_{\Bbbk}=0.2$ , $\beta\in\{10^{-2},10^{-3},10^{-4},10^{-5}\}$ and $n_{h}=1089(h=\frac{1}{2^{5}})$ .
| $n_{\xi}$ $n_{\tau}$ | # iter(t) 20 1 | # iter(t) 70 $m_{\xi}$ +1 | # iter(t) 84 $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ | 1 | $m_{\xi}$ +1 | $n_{A}$ |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Chebyshev-5+hGS-1 | | | | | | | | | |
| $\beta=10^{-2}$ | 52(10.0) | 32(7.9) | 30(13.5) | 60(272.3) | 32(149.2) | 30(363.0) | 56(596.9) | 32(338.9) | 30(626.9) |
| $\beta=10^{-3}$ | 52(10.3) | 34(8.3) | 34(15.5) | 60(273.8) | 35(161.7) | 34(411.4) | 56(594.0) | 34(360.7) | 34(709.7) |
| $\beta=10^{-4}$ | 43(11.7) | 37(7.6) | 35(15.6) | 62(279.7) | 37(172.1) | 35(299.2) | 58(629.0) | 37(393.4) | 35(731.8) |
| $\beta=10^{-5}$ | 41(11.2) | 38(7.9) | 37(17.0) | 60(270.6) | 38(175.9) | 37(317.9) | 56(589.0) | 38(413.3) | 38(795.5) |
| Chebyshev-10+hGS-1 | | | | | | | | | |
| $\beta=10^{-2}$ | 48(8.7) | 30(6.1) | 29(11.9) | 56(420.8) | 29(147.9) | 29(359.7) | 52(577.6) | 34(370.8) | 24(653.1) |
| $\beta=10^{-3}$ | 48(8.8) | 31(6.4) | 30(12.6) | 56(432.3) | 31(151.6) | 30(406.9) | 52(554.8) | 31(339.6) | 30(675.8) |
| $\beta=10^{-4}$ | 48(8.8) | 30(7.1) | 30(19.4) | 56(419.0) | 31(151.2) | 30(390.6) | 52(575.0) | 31(339.6) | 30(676.1) |
| $\beta=10^{-5}$ | 47(8.7) | 30(7.0) | 30(13.2) | 54(406.4) | 30(146.0) | 30(383.9) | 51(554.0) | 30(333.0) | 30(675.1) |
| Chol+hGS-1 | | | | | | | | | |
| $\beta=10^{-2}$ | 47(12.4) | 29(6.9) | 27(11.9) | 55(380.1) | 29(214.4) | 27(354.4) | 51(975.8) | 29(483.9) | 27(595.1) |
| $\beta=10^{-3}$ | 47(12.0) | 29(6.4) | 27(11.9) | 55(396.0) | 29(215.9) | 27(364.5) | 51(876.9) | 29(455.6) | 27(596.0) |
| $\beta=10^{-4}$ | 47(10.4) | 29(7.4) | 27(12.2) | 55(387.9) | 29(216.0) | 27(343.9) | 53(847.9) | 29(432.8) | 27(595.5) |
| $\beta=10^{-5}$ | 45(10.0) | 29(6.8) | 27(12.0) | 53(369.6) | 29(205.5) | 27(359.3) | 51(829.3) | 29(438.9) | 27(594.1) |
### 5.2 Time-dependent case
This subsection evaluates the all-at-once preconditioner (30) for time-dependent optimal control problems. The discretization results in KKT systems of dimension $n_{t}\times n_{\xi}\times n_{h}$ , where $n_{t}$ denotes the number of time steps. We investigate the scalability with respect to: (i) mesh refinement (Table 5), (ii) regularization parameter $\beta$ (Table 6), (iii) variance $\sigma_{\Bbbk}$ (Table 7), (iv) time-step size discretization $\tau$ (Table 8), and (v) stochastic dimension $(m_{\xi},p)$ (Table 9). Based on the steady-state findings, we employ the Chebyshev-5+hGS-1 configuration unless otherwise noted, reporting results for both stringent ( $10^{-6}$ ) and moderate ( $10^{-4}$ ) tolerances to illustrate practical convergence behavior. As in the steady-state case, we use $N=1$ in the Richardson iteration (Algorithm 5) for the (3,3)-block.
From our observations in the steady-state problem, a combination of a 5-step Chebyshev smoother with one step of our hGS method achieves a good balance between the computational cost of matrix operations and GMRES iterations; thus, we typically adopt this combination when testing time-dependent cases as well.
As indicated in Table 5, the 5-step Chebyshev smoother yields consistent iteration counts compared to either the 10-step smoother or direct Cholesky decomposition. As the spatial discretization is refined from $h=\frac{1}{2^{3}}$ to $h=\frac{1}{2^{6}}$ , representing a growth from 116,640 to 6,084,000 DoF, the iteration count for $n_{\tau}=m_{\xi}+1$ exhibits sub-linear growth consistent with the near mesh-independence predicted by the spectral theory in Section 4. The $n_{\tau}=m_{\xi}+1$ truncation achieves iteration counts comparable to the full expansion ( $n_{\tau}=n_{A}$ ) while reducing the cost of assembling and applying the (3,3)-block preconditioner—a trade-off that becomes increasingly favorable as the problem dimension grows.
Table 5: Results using hGS method with different truncation settings $n_{\tau}$ with the time-dependent model for different stopping tolerance and mesh size with $m_{\xi}$ =3, $p=3$ $\sigma_{\Bbbk}$ =0.2, $\beta=10^{-4}$ and number of time steps=8 ( $\tau=\frac{1}{2^{3}}$ ).
| $\frac{1}{2^{3}}$ $\frac{1}{2^{4}}$ $\frac{1}{2^{5}}$ | 116,640 416,160 1,568,160 | 57(11.9) 75(39.3) 83(150.9) | 39(8.8) 47(29.6) 53(103.3) | 39(15.7) 45(39.6) 51(158.6) | 41(6.8) 55(21.7) 66(90.6) | 31(5.7) 35(14.9) 43(58.8) | 29(8.8) 35(24.7) 39(92.1) |
| --- | --- | --- | --- | --- | --- | --- | --- |
| $\frac{1}{2^{6}}$ | 6,084,000 | 84(570.8) | 55(383.3) | 53(585.1) | 68(401.4) | 42(252.2) | 42(424.1) |
Table 6 examines four orders of magnitude for $\beta$ , ranging from $10^{-2}$ (control-dominant) to $10^{-8}$ (tracking-dominant). The mean-based preconditioner ( $n_{\tau}=1$ ) exhibits strong dependence on $\beta$ , with iteration counts decreasing as $\beta$ decreases (since smaller $\beta$ yields problems dominated by the simpler tracking term). In contrast, the $n_{\tau}=m_{\xi}+1$ truncation maintains stable iteration counts across all tested values, demonstrating that the hGS preconditioner automatically adapts to the problem structure without manual parameter tuning. This robustness confirms the applicability of the theoretical framework across diverse parameter regimes.
Table 6: Results using hGS method with different truncations setting $n_{\tau}$ with the time-dependent model for different stopping tolerance and $\beta$ , mesh size $h=\frac{1}{2^{5}}$ , $m_{\xi}$ =3, $p=3$ , $\sigma_{\Bbbk}$ =0.2, and number of time steps=8 ( $\tau=\frac{1}{2^{3}}$ ), which results in 1,568,160 DoF.
| $10^{-2}$ $10^{-3}$ $10^{-6}$ | 107(206.1) 116(164.4) 97(137.0) | 69(136.3) 74(109.4) 77(115.6) | 69(214.8) 74(174.6) 77(182.5) | 74(102.0) 66(90.6) 57(77.4) | 48(70.4) 42(60.7) 45(65.3) | 47(112.6) 42(100.0) 44(103.7) |
| --- | --- | --- | --- | --- | --- | --- |
Table 7 explores the range $\sigma_{\Bbbk}\in\{0.01,0.02,0.05,0.1,0.2,0.4\}$ , spanning from nearly deterministic to highly uncertain regimes. The mean-based preconditioner ( $n_{\tau}=1$ ) exhibits significant degradation as uncertainty increases, whereas the $n_{\tau}=m_{\xi}+1$ truncation maintains stable iteration counts across the entire range. Especially when $\sigma_{\Bbbk}$ increases from 0.2 to 0.4, the mean-based preconditioner performs poorly with a large number of iterations, but the hGS method maintains robust performance. This robustness confirms that the hierarchical preconditioner effectively captures the essential stochastic structure without requiring full expansion of all coupling coefficients.
Table 7: Results using hGS method with different truncations setting $n_{\tau}$ with the time-dependent model for different stopping tolerance and $\sigma_{\Bbbk}$ , $\beta=10^{-4}$ , mesh size $n_{h}=1089$ , $m_{\xi}$ =3, $p=3$ , and number of time steps=8 ( $\tau=\frac{1}{2^{3}}$ ), which results in 1,568,160 DoF.
| 0.01 0.02 0.05 | 51(77.7) 53(80.8) 59(90.2) | 51(85.3) 51(82.2) 51(82.1) | 51(133.9) 51(130.0) 51(130.9) | 43(68.6) 43(70.2) 47(78.3) | 41(72.7) 41(71.1) 41(72.6) | 41(112.4) 41(113.7) 41(113.1) |
| --- | --- | --- | --- | --- | --- | --- |
| 0.1 | 65(99.2) | 51(82.4) | 51(130.4) | 53(86.2) | 42(70.6) | 41(116.5) |
| 0.2 | 83(150.9) | 53(103.3) | 51(158.6) | 66(90.6) | 43(58.8) | 39(92.1) |
| 0.4 | 129(200.9) | 57(92.3) | 51(130.5) | 100(179.3) | 47(82.3) | 43(119.2) |
Table 8 investigates the all-at-once system scalability by varying $n_{t}$ from 4 to 256 (time steps $\tau\in\{1/4,1/16,1/64,1/256\}$ ), corresponding to total system sizes ranging from 784,080 to over 12.5 million DoF. As the temporal resolution increases, the coupled space-time-stochastic system grows proportionally, yet the $n_{\tau}=m_{\xi}+1$ truncation maintains sub-linear iteration growth relative to system size. The computational time scales approximately linearly with DoF, confirming the efficiency of the all-at-once preconditioner for massively coupled systems.
Table 8: Results using hGS method with different truncations setting $n_{\tau}$ with the time-dependent model for different stopping tolerance and $\sigma_{\Bbbk}$ , $\beta=10^{-4}$ , mesh size $n_{h}=1089$ , $m_{\xi}=3,\,p=3$ , and number of time steps=8 ( $\tau=\frac{1}{2^{3}}$ ).
| $1/2^{2}$ $1/2^{4}$ $1/2^{6}$ | 784,080 3,136,320 12,545,280 | 81(78.6) 85(238.2) 101(1029.8) | 51(51.1) 55(166.4) 67(733.3) | 51(71.1) 53(258.7) 65(1205.8) | 66(53.9) 68(174.2) 78(840.5) | 43(37.6) 44(124.3) 53(604.2) | 42(57.2) 43(203.7) 52(1028.6) |
| --- | --- | --- | --- | --- | --- | --- | --- |
Finally, Table 9 compares three gPC configurations: $(m_{\xi},p)\in\{(3,3),(4,4),(6,3)\}$ , yielding $n_{\xi}\in\{20,70,84\}$ basis functions. Table 9 varies the stochastic discretization parameters $(m_{\xi},p)$ , exploring both the number of random variables and polynomial order. As $n_{\xi}$ increases from 20 to 84, the iteration count for $n_{\tau}=m_{\xi}+1$ grows modestly, demonstrating near-independence from the stochastic discretization level. This behavior confirms the effectiveness of the hierarchical truncation strategy in maintaining spectral properties across varying gPC expansion settings.
Table 9: Results using hGS method with different truncations setting $n_{\tau}$ with the time-dependent model for different stopping tolerance and stochastic setting, $\beta=10^{-4}$ , mesh size $n_{h}=1089$ , $\sigma_{\Bbbk}$ =0.2, and number of time steps=8 ( $\tau=\frac{1}{2^{3}}$ ).
| 20 70 84 | 1,568,160 5,488,560 6,586,272 | 83(150.9) 95(2435.3) 93(4092.4) | 53(103.3) 53(1397.4) 56(2670.3) | 51(158.6) 52(2964.7) 52(5207.4) | 66(90.6) 74(1809.2) 70(3059.4) | 43(58.8) 43(1044.9) 43(2107.7) | 39(92.1) 43(2476.2) 43(4413.4) |
| --- | --- | --- | --- | --- | --- | --- | --- |
The time-dependent experiments establish that the proposed all-at-once preconditioner maintains robust performance across a wide range of problem parameters and discretization levels. Four key findings emerge from these results. First, regarding near mesh-independence, iteration growth remains sub-linear with spatial refinement (Table 5), consistent with the spectral bounds derived in Section 4. Second, concerning parameter robustness, the hGS method adapts automatically to varying $\beta$ (Table 6) and $\sigma_{\Bbbk}$ (Table 7) without manual tuning. Third, regarding truncation efficiency, the $n_{\tau}=m_{\xi}+1$ strategy consistently delivers performance comparable to the full expansion at significantly reduced cost. These results demonstrate that the hierarchical preconditioning framework extends seamlessly from steady-state to time-dependent problems, providing a practical and theoretically-grounded solution for large-scale stochastic optimal control problems.
## 6 Conclusions
In the paper, we designed, analyzed, and implemented a novel hierarchical preconditioning strategy for large-scale stochastic optimal control problems. Our approach leverages a truncated stochastic expansion within a block-structured preconditioner for the Karush-Kuhn-Tucker (KKT) system, striking an effective balance between computational cost and preconditioning quality. Numerical results confirm that the proposed hGS method consistently outperforms both standard mean-based preconditioner and computationally intensive full-expansion preconditioner across a wide range of problem parameters.
A key contribution of this work is the extension of the truncated gPC framework to time-dependent problems. We developed and tested a tailored hGS preconditioner within an all-at-once discretization scheme, demonstrating the versatility and effectiveness of our approach for these challenging, large-scale scenarios. Comprehensive numerical experiments on benchmark problems have validated the robustness and numerical efficiency of the proposed algorithms.
## References
## Appendix A Auxiliary Results
### A.1 Proof of Lemma 4.2
**Proof A.1**
*Consider the Rayleigh quotient:
$$
\frac{\mathbf{v}^{T}A\mathbf{v}}{\mathbf{v}^{T}B\mathbf{v}}=\frac{\mathbf{v}^{T}(B+E)\mathbf{v}}{\mathbf{v}^{T}B\mathbf{v}}=1+\frac{\mathbf{v}^{T}E\mathbf{v}}{\mathbf{v}^{T}B\mathbf{v}}.
$$
Applying the assumption $|\frac{\mathbf{v}^{T}E\mathbf{v}}{\mathbf{v}^{T}B\mathbf{v}}|\leq\delta$ directly yields the bounds.*
### A.2 Proof of Lemma 4.3
**Proof A.2**
*Let $(\lambda,\mathbf{x})$ be an eigenpair satisfying the generalized eigenvalue problem $C\mathbf{x}=\lambda D\mathbf{x}$ , with eigenvector $\mathbf{x}\neq\mathbf{0}$ . We perform a change of variables by setting $\mathbf{x}=Q\mathbf{y}$ . Since $Q$ is nonsingular, $\mathbf{x}\neq\mathbf{0}$ implies that the transformed vector $\mathbf{y}\neq\mathbf{0}$ . Substituting $\mathbf{x}=Q\mathbf{y}$ into the original problem gives:
$$
C(Q\mathbf{y})=\lambda D(Q\mathbf{y}).
$$
Multiplying from the left by $Q^{T}$ , we obtain:
$$
(Q^{T}CQ)\mathbf{y}=\lambda(Q^{T}DQ)\mathbf{y}.
$$
This final expression is the generalized eigenvalue problem for the pair $(Q^{T}CQ,Q^{T}DQ)$ , which is satisfied by the same eigenvalue $\lambda$ with the transformed eigenvector $\mathbf{y}$ . Therefore, the sets of eigenvalues for both pairs are identical.*
### A.3 Proof of Lemma 4.4
**Proof A.3**
*Let the SVD of $A$ be $A=U\Sigma V^{*}$ , where $U,V$ are unitary and $\Sigma=\operatorname{diag}(\sigma_{1},\ldots,\sigma_{n})$ with $\sigma_{1}\geq\cdots\geq\sigma_{n}=\sigma_{\min}(A)\geq 0$ . For any nonzero vector $\mathbf{v}$ , let $\mathbf{w}=V^{*}\mathbf{v}$ . Since $V$ is unitary, $\|\mathbf{w}\|_{2}=\|\mathbf{v}\|_{2}$ . Then
$$
\|A\mathbf{v}\|_{2}^{2}=\|U\Sigma V^{*}\mathbf{v}\|_{2}^{2}=\|\Sigma\mathbf{w}\|_{2}^{2}=\sum_{i=1}^{n}\sigma_{i}^{2}|w_{i}|^{2}\geq\sigma_{\min}^{2}\sum_{i=1}^{n}|w_{i}|^{2}=\sigma_{\min}^{2}\|\mathbf{w}\|_{2}^{2}=\sigma_{\min}^{2}\|\mathbf{v}\|_{2}^{2}.
$$
Taking square roots on both sides gives $\|A\mathbf{v}\|_{2}\geq\sigma_{\min}(A)\|\mathbf{v}\|_{2}$ .*
### A.4 Proof of Lemma 4.5
**Proof A.4**
*From the triangle inequality, for any vector $x$ , with $\|x\|_{2}=1$ , we have
$$
\|(A+B)x\|_{2}=\|Bx-(-A)x\|_{2}\geq\|Bx\|_{2}-\|Ax\|_{2}.
$$
Taking the minimum over all unit vectors $x$ on both sides of the inequality, we get
$$
\min_{\|x\|_{2}=1}\|(A+B)x\|_{2}\geq\min_{\|x\|_{2}=1}\left(\|Bx\|_{2}-\|Ax\|_{2}\right).
$$
Using the property that $\min(f-g)\geq\min(f)-\max(g)$ , we obtain
$$
\min_{\|x\|_{2}=1}\left(\|Bx\|_{2}-\|Ax\|_{2}\right)\geq\min_{\|x\|_{2}=1}\|Bx\|_{2}-\max_{\|x\|_{2}=1}\|Ax\|_{2}.
$$
By the definitions of the minimum singular value and the operator norm, the above expression is equivalent to $\sigma_{\min}(A+B)\geq\sigma_{\min}(B)-\|A\|_{2}$ .*