A multigrid solver for PDE-constrained optimization with uncertain inputs

TL;DR

This paper introduces a multigrid solver tailored for large saddle-point systems in PDE-constrained optimization under uncertainty, demonstrating optimal complexity and robustness across various problem types.

Contribution

The paper develops a novel collective multigrid algorithm with convergence analysis, enabling efficient solutions of large saddle-point systems in uncertain PDE-constrained optimization.

Findings

Optimal $O(N)$ complexity for reduced saddle-point systems

Effective as a solver and preconditioner in diverse PDE optimization problems

Robust performance across different problem settings

Abstract

In this manuscript, we present a collective multigrid algorithm to solve efficiently the large saddle-point systems of equations that typically arise in PDE-constrained optimization under uncertainty, and develop a novel convergence analysis of collective smoothers and collective two-level methods. The multigrid algorithm is based on a collective smoother that at each iteration sweeps over the nodes of the computational mesh, and solves a reduced saddle-point system whose size is proportional to the number $N$ of samples used to discretized the probability space. We show that this reduced system can be solved with optimal $O(N)$ complexity. The multigrid method is tested both as a stationary method and as a preconditioner for GMRES on three problems: a linear-quadratic problem, possibly with a local or a boundary control, for which the multigrid method is used to solve directly the linear optimality system; a nonsmooth problem with box constraints and $L^1$-norm penalization on the control, in which the multigrid scheme is used as an inner solver within a semismooth Newton iteration; a risk-averse problem with the smoothed CVaR risk measure where the multigrid method is called within a preconditioned Newton iteration. In all cases, the multigrid algorithm exhibits excellent performances and robustness with respect to the parameters of interest.

Tables

Table 1. Table 1: Number of V-cycle (left) and preconditioned GMRES (right) iterations to solve ( 10 ) for a linear quadratic problem on the L-shaped domain 𝒟 = ( 0 , 1 ) 2 ∖ ( 0.5 , 1 ) ¯ 2 𝒟 superscript 0 1 2 superscript ¯ 0.5 1 2 \mathcal{D}=(0,1)^{2}\setminus\overline{(0.5,1)}^{2} with a distributed control.

$\nu$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$
It.	18  $|$  11	19  $|$  13	19  $|$  15	19  $|$  15

Table 2. Table 2: Number of V-cycle (left) and preconditioned GMRES (right) iterations to solve ( 10 ) for a linear quadratic problem on the square domain 𝒟 = ( 0 , 1 ) 2 𝒟 superscript 0 1 2 \mathcal{D}=(0,1)^{2} with a local control acting on 𝒟 0 = ( 0.25 , 0.75 ) 2 subscript 𝒟 0 superscript 0.25 0.75 2 \mathcal{D}_{0}=(0.25,0.75)^{2} .

$\nu$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$
It.	17  $|$  11	20  $|$  13	26  $|$  16	26  $|$  18

Table 3. Table 3: Number of V-cycle (left) and preconditioned GMRES (right) iterations to solve ( 10 ) for a linear quadratic problem on the square domain 𝒟 = ( 0 , 1 ) 2 𝒟 superscript 0 1 2 \mathcal{D}=(0,1)^{2} with a boundary control acting on Γ = ( 0 , 1 ) × { 0 } Γ 0 1 0 \Gamma=(0,1)\times\left\{0\right\} .

$\nu$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$
It.	17  $|$  14	22  $|$  15	23  $|$  17	21  $|$  16

Table 4. Table 4: Number of semismooth Newton iterations (left), and average number of V-cycle (center) and preconditioned GMRES (right) iterations (in brackets).

${\sigma }^2$	0.1	0.5	1	1.5
It.	4   $|$  22.5  $|$  14	5   $|$  22.6  $|$  14.2	8  $|$  23.0  $|$  11.8	14.9  $|$  22.9  $|$  15.0

Table 5. Table 5: Number of semismooth Newton iterations, of V-cycle iterations and of preconditioned GMRES iterations (in brackets). In the second row, the semismooth Newton method starts from a warm-up initial guess obtained through continuation.

$\nu$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$
It.	2  $|$  23.0  $|$  14.5	5  $|$  22.7  $|$   14.2	17  $|$   25.6  $|$   15.0	50  $|$   41.4  $|$   17.2
It.	2  $|$  23.0  $|$   14.5	4  $|$   22.7  $|$   14.2	5  $|$  22.25  $|$  15.4	8  $|$   58.8  $|$   20.9

Table 6. Table 6: For each numerical experiment, we report from the left to the right: the number of outer preconditioned Newton iterations, the total number of inner Newton iterations, the averaged number of V-cycle iterations and the averaged number of preconditioned GMRES iterations.

$N_h$($N_L$)	161 (2)	705 (3)	2945 (4)
It.	5  $|$  62  $|$  23.0  $|$  13.9	6  $|$  79  $|$  28.0  $|$  15.5	6  $|$  79  $|$  26.2  $|$  14.8

Equations

$$\mathcal{S}\mathbf{x}=\mathbf{f},$$

$$\mathcal{S}=\left(\begin{array}{ccccccc}{C}_1& & & & A_1^{\top }& & \\ & \ddots & & & & \ddots & \\ & & C_N& & & & A_N^{\top }\\ & & & G& D_1& \dots & D_N\\ A_1& & & E_1& & & \\ & \ddots & & ⋮& & & \\ & & A_N& E_N& & & \end{array}\right),$$

$$\underset{u\in U}{\min }\mathcal{R}\left[Q(y(\omega ),u)\right]$$

$$\text{s.t. }y(\omega )\in V\text{ satisfies}$$

$$⟨e(y(\omega ),u,\omega ),v⟩=0\forall v\in V,\text{ a.e. }\omega \in \Omega ,$$

$$a_{\omega }(y,v)=⟨\mathcal{B}u,v⟩,\forall v\in V,\mathbb{P}\text{-a-e. }\omega \in \Omega ,$$

$$0<a_{\min }(\omega )\Vert v{\Vert }_V^2\le a_{\omega }(v,v)\le a_{\max }(\omega )\Vert v{\Vert }_V^2,\forall v\in V,\mathbb{P}\text{-a.e. }\omega \in \Omega ,$$

$$\underset{u\in U,y\in L^2(\Omega \text{mathchar}24635V)}{\min }\frac{1}{2}\mathbb{E}\left[\Vert \mathcal{I}{y}_{\omega }-y_d{\Vert }_{L^2(\mathcal{D})}^2\right]+\frac{\nu }{2}\Vert u{\Vert }_U^2,$$

$$a_{\omega }(y_{\omega },v)=⟨\mathcal{B}u+f,v⟩,\forall v\in V,\mathbb{P}\text{-a.e. }\omega \in \Omega ,$$

$$\underset{u\in U}{\min }\frac{1}{2}\mathbb{E}\left[\Vert \mathcal{I}S(\mathcal{B}u)-y_d{\Vert }_{L^2(\mathcal{D})}^2\right]+\frac{\nu }{2}\Vert u{\Vert }_U^2.$$

$$(\nu \overline{u}-{\Lambda }_U{\mathcal{B}}^{\star }{S}^{\star }{\mathcal{I}}^{\star }(y_d-S(\overline{u}+f)),v{)}_U=0,\forall v\in U,$$

$$a_{\omega }(v,p_{\omega })=⟨z(\omega ),v⟩,\forall v\in V,\mathbb{P}\text{-a-e. }\omega \in \Omega .$$

