# Relaxation approach for learning regularizers by neural networks for a class
of optimal control problems††thanks: This work is supported by…

Sébastien Court Department of Mathematics & Digital Science Center, University
of Innsbruck, Austria, email<EMAIL_ADDRESS>Markus Haltmeier
Department of Mathematics, University of Innsbruck, Austria, email:
<EMAIL_ADDRESS>
###### Abstract

The present paper deals with the data-driven design of regularizers in the
form of artificial nerual networks, for solving inverse problems formulated as
optimal control problems. These regularizers aim at improving accuracy,
wellposedness or compensating uncertainties for a class of optimal control
problems (inner-problems). Parameterized as neural networks, their weights are
chosen in order to reduce a misfit between data and observations of the state
solution of the inner- optimal control problems. Learning these weights
constitutes the outer-problem. Based on necessary first-order optimality
conditions for the inner-problems, a relaxation approach is proposed in order
to implement efficient solving of the inner-problems, namely the forward
operator of the outer-problem. Optimality conditions are derived for the
latter, and numerical illustrations show the feasibility of the relaxation
approach, first for rediscovering standard $L^{2}$-regularizers, and next for
designing regularizers that compensate unknown noise on the observed state of
the inner-problem.

Keywords: Optimal control theory, Regularization, Neural networks,
Identification problem, Data-driven optimization, Learning problems.

AMS subject classifications (2020): 49N45, 49K99, 68T07, 68T05, 68T20, 65J20.

###### Contents

  1. 1 Introduction
    1. 1.1 General learning problem
    2. 1.2 Motivation of data-driven regularization and state-of-the-art
    3. 1.3 Approach and strategy
    4. 1.4 Plan
  2. 2 Notation, assumptions and preliminaries
    1. 2.1 Functional setting
    2. 2.2 The regularizer as realization of a neural network
    3. 2.3 General assumptions
    4. 2.4 On the control-to-state mapping for the inner-problem
  3. 3 Necessary optimality conditions for the inner-problem
    1. 3.1 Differential properties of neural networks
    2. 3.2 Direct approach
    3. 3.3 Optimal control formalism for the inner-problem ($\ast$)
  4. 4 Relaxation approach for the outer-problem using necessary optimality conditions of the inner-problems
    1. 4.1 Relaxed formulation of the main problem
    2. 4.2 Existence of minimizers for the relaxed problem
    3. 4.3 Derivatives of second-order for the neural network
    4. 4.4 Properties of the adjoint operator
      1. 4.4.1 When the adjoint operator is invertible
      2. 4.4.2 When the adjoint operator is not invertible
    5. 4.5 Necessary optimality conditions for the relaxed problem
  5. 5 A gradient-based algorithm
    1. 5.1 Nesterov gradient descent for solving the inner-problems
    2. 5.2 Barzilai-Borwein gradient rule for the relaxed outer-problem
  6. 6 Illustration: Identification of coefficients in an elliptic equation
    1. 6.1 Presentation of the problem
    2. 6.2 Rediscovering the $L^{2}$-norm, or not
    3. 6.3 Regularizers for compensating noise on the data
  7. 7 Conclusion
  8. A Appendix
    1. A.1 On the optimal control approach for solving an inverse problem
    2. A.2 Results related to the inverse conductivity problem

##  1 Introduction

Designing regularizers for solving inverse problems is an important data-
science problem for which data-driven approaches have been recently developed.
In particular, the design of regularizers as artificial neural networks
demonstrate interesting capacities for improving characteristics of inverse
problems. Among these characteristics, one can aim at recovering uniqueness of
solutions for initially ill-posed inverse problems, facilitating their
solving, obtaining specific qualitative properties like sparsity, or even
compensating unknown real-world noise that leads to uncertainties on the
measured data. The design of such regularizers can be achieved towards an
optimal control approach, consisting in optimizing parameters of a machine,
namely an input-output map, in order to minimize a distance between real-world
data and outputs of this machine.

###  1.1 General learning problem

Consider a series of $K$ problems that belong to a given class of optimal
control problems indexed by $1\leq k\leq K$, and assume that for each of them
we know a solution. More precisely, we consider a set of tasks $\hat{z}_{k}$
that we know to be achieved by the means of a control $\hat{u}_{k}$. The $K$
pairs $\left\\{(\hat{z}_{k},\hat{u}_{k}),1\leq k\leq K\right\\}$ constitutes
our data set, and our main goal is, from this data set, to train a regularizer
that improves the solving of this class of control problem. Our supervised-
learning problem can be formulated as follows:

$\left\\{\begin{array}[]{l}\displaystyle\min_{r\in\mathcal{R}}\quad\frac{1}{2K}\sum_{k=1}^{K}\|u_{k}-\hat{u}_{k}\|_{\mathcal{U}}^{2},\\\\[15.0pt]
\text{subject to: $u_{k}$ is solution of Problem~{}\eqref{pbstandard} with
$(\hat{z}_{k},r)$, for all $1\leq k\leq K$,}\end{array}\right.$ ($\ast\ast$)

where Problem $\eqref{pbstandard}$ with $(\hat{z},r)=(\hat{z}_{k},r)$ is given
by

$\left\\{\begin{array}[]{l}\displaystyle\min_{(y,u)\in\mathcal{Y}\times\mathcal{U}}\big{(}c(y,\hat{z})+\gamma\circ
r(u)\big{)},\\\\[10.0pt] \text{subject to: }\varphi(y,u)=0.\end{array}\right.$
($\ast$)

The set $\mathcal{R}$ describes the class of regularizers in which we design a
– possible – optimal function $r$. The state is denoted by $y$, and the state
equation is given by $\varphi(y,u)=0$. The cost $c(y,\hat{z})$ corresponds to
a misfit between the state $y$ and the task $\hat{z}$ to achieve. The mapping
$\gamma$ is given, with values in $\mathbb{R}$, and can correspond for
instance to a norm.

###  1.2 Motivation of data-driven regularization and state-of-the-art

In the formulation of the inner-problem ($\ast$ ‣ 1.1), the additional term
$r(u)$ corresponds to a so-called regularizer. Its role consists in forcing
the control $u$ to present certain properties, like regularity, boundedness,
or sparsity for instance. Besides the qualitative properties aforementioned,
designing a regularizer for a class of identification problems can address
several issues. For instance, inverse problems are well-known to present non-
uniqueness issues. A regularizer designed with a data-driven approach could
help to select a relevant solution. An other issue could be due to uncertainty
on measurements: Data can be subject to a noise phenomenon whose nature is
unknown, and that would lead to inaccuracy on the solution. An appropriate
regularizer, implementing implicitly the effects of this noise, could help to
compensate this lack of accuracy. We illustrate this last point in our
numerical realizations (see section 6).

In order to design such a regularizer, one need to parameterize the set
$\mathcal{R}$, that we choose as a set of artificial neural networks (ANNs)
with a prescribed number of layers. A regularizer $(u\mapsto
r(w,u))\in\mathcal{R}$ is then a mapping parameterized by a set of weights
$w$. We refer to section 3.1 for further details on the definitions and
notation related to the ANNs we will consider. Beside the theoretical
approximation capacities of ANNs (see [GPEB19] for instance), data-driven
methods training ANNs have lead to the design of regularizers that outperform
classical solving methods for inverse problems.

For example, the data-driven modeling of regularizers with NETT (Network
Tiknonov) was investigated in [LSAH20a], in order to solve inverse problems
subject to unknown noise, with applications to Computer Tomography. Let us
also cite [OSH19] in the context of trained decoder networks for sparse
reconstruction, and refer to [LSAH20b] for a theoretical analysis of such an
approach. Further developments were studied in [ONSH19, ONSH20], considering
aNETT (augmented Network Thikonov), designing a regularizer and also a penalty
(misfit) function in the form of a neural network. Still on applications in
medical imaging, the use of adversarial neural networks as regularizers was
studied in [LOS18]. The use of deep neural networks for solving ill-posed
inverse problems related to medical imaging was also studied in [AO17]. The
special case of input-convex neural networks (ICNNs) for modeling regularizers
with a data-driven approach was considered in [MDS+20]. More generally, we
orient the reader to [HN20] for a comprehensive review on techniques related
to this topic.

###  1.3 Approach and strategy

In [HKB18] the authors investigated the problem of learning the Thikhonov
regularization parameter for a class of inverse problems, with a bi-level
approach. This approach consists in adopting a hierarchy in the formulation of
the regularizer learning problem: The identification problems ($\ast$ ‣ 1.1)
for which we want to design a regularizer are called inner-problems, and the
main problem ($\ast\ast$ ‣ 1.1) consists in optimizing the parameters of the
regularizer for a given set of identification problems for which we know the
ground-truth (the data set). The bi-level approach was initially developed for
image denoising in [KP13].

In the present paper, we adopt this bi-level approach, where the inner-problem
($\ast$ ‣ 1.1) is considered as a constraint of the outer-problem. Compared
with [HKB18], in our case the range of possible regularizers is a priori wider
than a Thikonov regularizer, as the regularizers are chosen in the form of
artificial neural networks parameterized with weights. The price is to pay is
a theoretical analysis (existence and uniqueness of minimizers) that becomes
out-of-reach in the general case, due in particular to the purely nonlinear
aspects of neural networks. Therefore, our main contribution is the
development of a relaxation approach, consisting in replacing the optimality
constraint of the inner-problems by the corresponding necessary Karush-Kuhn-
Tucker conditions. The latter are equality constraints (vanishing gradient
equalities) that transform the outer-problem into a standard optimal control
problem of state constraint type.

This relaxation approach is inspired by optimal control problems with partial
differential equations constraints. Indeed, partial differential equations
modeling physical dynamical systems can often be derived from the least action
principle. They are merely first-order optimality conditions of a minimization
problem. Wellposedness questions for these partial differential equations are
then related to existence and uniqueness of minimizers for functionals called
actions.

The main theoretical result is given in Theorem 2, where we provide the
sensitivity of the objective function of Problem ($\ast\ast$ ‣ 1.1)111More
precisely the sensitivity of the objective function of the relaxed problem
($\tilde{\ast\ast}$ ‣ 4.1) introduced in section 4, but the expression of this
objective function is the same as the one of Problem ($\ast\ast$ ‣ 1.1). with
respect to the weights of the neural network. This allows us to implement the
determination of these weights in Algorithm 1, with gradient rules.
Illustrations are proposed by considering as inner-problems the classical
identification of coefficients in elliptic partial differential equations.
Regularizers are designed, in particular for compensating the effects of a
Gaussian noise introduced on the state variable.

###  1.4 Plan

The paper is organized as follows: Notation and assumptions are given in
section 2. In section 3 we study the inner-problems with ANNs as regularizers.
The latter are described in section 3.1. The necessary optimality conditions
related to the inner-problems are investigated in sections 3.2 and 3.3. The
relaxation approach is presented in section 4.1. We prove existence of
minimizers for the relaxed outer-problem in section 4.2, and derive optimality
conditions in section 4.5 for the latter. In section 5 we present algorithms
for implementing the relaxation approach. Section 6 is devoted to numerical
illustrations considering the identification of coefficients in an elliptic
partial differential equation. Concluding remarks are given in section 7. In
Appendix A.1 we present advantages that the optimal control approach can offer
when dealing with state-constrained identification problems. Appendix A.2
recalls classical results related to the inverse conductivity problem.

##  2 Notation, assumptions and preliminaries

###  2.1 Functional setting

We denote by $\mathcal{Y}$, $\mathcal{U}$, $\mathcal{Q}$ real Banach spaces,
and by $\varphi:\mathcal{Y}\times\mathcal{U}\rightarrow\mathcal{Q}^{\prime}$
the state mapping. We assume that $\mathcal{Q}$ is reflexive. The control
space $\mathcal{U}$ is assumed to be finite-dimensional, of dimension
$n_{\mathrm{in}}$, so that $\mathcal{U}\simeq\mathbb{R}^{n_{\mathrm{in}}}$.
Therefore in particular $\mathcal{U}$ is a Hilbert space, and reflexive. The
range space $\mathcal{V}$ of the neural networks is finite-dimensional, of
dimension $n_{\mathrm{out}}$, and
$\mathcal{V}\simeq\mathbb{R}^{n_{\mathrm{out}}}$. The space of weights is
denoted by $\mathcal{W}$, and described in section 3.1. Given weights
$w\in\mathcal{W}$, we denote by $r(w,\cdot):\mathcal{U}\rightarrow\mathcal{V}$
the associated neural network of $\mathcal{R}$, and by
$\gamma:\mathcal{V}\rightarrow\mathbb{R}$ a smooth given mapping, possibly
chosen as a norm on $\mathbb{R}^{n_{\mathrm{out}}}$.

The data $\left\\{(\hat{z}_{k},\hat{u}_{k}),1\leq k\leq K\right\\}$ are of
finite number. The controls $(\hat{u}_{k})_{1\leq k\leq K}$ are assumed to lie
in $\mathcal{U}$, while the tasks $(\hat{z}_{k})_{1\leq k\leq K}$ lie in a
finite-dimensional subspace of $\mathcal{Y}$.

###  2.2 The regularizer as realization of a neural network

Given $L\geq 1$, we denote by $w=(w_{1},w_{2},\dots,w_{L})$ a vector of affine
functions. The set of regularizers, denoted by $\mathcal{R}$, is made of
feedforward neural networks of $L$ layers, for which the dimensions $n_{\ell}$
($1\leq\ell\leq L$) of the hidden layers are prescribed:

$\displaystyle\mathcal{R}$ $\displaystyle=$
$\displaystyle\left\\{\mathcal{U}\ni u\mapsto
w_{L}\rho\left(w_{L-1}\rho\left(\left(\dots
w_{2}\rho(w_{1}u)\right)\right)\right)\in\mathcal{V}\mid\
w=(w_{1},\dots,w_{L})\in\prod_{\ell=1}^{L}\mathrm{Aff}(n_{\ell},n_{\ell+1},\mathbb{R})\right\\}$

We denoted by $\mathrm{Aff}(n_{\ell},n_{\ell+1},\mathbb{R})$ the set of affine
functions from $\mathbb{R}^{n_{\ell}}$ to $\mathbb{R}^{n_{\ell+1}}$. For
$1\leq\ell\leq L$, its components write

$w_{\ell}:u\mapsto A_{\ell}u+b_{\ell},$

where $A_{\ell}\in\mathbb{R}^{n_{\ell+1}\times n_{\ell}}$ and
$b_{\ell}\in\mathbb{R}^{n_{\ell+1}}$, with
$n_{1}=n_{\mathrm{in}}=\mathrm{dim}(\mathcal{U})$ and
$n_{L+1}=n_{\mathrm{out}}=\mathrm{dim}(\mathcal{V})$. For the sake of
concision, we will denote throughout the paper

$\mathcal{W}:=\prod_{\ell=1}^{L}\mathrm{Aff}(n_{\ell},n_{\ell+1},\mathbb{R}).$

We introduce
$\mathcal{L}_{\mathcal{W}}:=\displaystyle\sum_{\ell=1}^{L-1}n_{\ell}$. We
endow the space $\mathcal{W}$ with the Euclidean norm given by

$\|w\|^{2}_{\mathcal{W}}:=\sum_{\ell=1}^{L}\left(|A_{\ell}|^{2}_{\mathbb{R}^{n_{\ell+1}\times
n_{\ell}}}+|b_{\ell}|^{2}_{\mathbb{R}^{n_{\ell+1}}}\right),$

where $|A_{\ell}|^{2}_{\mathbb{R}^{n_{\ell+1}\times
n_{\ell}}}=\mathrm{trace}(A_{\ell}^{T}A_{\ell})$ and $|b_{\ell}|$ is the
Euclidean norm of $\mathbb{R}^{n_{\ell+1}}$. Given a smooth activation
function $\rho:\mathbb{R}\rightarrow\mathbb{R}$ and a vector $w$ of $L$
weights, a neural network $r(w,\cdot)\in\mathcal{R}$ of $L$ layers writes as
follows:

$\displaystyle r(w,u)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}w_{2}(\rho(w_{1}(u)))&\text{if
}L=2,\\\\[5.0pt] w_{L}(\rho(w_{L-1}(\rho(\dots\rho(w_{1}(u))))))&\text{if
}L\geq 3.\end{array}\right.$

In the expression above, since the functions $w_{\ell}$ have values in
$\mathbb{R}^{n_{\ell+1}}$, the compositions by $\rho$ have to be understood
coordinate-wise as
$\rho(w_{\ell}(\mathrm{z}))_{i}=\rho((w_{\ell}(\mathrm{z}))_{i})$, for
$\mathrm{z}\in\mathbb{R}^{n_{\ell}}$ and $i\in\\{1,\dots,n_{\ell+1}\\}$. Thus,
the regularizer space $\mathcal{R}$ is parameterized by the weights lying in
$\mathcal{W}$, and the outer-problem ($\ast\ast$ ‣ 1.1) is formulated as
follows:

$\left\\{\begin{array}[]{l}\displaystyle\min_{w\in\mathcal{W}}\quad\frac{1}{2}\sum_{k=1}^{K}\|\hat{u}_{k}-u_{k}\|_{\mathcal{U}}^{2},\\\\[15.0pt]
\text{subject to: $u_{k}$ is solution of Problem~{}\eqref{pbstandard} with
$(\hat{z},r)=\big{(}\hat{z}_{k},r(\cdot,w)\big{)}$, for all $1\leq k\leq
K$.}\end{array}\right.$ ($\ast\ast$)

###  2.3 General assumptions

Throughout the paper, we will assume that the following set of assumptions
holds:

$(\mathbf{A1})$

The mapping $\varphi$ is of class $C^{2}$ over $\mathcal{Y}\times\mathcal{U}$.
For all $u\in\mathcal{U}$, there exists $y\in\mathcal{Y}$ such that
$\varphi(y,u)=0$.

$(\mathbf{A2})$

For all $(\overline{y},\overline{u})\in\mathcal{Y}\times\mathcal{U}$, the
linear mapping
$\varphi^{\prime}_{y}(\overline{y},\overline{u}):\mathcal{Y}\rightarrow\mathcal{Q}^{\prime}$
is surjective. Moreover, there exists a constant
$C(\overline{y},\overline{u})>0$, depending only on
$(\overline{y},\overline{u})$ such that for all $y\in\mathcal{Y}$ the
following estimate holds:

$\displaystyle\|y\|_{\mathcal{Y}}$ $\displaystyle\leq$ $\displaystyle
C\left\|\varphi^{\prime}_{y}(\overline{y},\overline{u}).y\right\|_{\mathcal{Q}^{\prime}}.$

Instead of $(\mathbf{A2})$, we may consider a weaker assumption, namely
$(\mathbf{A2})^{\prime}$:

$(\mathbf{A2})^{\prime}$

For all $(\overline{y},\overline{u})\in\mathcal{Y}\times\mathcal{U}$, the
linear mapping
$\varphi^{\prime}_{y}(\overline{y},\overline{u}):\mathcal{Y}\rightarrow\mathcal{Q}^{\prime}$
is almost surjective, namely its range is dense in $\mathcal{Q}^{\prime}$.

$(\mathbf{A3})$

When $r\equiv 0$, Problem ($\ast$ ‣ 1.1) admits at least one global minimizer.

$(\mathbf{A4})$

The cost function $c:\mathcal{Y}\times\mathcal{U}$ is twice continuously
Fréchet-differentiable.

$(\mathbf{A5})$

The activation function $\rho$ is of class $C^{2}$ on $\mathbb{R}$.

Assumption $(\mathbf{A1})$ combined with $(\mathbf{A2})$ imply the existence
of a control-to-state ammping for the inner-problem (see section 2.4).
Assumption $(\mathbf{A3})$ is natural, as we aim at designing regularizers for
(inner-)problems that are solvable. From $(\mathbf{A3})$ we deduce easily
existence of minimizers when the regularizer is non-trivial (Proposition 2).
While addressing the inner-problems only requires first-order
differentiability (see section 3), the derivation of optimality conditions for
the outer-problem (see section 4) necessitates second-order derivatives, in
particular for the objective function of the inner-problem, including the cost
function $c$ (Assumption $(\mathbf{A4})$) and the neural network (Assumption
$(\mathbf{A5})$). Throughout the paper these assumptions will be implicitly
assumed.

###  2.4 On the control-to-state mapping for the inner-problem

The state constraint is given by the implicit function
$\varphi:\mathcal{Y}\times\mathcal{U}\rightarrow\mathcal{Q}^{\prime}$. Under
assumptions $(\mathbf{A1})$-$(\mathbf{A2})$, we are able to define the
control-to-state mapping that we denote by $\mathbb{S}$:

######  Proposition 1.

Assume $(\mathbf{A1})$-$(\mathbf{A2})$. There exists a mapping
$\mathbb{S}:\mathcal{U}\rightarrow\mathcal{Y}$ of class $C^{2}$ such that for
all $(y,u)\in\mathcal{Y}\times\mathcal{U}$ we have

$\displaystyle\varphi(y,u)=0$ $\displaystyle\Leftrightarrow$ $\displaystyle
y=\mathbb{S}(u).$

Furthermore, $\mathbb{S}$ is of class $C^{2}$ over $\mathcal{U}$.

###### Proof.

The result is due to the implicit function theorem. ∎

Depending on the context, for the sake of convenience we may consider the
original equality constraint $\varphi(y,u)=0$, or substitute the variable $y$
by $\mathbb{S}(u)$. The lack of control-to-state mapping $\mathbb{S}$ for the
inner-problem would introduce further difficulties that we choose to not
address in this paper.

##  3 Necessary optimality conditions for the inner-problem

Given $w\in\mathcal{W}$, when the regularizer is chosen in the form of a
neural network as described in section 3.1, the inner-problems correspond to
the following class of optimal control problems:

$\left\\{\begin{array}[]{l}\displaystyle\min_{(y,u)\in\mathcal{Y}\times\mathcal{U}}\Big{(}J(y,u,w):=c(y,\hat{z})+\gamma\circ
r(w,u)\Big{)},\\\\[10.0pt] \text{subject to: }\varphi(y,u)=0\text{ in
}\mathcal{Q}^{\prime}.\end{array}\right.$ ($\ast$)

The aim of this section is to give necessary optimality conditions for the
inner-problem ($\ast$ ‣ 1.1). For that purpose, we first need to investigate
the differential properties owned by neural networks of $\mathcal{W}$.

###  3.1 Differential properties of neural networks

Consider the partial neural network with $\ell$ layers ($1\leq\ell\leq L$):

$\displaystyle r^{(\ell)}(w^{(\ell)},u)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}w_{1}(u)&\text{if }\ell=1,\\\
w_{2}(\rho(w_{1}(u)))&\text{if }\ell=2,\\\\[5.0pt]
w_{\ell}(\rho(w_{\ell-1}(\rho(\dots\rho(w_{1}(u))))))&\text{if }\ell\geq
3.\end{array}\right.$

Recall that the weights $(w_{\ell})$ are affine functions of
$\mathbb{R}^{n_{\ell}}$, of the form
$w_{\ell}(\mathrm{z})=A_{\ell}\mathrm{z}+b_{\ell}$ with
$A_{\ell}\in\mathbb{R}^{n_{\ell+1}\times n_{\ell}}$ and
$b_{\ell}\in\mathbb{R}^{n_{\ell+1}}$. We have denoted

$w^{(\ell)}=(w_{1},\dots,w_{\ell})\in\prod_{k=1}^{\ell}\mathrm{Aff}(\mathbb{R}^{n_{k}},\mathbb{R}^{n_{k+1}}).$

Observe that for $1\leq\ell\leq L-1$ the following formula holds:

$\displaystyle r^{(\ell+1)}(w^{(\ell+1)},u)$ $\displaystyle=$ $\displaystyle
w_{\ell+1}\rho\big{(}r^{(\ell)}(w^{(\ell)},u)\big{)}.$ (3)

From (3), by induction we can verify that the sensitivity of $r^{(L)}$ with
respect to $u$ satisfies the identity

$\frac{\partial r^{(L)}}{\partial
u}(w^{(L)},u)=A_{L}\mathrlap{\prod_{\ell=1}^{L-1}}{\longrightarrow}\rho^{\prime}(r^{(\ell)}(w^{(\ell)},u))A_{\ell}.$
(4)

The vector/matrix product $\rho^{\prime}(r^{(\ell)}(w^{(\ell)},u))A_{\ell}$ is
a matrix, and has to be understood as follows:

$(\rho^{\prime}\,A)_{ij}=\rho^{\prime}_{i}\,A_{ij}.$ (5)

The symbol $\displaystyle\mathrlap{\prod_{\ell=1}^{L-1}}{\longrightarrow}$ is
used for calculating the non-commutative product of matrices, by multiplying
them iteratively to the left when the index $\ell$ increases:
$\displaystyle\mathrlap{\prod_{\ell=1}^{L}}{\longrightarrow}B_{\ell}=B_{L}B_{L-1}\dots
B_{2}B_{1}$. Without ambiguity, we will simply denote $r=r^{(L)}$ and
$w=w^{(L)}$. Using the mean-value theorem, by induction we can estimate

$\displaystyle\|r(w,u)-r(w,0)\|_{\mathcal{V}}$ $\displaystyle\leq$
$\displaystyle\left(\prod_{\ell=1}^{L}|A_{\ell}|_{\mathbb{R}^{n_{\ell+1}\times
n_{\ell}}}\right)\|\rho^{\prime}\|_{\infty}^{\mathcal{L}_{\mathcal{W}}}\|u\|_{\mathcal{U}},$
$\displaystyle\|r(w,u_{1})-r(w,u_{2})\|_{\mathcal{V}}$ $\displaystyle\leq$
$\displaystyle\left(\prod_{\ell=1}^{L}|A_{\ell}|_{\mathbb{R}^{n_{\ell+1}\times
n_{\ell}}}\right)\|\rho^{\prime}\|_{\infty}^{\mathcal{L}_{\mathcal{W}}}\|u_{1}-u_{2}\|_{\mathcal{U}},$

where we recall
$\mathcal{L}_{\mathcal{W}}=\displaystyle\sum_{\ell=1}^{L-1}n_{\ell}$. The
differentiation of $r$ with respect to the weights $w_{\ell}$ is given by the
following formula, for all
$\tilde{w}\in\mathrm{Aff}(\mathbb{R}^{n_{\ell}},\mathbb{R}^{n_{\ell+1}})$:

$\frac{\partial r}{\partial
w_{\ell}}(w,u).\tilde{w}=A_{L}\left(\mathrlap{\prod_{k=\ell}^{L-1}}{\hskip
1.0pt\longrightarrow}\rho^{\prime}(r^{(k)}(w^{(k)},u))A_{k}\right)\left(\tilde{w}(\rho(r^{(\ell-1)}(w^{(\ell-1)},u)))\right).$

We can obtain this formula also by induction, using (3). Here the product has
to be understood as follows for $1\leq i\leq n_{\mathrm{out}}$:

$\displaystyle\left(\frac{\partial r}{\partial
w_{\ell}}(w,u).\tilde{w}\right)_{i}$ $\displaystyle=$
$\displaystyle\left(A_{L}\right)_{ij}\left(\mathrlap{\prod_{k=\ell}^{L-1}}{\longrightarrow}\left(\rho^{\prime}(r^{(k)}(w^{(k)},u))A_{k}\right)\right)_{jh}\left(\tilde{w}(\rho\left(r^{(\ell-1)}(w^{(\ell-1)},u)\right))\right)_{h}.$

These expressions for the derivatives of the neural network will be used in
practice for the implementation of the gradient associated with the inner-
problem. Let us now derive expressions for the latter.

###  3.2 Direct approach

The control-to-state mapping $\mathbb{S}$ obtained in Proposition 1 enables us
to derive classically necessary optimality conditions for the inner-problems
directly by omitting the state constraint. Indeed, using $y=\mathbb{S}(u)$,
and given $w\in\mathcal{W}$, we can rewrite Problem ($\ast$ ‣ 1.1) as

$\min_{u\in\mathcal{U}}\Big{(}\tilde{J}(u,w):=J(\mathbb{S}(u),u,w)=c(\mathbb{S}(u),\hat{z})+\gamma(r(w,u))\Big{)}.$

Using Assumption $(\mathbf{A3})$ and the formulation above, we obtain
existence of minimizers for the inner-problem ($\ast$ ‣ 3.3):

######  Proposition 2.

Under Assumption $(\mathbf{A3})$, for all $w\in\mathcal{W}$ Problem ($\ast$ ‣
3.3) admits a global minimizer.

###### Proof.

From Assumption $(\mathbf{A3})$ and Proposition 1, the mapping $u\mapsto
c(\mathbb{S}(u),\hat{z})$ admits global minimizers, and due to the regularity
of $u\mapsto\gamma(r(w,u))$, the result follows straightforwardly. ∎

We next obtain necessary optimality conditions involving the derivative of
$\mathbb{S}$:

######  Proposition 3.

Let be $w\in\mathcal{W}$. If
$(y=\mathbb{S}(u),u)\in\mathcal{Y}\times\mathcal{U}$ is solution of Problem
($\ast$ ‣ 1.1), then the gradient

$G(u,w,\hat{z}):=\frac{\partial\tilde{J}}{\partial
u}(u,w)=\mathbb{S}^{\prime}(u)^{\ast}c_{y}(\mathbb{S}(u),\hat{z})+r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))$
(6)