$$a_{\omega }(y_{\omega },v)=⟨\mathcal{B}\overline{u}+f,v⟩,\forall v\in V,\mathbb{P}\text{-a-e. }\omega \in \Omega ,$$

$$a_{\omega }(v,p_{\omega })=⟨{\mathcal{I}}^{\star }(y_d-y_{\omega }),v⟩,\forall v\in V,\mathbb{P}\text{-a-e. }\omega \in \Omega ,$$

$$(\nu \overline{u}-{\Lambda }_U{\mathcal{B}}^{\star }\mathbb{E}\left[p_{\omega }\right],v{)}_U=0,\forall v\in U.$$

$$\mathbb{E}\left[X\right]\approx \hat{\mathbb{E}}\left[X\right]:=\sum _{i=1}^N{\zeta }_{i}X({\omega }_i),\text{with}\sum _{i=1}^N{\zeta }_i=1.$$

$$\left(\begin{array}{ccccccc}M& & & & A_1^{\top }& & \\ & \ddots & & & & \ddots & \\ & & M& & & & A_N^{\top }\\ & & & \nu M_U& -{\zeta }_1{B}^{\top }& \dots & -{\zeta }_N{B}^{\top }\\ A_1& & & -B& & & \\ & \ddots & & ⋮& & & \\ & & A_N& -B& & & \end{array}\right)\left(\begin{array}{c}{\mathbf{y}}_1\\ ⋮\\ {\mathbf{y}}_N\\ \mathbf{u}\\ {\mathbf{p}}_1\\ ⋮\\ {\mathbf{p}}_N\end{array}\right)=\left(\begin{array}{c}M{\mathbf{y}}_d\\ ⋮\\ M{\mathbf{y}}_d\\ 0\\ M\mathbf{f}\\ ⋮\\ M\mathbf{f}\end{array}\right),$$

$${\tilde{\mathbf{y}}}_i=\left(\begin{array}{c}({\mathbf{y}}_1{)}_i\\ ⋮\\ ({\mathbf{y}}_N{)}_i\end{array}\right)\in {\mathbb{R}}^N,{\tilde{\mathbf{p}}}_i=\left(\begin{array}{c}({\mathbf{p}}_1{)}_i\\ ⋮\\ ({\mathbf{p}}_N{)}_i\end{array}\right)\in {\mathbb{R}}^N,i=1,\dots ,N_h,$$

$$R_i\left(\begin{array}{c}\mathbf{y}\\ \mathbf{u}\\ \mathbf{p}\end{array}\right)=\left(\begin{array}{c}{\tilde{\mathbf{y}}}_i\\ u_i\\ {\tilde{\mathbf{p}}}_i\end{array}\right)=:{\mathbf{x}}_i.$$

$${\tilde{S}}_i=\left(\begin{array}{ccc}\text{diag}({\mathbf{c}}_i)& 0& \text{diag}({\mathbf{a}}_i)\\ 0& (G{)}_{i,i}& {\mathbf{d}}_i^{\top }\\ \text{diag}({\mathbf{a}}_i)& {\mathbf{e}}_i& 0\end{array}\right)$$

$${\mathbf{x}}^n={\mathbf{x}}^{n-1}+\theta \sum _{i=1}^{N_h}{P}_i{\tilde{S}}_i^{-1}{R}_i\left(\mathbf{f}-S{\mathbf{x}}^{n-1}\right),$$

$$u_i$$

$${\tilde{\mathbf{y}}}_i$$

$${\tilde{\mathbf{p}}}_i$$

$${\mathbf{d}}_i^{\top }\text{diag}({\mathbf{a}}_i{)}^{-1}\text{diag}({\mathbf{c}}_i)\text{diag}({\mathbf{a}}_i{)}^{-1}{\mathbf{e}}_i=(M{)}_{i,i}^3\sum _{j=1}^N{\zeta }_j(A_j{)}_{i,i}^{-2}>0,$$

$${\tilde{\mathbf{y}}}_i$$

$${\tilde{\mathbf{p}}}_i$$

$$a_{\omega }(y_{\omega },v)={\int }_{\mathcal{D}}\kappa (x,\omega )\nabla y(x,\omega )\cdot \nabla v(x)dx={\int }_{\mathcal{D}}u(x)v(x)dx,\forall v\in V,\mathbb{P}\text{-a.e. }\omega \in \Omega ,$$

$$\kappa (x,\omega )=e^{\sigma {\sum }_{j=1}^M\sqrt{{\lambda }_j}{b}_j(x)N_j(\omega )}\approx e^{g(x,\omega )},$$

$$\underset{u\in U_{ad}}{\min }\frac{1}{2}\mathbb{E}\left[\Vert y_{\omega }(u)-y_d{\Vert }_{L^2(\mathcal{D})}^2\right]+\frac{\nu }{2}\Vert u{\Vert }_{L^2(\mathcal{D})}^2+\beta \Vert u{\Vert }_{L^1(\mathcal{D})},$$

$$a_{\omega }(y_{\omega }(u),v)=(u+f,v{)}_{L^2(\mathcal{D})},\forall v\in V,\mathbb{P}\text{-a-e. }\omega \in \Omega ,$$

$$U_{ad}:=\left\{v\in L^2(\mathcal{D}):a\le u\le b\text{almost everywhere in }\mathcal{D}\right\},$$

$$(\nu \overline{u}-S^{\star }(y_d-S(\overline{u}+f)),\overline{u}-v)+\phi (\overline{u})-\phi (v)\ge 0,\forall v\in U_{ad}.$$

$$\mathcal{F}(u):=u-\frac{1}{\nu }\left(\max (0,\mathcal{T}u-\beta )+\min (0,\mathcal{T}u+\beta )-\max (0,\mathcal{T}u-\beta -\nu b)-\min (0,\mathcal{T}u+\beta -\nu a)\right),$$

$$u^{k+1}=u^k+d{u}^k,\text{with}\mathcal{G}(u^k)d{u}^k=-\mathcal{F}(u^k),$$

Full text

We test the multigrid method both as a stationary method and as a preconditioner for GMRES on three problems: a linear-quadratic problem, possibly with a local or a boundary control, for which the multigrid method is used to solve directly the linear optimality system; a nonsmooth problem with box constraints and $L^{1}$-norm penalization on the control, in which the multigrid scheme is used as an inner solver within a semismooth Newton iteration; a risk-adverse problem with the smoothed CVaR risk measure where the multigrid method is called within a preconditioned Newton iteration. In all cases, the multigrid algorithm exhibits excellent performances and robustness with respect to the parameters of interest.

MSC Classification]65M55 - 65F10 - 65K10 - 49J55

1 Introduction

In this work, we present a multigrid method to solve the saddle point system

[TABLE]

where $\mathbf{x}=(\mathbf{y},\mathbf{u},\mathbf{p})=(\mathbf{y}_{1},\dots,\mathbf{y}_{N},\mathbf{u},\mathbf{p}_{1},\dots,\mathbf{p}_{N})^{\top}$, $\mathcal{S}$ has the block structure

[TABLE]

and all submatrices involved represent the discretization of some differential operators. More details on each block are provided in Section 2. Matrices such as (2) are often encountered while solving PDE-constrained optimization problems under uncertainty of the form

[TABLE]

where $u$ is the unknown deterministic control, $y(\omega)$ is the state variable which satisfies a random PDE constraint expressed by $e(\cdot,\cdot,\omega)$ for almost every realization $\omega$ of the randomness, $Q$ is a real-valued quantity of interest (cost functional) and $\mathcal{R}$ is a risk measure. The vectors $\left\{\mathbf{y}_{j}\right\}_{j=1}^{N}$ and $\left\{\mathbf{p}_{j}\right\}_{j=1}^{N}$ are the discretizations of the state and adjoint variables $y(\omega)$ and $p(\omega)$ at the $N$ samples in which the random PDE constraint is collocated. The vector $\mathbf{u}$ is the discretization of the deterministic control $u$. Problems of the form (3) are increasingly employed in applications. The PDE constraints typically represent some underlying physical model whose behaviour should be optimally controlled, and the randomness in the PDE allows one to take into account the intrisinc variability or lack of knowledge on some parameters entering the model. The introduction of a risk measure in (3) allows one to construct robust controls that take into account the distribution of the cost over all possible realizations of the random parameters. Therefore, the topic has received a lot of attention in the last years, see, e.g. kouri2018optimization ; Kouri_Cvar ; martinez2018optimal ; doi:10.1137/19M1294952 ; geiersbach2020stochastic ; antil2021ttrisk ; nobile_vanzan2 ; eigel2018risk ; ASADPOURE20111131 .

However, few works have focused on efficient solvers for the optimality systems (1). A popular approach is to perform a Schur complement on $\mathbf{u}$ and solve the reduced system with a Krylov method (possibly with Conjugate Gradient), despite each iteration would then require the solution of $2N$ PDEs, with $A_{j}$ and $A_{j}^{\top}$ for $j=1,\dots,N$ Kourisparse . For a full-space formulation, block diagonal preconditioners have been proposed in Kouri2018 and analyzed in nobile_vanzan , using both an algebraic approach based on Schur complement approximations and an operator preconditioning framework.

In this manuscript, we present a multigrid method to solve general problems of the form (1) and show how this strategy can be used for the efficient solution of three different Optimal Control Problems Under Uncertainty (OCPUU). First, we consider a linear-quadratic OCPUU and use the multigrid algorithm directly to solve the linear optimality system. Second, we consider a nonsmooth OCPUU with box constraints and $L^{1}$ regularization on the control. To solve such problem, we use the collective multigrid method as an inner solver within an outer semismooth Newton iteration. Incidentally, we show that the theory developed for the deterministic OCPs with $L^{1}$ regularization can be naturally extended to the class of OCPUU considered here. Third, we study a risk-adverse OCPUU involving the smoothed Conditional Value at Risk (CVaR) and test the performance of the multigrid scheme in the context of a nonlinear preconditioned Newton method.

The multigrid algorithm is based on a collective smoother borzi2005multigrid ; borzi2009multigrid ; takacs2011convergence that, at each iteration, loops over all nodes of the computational mesh (possibly in parallel), collects all the degrees of freedom related to a node, and updates them collectively by solving a reduced saddle-point problem. For classical (deterministic) PDE-constrained optimization problems with a distributed control, this reduced system has size $3\times 3$, thus its solution is immediate borzi2009multigrid . In our context, the reduced problem has size $(2N+1)\times(2N+1)$, which can be large when dealing with a large number of samples. Fortunately, we show that it can be solved with optimal $O(N)$ complexity. Let us remark that collective multigrid strategies have been applied to OCPUU in Borzi ; Borzi2 and in rosseel2012optimal . This manuscript differs from the mentioned works since, on the one hand, Borzi ; Borzi2 considers a stochastic control $u$, therefore for (almost) every realization of the random parameters a different control $u(\omega)$ is computed through the solution of a standard deterministic OCP. On the other hand, rosseel2012optimal considers a stochastic Galerkin discretization, and hence the correspoding optimality system has a structure which is very different from (2).

The multigrid algorithm presented here assumes that all state and adjoint variables are discretized on the same finite element mesh. The control can instead live on a subregion of the computational mesh, so that the algorithm is applicable also to optimization problems with local or boundary controls.

Finally, we remark that the multigrid solver proposed is based on a hierarchy of spatial discretizations corresponding to different levels of approximation, but the discretization of the probability space remains fixed, that is, the number of samples remains constant across the multigrid hierarchy. The extension of the multigrid algorithm to coarsening procedures also in the probability space will be the subject of future endeavours. We hint at possible approaches and challenges in Section 3 (see Remark 1). Nevertheless, we stress that the multigrid algorithm can already be incorporated within outer optimization routines that take advantange of different levels of approximations of the probability space, see, e.g., Kourisparse ; kouri2014multilevel ; nobile_vanzan2 .

The rest of the manuscript is organized as follows. In Section 2 we introduce the notation, a classical linear-quadratic OCPUU, and interpret (2) as the matrix associated to the optimality system of a discretized OCPUU. In Section 3, we present the collective multigrid algorithm, discuss implementation details and show its performance to solve the linear-quadratic OCPUU. In Section 4, we consider a nonsmooth OCPUU with box constraints and a $L^{1}$ regularization on the control. Section 5 deals with a risk-adverse OCPUU. For each of these cases, we first show how the multigrid approach can be integrated into the solution process, by detailing concrete algorithms, and then we present extensive numerical experiments to show the efficiency of the proposed framework. Finally, we draw our conclusions in Section 6.

2 A linear-quadratic optimal control problem under uncertainty

Let $\mathcal{D}\subset\mathbb{R}^{d}$ be a Lipschitz bounded domain, $V\subset L^{2}(\mathcal{D})$ a Sobolev space (e.g. $H^{1}(\mathcal{D})$ equipped with suitable boundary conditions), and $(\Omega,\mathcal{F},\mathbb{P})$ a complete probability space. Given a function $u$ belonging to a Hilbert space $U$, we consider the linear elliptic random PDE

[TABLE]

where $a_{\omega}(\cdot,\cdot):V\times V\rightarrow\mathbb{R}$ is a bilinear form and $\langle\cdot,\cdot\rangle$ denotes the duality between $V$ and $V^{\prime}$. $\mathcal{B}:U\rightarrow V^{\prime}$ is a continuous control operator allowing possibly for a local control (i.e. a control acting only on a subset $\mathcal{D}_{0}\subset\mathcal{D})$ or a boundary control (i.e. a control acting as Neumann condition on a subset of $\partial\mathcal{D}$). To assure uniqueness and sufficient integrability of the solution of (4), we make the following additional assumption.

Assumption 1**.**

There exist two random variables $a_{\min}(\omega)$ and $a_{\max}(\omega)$ such that

[TABLE]

and further $a^{-1}_{\min}$ and $a_{\max}$ are in $L^{p}(\Omega)$ for some $p\geq 4$.

Under Assumption 1, it is well-known (see, e.g., lord_powell_shardlow_2014 ; Scheichl ) that (4) admits a solution in $V$ for $\mathbb{P}\text{-a.e. }\omega$, and the solution $y$, interpreted as a function-valued random variable $y:\omega\in\Omega\mapsto y(\omega)\in V$, lies in the Bochner space $L^{q}(\Omega\mathchar 24635\relax\;V)$, $q\leq p$, cohn2013measure . We often use the shorthand notation $y_{\omega}=y(\cdot,\omega)$ when the dependence on $x$ is not needed, or $y_{\omega}(u)$ if we wish to highlight the dependence on the control function $u$.