vanishes.

###### Proof.

Following the previous comments and Proposition 1, the result follows from the
Karush-Kuhn-Tucker conditions. ∎

In certain problems, like shape optimization for example, evaluating
efficiently values of the control-to-state mapping $\mathbb{S}(u)$, and in
particular $\mathbb{S}^{\prime}(u)^{\ast}$, can be challenging numerically, or
expensive computationally. Therefore we prefer sometimes to adopt an optimal
control approach, also called the adjoint approach, that we develop in our
context in section 3.3 below. See remarks in [HPUU09, section 1.6.1, page 59]
for further explanations. Further comments on the possible advantages of the
optimal control approach are given in Appendix A.1.

###  3.3 Optimal control formalism for the inner-problem ($\ast$ ‣ 1.1)

Let us recall the initial description of the inner-problem, when we keep the
original state constraint:

$\left\\{\begin{array}[]{l}\displaystyle\min_{(y,u)\in\mathcal{Y}\times\mathcal{U}}\Big{(}J(y,u,w):=c(y,\hat{z})+\gamma\circ
r(w,u)\Big{)},\\\\[10.0pt] \text{subject to: }\varphi(y,u)=0\text{ in
}\mathcal{Q}^{\prime}.\end{array}\right.$ ($\ast$)

Unlike the direct approach which utilizes $\varphi(y,u)=0\Leftrightarrow
y=\mathbb{S}(u)$, the optimal control approach – also called the adjoint
approach – consists in taking into account the state constraint by introducing
a Lagrange multiplier. We define the Lagrangian of problem ($\ast$ ‣ 3.3) as
follows

$\displaystyle\begin{array}[]{rrcl}\mathscr{L}:&(\mathcal{Y}\times\mathcal{U})\times\mathcal{Q}&\rightarrow&\mathbb{R}\\\
&((y,u),p)&\mapsto&c(y,\hat{z})+\displaystyle\gamma\circ r(w,u)+\left\langle
p\,,\varphi(y,u)\right\rangle_{\mathcal{Q},\mathcal{Q}^{\prime}}.\end{array}$

Problem ($\ast$ ‣ 3.3) can be solved by finding a saddle-point to the
Lagrangian $\mathscr{L}$, minimized with respect to variables $(y,u)$ and
maximizing with respect to $p$. The sensitivity of $\mathscr{L}$ with respect
to variable $y$ yields the adjoint equation, whose unknown is the adjoint
state $p\in\mathcal{Q}$. It is given in $\mathcal{Y}^{\prime}$ by

$\displaystyle\varphi^{\prime}_{y}(y,u)^{\ast}.p+c^{\prime}_{y}(y,\hat{z})=0.$
(8)

We define a solution of system (8) by transposition.

######  Definition 1.

We say that $p$ is a solution of system (8) if for all
$f\in\mathcal{Q}^{\prime}$ we have

$\displaystyle\langle
p\,,f\rangle_{\mathcal{Q},\mathcal{Q}^{\prime}}+c^{\prime}_{y}(y,\hat{z}).\tilde{y}=0,$
(9)

where $\tilde{y}\in\mathcal{Y}$ denotes the solution of the following linear
system in $\mathcal{Q}^{\prime}$:

$\displaystyle\varphi^{\prime}_{y}(y,u).\tilde{y}=f.$ (10)

From assumption $(\mathbf{A2})$, system (10) is well-posed. The existence of
solutions for adjoint system (8), in the sense of Definition 1, is stated as
follows:

######  Proposition 4.

Let be $(y,u)\in\mathcal{Y}\times\mathcal{U}$. There exists a unique solution
$p\in\mathcal{Q}$ to system (8), in the sense of Definition 1, and there
exists a constant $C(y,u)>0$ (depending only on $(y,u)$) such that

$\displaystyle\|p\|_{\mathcal{Q}}$ $\displaystyle\leq$ $\displaystyle
C(y,u)\|c^{\prime}_{y}(y,\hat{z})\|_{\mathcal{Y}^{\prime}}.$

###### Proof.

From assumption $(\mathbf{A2})^{\prime}$, it is easy to verify that the
mapping $\varphi_{y}(y,u)^{\ast}$ is injective, and so uniqueness holds. Let
us prove existence. Denote by
$\Lambda(y,u):\mathcal{Q}^{\prime}\rightarrow\mathcal{Y}$ the linear operator
which maps $f$ to $\tilde{y}$, where $\tilde{y}$ is the solution of (10). From
Assumption $(\mathbf{A2})$, the operator $\Lambda(y,u)$ is well-defined and
bounded, and consequently
$\Lambda(y,u)^{\ast}:\mathcal{Y}^{\prime}\rightarrow\mathcal{Q}^{\prime\prime}\simeq\mathcal{Q}$
too. Now define $p=-\Lambda(y,u)^{\ast}c^{\prime}_{y}(y,\hat{z})$. Then, for
all $f\in\mathcal{Q}^{\prime}$, we verify that

$\displaystyle\langle p\,,f\rangle_{\mathcal{Q}^{\prime},\mathcal{Q}}$
$\displaystyle=$ $\displaystyle-\langle
c^{\prime}_{y}(y,\hat{z})\,,\Lambda(y,u)f\rangle_{\mathcal{Y}^{\prime},\mathcal{Y}}$
$\displaystyle=$ $\displaystyle-\langle
c^{\prime}_{y}(y,\hat{z})\,,z\rangle_{\mathcal{Y}^{\prime},\mathcal{Y}}=-c^{\prime}_{y}(y,\hat{z}).z,$

and the announced estimate follows from boundedness of $\Lambda(y,u)^{\ast}$,
which completes the proof, as uniqueness is straightforward.∎

Consequently, the adjoint approach provides another expression for the
gradient associated with the inner-problem, namely the quantity
$G(u,w,\hat{z})$ introduced in (6).

######  Proposition 5.

Let be $w\in\mathcal{W}$. If $u$ is solution to Problem ($\ast$ ‣ 1.1), then
the gradient $G(u,w,\hat{z})$ introduced in Proposition 3 writes

$G(u,w,\hat{z})=r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))+\varphi^{\prime}_{u}(y,u)^{\ast}.p,$
(11)