In this manuscript, we consider the minimization of functionals constrained by (4). Let us first focus on the linear-quadratic problem

[TABLE]

where $y_{d}\in L^{2}(\mathcal{D})$ is a target state, $f\in V^{\prime}$, $\mathbb{E}:L^{1}(\Omega)\rightarrow\mathbb{R}$ is the expectation operator, $\nu>0$, and $\mathcal{I}$ is the embedding operator from $V$ to $L^{2}(\mathcal{D})$. Introducing the linear control-to-state map $S:g\in V^{\prime}\rightarrow y_{\omega}(g)\in L^{2}(\Omega\mathchar 24635\relax\;V)$, the reduced formulation of (5) is

[TABLE]

Existence and uniqueness of the minimizer of (6) follows directly from standard variational arguments lions1971optimal ; hinze2008optimization ; troltzsch2010optimal ; kouri2018optimization . Furthermore, due to Assumption 1, the optimal control $\overline{u}$ satisfies the variational equality

[TABLE]

where $\Lambda_{U}$ is the Riesz operator of $U$. The adjoint operator $S^{\star}:L^{2}(\Omega\mathchar 24635\relax\;V^{\prime})\rightarrow V$ is characterized by $S^{\star}z=\mathbb{E}\left[p\right]$ where $p=p_{\omega}(x)$ is the solution of the adjoint equation

[TABLE]

The optimality condition (7) can thus be formulated as the optimality system

[TABLE]

To solve numerically (5), we replace the exact expectation operator $\mathbb{E}$ of the objective functional by a quadrature formula $\widehat{\mathbb{E}}$ with $N$ nodes $\left\{\omega_{i}\right\}_{i=1}^{N}$ and positive weights $\left\{\zeta_{i}\right\}_{i=1}^{N}$, namely

[TABLE]

Common quadrature formulae are Monte Carlo, Quasi-Monte Carlo and Gaussian formulae. The latter requires that the probability space can be parametrized by a sequence (finite or countable) of random variables $\left\{\chi_{j}\right\}_{j}$, each with distribution $\mu_{j}$, and the existence of a complete basis of tensorized $L^{2}_{\mu_{j}}$-orthonormal polynomials. Hence for the semi-discrete OCP, the $\mathbb{P}$-a.e. PDE-constraint is naturally collocated onto the nodes of the quadrature formula.

Concerning the space domain, we consider a family of regular triangulations $\left\{\mathcal{T}_{h}\right\}_{h>0}$ of $\mathcal{D}$, and a Galerkin projection onto a conforming finite element space $V^{h}\subset V$ of continuous piecewise polynomial functions of degree $r$ over $\mathcal{T}_{h}$. $N_{h}$ is the dimension of $V^{h}$ and $\left\{\phi_{i}\right\}_{i=1}^{N_{h}}$ is a nodal Lagrangian basis. We discretize the state and adjoint variables on the same finite element space. The control variable is discretized on the finite element space $U_{h}=\text{span}\left\{\phi_{i}\right\}_{i=1}^{N_{u}}$, where $N_{u}$ is possibly strictly smaller than $N_{h}$ in case of a local or a boundary control.

Once fully discretized, (9) can be expressed as

[TABLE]

where $A_{j}$ are the stiffness matrices associated to the bilinear forms $a_{\omega_{j}}(\cdot,\cdot)$, $M$ and $M_{U}$ are mass matrices corresponding to the finite element spaces $V_{h}$ and $U_{h}$, $B$ is the discretization of the control operator, $\mathbf{y}_{d}$ and $\mathbf{f}$ are the finite element discretizations of $y_{d}$ and $f$ respectively, while $\mathbf{y}_{j}$ and $\mathbf{p}_{j}$ are the discretizations of $y_{\omega_{j}}$ and $p_{\omega_{j}}$. Notice that the matrix in (10) could be symmetrized by multipling the first and the last $N$ rows by the weights $\left\{\zeta_{i}\right\}_{i=1}^{N}$. This would be also consistent with the theoretical interpretation of the blocks of the saddle point system as discretizations of continuous inner products. From the numerical point of view, we have not observed relevant advantanges in maintaining the weights. Since for more general problems (see, e.g., Sec. 4) the symmetry of the saddle point system cannot be recovered by multiplying some equations by the quadrature weights, we omit them in our discussion.

3 Collective multigrid scheme

In this section, we describe the multigrid algorithm to solve the full space optimality system (10). First, we consider a distributed control, so that $u$ lives on the whole computational mesh and $B=M$. Local and boundary controls are discussed at the end of the section. Second, for the sake of generality, we consider the more general matrix (2), so that our discussion covers also the different saddle-point matrices obtained in Sections 4 and 5.

For each node of the triangulation, let us introduce the vectors $\widetilde{\mathbf{y}}_{i}$ and $\widetilde{\mathbf{p}}_{i}$,

[TABLE]

which collect the degrees of freedom associated to the $i$-th node, the scalar $u_{i}=(\mathbf{u})_{i}$, and the restriction operators $R_{i}\in\mathbb{R}^{(2N+1)N_{h}\times 2N+1}$ such that

[TABLE]

The prolongation operators are $P_{i}:=R_{i}^{\top}$, while the reduced matrices $\widetilde{S}_{i}:=R_{i}SP_{i}\in\mathbb{R}^{2N+1\times 2N+1}$ represent a condensed saddle-point matrix on the $i$-th node, and satisfy

[TABLE]

with $\mathbf{c}_{i}:=((C_{1})_{i,i},\dots,(C_{N})_{i,i})^{\top}$, $\mathbf{a}_{i}:=((A_{1})_{i,i},\dots,(A_{N})_{i,i})^{\top}$, $\mathbf{e}_{i}=((E_{1})_{i,i},\dots,(E_{N})_{i,i})^{\top}$, $\mathbf{d}_{i}=((D_{1})_{i,i},\dots,(D_{N})_{i,i})^{\top}$, where $\text{diag}(\mathbf{v})$ denotes a diagonal matrix with the components of $\mathbf{v}$ on the main diagonal.

Given an initial vector $\mathbf{x}^{0}$, a Jacobi-type collective smoothing iteration computes for $n=1,\dots,n_{1}$,

[TABLE]

where $\theta\in(0,1]$ is a damping parameter. Gauss-Seidel variants can straightforwardly be defined. Next, we consider a sequence of meshes $\left\{\mathcal{T}_{h_{\ell}}\right\}_{\ell=\ell_{\min}}^{\ell_{\max}}$, which we assume for simplicity to be nested, and restriction and prolongator operators $R_{\ell-1}^{\ell}$, $P_{\ell-1}^{\ell}$ which map between grids $\mathcal{T}_{h_{\ell-1}}$ and $\mathcal{T}_{h_{\ell}}$. In the numerical experiments, the coarse matrices are defined recursively in a Galerkin fashion starting from the finest one, namely $S_{\ell}:=R^{\ell+1}_{\ell}S_{\ell+1}P^{\ell+1}_{\ell}$ for $\ell\in\left\{1,\dots,\ell_{\max}-1\right\}$. Nevertheless it is obviously possible to define $S_{\ell}$ as the discretization of the continuous saddle-point system onto the mesh $\mathcal{T}_{h_{\ell}}$. With this notation, the V-cycle collective multigrid is described by Algorithm 1, which can be repeated until a certain stopping criterium is satisfied. We used the notation Collective_Smoothing$(\cdot,\cdot,\cdot)$ to denote possible variants of (12) (e.g. Gauss-Seidel).

Notice that (12) requires to invert the matrices $S_{i}$ for each computational node. We now show that this can be done with optimal $O(N)$ complexity. Indeed, performing a Schur complement on $u_{i}$, the system $\widetilde{S}_{i}\mathbf{x}_{i}=\mathbf{f}_{i}$, with $\mathbf{f}_{i}=(\mathbf{f}_{p_{i}},b_{u_{i}},\mathbf{f}_{y_{i}})^{\top}$ can be solved exclusively computing inverses of diagonal matrices and scalar products between vectors through

[TABLE]

Notice that we should guaranteee that $\text{diag}(\mathbf{a}_{i})$ admits an inverse and that $(G)_{i,i}+\mathbf{d}_{i}^{\top}\text{diag}(\mathbf{a}_{i})^{-1}\text{diag}(\mathbf{c}_{i})\text{diag}(\mathbf{a}_{i})^{-1}\mathbf{e}_{i}\neq 0$. This has to be verified case by case, so we now focus on the specific matrix (10). On the one hand, the vectors $\mathbf{a}_{i}$ are strictly positive componentwise, since $(\mathbf{a}_{i})_{j}=a_{\omega_{j}}(\phi_{i},\phi_{i})>0$ $\forall i=1,\dots,N_{h}$, $j=1,\dots,N$ (due to Assumption 1). On the other hand, $(G)_{i,i}=\int_{\mathcal{D}}\phi^{2}_{i}(x)\ dx>0$, while a direct calculation shows that

[TABLE]

which implies that the denominator in the first equation of (13) is strictly positive.

The collective smoother can be easily adjusted to accomodate local or boundary controls as discussed in borzi2003multigrid for deterministic OCPs. For all nodes $i$ for which a control basis function is present, the smoothing procedure remains that of (13). For all others computational nodes for which there is not a control basis function associated, the smoothing procedure becomes

[TABLE]

which is consistently obtained from (13) setting $u_{i}=0$.

To conclude this section we remark that, if the V-cycle algorithm requires a constant number of iterations to reach a given tolerance as both $N_{h}$, $N$ and the number of levels increase (so that the size of the coarse problem becomes neglegible in the limit), Algorithm 1 exhibits an overall optimal linear complexity $O(N_{h}N)$ since asymptotically the cost of each iteration is dominated by the $N_{h}$ relaxation steps, each of a cost proportional to $N$ due to the above discussion. In the next numerical experiments sections, we show indeed that the number of iterations remain constant for several test cases.

Remark 1** (Extension to a hierarchy of samples).**

*The multigrid algorithm presented is based on a hierarchy of spatial discretizations. However, the sample to discretize the probability space remains fixed among the levels. If one relies on the stochastic collocation method to discretize the probability space, it is possible to envisage a multigrid algorithm that also involves a coarsening of the sample size, since for each sample one could consider the associated stable interpolator which can then be evaluated onto a coarser or finer set of samples. Nevertheless, it is not clear at the moment the interplay between the smoothing and coarsening procedures, which is key for the efficient behaviour of a multigrid scheme. Future endeavours will investigate this interesting direction. For the rest of the manuscript we restrict oursevels to a hierarchy of spatial discretizations since on the one hand, the multigrid algorithm can already be embedded in other outer optimization algorithms that involve a hierarchy of samples Kourisparse ; kouri2014multilevel ; nobile_vanzan2 ; van2019robust . On the other hand, the reduced system can be solved with optimal $O(N)$ linear complexity, so that a coarsening in the number of samples may be superfluous. *

3.1 Numerical experiments

We now show the performance of Alg. 1 and its robustness with respect to several parameters for the solution of (10). We first consider the state equation

[TABLE]

in the L-shaped domain $\mathcal{D}=(0,1)^{2}\setminus\overline{(0.5,1)}^{2}$ discretized with a regular mesh of squares of edge $h_{\ell}=2^{-\ell}$, which are then decomposed into two right triangles. We choose $\kappa(x,\omega)$ as an approximated log-normal diffusion field

[TABLE]

where $g(x,\omega)$ is a mean zero Gaussian field with Covariance function $Cov_{g}(x,y)=\sigma^{2}e^{\frac{-\|x-y\|_{2}^{2}}{L^{2}}}$. The parameter $\sigma^{2}$ tunes the variance of the random field, while $L$ denotes the correlation length. The pairs $(b_{j}(x),\sigma^{2}\lambda_{j})$ are the eigenpairs of $T:L^{2}(\mathcal{D})\rightarrow L^{2}(\mathcal{D})$, $(Tf)(x)=\int_{\mathcal{D}}Cov_{g}(x,y)f(y)\ dy$, and $N_{j}\sim\mathcal{N}(0,1)$. Assumption 1 is satisfied since $a_{\min}(\omega)=\left(\text{ess}\inf_{x\in\mathcal{D}}\kappa(x,\omega)\right)^{-1}$ and $a_{\max}(\omega)=\|\kappa(\cdot,\omega)\|_{L^{\infty}(\mathcal{D})}$ are in $L^{p}(\Omega)$ for every $p<\infty$ Charrier . The target state is $y_{d}=e^{y^{2}}\sin(2\pi x)\sin(2\pi y)$.

Table 1 shows the number of V-cycle iterations (Alg. 1) and of GMRES iterations preconditioned by the V-cycle to solve (10) up to a tolerance of $10^{-9}$ on the relative (unpreconditioned) residual. Inside the V-cycle algorithm, we use $n_{1}=n_{2}=2$ pre- and post-smoothing iterations based on the Jacobi relaxation (12) with a damping parameter $\theta=0.5$ (the same value will be used for all numerical experiments in this manuscript). The number of levels of the V-cycle hierachy is denoted with $N_{L}$.

The first four subtables are based on a discretization of the probability space using the Stochastic Collocation method babuvska2010stochastic on Gauss-Hermite tensorized quadrature nodes, since for $L^{2}=0.5$, setting $M=3$ into (16) is enough to preserve $99\%$ of the variance. Remark that the multigrid algorithm has a perfect robustness with respect to all parameters considered. In particular, the convergence does not deteriorate as $\nu\rightarrow 0$, which is a well-known troublesome limit for several preconditioners (see, e.g., nobile_vanzan ; rees2010optimal ; pearson2012new ; ciaramella2022iterative and references therein). Further, the third sub-table shows the robustness of the algorithm with respect to the number of levels as the fine grid is refined. The fourth sub-table instead shows that the number of iteration remains constant as the discretization of the probability space is refined. As a complementary experiment, we set $L^{2}=0.1$ and use the Monte Carlo method, since we need $M=15$ to preserve $99\%$ of the variance of the random field, and the Stochastic Collocation method suffers from the curse of dimensionality. The fifth subtable confirms the robustness of the algorithm even for much larger numbers of samples.

Next, we consider the unit square domain $\mathcal{D}=(0,1)^{2}$ and a local control acting on the subset $\mathcal{D}_{0}=(0.25,0.75)^{2}\subset\mathcal{D}$, and a Neumann boundary control acting on $\Gamma=(0,1)\times\left\{0\right\}\subset\partial\mathcal{D}$. Tables 2 and 3 report the performances of the multigrid algorithm for these two cases. We stress once more the excellent robustness and efficiency of the multigrid algorithm in all regimes.