where $y=\mathbb{S}(u)$, $p$ is the solution of (8) corresponding to $(y,u)$,
and necessarily $G(u,w)=0$.

###### Proof.

From Proposition 1, Problem ($\ast$ ‣ 1.1) reduces to
$\displaystyle\min_{u\in\mathcal{U}}J(\mathbb{S}(u),u,w)$. Following the
Karush-Kuhn-Tucker conditions, if $u$ is optimal to Problem ($\ast$ ‣ 1.1),
then $\frac{\partial}{\partial u}\left(J(\mathbb{S}(u),u,w)\right)=0$, with
for all $\tilde{u}\in U$

$\displaystyle\frac{\partial}{\partial
u}\left(J(\mathbb{S}(u),u,w)\right).\tilde{u}$ $\displaystyle=$ $\displaystyle
J^{\prime}_{y}(\mathbb{S}(u),u,w).(\mathbb{S}^{\prime}(u).\tilde{u})+J^{\prime}_{u}(\mathbb{S}(u),u,w).\tilde{u}$
(12) $\displaystyle=$ $\displaystyle\langle
c^{\prime}_{y}(y,\hat{z})\,,\mathbb{S}^{\prime}(u).\tilde{u}\rangle_{\mathcal{U}}+\langle\gamma^{\prime}(r(w,u))\,,r^{\prime}_{u}(w,u).\tilde{u}\rangle_{\mathcal{V}}$
$\displaystyle=$
$\displaystyle\langle\mathbb{S}^{\prime}(u)^{\ast}c^{\prime}_{y}(y,\hat{z})+r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))^{\ast}\,,\tilde{u}\rangle_{\mathcal{U}}$
$\displaystyle=$
$\displaystyle\langle-\mathbb{S}^{\prime}(u)^{\ast}\varphi^{\prime}_{y}(y,u)^{\ast}.p+r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))\,,\tilde{u}\rangle_{\mathcal{U}},$

where $y=\mathbb{S}(u)$ and $p$ satisfies (8). By differentiating the identity
$\varphi(\mathbb{S}(u),u)=0$, for all $\tilde{u}\in\mathcal{U}$ we obtain

$\varphi^{\prime}_{u}(y,u).\tilde{u}+\varphi_{y}^{\prime}(y,u).(\mathbb{S}^{\prime}(u).\tilde{u})=0,$
(13)

which implies that

$-\mathbb{S}^{\prime}(u)^{\ast}\varphi^{\prime}_{y}(y,u)^{\ast}=\varphi^{\prime}_{u}(y,u)^{\ast},$

and thus (12) reduces to

$\frac{\partial}{\partial
u}\left(J(\mathbb{S}(u),u,w)\right).\tilde{u}=\langle\varphi^{\prime}_{u}(y,u)^{\ast}.p+r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))\,,\tilde{u}\rangle_{\mathcal{U}^{\prime},\mathcal{U}}=0$

for all $\tilde{u}\in\mathcal{U}$, which completes the proof. ∎

######  Remark 1.

One could also have obtained the necessary conditions of Proposition 5 by
considering a stationary point of the Lagrangian functional $\mathscr{L}$.
Indeed, vanishing the derivative of $\mathscr{L}$ with respect to $p$ yields
$\varphi(y,u)=0\Leftrightarrow y=\mathbb{S}(u)$, in virtue of Proposition 1.
Vanishing the derivative of $\mathscr{L}$ with respect to $y$ leads to the
adjoint system (8). With the chain rule, we calculate for all
$\tilde{u}\in\mathcal{U}$:

$\displaystyle\frac{\partial\mathscr{L}}{\partial u}(w,u).\tilde{u}$
$\displaystyle=$
$\displaystyle\langle\gamma^{\prime}(r(w,u))\,,r_{u}^{\prime}(w,u).\tilde{u}\rangle_{\mathcal{V}}+c^{\prime}_{y}(\mathbb{S}(u),\hat{z}).(\mathbb{S}^{\prime}(u).\tilde{u})$
(14) $\displaystyle=$ $\displaystyle\langle
r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))\,,\tilde{u}\rangle_{\mathcal{U}}-\langle\varphi^{\prime}_{y}(\mathbb{S}(u),u)^{\ast}p\,,\mathbb{S}^{\prime}(u).\tilde{u}\rangle_{\mathcal{Y}^{\prime},\mathcal{Y}}$
$\displaystyle=$ $\displaystyle\langle
r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))\,,\tilde{u}\rangle_{\mathcal{U}}-\langle
p\,,\varphi^{\prime}_{y}(\mathbb{S}(u),u).(\mathbb{S}^{\prime}(u).\tilde{u})\rangle_{\mathcal{U}},$

where $p$ is solution of (8) with $y=\mathbb{S}(u)$. Using (13), we see that
(14) reduces to

$\frac{\partial\mathscr{L}}{\partial u}(w,u)=\langle
r^{\prime}_{u}(w,u)^{\ast}\gamma^{\prime}(r(w,u))\,,\tilde{u}\rangle_{\mathcal{U}}+\langle
p\,,\varphi^{\prime}_{u}(y,u).\tilde{u}\rangle_{\mathcal{U}^{\prime},\mathcal{U}}=\langle
G(u,w,\hat{z}),\tilde{u}\rangle_{\mathcal{U}^{\prime},\mathcal{U}},$

which leads to the same result as Proposition 5. Therefore the identities
given in Proposition 5 can be summarized as

$\left(\frac{\partial\mathscr{L}}{\partial
y}(y,u,p),\frac{\partial\mathscr{L}}{\partial
u}(y,u,p),\frac{\partial\mathscr{L}}{\partial p}(y,u,p)\right)=0.$ (15)

The necessary optimality conditions given in Proposition 3 are self-
consistent, and thus more appropriate for theoretical analysis, while those
provided by Proposition 5 are more convenient more numerical realization.
Following these two results, the constraint of Problem ($\ast\ast$ ‣ 1.1) is
then relaxed, by being simply replaced by the identity $G(u,w,\hat{z})=0$.

##  4 Relaxation approach for the outer-problem using necessary optimality
conditions of the inner-problems

###  4.1 Relaxed formulation of the main problem

The aim of the present section is to translate the abstract constraint of
Problem ($\ast\ast$ ‣ 1.1), namely ”$u_{k}$ is solution of Problem ($\ast$ ‣
1.1)[…]”, into an equality constraint. For that purpose, we use the necessary
vanishing gradient condition given in Propositions 3 and 5. Note that a
solution to Problem ($\ast$ ‣ 1.1) with $(\hat{z},r)=(\hat{z}_{k},r(w,\cdot))$
satisfies necessarily $G(u_{k},w,\hat{z}_{k})=0$, but if it satisfies this
identity only, then in general this is not necessarily a solution of Problem
($\ast$ ‣ 1.1). Therefore we consider the relaxed version of ($\ast\ast$ ‣
1.1) below:

$\left\\{\begin{array}[]{l}\displaystyle\min_{w\in\mathcal{W}}\quad\frac{1}{2K}\sum_{k=1}^{K}\|u_{k}-\hat{u}_{k}\|_{\mathcal{U}}^{2}+\frac{\nu}{2}\|w\|^{2}_{\mathcal{W}},\\\\[15.0pt]
\text{subject to: $G(u_{k},w,\hat{z}_{k})=0$ for all $1\leq k\leq
K$.}\end{array}\right.$ ($\tilde{\ast\ast}$)

The parameter $\nu\geq 0$ is given, and introduced for theoretical purpose
mainly (see Theorem 1). This regularization term for the set of weights $w$
could also facilitate the numerical solving of Problem ($\tilde{\ast\ast}$ ‣
4.1), by forcing the weights to be bounded. Note that a parameter-to-solution
mapping $w\mapsto u$ characterizing $G(u,w,\hat{z})=0$ is not necessarily
available. Even if such a mapping exists, its derivative is difficult to
describe explicitly. See section 4.4 for more details. For taking into account
the equality constraints of problem ($\tilde{\ast\ast}$ ‣ 4.1), we adopt an
adjoint approach and introduce the multipliers $\mu\in\mathcal{U}$, with the
following Lagrangian mapping:

$\displaystyle\begin{array}[]{rccl}\mathbf{L}:&\mathcal{U}\times\mathcal{W}\times\mathcal{U}&\rightarrow&\mathbb{R}\\\\[5.0pt]
&(u,w,\mu)&\mapsto&\displaystyle\frac{1}{2}\|u-\hat{u}\|_{\mathcal{U}}^{2}+\frac{\nu}{2}\|w\|^{2}_{\mathcal{W}}+\left\langle\mu\,,G(u,w,\hat{z})\right\rangle_{\mathcal{U},\mathcal{U}^{\prime}}\end{array}$
(18)

Like in section 3.3 for the inner-problems (see Remark 1), a solution of
Problem ($\tilde{\ast\ast}$ ‣ 4.1) can be sought – when $K=1$ – as a saddle-
point of $\mathbf{L}$ with respect to the variables $(u,w)$ and $\mu$. In
practice, we shall consider $(u_{k})_{1\leq k\leq K}$, $\mu=(\mu_{k})_{1\leq
k\leq K}$ and

$\displaystyle\begin{array}[]{rccl}\mathbf{L}^{(K)}:&\mathcal{U}^{K}\times\mathcal{W}\times\mathcal{U}^{K}&\rightarrow&\mathbb{R}\\\\[5.0pt]
&(u,w,\mu)&\mapsto&\displaystyle\frac{1}{2K}\sum_{k=1}^{K}\|u_{k}-\hat{u}_{k}\|_{\mathcal{U}}^{2}+\frac{\nu}{2}\|w\|^{2}_{\mathcal{W}}+\frac{1}{K}\sum_{k=1}^{K}\left\langle\mu_{k}\,,G(u_{k},w,\hat{z}_{k})\right\rangle_{\mathcal{U},\mathcal{U}^{\prime}}.\end{array}$
(21)

Up to a change of the functional framework, for the theoretical questions we
can simply consider the functional introduced in (18), by considering the
vector version of the different variables, namely $u=(u_{k})_{1\leq k\leq K}$,
$\hat{u}=(\hat{u}_{k})_{1\leq k\leq K}$ and $\mu=(\mu_{k})_{1\leq k\leq K}$.
Before deriving conditions of stationarity for $\mathbf{L}$, let us discuss
about existence of solutions for Problem ($\tilde{\ast\ast}$ ‣ 4.1).

###  4.2 Existence of minimizers for the relaxed problem

Existence of solutions for Problem ($\tilde{\ast\ast}$ ‣ 4.1) follows from
standard arguments:

######  Theorem 1.

Assume that $\nu>0$. Then Problem ($\tilde{\ast\ast}$ ‣ 4.1) admits a global
minimizer.

###### Proof.

The existence of feasible solutions for Problem ($\tilde{\ast\ast}$ ‣ 4.1) is
due to Proposition 2, combined with Proposition 3 (or also Proposition 5).
Since the functional

$\mathcal{U}\times\mathcal{W}\ni(u,w)\mapsto\frac{1}{2}\|u-\hat{u}\|^{2}_{\mathcal{U}}+\frac{\nu}{2}\|w\|_{\mathcal{W}}^{2}$
(22)

is bounded from below, Problem ($\tilde{\ast\ast}$ ‣ 4.1) admits a minimizing
sequence $(\mathrm{u}^{(n)},w^{(n)})_{n}$ of feasible solutions. Further, from
(22), the sequences $(u^{(n)})_{n}$ and $(w^{(n)})_{n}$ are bounded in
$\mathcal{U}$ and $\mathcal{W}$, respectively. Remind that the spaces
$\mathcal{U}$ and $\mathcal{W}$ are finite-dimensional. Consequently, up to
extraction, there exists $(u,w)\in\mathcal{U}\times\mathcal{W}$ such that
$(u^{(n)},w^{(n)})_{n}$ converges strongly towards $(u,w)$. It is clear that
$\mathcal{U}\times\mathcal{W}\ni(u,w)\mapsto
G(u,w,\hat{z})\in\mathcal{U}^{\prime}$ is continuous, so that we can pass to
the limit in $G(u^{(n)},w^{(n)},\hat{z})=0$ for obtaining $G(u,w,\hat{z})=0$,
which concludes the proof, as $(u,w)$ is then solution of Problem
($\tilde{\ast\ast}$ ‣ 4.1). ∎

######  Remark 2.

If we consider the case $\nu=0$, the boundedness of the sequence $(w_{n})_{n}$
in the proof above is a priori not guaranteed. For proceeding like we did, in
such a case we may need to assume further hypotheses on the neural network.
For example, we could exploit the equality $G(u^{(n)},w^{(n)})=0$, that is

$\mathbb{S}^{\prime}(u^{(n)})^{\ast}c_{y}(\mathbb{S}(u^{(n)}),\hat{z})+r^{\prime}_{u}(w^{(n)},u^{(n)})^{\ast}\gamma^{\prime}(r(w^{(n)},u^{(n)}))=0,$

due to the identity (6). The sequence $(u^{(n)})_{n}$ still converges, which
implies that the quantity

$r^{\prime}_{u}(w^{(n)},u^{(n)})^{\ast}\gamma^{\prime}(r(w^{(n)},u^{(n)}))$

is uniformly bounded with respect to $n$. In the simple one-dimensional case
where $n_{\mathrm{in}}=n_{\mathrm{out}}=1$, and $\gamma\equiv\mathrm{Id}$,
this would lead us to consider that the boundedness of the gradient of the
neural network implies the boundedness of its weights, which is not true in
general. Therefore considering $\nu>0$ seems to be reasonable, as this
provides a simple proof to the existence of minimizers.

###  4.3 Derivatives of second-order for the neural network

In section 4.5, we will need expressions of the second-order derivatives of
the neural network. From the formula (4) obtained in section 3.1, we derive

$r^{\prime\prime}_{uu}(w,u)=A_{L}\sum_{\ell=1}^{L-1}\left(\left(\mathrlap{\prod_{k=1}^{\ell-1}}{\longrightarrow}\rho^{\prime}(r^{(k)})A_{k}\right)\rho^{\prime\prime}(r^{(\ell)})\frac{\partial
r^{(\ell)}}{\partial u}A_{\ell}\left(\mathrlap{\prod_{k=\ell+1}^{L-1}}{\hskip
4.5pt\longrightarrow}\rho^{\prime}(r^{(k)})A_{k}\right)\right),$ (23)

and

$r^{\prime\prime}_{u{w_{s}}}(w,u).\tilde{w}=A_{L}\sum_{\ell=1}^{L-1}\left(\left(\mathrlap{\prod_{k=1}^{\ell-1}}{\longrightarrow}\rho^{\prime}(r^{(k)})A_{k}\right)\frac{\partial}{\partial
w_{s}}\left(\rho^{\prime}(r^{(\ell)}(w^{(\ell)},u))A_{\ell}\right).\tilde{w}\left(\mathrlap{\prod_{k=\ell+1}^{L-1}}{\hskip
4.5pt\longrightarrow}\rho^{\prime}(r^{(k)})A_{k}\right)\right)$ (24)

for all $1\leq s\leq L-1$, where we have simply denoted
$r^{(k)}=r^{(k)}(w^{(k)},u)$. We recall

$\displaystyle\frac{\partial r^{(\ell)}}{\partial
u}(w^{(\ell)},u)=A_{L}\mathrlap{\prod_{k=1}^{\ell-1}}{\longrightarrow}\rho^{\prime}(r^{(k)})A_{k},$
$\displaystyle\frac{\partial r^{(\ell)}}{\partial
w_{s}}(w^{(\ell)},u).\tilde{w}=A_{L}\left(\mathrlap{\prod_{k=s}^{\ell-1}}{\longrightarrow}\rho^{\prime}(r^{(k)})A_{k}\right)\left(\tilde{w}(\rho(r^{(s-1)}))\right),$

and for all $\tilde{w}:\mathrm{z}\mapsto\tilde{A}\mathrm{z}+\tilde{b}$

$\frac{\partial}{\partial
w_{s}}\left(\rho^{\prime}(r^{(\ell)}(w^{(\ell)},u))A_{\ell}\right).\tilde{w}=\displaystyle\left\\{\begin{array}[]{cl}\displaystyle
0&\text{if }s>\ell,\\\\[10.0pt]
\displaystyle\rho^{\prime}(r^{(\ell)}(w^{(\ell)},u))\tilde{A}+\rho^{\prime\prime}(r^{(\ell)})\left(\tilde{w}\rho(r^{(\ell-1)})\right)A_{\ell}&\text{if
}s=\ell,\\\\[10.0pt]
\displaystyle\rho^{\prime\prime}(r^{(\ell)})\left(\frac{\partial
r^{(\ell)}}{\partial w_{s}}.\tilde{w}\right)A_{\ell}&\text{if
}s<\ell.\end{array}\right.$

Like for (5), the products of type
$(\rho^{\prime\prime}\tilde{r})_{i}=\rho^{\prime\prime}_{i}\tilde{r}_{i}$ are
defined componentwise. For $s=L$, we have for all
$\tilde{w}:\mathrm{z}\mapsto\tilde{A}\mathrm{z}+\tilde{b}$

$\displaystyle
r^{\prime\prime}_{u{w_{L}}}(w,u).\tilde{w}=\tilde{A}\mathrlap{\prod_{k=1}^{\ell-1}}{\longrightarrow}\rho^{\prime}(r^{(k)})A_{k}.$

These expressions in themselves are not essential for what follows in the
present section, but they can help the reader to reproduce the numerical
results we obtained in section 6.

###  4.4 Properties of the adjoint operator

The adjoint equation for Problem ($\tilde{\ast\ast}$ ‣ 4.1) is obtained by
differentiating the Lagrangian $\mathbf{L}$ with respect to the variable $u$,
which gives

$G^{\prime}_{u}(u,w,\hat{z})^{\ast}.\mu=-(u-\hat{u}).$ (25)

Let us take a look at the operator
$G^{\prime}_{u}(u,w,\hat{z})^{\ast}\in\mathcal{L}(\mathcal{U},\mathcal{U}^{\prime})$
which arises in this equation. Since
$\mathcal{U}\simeq\mathbb{R}^{n_{\mathrm{in}}}$ is assumed to be of finite
dimension, for each $(u,w)\in\mathcal{U}\times\mathcal{W}$ the operator
$G^{\prime}_{u}(u,w,\hat{z})$ can be represented by a matrix of
$\mathbb{R}^{n_{\mathrm{in}}\times n_{\mathrm{in}}}$. Further, reminding that
$G(u,w,\hat{z})=\frac{\partial}{\partial u}J(\mathbb{S}(u),u,w)$, it yields
that $G^{\prime}_{u}(u,w,\hat{z})$ is an Hessian matrix, and so it is
represented by a symmetric matrix. Let us give the expression of
$G^{\prime}_{u}(u,w,\hat{z})$, by denoting $y=\mathbb{S}(u)$ :

$G^{\prime}_{u}(u,w,\hat{z})=\mathbb{S}^{\prime}(u)^{\ast}c^{\prime\prime}_{yy}(y,\hat{z})\mathbb{S}^{\prime}(u)+c^{\prime}_{y}(y,\hat{z})\mathbb{S}^{\prime\prime}(u)+\gamma^{\prime}(r(w,u))r_{uu}^{\prime\prime}(w,u)+r^{\prime}(w,u)^{\ast}\gamma^{\prime\prime}(r(w,u))r^{\prime}(w,u).$

More specifically, for all $v_{1}$, $v_{2}\in\mathcal{U}$ we have:

$\displaystyle\left\langle
G^{\prime}_{u}(u,w,\hat{z}).v_{1};v_{2}\right\rangle_{\mathcal{U}^{\prime},\mathcal{U}}$
$\displaystyle=$ $\displaystyle
c^{\prime\prime}_{yy}(y,\hat{z}).(\mathbb{S}^{\prime}(u).v_{1},\mathbb{S}^{\prime}(u).v_{2})+c^{\prime}_{y}(y,\hat{z}).(\mathbb{S}^{\prime\prime}(u).(v_{1},v_{2}))$
$\displaystyle+\gamma^{\prime\prime}(r(w,u)).(r^{\prime}(w,u).v_{1},r^{\prime}(w,u).v_{2})+\gamma^{\prime}(r(w,u)).(r_{uu}^{\prime\prime}(w,u).(v_{1},v_{2})).$

Solving equation (25) determines the value of the variable $\mu$, and thus the
question arises whether the operator
$G^{\prime}_{u}(u,w,\hat{z})^{\ast}=G^{\prime}_{u}(u,w,\hat{z})$ is
invertible. Let us comment on the general case where
$G^{\prime}_{u}(u,w,\hat{z})$ is not necessarily invertible.

####  4.4.1 When the adjoint operator is invertible

Consider any $(u,w)\in\mathcal{U}\times\mathcal{W}$ such that
$G^{\prime}_{u}(u,w,\hat{z})$ is an invertible matrix. Analogously to the
inner-problems, such a case corresponds to the assumption $(\mathbf{A2})$, and
following section 2.4, we can prove via the implicit function theorem the
existence of a mapping $w\mapsto u$ characterizing the equality
$G(u,w,\hat{z})=0$. This is a mapping of type control-to-state, since for the
outer-problem the variable $u$ plays the role of a state and $w$ is the
command. This case corresponds to variables $w\in\mathcal{W}$ such that the
associated inner-problem admits a unique solution $u\in\mathcal{W}$. Since
$G^{\prime}_{u}(u,w,\hat{z})$ is a square matrix, the matrix
$G^{\prime}_{u}(u,w,\hat{z})^{\ast}$ is also invertible, and so the adjoint
equation (25), namely $G^{\prime}_{u}(u,w,\hat{z})^{\ast}\mu+(u-\hat{u})=0$ in
the generic case, admits the unique solution

$\mu=-G^{\prime}_{u}(u,w,\hat{z})^{-T}(u-\hat{u}).$

Therefore in this case we are able to derive the necessary Karush-Kuhn-Tucker
conditions for Problem ($\tilde{\ast\ast}$ ‣ 4.1), as the associated
qualification constraints are satisfied without ambiguity (see Theorem 2).

####  4.4.2 When the adjoint operator is not invertible

In the case where $G^{\prime}_{u}(u,w,\hat{z})$ is not invertible, we
introduce a perturbation of the latter, based on the following basic result:

######  Lemma 1.

Let $A$ be a matrix. There exists $\varepsilon>0$ such that
$A-\varepsilon\mathrm{I}$ is invertible.

###### Proof.

Denote by $\chi(\lambda):=\mathrm{det}(A-\lambda\mathrm{I})$ the
characteristic polynomial of $A$. If $\mathrm{Sp(A)}=\\{0\\}$, any
$\varepsilon\neq 0$ is suitable. Otherwise, choose any $\varepsilon$ such that
$0<\varepsilon<\inf\left\\{|\lambda|,\
\lambda\in\mathrm{Sp}(A)\setminus\\{0\\}\right\\}$. Then we have
$\mathrm{det}(A-\varepsilon\mathrm{I})=\chi(\varepsilon)\neq 0$, which
concludes the proof. ∎

Therefore, when $G^{\prime}_{u}(u,w,\hat{z})$ is singular, we can find
$\varepsilon(u,w,\hat{z})>0$ small enough such that
$G^{\prime}_{u}(u,w,\hat{z})-\varepsilon(u,w,\hat{z})\mathrm{I}$ is
invertible. Note that, up to adapting the proof above, we can also choose
$\varepsilon(u,w,\hat{z})>0$ small enough such that
$G^{\prime}_{u}(u,w,\hat{z})+\varepsilon(u,w,\hat{z})\mathrm{I}$ is
invertible. Thus we replace $G^{\prime}_{u}(u,w,\hat{z})$ by
$G^{\prime}_{u}(u,w,\hat{z})\pm\varepsilon(u,w,\hat{z})\mathrm{I}$. This is
equivalent to add the small $L^{2}$-type regularizer
$u\mapsto\mp\frac{\varepsilon}{2}\|u\|_{\mathcal{U}}^{2}$ to the functional of
the inner-problem ($\ast$ ‣ 1.1). Note that the data $\hat{z}_{k}$ are of
finite number, and thus the corresponding solutions $u_{k}$ of the inner-
problems are of finite number too. Therefore $\varepsilon>0$ can be chosen
independent of $u$ and $\hat{z}$. In practice, we can choose
$\varepsilon=\varepsilon(w)>0$ as small as desired, as long as
$G^{\prime}_{u}(u,w,\hat{z})\pm\varepsilon(w)\mathrm{I}$ is invertible.

###  4.5 Necessary optimality conditions for the relaxed problem

The following result gives necessary conditions that an optimal solution of
Problem ($\tilde{\ast\ast}$ ‣ 4.1) has to fulfill.

######  Theorem 2.

Assume that $w\in\mathcal{W}$ is a solution of Problem ($\tilde{\ast\ast}$ ‣
4.1). Then for all $1\leq k\leq K$ we have $G(u_{k},w,\hat{z}_{k})=0$, and
necessarily

$\mathbf{G}_{\ell}(w):=\nu
w_{\ell}+\frac{1}{K}\sum_{k=1}^{K}G^{\prime}_{w_{\ell}}(u_{k},w,\hat{z}_{k})^{\ast}.\mu_{k}=0$
(26)

for all $1\leq\ell\leq L$, where $\mu_{k}$ is the solution of

$G^{\prime}_{u}(u_{k},w,\hat{z}_{k})^{\ast}.\mu_{k}+(u_{k}-\hat{u}_{k})=0$
(27)

for each $1\leq k\leq K$.

###### Proof.

Since the linear mappings
$G^{\prime}_{u}(u_{k},w,\hat{z}_{k})\in\mathcal{L}(\mathcal{U},\mathcal{U}^{\prime})$
are invertible for each $1\leq k\leq K$, or can be perturbed such that they
become invertible (see section 4.4.2), the qualification constraints of the
Karush-Kuhn-Tucker conditions are satisfied for the Lagrangian (21). Therefore
a solution of Problem ($\tilde{\ast\ast}$ ‣ 4.1) is necessary a critical point
of $\mathbf{L}^{(K)}$, which leads to the announced identities. ∎

In practice, we solve Problem ($\tilde{\ast\ast}$ ‣ 4.1) iteratively, by
updating the weights $w$ such that the norm of (26) converges to zero. Given a
set of weights $w$, we first compute each $u_{k}$ solution of ($\ast$ ‣ 1.1)
corresponding to $(\hat{z},r)=(\hat{z}_{k},r(w,\cdot))$. Note that the
variable $u=(u_{k})_{1\leq k\leq K}$ is a state variable in
($\tilde{\ast\ast}$ ‣ 4.1). Next we compute $\mu=(\mu_{k})_{1\leq k\leq K}$,
where each $\mu_{k}$ is solution of (27). Finally we are able to determine the
values of the left-hand-sides $\mathbf{G}_{\ell}(w)$ in (26), that we use in a
gradient rule.

##  5 A gradient-based algorithm

This section is devoted to algorithmic implementation of the relaxed approach
for Problem ($\ast\ast$ ‣ 1.1). Let us first explain how the inner-problems
($\ast$ ‣ 1.1) can be solved in practice.

###  5.1 Nesterov gradient descent for solving the inner-problems

We solve the inner-problems with Nesterov algorithm [Nes83], also called the
accelerated gradient method. It offers a nice compromise between robustness
and rapidity of convergence, as inner-problems need to be solved many times,
when solving the outer-problem. The corresponding method is given in Algorithm
‣ 5.1 below.

Initialization:

-

Initialize $u^{(0)}$ as $\displaystyle\frac{1}{K}\sum_{k=1}^{K}\hat{u}_{k}$.

Initial gradient:

From $u^{(0)}$, compute the (initial) gradient as follows:
$\left.\begin{minipage}{327.20668pt} \begin{description}
\par\description@item@[-] Compute the state $y$ corresponding to $u^{(0)}$, by
solving $\varphi(y,u^{(0)})=0$. \par\description@item@[-] Compute the costate
$p$ as solution of
$\varphi^{\prime}_{y}(y,u)^{\ast}.p=-c^{\prime}_{y}(y,\hat{z})$.
\par\description@item@[-] Compute the gradient $G^{(0)}=G(u,w,\hat{z})$.
\par\end{description} \end{minipage}\right\\}$ Compute the gradient in 3 steps
Store $u^{(0)}$ and $G^{(0)}$.

Armijo rule:

Choose $\alpha=0.5$.

-

Find the smallest $n\in\mathbb{N}$ such that
$J(u^{(0)}-\alpha^{n}G^{(0)})<J(u^{(0)})$.

-

Define $u^{(1)}=u^{(0)}-\alpha^{n}G^{(0)}$.

-

Compute the gradient $G^{(1)}$ as above, corresponding to $u^{(1)}$.

-

Store $u^{(1)}$ and $G^{(1)}$.

Nesterov gradient steps:

Initialization with $(u^{(0)},G^{(0)})$ and $(w^{(1)},G^{(1)})$.
Compute iteratively $u^{(n)}$ ($n\geq 2$) with the Nesterov steps.
While $||G(u^{(n)},w,\hat{z})||_{\mathcal{U}}>1.e^{-10}$, do gradient steps.

End:

Obtain $u$, approximated solution of ($\ast$ ‣ 1.1) with
$(\hat{z},r)=(\hat{z},r(w,\cdot))$.

Algorithm 0 Solving the inner-problem for a given regularizer $u\mapsto
r(u,w)$ via first-order necessary optimality conditions.

Note that for evaluating (26), Algorithm  ‣ 5.1 needs to be performed $K$
times, since we need to determine each $u_{k}$ corresponding to $\hat{z}_{k}$.

###  5.2 Barzilai-Borwein gradient rule for the relaxed outer-problem

We solve the outer-problem with the Barzilai-Borwein algorithm [Ray97]. The
corresponding method is explained in Algorithm 1, choosing $\nu=0$.

Initialization:

-

Load/create the data set.

-

Choose random weights as initial weights $w^{(0)}$ for the neural network.

Initial gradient:

From $w^{(0)}$, compute the (initial) gradient. For each $1\leq k\leq K$:
$\left.\begin{minipage}{327.20668pt} \begin{description}
\par\description@item@[-] With Alogirthm~{}\ref{algo-inner}, find the
approximated solution $u_{k}$ of the inner-problem~{}\eqref{pbstandard}
corresponding to $(\hat{z},r)=(\hat{z}_{k},r(\cdot,w^{(0)}))$.
\par\description@item@[-] Compute the costate $\mu_{k}$ as solution of
$G^{\prime}_{u}(u_{k},w^{(0)},\hat{z}_{k})^{\ast}.\mu_{k}=-(u_{k}-\hat{u}_{k})$.
\par\description@item@[-] Compute the gradient
$\mathbf{G}_{\ell}^{(0)}=G^{\prime}_{w_{\ell}}(u_{k},w^{(0)},\hat{z}_{k})^{\ast}.\mu_{k}$.
\par\end{description} \end{minipage}\right\\}$ Compute the gradient in 3 steps
Store $w^{(0)}$ and
$\mathbf{G}^{(0)}=(\mathbf{G}_{1}^{(0)},\dots,\mathbf{G}_{L}^{(0)})$.