4 An optimal control problem under uncertainty with box-constraints and $L^{1}$ penalization

In this section, we consider the nonsmooth OCPUU111To keep a light notation, we omitted the continuous embedding operator from $L^{2}(\Omega\mathchar 24635\relax\;V)$ to $L^{2}(\Omega\mathchar 24635\relax\;L^{2}(\mathcal{D}))$ and from $L^{2}(\mathcal{D})$ to $V^{\prime}$, see Section 2 and, e.g., (troltzsch2010optimal, , Section 2.13).

[TABLE]

with $a<0<b$ and $\nu,\beta>0$. Deterministic OCPs with a $L^{1}$ penalization lead to optimal controls which are sparse, i.e. they are nonzero only on certain regions of the domain $\mathcal{D}$ stadler2009elliptic ; casas2017review . Sparse controls can be of great interest in applications, because it is often not desirable, or even impossible, to control the system over the whole domain $\mathcal{D}$. For sparse OCPUU, we mention li2019sparse where the authors considered both a simplified version of (17) in which the randomness enters linearly into the state equation as a force term, and a different optimization problem whose goal is to find a stochastic control $u(\omega)$ which has a similar sparsity pattern regardless of the realization $\omega$. Note further that the assumption $\nu>0$ does not eliminate the nonsmoothness of the objective functional, but it regularizes the optimal solution $u$, and is needed to use the fast optimation algorithm described in the following.

The well-posedness of (17) follows directly from standard variational arguments troltzsch2010optimal ; hinze2008optimization , being $U_{ad}$ a convex set, $\varphi(u):=\beta\|u\|_{L^{1}(\mathcal{D})}$ a convex function and the objective functional coercive. In particular, the optimal solution $\overline{u}$ satisfies the variational inequality ((ekeland1999convex, , Proposition 2.2))

[TABLE]

Through a pointwise discussion of the box constraints and an analysis of a Lagrange multiplier belonging to the subdifferential of $\varphi$ in $\overline{u}$, stadler2009elliptic showed that (18) can be equivalently formulated as the nonlinear equation $\mathcal{F}(\overline{u})=0$, with $\mathcal{F}:L^{2}(\mathcal{D})\rightarrow L^{2}(\mathcal{D})$ defined as

[TABLE]

where $\mathcal{T}:L^{2}(\mathcal{D})\ni u\rightarrow-S^{\star}(Su)+S^{\star}(y_{d}-Sf)\in L^{2}(\mathcal{D})$. Notice that $\mathcal{F}$ is nonsmooth due to the presence of the Lipschitz functions $\max(\cdot)$ and $\min(\cdot)$. Nevertheless, $\mathcal{F}$ can be shown to be semismooth hinze2008optimization , provided that $\mathcal{T}$ is continuously Fréchet differentiable, and further Lipschitz continuous interpreted as map from $\mathcal{T}:L^{2}(\mathcal{D})\rightarrow L^{r}(\mathcal{D})$, with $r>2$ doi:10.1137/1.9781611970692 ; hinze2008optimization . These conditions are satisfied also in our settings since $\mathcal{T}$ is affine and further the adjoint variable $p_{\omega}$, solution of (8) with $z=y_{d}-S(u+f)$, lies in $L^{2}(\Omega,H^{1}_{0}(\mathcal{D}))$ so that $\mathcal{T}u=\mathbb{E}\left[p_{\omega}\right]\in H^{1}_{0}(\mathcal{D})\subset L^{r}(\mathcal{D})$, where $r>2$ follows from Sobolev embeddings.

Hence, to solve (19) we use the semismooth Newton method whose iteration reads for $k=1,2,\dots$ until convergence,

[TABLE]

$\mathcal{G}(u):L^{2}(\mathcal{D})\rightarrow L^{2}(\mathcal{D})$ being the generalized derivative of $\mathcal{F}$. Using the linearity of $\mathcal{T}$ and considering the supports of the weak derivatives of $\max(0,x)$ and $\min(0,x)$, we obtain that

[TABLE]

where $\chi$ is the charateristic function of the union of the disjoint sets

[TABLE]

It is possible to show that the generalized derivative $\mathcal{G}(u)$ is invertible with bounded inverse for all $u$, the proof being identical to the deterministic case treated in Stadler2 . This further implies that the semismooth Newton method (20) converges locally superlinearly doi:10.1137/1.9781611970692 . We briefly summarize these results in the following proposition.

Proposition 2**.**

Let the inizialization $u^{0}$ be sufficiently close to the solution $\overline{u}$ of (17). Then the iterates $u^{k}$ generated by (20) converge superlinearly to $\overline{u}\in L^{2}(\mathcal{D})$.

Introducing the supporting variables $dy^{k}_{\omega}$ and $dp^{k}_{w}$ in $L^{2}(\Omega\mathchar 24635\relax\;H^{1}_{0}(\mathcal{D}))$, the semismooth Newton equation $\mathcal{G}(u^{k})du^{k}=-\mathcal{F}(u^{k})$ may be rewritten as the equivalent saddle point system

[TABLE]

Further, if we set $y^{0}=S(f+u^{0})$ and $p^{0}=S^{\star}(y_{d}-y^{0})$, due to the linearity of $S$ and $S^{\star}$, it holds $y^{k+1}=S(u^{k+1})=y^{k}+dy^{k}$ and similarly $p^{k+1}=p^{k}+dp^{k}$. Once fully discretized and using the notation $\widehat{\mathbb{E}}\left[p_{\omega}\right]=\sum_{j=1}^{N}\zeta_{j}\mathbf{p}_{j}$, the optimality condition (19) can be expressed through the nonlinear finite-dimensional map $\mathbf{F}:\mathbb{R}^{N_{h}}\rightarrow\mathbb{R}^{N_{h}}$,

[TABLE]

where the $\max(\cdot)$ and $\min(\cdot)$ functions act componentwise. Equation (22) leads to the saddle point system

[TABLE]

where $H^{k}\in\mathbb{R}^{N_{h}\times N_{h}}$ is a diagonal matrix representing the charateristic function $\chi_{I_{k}^{+}\cup I_{k}^{-}}$, namely

[TABLE]

with

[TABLE]

To derive the expression of $H$, we assumed that a Lagrangian basis is used for the finite element space. Notice that (24) fits into the general form (2), and thus we use the collective multigrid algorithm to solve it. Further, with the notation of (2), it holds

[TABLE]

if $i\in I^{+}\cup I^{-}$, and

[TABLE]

if $i\notin I^{+}\cup I^{-}$. The collective multigrid iteration is then well-defined.

The overall semismooth Newton Algorithm is summarized in Algorithm 2. At each iteration we solve (24) using the collective multigrid algorithm (line 4) and update the active sets given the new iteration (line 10). Notice that in order to globalize the convergence, we consider a line-search step (lines 6-8) performed on the merit function $\phi(\mathbf{u})=\sqrt{\mathbf{F}(\mathbf{u})^{\top}M\mathbf{F}(\mathbf{u})}$ martinez1995inexact .

4.1 Numerical experiments

In this section we test the semismooth Newton algorithm for the solution of (19) and the collective multigrid algorithm to solve the related optimality system (24). We consider the random PDE-constraint (15) with the random diffusion coefficient (16) set on the L-squared domain. The semismooth iteration is stopped when $\phi(\mathbf{u}^{k})<10^{-9}$. The inner linear solvers are stopped when the relative (unpreconditioned) residual is smaller than $10^{-11}$.

Table 4 reports the number of semismooth Newton iterations and in brackets the averaged number of iterations of the V-cycle algorithm used as a solver (left) or as preconditioner for GMRES (right). Table 4 confirms the effectiveness of the multigrid algorithm, which requires essentially the same computational effort as in the linear-quadratic case.

More challenging is the limit $\nu\rightarrow 0$ reported in Table 5. The performance of both the (globalized) semismooth Newton iteration and the inner multigrid solver deteriorates. The convergence of the outer nonlinear algorithm can be improved by performing a continuation method, namely we consider a sequence of $\nu=10^{-j}$, $j=2,\dots,8$, and we start the $j$-th problem using as initial condition the optimal solution computed for $\nu=10^{-j+1}$. Concerning the inner solver, the stand-alone multigrid algorithm struggles since for small values of $\nu$ the optimal control is of bang-bang type, that is satisfies $u=a$, $u=b$ or $u=0$ for almost every point of the mesh (for $\nu=10^{-8}$ only five nodes are nonactive at the optimum). The matrices $H^{k}$ are then close to zero, and the multigrid hierarchy struggles to capture changes at such small scale. Nevertheless, the multigrid algorithm remains a very efficient preconditioner for GMRES even in this challenging limit.

Fig. 1 shows a sequence of optimal controls for different values of $\beta$ with and without box-constraints. The optimal control for $\beta=0$ and without box-constraints corresponds to the minimizer of the linear-quadratic OCP (5). We observe that $L^{1}$ penalization indeed induces sparsity, since the optimal controls are more and more localized as $\beta$ increases. Numerically we have verified that for sufficiently large $\beta$, the optimal control is identically equal to zero, a property shown in stadler2009elliptic .

5 A risk-adverse optimal control problem under uncertainty

In this section we consider an instance of risk-adverse OCPUU. This class of problems has recently drawn lot of attention since in engineering applications it is important to compute a control that minimizes the quantity of interest even in rare, but often troublesome, scenarios Kouri_Cvar ; Kouri_ex ; antil2021ttrisk ; kouri2022primal . As a risk-measure shapiro2014lectures , we use the Conditional Value-At-Risk (CVaR) of confidence level $\lambda\in(0,1)$,

[TABLE]

that is, the expected value of a quantity of interest $X$ given that the latter is greater than or equal to its $\lambda$-quantile, here denoted by $\text{VaR}_{\lambda}\left(X\right)$. Rockafellar and Uryasev rockafellar2000optimization proved that $\text{CVaR}_{\lambda}\left(X\right)$ admits the equivalent formulation

[TABLE]

where $(\cdot)^{+}:=\max(0,\cdot)$, if the distribution of $X$ does not have an atom at $\text{VaR}_{\lambda}\left(X\right)$. In order to use tools from smooth optimization, we rely on a smoothing approach proposed in Kouri_Cvar , which consists in replacing $(\cdot)^{+}$ with a smooth function $g_{\varepsilon}$, $\varepsilon\in\mathbb{R}^{+}$, such that $g_{\varepsilon}\rightarrow(\cdot)^{+}$ in some functional norm as $\varepsilon\rightarrow 0$. Specifically, we choose the $C^{2}$-differentiable approximation

[TABLE]

Then, the smoothed risk-adverse OCPUU is

[TABLE]

where $\nu\in\mathbb{R}^{+}$ and $\lambda\in[0,1)$. The well-posedness of (28), the differentiability of its objective functional, as well as bounds for the error introduced by replacing $(\cdot)^{+}$ with $g_{\varepsilon}(\cdot)$, have been analyzed in Kouri_Cvar . Further, defining $Q_{\omega}=\frac{1}{2}\|y_{\omega}-y_{d}\|^{2}_{L^{2}(\mathcal{D})}-t$, the optimality conditions form the nonlinear system,

[TABLE]

Approximating $V$ and $\mathbb{E}$ with $V_{h}$ and $\widehat{\mathbb{E}}$, and letting $\widetilde{\mathbf{x}}=(\mathbf{y},\mathbf{u},\mathbf{p},t)$, the finite-dimensional discretization of (29) correponds to the nonlinear system $\widetilde{\mathbf{F}}(\widetilde{\mathbf{x}})=\mathbf{0}$, where $\widetilde{\mathbf{F}}:\mathbb{R}^{(2N+1)N_{h}+1}\rightarrow\mathbb{R}^{(2N+1)N_{h}+1}$,

[TABLE]

with $A=\text{diag}(A_{1},\dots,A_{N})$, $I=[I_{Nh},\dots,I_{Nh}]\in\mathbb{R}^{N_{h}\times N_{h}N}$, $I_{h}$ being the identity matrix, $\mathbf{y}_{d}$ is the discretization of $y_{d}$, and

[TABLE]

for $j=1,\dots,N$.

A possible approach to solve (30) is to use a Newton method, which given $\mathbf{x}^{k}=(\mathbf{y}^{k},\mathbf{u}^{k},\mathbf{p}^{k},t^{k})$ computes the corrections $\widetilde{\mathbf{dx}}^{k}=(\mathbf{dy}^{k},\mathbf{du}^{k},\mathbf{dp}^{k},dt^{k})$ solution of $\widetilde{\mathbf{J}}^{k}\widetilde{\mathbf{dx}}^{k}=-\widetilde{\mathbf{F}}(\widetilde{\mathbf{x}}^{k})$, where

[TABLE]

with

[TABLE]

for $i=1,\dots,N$. Unfortunately, $\widetilde{\mathbf{J}}^{k}$ can be singular away from the optimum, in particular whenever $\widehat{\mathbb{E}}\left[g_{\varepsilon}^{\prime\prime}(Q_{\omega}^{k})\right]=0$ which implies

[TABLE]

which is not unlikely for small $\varepsilon$ since $\text{supp}(g_{\varepsilon}^{\prime\prime})=(-\frac{\varepsilon}{2},\frac{\varepsilon}{2})$. Splitting strategies have been proposed (e.g. Markowski2022 in a reduced approach), in which whenever (35) is satisfied, an intermediate value of $t$ is computed by solving $\widetilde{F}_{4}(t\mathchar 24635\relax\;\mathbf{y},\mathbf{u},\mathbf{p})=0$ so to violate (35). In the next section, we discuss a similar splitting approach. To speed up the convergence of the outer nonlinear algorithm, we use a preconditioned Newton method based on nonlinear elimination doi:10.1137/S106482759325154X . At each iteration we will need to invert saddle-point matrices like (2), possibly several times. To do so, we rely on the collective multigrid algorithm.

5.1 Nonlinear preconditioned Newton method

Nonlinear elimination is a nonlinear preconditioning technique based on the identification of variables and equations of $\mathbf{F}$ (e.g. strong nonlinearities) that slow down the convergence of Newton method. These components are then eliminated through the solution of a local nonlinear problem at every step of an outer Newton. This elimination step provides a better initial guess for the outer iteration, so that a faster convergence is achieved doi:10.1137/S106482759325154X ; doi:10.1137/15M104075X .

In light of the possible singularity of $\widetilde{\mathbf{J}}$, we split the discretized variables $\widetilde{\mathbf{x}}$ into $\widetilde{\mathbf{x}}=(\mathbf{x},t)$, and we aim to eliminate the variables $\mathbf{x}$ to obtain a scalar nonlinear equation only for $t$. To do so, we partition (29) as