Armijo rule:

Choose $\beta=0.5$.

-

Find the smallest $n\in\mathbb{N}$ such that
$\mathcal{J}(w^{(0)}-\beta^{n}\mathbf{G}^{(0)})<\mathcal{J}(w^{(0)})$.

-

Define $w^{(1)}=w^{(0)}-\beta^{n}\mathbf{G}^{(0)}$.

-

Compute the gradient $\mathbf{G}^{(1)}$ as above, corresponding to $w^{(1)}$.

-

Store $w^{(1)}$ and $\mathbf{G}^{(1)}$.

Barzilai-Borwein steps:

Initialization with $(w^{(0)},\mathbf{G}^{(0)})$ and
$(w^{(1)},\mathbf{G}^{(1)})$.
Compute iteratively $w^{(n)}$ ($n\geq 2$) with the Barzilai-Borwein steps.
While $\displaystyle||\mathbf{G}(w^{(n)})||_{\mathcal{W}}>1.e^{-8}$, do
gradient steps.

End:

Obtain $w\in\mathcal{W}$, approximated solution of ($\tilde{\ast\ast}$ ‣ 4.1).

Algorithm 1 Solving the first-order optimality conditions for the outer-
problem.

The Barzilai-Borwein is not a monotonous method, but presents second-order
convergence behavior, as the expressions of its gradient steps can also be
obtained via a quasi-Newton method, more precisely a BFGS method (see [NW06,
Chapter 6] for more details).

##  6 Illustration: Identification of coefficients in an elliptic equation

In this section we propose to address the case of coefficient identification
in an elliptic partial differential equations.

###  6.1 Presentation of the problem

We propose to illustrate the relaxation approach by considering the elliptic
equation identification problem in a smooth domain $\Omega$ of
$\mathbb{R}^{d}$ ($d\geq 1$). Denote by $\mathrm{L}^{\infty}_{F}(\Omega)$ any
finite-dimensional subspace of $\mathrm{L}^{\infty}(\Omega)$, with
$\dim(\mathrm{L}^{\infty}_{F}(\Omega))=n_{\mathrm{in}}$. For instance,
$\mathrm{L}^{\infty}_{F}(\Omega)$ can be made of piecewise constant functions
on a given grid of $\Omega$. Define for $M>m>0$ the following set:

$\mathcal{U}_{m,M}=\left\\{u\in\mathrm{L}_{F}^{\infty}(\Omega)\mid m\leq
u(x)\leq M,\text{ for a.e. }x\in\Omega\right\\}.$

The control-to-state mapping is given by
$\mathbb{S}:\mathcal{U}:=\mathcal{U}_{m,M}\ni u\mapsto
y\in\mathcal{Y}:=\mathrm{H}^{1}(\Omega)$, where $y$ satisfies the following
system in $\mathcal{Q}^{\prime}=\mathrm{H}^{-1}(\Omega)$, with
$\mathcal{Q}:=\mathrm{H}^{1}_{0}(\Omega)$:

$\left\\{\begin{array}[]{rl}-\operatorname{div}(u\nabla y)=f,&\text{in
}\Omega,\\\\[5.0pt] y_{|\partial\Omega}=g,&\text{on
}\partial\Omega.\end{array}\right.$ (28)

In system (28) the right-hand-sides $f$ and $g$ are given. Standard results
related to system (28) are presented in Appendix A.2.

#### On the inner-problem

Given $w\in\mathcal{W}$, as inner-problem we consider the following one:

$\left\\{\begin{array}[]{l}\displaystyle\min_{(y,u)\in\mathcal{Y}\times\mathcal{U}}\quad\frac{1}{2}\|y-\hat{z}\|_{\mathrm{L}^{2}(\Omega)}^{2}+r(w,u),\\\\[15.0pt]
\text{subject to~{}\eqref{sys-a}.}\end{array}\right.$ ($\ast$)

Note that we have chosen $\gamma\equiv\mathrm{Id}$, and thus necessarily
$n_{\mathrm{out}}=1$. The gradient for problem ($\ast$ ‣ 6.1) is given by

$G(u,w,\hat{z}).\tilde{u}=r^{\prime}_{u}(w,u).\tilde{u}+\int_{\Omega}\tilde{u}\nabla
y\cdot\nabla p\,\mathrm{d}\Omega.$

Following Proposition 5, the adjoint state
$p\in\mathcal{Q}=\mathrm{H}^{1}_{0}(\Omega)$ is given as the solution of
system (8), which corresponds to the following one:

$\left\\{\begin{array}[]{rl}-\operatorname{div}(u\nabla p)=-(y-y_{d})&\text{in
}\Omega,\\\\[5.0pt] p_{|\partial\Omega}=0&\text{on
}\partial\Omega.\end{array}\right.$ (29)

Proving existence and uniqueness of solutions for system (29) is
straightforward.

#### On the relaxed outer-problem

The relaxed outer-problem is formulated as follows:

$\left\\{\begin{array}[]{l}\displaystyle\min_{w\in\mathcal{W}}\quad\frac{1}{2K}\sum_{k=1}^{K}\|u_{k}-\hat{u}_{k}\|_{\mathrm{L}^{\infty}_{F}(\Omega)}^{2}+\frac{\nu}{2}\|w\|^{2}_{\mathcal{W}},\\\\[15.0pt]
\text{subject to: $G(u_{k},w,\hat{z}_{k})=0$ for all $1\leq k\leq
K$.}\end{array}\right.$ ($\tilde{\ast\ast}$)

Instead of detailing the different steps of the theoretical findings for
solving ($\tilde{\ast\ast}$ ‣ 6.1) in this particular example, we leave this
task to the reader, and prefer to focus now on the numerical results.

In what follows, we consider $\Omega=(0,1)$ (and thus $d=1$), and discretize
the different elliptic equations with P1-finite elements on a uniform mesh
made of 100 subdivisions.

###  6.2 Rediscovering the $L^{2}$-norm, or not

Data $\\{(u_{k},y_{k});1\leq k\leq K\\}$ are generated artificially by solving
the inner-problems corresponding to the $L^{2}$ regularizer
$u\mapsto\frac{3}{2}\|u\|_{\mathcal{U}}^{2}$, with $\mathcal{U}=\mathbb{R}$.
Next a feed-forward neural network is trained with this data set. The purpose
of this test is to compare the trained neural network with the $L^{2}$
function that represents the ground-truth regularizer to (re-)discover.
Convergence of the objective function is presented in Figure 1, with values
given in Table 1. Figure 2 represents the graph of the trained ANN compared
with the one of the regularizer with which we generated the artificial data,
namely $u\mapsto\frac{3}{2}|u|_{\mathbb{R}}^{2}$.

linear scale

log scale

Figure 1: Misfit evolution during gradient steps for training a neural network
from data generated with a $\mathrm{L}^{2}$-regularizer, in linear scale
(left) and log scale (right).

Figure 2: Graph of the trained neural network in blue, and of the ground-truth
$L^{2}$ regularizer in red. (For interpretation of the references to color in
this figure legend, the reader is referred to the digital version of this
article.)

In Figure 1 and Table 1 we see that the relative misfit decreased from $54$%
to $0.27$% after a few iterations. The curves of Figure 2 suggest non-
uniqueness for Problem ($\tilde{\ast\ast}$ ‣ 4.1), as the residual misfit is
of the same range as the approximation error made when discretizing (28).

$\displaystyle\begin{array}[]{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}\hline\cr\text{steps}&1&2&3&4&5&6&7&8&9&10&11&12&13&14&15\\\
\hline\cr\text{misfit
(\%)}&54.6&31.1&18.4&5.62&2.00&1.25&1.14&7.09&10.2&2.83&2.00&1.17&1.16&1.15&1.14\\\
\hline\cr\end{array}$
$\displaystyle\begin{array}[]{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}\hline\cr\text{steps}&16&17&18&19&20&21&22&23&24&25&26&27&28&29\\\
\hline\cr\text{misfit
(\%)}&1.46&1.19&1.18&1.19&2.68&1.09&0.33&0.30&0.30&0.29&0.30&0.32&0.27&0.27\\\
\hline\cr\end{array}$
Table 1: Misfit vales (in %) through the different gradient steps.

###  6.3 Regularizers for compensating noise on the data

Given a set of $K=10$ piecewise constant coefficients
$u_{\mathrm{true}}=(u_{\mathrm{true},k})_{1\leq k\leq K}$, we compute the
corresponding states $(y_{\mathrm{true},k})_{1\leq k\leq K}$ solutions of
(28). We then introduce artificially a noise following a Gaussian distribution
that transforms each $y_{\mathrm{true},k}$ into $y_{\mathrm{noisy},k}$. The
pairs
$\\{(\hat{z}_{k},\hat{u}_{k}):=(y_{\mathrm{noisy},k},u_{\mathrm{true},k}),1\leq
k\leq K\\}$ constitute our data set on which we will train the neural network
as regularizer for the inner-problem. The aim of such a regularizer consists
in compensating the effects of the noise. Problem ($\ast$ ‣ 1.1) equipped with
the so-trained regularizer will then enable us to determine the desired
coefficients $u$, in spite of the noisy state $\hat{z}$ on which we can rely
only.

linear scale

log scale

Figure 3: Misfit evolution during gradient steps for training a noise-
compensating neural network, in linear scale (left) and log scale (right).

$\displaystyle\begin{array}[]{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}\hline\cr\text{steps}&1&2&3&4&5&6&7&8&9&10&11&12&13&14&15&16&17&18\\\
\hline\cr\text{misfit
(\%)}&54.2&49.4&31.1&50.0&52.8&5.37&13.9&4.25&4.24&3.58&3.45&4.63&3.74&3.33&3.67&3.60&3.65&3.78\\\
\hline\cr\end{array}$
$\displaystyle\begin{array}[]{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}\hline\cr\text{steps}&19&20&21&22&23&24&25&26&27&28&29&30&31&32&33&34&35&36\\\
\hline\cr\text{misfit
(\%)}&3.79&3.56&6.80&3.70&3.90&3.50&3.79&4.17&3.79&3.50&3.79&3.91&3.91&3.79&3.79&3.91&3.91&3.79\\\
\hline\cr\end{array}$
Table 2: Misfit vales (in %) through the different gradient steps.

In the illustration presented in Figure 3 and Table 2, the initial misfit –
corresponding to a neural network close to zero – is about 54%. We see that
few iterations are required for obtaining a misfit lower than 5%. The type of
neural network we consider does not allow us to reduce the misfit under 3,32%.

##  7 Conclusion

We have proposed a method for addressing in practice the data-driven design of
regularizers for improving a class of inverse problems formulated as optimal
control problems. The class of regularizers is restricted to a set of
artificial neural networks, which goes in the sense of the purely nonlinear
nature of such a problem. Our approach is bi-level, and based on a relaxation
exploiting the necessary optimality conditions for these problems. We saw that
such a type of approach enables practical and efficient implementation,
providing that the different objective functions and activation functions of
the neural network are smooth. Developing a similar approach when quantities
involved are less smooth could present an interesting challenge. Further,
investigating applicability of our approach for more complex problems, like
image segmentation for instance, would demonstrate even more its relevance. In
such a context, one may need to take into account neural networks
architectures that allow more sophisticated realizations. Finally, noticing
that our approach necessitates the evaluation of second-order derivatives of
the control-to-state mapping for the inner-problem, developing higher order
adjoint-based methods for deriving the optimality conditions for the outer-
problem could potentially improve the implementation and the applicability to
large-scale problems.

##  Appendix A Appendix

###  A.1 On the optimal control approach for solving an inverse problem

Let us fix $w$ for this subsection. A classical approach for solving inverse
problems consists in taking into account the equality constraint
$\varphi(y,u)=0$ by incorporating it in the functional to minimize, as
follows:

$\displaystyle\min_{u\in\mathcal{U},y\in\mathcal{Y}}\left(J(y,u,w):=\alpha\|\varphi(y,u)\|_{\mathcal{Q}}^{2}+c(y,\hat{z})+\gamma\circ
r(w,u)\right),$

where $\alpha>0$ is a given parameter and $w$ is also given. For example,
consider the problem of identifying a conductivity represented by the variable
$u$ in the following elliptic partial differential equation:

$\displaystyle\left\\{\begin{array}[]{rl}-\operatorname{div}(u\nabla
y)=f,&\text{in }\Omega,\\\ y_{|\partial\Omega}=g,&\text{on
}\partial\Omega.\end{array}\right.$

Here we set $\mathcal{Y}=\mathrm{H}^{1}(\Omega)$,
$\mathcal{Q}=\mathrm{H}^{1}_{0}(\Omega)$ and
$\mathcal{Q}^{\prime}=\mathrm{H}^{-1}(\Omega)$. Right-hand-sides are given by
$f\in\mathrm{H}^{-1}(\Omega)$ and $g\in\mathrm{H}^{1/2}(\partial\Omega)$. In
that case, the contribution of the state equation in the functional to
minimize is $\|-\operatorname{div}(u\nabla y)-f\|_{\mathcal{Q}}^{2}$. This
approach would imply differentiation of the data, and thus could lead to
approximation biases, especially when the set of data is incomplete. Instead
of that, we consider the state equation as a constraint, that we impose by
introducing a Lagrange multiplier denoted by
$p\in\mathcal{Q}\simeq\mathcal{Q}^{\prime\prime}$, namely the costate or
adjoint variable, leading to the search of a saddle-point of a Lagrangian
functional, given as:

$\mathscr{L}(y,u,p)=c(y,\hat{z})+\gamma\circ r(w,u)+\langle
p\,,\varphi(y,u)\rangle_{\mathcal{Q},\mathcal{Q}^{\prime}}.$

A saddle-point of $\mathscr{L}$ is obtained as a critical point with respect
to the variables $(y,u)$ and $p$. The duality product $\langle
p\,,\varphi(y,u)\rangle_{\mathcal{Q},\mathcal{Q}^{\prime}}$ corresponds to a
weak variational formulation, namely

$\displaystyle\langle
p\,,\varphi(y,u)\rangle_{\mathcal{Q},\mathcal{Q}^{\prime}}$ $\displaystyle=$
$\displaystyle\int_{\Omega}p\cdot\operatorname{div}(u\nabla
y)\,\mathrm{d}\Omega-\int_{\Omega}p\cdot f\,\mathrm{d}\Omega$ $\displaystyle=$
$\displaystyle-\int_{\Omega}\nabla p:(u\nabla
y)\,\mathrm{d}\Omega+\int_{\partial\Omega}p\cdot\left(u\frac{\partial
y}{\partial n}\right)\,\mathrm{d}\partial\Omega-\int_{\Omega}\nabla p:(u\nabla
y)\,\mathrm{d}\Omega,$

and thus we see that the data represented by the control function $u$ does not
need to be differentiated.

###  A.2 Results related to the inverse conductivity problem

Recall that
$\mathrm{H}^{-1}(\Omega)=\left(\mathrm{H}^{1}_{0}(\Omega)\right)^{\prime}$. It
is well-known that for all $u\in\mathcal{U}_{m,M}$,
$f\in\mathrm{H}^{-1}(\Omega)$ and $g\in\mathrm{H}^{1/2}(\partial\Omega)$
system (28) admits a unique solution $y\in\mathrm{H}^{1}(\Omega)$. Moreover,
there exists a constant $C>0$ depending only on $\Omega$ such that the
following estimate holds:

$\|y\|_{\mathrm{H}^{1}(\Omega)}\leq\frac{1}{m}\left(\|f\|_{\mathrm{H}^{-1}(\Omega)}+M\|g\|_{\mathrm{H}^{1/2}(\partial\Omega)}\right).$
(35)

Furthermore, the following result show that the control-to-state mapping
$\mathbb{S}$ is Lipschitz:

######  Lemma 2.

Let be $M>m>0$ and assume that $u_{1},u_{2}\in\mathcal{U}_{m,M}$. Then the
solutions $y_{1}$ and $y_{2}$ of system (28) corresponding to $u=u_{1}$ and
$u=u_{2}$ respectively satisfy

$\|y_{1}-y_{2}\|_{\mathrm{H}_{0}^{1}(\Omega)}\leq\frac{1}{m^{2}}\|u_{1}-u_{2}\|_{\mathrm{L}^{\infty}(\Omega)}\left(\|f\|_{\mathrm{H}^{-1}(\Omega)}+M\|g\|_{\mathrm{H}^{1/2}(\partial\Omega)}\right),$
(36)

where
$\|\varphi\|_{\mathrm{H}_{0}^{1}(\Omega)}:=\|\nabla\varphi\|_{\mathrm{L}^{2}(\Omega)}$.

###### Proof.

Substituting the first equations of (28) satisfied by $y_{1}$ and $y_{2}$
leads us to

$-\operatorname{div}(u_{1}\nabla y_{1})+\operatorname{div}(u_{2}\nabla
y_{2})=0\Rightarrow-\operatorname{div}(u_{1}\nabla(y_{1}-y_{2}))+\operatorname{div}((u_{2}-u_{1})\nabla
y_{2})=0.$

Note that $y_{1}-y_{2}\in\mathrm{H}^{1}_{0}(\Omega)$. Then, taking the inner
product of the last identity above by $y_{1}-y_{2}$ and integrating by parts
yields

$\displaystyle\int_{\Omega}u_{1}|y_{1}-y_{2}|^{2}\,\mathrm{d}\Omega=\int_{\Omega}(u_{2}-u_{1})\nabla
y_{2}\cdot\nabla(y_{1}-y_{2})\,\mathrm{d}\Omega$ $\displaystyle\Rightarrow$
$\displaystyle
m\|\nabla(y_{1}-y_{2})\|^{2}_{\mathrm{L}^{2}(\Omega)}\leq\|u_{1}-u_{2}\|_{\mathrm{L}^{\infty}(\Omega)}\|\nabla
y_{2}\|_{\mathrm{L}^{2}(\Omega)}\|\nabla(y_{1}-y_{2})\|_{\mathrm{L}^{2}(\Omega)}$
$\displaystyle\Rightarrow$ $\displaystyle
m\|y_{1}-y_{2}\|_{\mathrm{H}^{1}_{0}(\Omega)}\leq\|u_{1}-u_{2}\|_{\mathrm{L}^{\infty}(\Omega)}\|y_{2}\|_{\mathrm{H}^{1}_{0}(\Omega)}.$

Using estimate (35) satisfied by $y_{2}$, the previous one implies

$m\|y_{1}-y_{2}\|_{\mathrm{H}^{1}_{0}(\Omega)}\leq\|u_{1}-u_{2}\|_{\mathrm{L}^{\infty}(\Omega)}\frac{1}{m}\left(\|f\|_{\mathrm{H}^{-1}(\Omega)}+M\|g\|_{\mathrm{H}^{1/2}(\partial\Omega)}\right).$

Thus estimate (36) follows. ∎

#### On the inner-problem

Let be $w\in\mathcal{W}$. Considering $y_{d}\in\mathrm{L}^{2}(\Omega)$, the
inner-problem is formulated as follows:

$\left\\{\begin{array}[]{l}\displaystyle\min_{(y,u)\in\mathrm{H}^{1}(\Omega)\times\mathcal{U}_{m,M}}J(y,u,w)=\frac{1}{2}\|y-y_{d}\|_{\mathrm{L}^{2}(\Omega)}^{2}+\gamma\circ
r(w,u)\\\\[10.0pt] \text{subject to \eqref{sys-a}.}\end{array}\right.$ (37)

The result below establishes existence of minimizers for Problem (37).

######  Lemma 3.

Given $w\in\mathcal{W}$, Problem (37) admits global minimizers.

###### Proof.

Since the function $J$ of Problem (37) is bounded from below, and since system
(28) is well-posed, $J$ admits a minimizing sequence $(y_{n},u_{n})_{n}$ of
feasible solutions. From the definition of $\mathcal{U}_{m,M}$, the sequence
$(u_{n})_{n}$ is uniformly bounded in $\mathrm{L}^{\infty}_{F}(\Omega)$, which
is of finite dimension. Therefore, up to extraction, the sequence
$(u_{n})_{n}$ converges strongly in $\mathrm{L}^{\infty}(\Omega)$ towards
$u\in\mathcal{U}_{m,M}$. From estimate (35), the sequence $(y_{n})_{n}$ is
bounded in $\mathrm{H}^{1}(\Omega)$. In virtue of the Banach-Alaoglu theorem,
up to extraction it converges weakly towards $y\in\mathrm{H}^{1}(\Omega)$. Let
us show that $y$ in solution of system (28) corresponding to $u$. For all
$\varphi\in\mathrm{H}^{1}_{0}(\Omega)$, using the Green formula, we have

$\displaystyle\langle-\operatorname{div}(u\nabla
y)-f;\varphi\rangle_{\mathrm{H}^{-1}(\Omega),\mathrm{H}^{1}_{0}(\Omega)}$
$\displaystyle=$ $\displaystyle\langle-\operatorname{div}(u\nabla
y)+\operatorname{div}(u_{n}\nabla
y_{n});\varphi\rangle_{\mathrm{H}^{-1}(\Omega),\mathrm{H}^{1}_{0}(\Omega)}$
$\displaystyle=$ $\displaystyle\langle u\nabla y-u_{n}\nabla
y_{n};\nabla\varphi\rangle_{[\mathrm{L}^{2}(\Omega)]^{d}}$ $\displaystyle=$
$\displaystyle\langle(u-u_{n})\nabla
y_{n};\nabla\varphi\rangle_{[\mathrm{L}^{2}(\Omega)]^{d}}+\langle
u\nabla(y-y_{n});\nabla\varphi\rangle_{[\mathrm{L}^{2}(\Omega)]^{d}}$
$\displaystyle=$ $\displaystyle\langle(u-u_{n})\nabla
y_{n};\nabla\varphi\rangle_{[\mathrm{L}^{2}(\Omega)]^{d}}-\langle\operatorname{div}(u\nabla\varphi);y-y_{n}\rangle_{\mathrm{H}^{-1}(\Omega),\mathrm{H}^{1}_{0}(\Omega)},$

where we have used that $y-y_{n}\in\mathrm{H}^{1}_{0}(\Omega)$. This yields

$\displaystyle\left|\langle-\operatorname{div}(u\nabla
y)-f;\varphi\rangle_{\mathrm{H}^{-1}(\Omega),\mathrm{H}^{1}_{0}(\Omega)}\right|$
$\displaystyle\leq$
$\displaystyle\|u-u_{n}\|_{\mathrm{L}^{\infty}(\Omega)}\|\nabla
y_{n}\|_{[\mathrm{L}^{2}(\Omega)]^{d}}\|\nabla\varphi\|_{[\mathrm{L}^{2}(\Omega)]^{d}}$
$\displaystyle+\left|\langle\operatorname{div}(u\nabla\varphi);y-y_{n}\rangle_{\mathrm{H}^{-1}(\Omega),\mathrm{H}^{1}_{0}(\Omega)}\right|.$

Passing to the limit, we deduce that $-\operatorname{div}(u\nabla y)=f$ in
$\mathrm{H}^{-1}(\Omega)$, and thus that $(y,u)$ is a solution of Problem
(37). ∎

Recall finally the following Lipschitz estimate for the Calderon problem (see
for instance [Uhl14, Theorem 2.6], refering to [Ale88]):

######  Lemma 4.

Assume that $\Omega$ is a bounded open set with smooth boundary, with $d\geq
3$. For $j=1,2$ and $s>d/2$, denote by $u_{j}\in\mathrm{H}^{s+2}(\Omega)\cap
A_{m,M}$ satisfying

$\|u_{j}\|_{\mathrm{H}^{s+2}(\Omega)}\leq M.$

There exists a constant $C=C(\Omega,d,M,m,s)>0$ such that

$\|u_{1}-u_{2}\|_{\mathrm{L}^{\infty}(\Omega)}\leq
C\omega\left(\|\Lambda_{u_{1}}-\Lambda_{u_{2}}\|_{\mathcal{L}(\mathrm{H}^{1/2}(\partial\Omega),\mathrm{H}^{-1/2}(\partial\Omega))}\right),$

where $\omega:x\mapsto|\log(x)|^{-\theta}+x$, for some $0<\theta<1$.

## References

  * [Ale88] Giovanni Alessandrini.  Stable determination of conductivity by boundary measurements.  Appl. Anal., 27(1-3):153–172, 1988.
  * [AO17] Jonas Adler and Ozan Öktem.  Solving ill-posed inverse problems using iterative deep neural networks.  Inverse Problems, 33(12):124007, 24, 2017.
  * [GPEB19] Philipp Grohs, Dmytro Perekrestenko, Dennis Elbrächter, and Helmut Bölcskei.  Deep neural network approximation theory.  ArXiv, abs/1901.02220, 2019.
  * [HKB18] Gernot Holler, Karl Kunisch, and Richard C. Barnard.  A bilevel approach for parameter learning in inverse problems.  Inverse Problems, 34(11):115012, 28, 2018.
  * [HN20] M. Haltmeier and Linh V. Nguyen.  Regularization of inverse problems by neural networks.  ArXiv, abs/2006.03972, 2020.
  * [HPUU09] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.  Optimization with PDE constraints, volume 23 of Mathematical Modelling: Theory and Applications.  Springer, New York, 2009.
  * [KP13] Karl Kunisch and Thomas Pock.  A bilevel optimization approach for parameter learning in variational models.  SIAM J. Imaging Sci., 6(2):938–983, 2013.
  * [LOS18] Sebastian Lunz, Ozan Öktem, and Carola-Bibiane Schönlieb.  Adversarial regularizers in inverse problems.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 8507–8516. Curran Associates, Inc., 2018.
  * [LSAH20a] Housen Li, Johannes Schwab, Stephan Antholzer, and Markus Haltmeier.  NETT: solving inverse problems with deep neural networks.  Inverse Problems, 36(6):065005, 23, 2020.
  * [LSAH20b] Housen Li, Johannes Schwab, Stephan Antholzer, and Markus Haltmeier.  NETT: solving inverse problems with deep neural networks.  Inverse Problems, 36(6):065005, jun 2020.
  * [MDS+20] Subhadip Mukherjee, Sören Dittmer, Zakhar Shumaylov, Sebastian Lunz, Ozan Öktem, and Carola-Bibiane Schönlieb.  Learned convex regularizers for inverse problems, 2020.
  * [Nes83] Yu. E. Nesterov.  A method for solving the convex programming problem with convergence rate $O(1/k^{2})$.  Dokl. Akad. Nauk SSSR, 269(3):543–547, 1983.
  * [NW06] Jorge Nocedal and Stephen J. Wright.  Numerical optimization.  Springer Series in Operations Research and Financial Engineering. Springer, New York, second edition, 2006.
  * [ONSH19] Daniel Obmann, Linh Nguyen, Johannes Schwab, and Markus Haltmeier.  Augmented nett regularization of inverse problems, 2019.
  * [ONSH20] D. Obmann, L. Nguyen, J. Schwab, and M. Haltmeier.  Sparse anett for solving inverse problems with deep learning.  In 2020 IEEE 17th International Symposium on Biomedical Imaging Workshops (ISBI Workshops), pages 1–4, 2020.
  * [OSH19] D. Obmann, J. Schwab, and M. Haltmeier.  Sparse synthesis regularization with deep neural networks.  In 2019 13th International conference on Sampling Theory and Applications (SampTA), pages 1–5, 2019.
  * [Ray97] M. Raydan.  The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem.  SIAM J. Optim., 7(1):26–33, 1997.
  * [Uhl14] Gunther Uhlmann.  30 years of Calderón’s problem.  In Séminaire Laurent Schwartz—Équations aux dérivées partielles et applications. Année 2012–2013, Sémin. Équ. Dériv. Partielles, pages Exp. No. XIII, 25. École Polytech., Palaiseau, 2014.