[TABLE]

where $\mathbf{F}_{1}=(\widetilde{\mathbf{F}}_{1}(\mathbf{x},t),\widetilde{\mathbf{F}}_{2}(\mathbf{x},t),\widetilde{\mathbf{F}}_{3}(\mathbf{x},t))$ and $F_{2}(\mathbf{x},t)=\widetilde{F}_{4}(\mathbf{x},t)$. Similarly, $\widetilde{\mathbf{J}}$ is partitioned into

[TABLE]

whose blocks have dimensions $\mathbf{J}_{1,1}\in\mathbb{R}^{(2N+1)N_{h}\times(2N+1)N_{h}}$, $\mathbf{J}_{1,2}\in\mathbb{R}^{(2N+1)N_{h}\times 1}$, $\mathbf{J}_{2,1}\in\mathbb{R}^{1\times(2N+1)N_{h}}$, and $\mathbf{J}_{2,2}\in\mathbb{R}$. Notice that $\mathbf{J}_{1,1}$ is always nonsingular, while $\mathbf{J}_{2,1}$, $\mathbf{J}_{1,2}$ and $J_{2,2}$ are identically zero if (35) is verified. Thus $\mathbf{F}_{1}$ allows us to define an implicit map $h:\mathbb{R}\rightarrow\mathbb{R}^{(2N+1)N_{h}}$, such that $\mathbf{F}_{1}(h(t),t)=0$, so that the first set of nonlinear equations in (36) are satisfied. We are then left to solve the nonlinear scalar equation

[TABLE]

To do so using the Newton method, we need the derivative of $F(t)$ evaluated at $t=t^{k}$ which, using implicit differentiation, can be computed as

[TABLE]

The nonlinear preconditioned Newton method is described in Alg. 3, and consists in solving (37) with Newton method. However, to overcome the possible singularity of $J^{k}_{2,2}$, $\mathbf{J}^{k}_{1,2}$ and $\mathbf{J}^{k}_{2,1}$, we check at each iteration $k$ if (35) is satisfied, and in the affermative case we update $\mathbf{x}^{k}$ by solving $\mathbf{F}_{1}(\mathbf{x}^{k+1},t^{k})=0$ using Newton method, and update $t^{k}$ by solving $F_{2}(\mathbf{x}^{k},t^{k+1})=0$. Notice further, that each iteration of the backtracking line-search requires to solve $F_{1}(h(t),t)=0$ using Newton method, thus additional linear systems with matrix $\mathbf{J}_{1,1}$ must be solved.

We report that we also tried to eliminate $t$ by computing the map $l$ such that $F_{2}(\mathbf{x},l(\mathbf{x}))=0$, while iterating on the variable $\mathbf{x}$. This has the advantage that $l$ can be evaluted very cheaply, being a scalar equation. However, we needed many more iterations both of the outer Newton method, and consequently of the inner linear solver. Thus, according to our experience, this second approach was less efficient and appealing.

5.2 Numerical experiments

In this section we report numerical tests to asses the performance of the preconditioned Newton algorithm to solve (37), and of the collective multigrid algorithm to invert the matrix $\mathbf{J}_{1,1}$. We consider the random PDE-constraint (15) with the random diffusion coefficient (16). Table 6 reports the number of outer and inner Newton iterations, and the average number of V-cycle iterations and of preconditioned GMRES iterations to solve the linear systems at each (inner/outer) Newton iterations. The outer Newton iteration is stopped when $|F(t^{k})|\leq 10^{-6}$, the inner Newton method to compute $h(\cdot)$ is stopped when $\max\left(\|\mathbf{F}_{1,1}(\mathbf{x}^{k}\mathchar 24635\relax\;t)\|_{2}/\|\mathbf{F}_{1,1}(\mathbf{x}^{0}\mathchar 24635\relax\;t)\|_{2},\|\mathbf{F}_{1,1}(\mathbf{x}^{k}\mathchar 24635\relax\;t)\|_{2}\right)\leq 10^{-8}$, and the linear solvers are stopped when the relative (unpreconditioned) residual is smaller than $10^{-9}$.

In Table 6, the number of outer Newton iterations is stable, while the number of inner Newton iterations varies between five and fifteen iterations per outer iteration. This is essentially due to how difficult it is to compute the nonlinear map $h(t)$ by solving $\mathbf{F}_{1}(\mathbf{x}\mathchar 24635\relax\;t)=0$ in line (5), (8) and (11) of Alg. 3. The average number of inner linear solver iterations is quite stable across all experiments. The most challenging case is the limit $\varepsilon\rightarrow 0$ in which we used the solution to the optimization problem as a warmed-up initial guess for the next smaller value of $\varepsilon$. Further, we emphasize that the top left blocks of $\mathbf{J}_{1,1}$ involve the matrices $C_{i}(\mathbf{y}^{k}_{i},t^{k})$ (see (33)) which contain a dense low-rank term if $g_{\varepsilon}^{\prime\prime}(Q^{k}_{\omega_{i}})\neq 0$. As $\varepsilon\rightarrow 0$, $g_{\varepsilon}^{\prime\prime}(\cdot)$ tends to a Dirac delta, so the dense term become dominant. Multigrid methods based on a pointwise relaxations are expected to be not very efficient for these matrices which may not be diagonally dominant. The standard V-cycle algorithm indeed suffers, however the Krylov acceleration performs better as it handles these low-rank perturbation with smaller effort. For $\varepsilon=10^{-4}$, we sometimes noticed that the GMRES residual stagnates after 20/30 iterations around $10^{-7}/10^{-8}$, due to a loss of orthogonality in the Krylov subspace, and thus resulting in higher number of iterations. We allowed a maximum number of 80 iterations per linear system.

Figure 2 compares the two optimal controls obtained minimizing either $\mathbb{E}\left[Q(y_{\omega})\right]$ or $\text{CVaR}_{0.99}\left[Q(y_{\omega})\right]$, and the cumulative distribution functions of $Q(y_{\omega_{j}})$ computed on 8000 out-of-sample realizations. The risk-adverse control indeed minimizes the risk of having large values of $Q(y_{\omega})$. The CVaR of level $\lambda=0.99$ is respectively $\text{CVaR}_{0.99}\left(Q(y_{\omega})\right)=2.79$ for the risk-neutral control and $\text{CVaR}_{0.99}\left(Q(y_{\omega})\right)=0.90$ for the risk-adverse control.

6 Conclusion

We have presented a multigrid method to solve the large saddle point linear systems that typically arise in full-space approaches to solve OCPUU. We tested the algorithm as a iterative solver and as a preconditioner on three test cases: a linear-quadratic OCPUU, a nonsmooth OCPUU, and a risk-adverse nonlinear OCPUU. Overall, the multigrid method shows very good performances and robustness with respect to the several parameters of the problems considered.

7 Declarations

Ethical Approval: Not applicable.

Data Availability: The codes used in this study are available from the corresponding author on reasonable request.

Competing interests: The authors declare that they have no conflict of interest.

Funding: Not applicable.

Authors’ contributions: All authors contributed equally to the manuscript.

Acknowledgements: Gabriele Ciaramella and Tommaso Vanzan are members of GNCS (Gruppo Nazionale per il Calcolo Scientifico) of INdAM. The present research is part of the activities of “Dipartimento di Eccellenza 2023-2027".