# To be or not to be stable, that is the question: understanding neural
networks for inverse problems

D. Evangelista, J. Nagy, E. Morotti, E. Loli Piccolomini

###### Abstract

The solution of linear inverse problems arising, for example, in signal and
image processing is a challenging problem, since the ill-conditioning
amplifies the noise on the data. Recently introduced deep-learning based
algorithms overwhelm the more traditional model-based approaches but they
typically suffer from instability with respect to data perturbation. In this
paper, we theoretically analyse the trade-off between neural networks
stability and accuracy in the solution of linear inverse problems. Moreover,
we propose different supervised and unsupervised solutions, to increase
network stability by maintaining good accuracy, by inheriting, in the network
training, regularization from a model-based iterative scheme. Extensive
numerical experiments on image deblurring confirm the theoretical results and
the effectiveness of the proposed networks in solving inverse problems with
stability with respect to noise.

Keywords— Neural Networks Stability, Inverse Problems, Stability, Image
Deblurring

##  1 Introduction

Linear inverse problems of the form:

$\displaystyle y=Ax^{gt}$ (1.1)

where $A\in\mathbb{R}^{m\times n}$ is the discretization of a linear operator,
$x^{gt}\in\mathbb{R}^{n}$ and $y\in\mathbb{R}^{m}$ arise in various
applications of signal or image processing, such as deblur, super-resolution
or tomographic reconstruction [13, 19]. It is well-know that such inverse
problems are ill-posed and their solution is very challenging when noise
affects the data. We consider here the case of data corrupted by Gaussian
noise, i.e.

$y^{\delta}=Ax^{gt}+e\qquad e\sim\mathcal{N}(0,\delta^{2}I)$ (1.2)

where $\delta$ denotes the standard deviation of the white additive Gaussian
noise.

Traditional approaches model the problem as a minimization of an objective
function containing a data-fit component and a regularization prior, with
possible further constraints on the solution [5, 8]. They compute accurate
solutions, but in general the computational time required is high and many
parameters must be set for any input data.

In the last few years neural networks have been introduced with great success
for the solution of problem (1.1) [2, 17]. Deep-learning based algorithms are
very fast, do not require parameter setting, but their accuracy is obtained at
the expense of stability. Hence, they are not efficient in dealing with noisy
data. In particular, they often produce poor results when applied to data
corrupted by noise different from the one learnt in the training phase. This
is usually referred to as network instability. Recent studies [7, 11]
demonstrate that it is not possible to achieve optimal stability and accuracy
at the same time. However, the stability of neural networks solving inverse
problems is not widely studied in literature up to now.

#### Contributions.

In this work, we theoretically analyse accuracy and stability of neural
newtorks solving linear inverse problems, by proving a mathematical relation
showing that it is not possible to increase stability without decreasing
accuracy. Then we propose three possible solutions to improve stability while
preserving good accuracy. First, starting from the theory developed before, we
introduce regularization in the neural network from the traditional
regularized variational approach, by training the networks on the solutions
computed with the regularized models. We call this approach as REgularized
Neural Network (ReNN). Second, we act on the network input to stabilize the
process by pre-processing the noisy data through a regularized model-based
approach. We call this approach as STabilized Neural Network (StNN). A
combination of ReNN and StNN (denoted as StReNN) is also proposed as a good
trade-off between stability and accuracy.

#### Notations.

In order to improve the readability of the work, we will introduce some
notations. In particular, we will always consider our data $x^{gt}$ to lie on
a subset $\chi$ of $\mathbb{R}^{n}$, the set of _admissible data_. We will
assume $\chi$ to be compact and locally connected, since we want that for any
admissible datum $x^{gt}\in\chi$, there exists a neighborhood $V$ of $x^{gt}$
such that any $x\in V$ is admissible. We will also consider
$\mathcal{Y}=Rg(A;\chi)$, the range of $A$ over $\chi$, where $A$ is a
continuous linear operator. By continuity of $A$, $\mathcal{Y}$ will also be
compact and locally connected. The compactness of $\mathcal{Y}$ implies that
$(C^{0}(\mathcal{Y},\mathcal{X}),||\cdot||_{\infty})$ is a Banach space. Thus,
when $f:\mathcal{Y}\to\mathcal{X}$ is continuous, we will always assume
$||f||$ to be the $||\cdot||_{\infty}$ norm. On the contrary, when
$x\in\mathcal{X}\subseteq\mathbb{R}^{n}$ or
$y\in\mathcal{Y}\subseteq\mathbb{R}^{m}$, then $||x||$ and $||y||$ will be
Euclidian norms. For any $\epsilon>0$, we also define
$\mathcal{Y}^{\epsilon}=\\{y+e;y\in\mathcal{Y},||e||\leq\epsilon\\}$. Since
$\mathcal{Y}^{\epsilon}$ is compact and locally connected if $\mathcal{Y}$ is,
then also $(C^{0}(\mathcal{Y}^{\epsilon},\mathcal{X}),||\cdot||_{\infty})$ is
a Banach space and the same discussion we did for $\mathcal{Y}$ applies.

#### Structure of the article.

The article is structured as follows. In Section 2 we introduce the concept of
reconstructor together with its accuracy and stability; in Section 3, after
introducing as reconstructors Tikhonov and neural networks, we propose a new
regularized reconstructor, named ReNN. A new technique to obtain two
stabilized reconstructors, StNN and StReNN, is proposed in Section 4. In
Section 5 and Section 6 we describe the experiments performed with the
proposed new reconstructors. Finally, in Section 7 we report some conclusions.

##  2 Reconstructors for the solution of Linear Inverse Problems

Given a problem like (1.2) where $m\geq n$, one aims at reconstructing a good
approximation of $x^{gt}$ as possible, using $y^{\delta}$ and the information
in $A$. To this extent, we formalize the concept of reconstructors.

######  Definition 2.1.

Any continuous function $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$, mapping $y$ to
$x=\Psi(y)$, is called a reconstructor.

We are mainly interested in reconstructors returning a good approximation of
the original signal $x^{gt}$ when applied to the noiseless datum $y$ of
Equation (1.1).

######  Definition 2.2.

A reconstructor $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$ is said to be
$\eta^{-1}$-accurate, with $\eta>0$, if

$\displaystyle||\Psi(Ax^{gt})-x^{gt}||\leq\eta\qquad\forall
x^{gt}\in\mathcal{X}.$

The smallest $\eta>0$ for which the above relationship is verified for any
$x^{gt}\in\mathcal{X}$ is called reconstruction error (or error parameter) of
$\Psi$.

Given the inverse problem statement where $A$ is a column full-rank matrix, it
is possible to build an accurate reconstructor, by exploiting the pseudo-
inverse matrix $A^{\dagger}=(A^{T}A)^{-1}A^{T}$. The resulting pseudo-inverse
reconstructor is defined by

$\Psi^{\dagger}(y)=A^{\dagger}y$ (2.1)

and the following proposition holds.

######  Proposition 2.1.

If $A$ is column full-rank in the inverse problem (1.2), the reconstructor
$\Psi^{\dagger}:\mathbb{R}^{m}\to\mathbb{R}^{n}$ is $\infty$-accurate (i.e. it
has reconstruction error $\eta=0$).

[1em]0em

###### Proof.

The proof is a simple computation.

$||\Psi^{\dagger}(Ax^{gt})-x^{gt}||=||(A^{T}A)^{-1}A^{T}(Ax^{gt})-x^{gt}||=||(A^{T}A)^{-1}(A^{T}A)x^{gt}-x^{gt}||=0$

which concludes the proof. $\blacksquare$

However, reconstructors are rarely applied to noise-free data, hence a study
on the robustness of reconstructors is necessary.

######  Definition 2.3.

Let $\epsilon>0$ and $\Psi$ be an $\eta^{-1}$-accurate reconstructor applied
to problem (1.2). We define the $\epsilon$-stability constant
$C_{\Psi}^{\epsilon}$ of $\Psi$ as

$C^{\epsilon}_{\Psi}=\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
e:||e||\leq\epsilon\end{subarray}}\frac{||\Psi(Ax^{gt}+e)-x^{gt}||-\eta}{||e||}$
(2.2)

######  Definition 2.4.

The reconstructor $\Psi$ is said to be $\epsilon$-stable for a given
$\epsilon$ if $C_{\Psi}^{\epsilon}\in[0,1)$.

An $\epsilon$-stable reconstructor $\Psi$ does not amplifies corruptions which
have norm less than $\epsilon$, since from (2.2) we derive

$||\Psi(Ax^{gt}+e)-x^{gt}||\leq\eta+C^{\epsilon}_{\Psi}||e||\qquad\forall
x^{gt}\in\mathcal{X},\forall||e||\leq\epsilon.$

######  Definition 2.5.

We define the stability radius $\bar{\epsilon}$ of $\Psi$ as

$\bar{\epsilon}=\sup\\{\epsilon>0;C^{\epsilon}_{\Psi}\in[0,1)\\}$ (2.3)

An example of reconstructor with infinite stability radius is the following.
Given $\epsilon>0$, if $\mu$ is a probability distribution over $\mathcal{X}$
(for example, $\mu$ is the normalized Lebesgue distribution over
$\mathcal{X}$), the reconstructor defined as

$\Psi^{\mathcal{X},\epsilon}(y^{\delta})=\int_{\mathcal{X}}x\mu(dx)\qquad\forall
y^{\delta}\in\mathcal{Y}^{\epsilon}$

is $\epsilon$-stable independently from the value of $\epsilon>0$. Indeed,

$||\Psi^{\mathcal{X},\epsilon}(Ax^{gt}+e)-x^{gt}||=||\int_{\mathcal{X}}x\mu(dx)-x^{gt}||\leq\rho(\mathcal{X})$

where $\rho(\mathcal{X})$ is the radius of $\mathcal{X}$, defined as
$\rho(\mathcal{X})=\inf\\{r>0:\mathcal{X}\subseteq\mathcal{B}(\int_{\mathcal{X}}x\mu(dx);r)\\}$.
As consequence, the stability constant is zero regardless $\epsilon$ and
$\Psi^{\mathcal{X}(y),\epsilon}$ has accuracy $\rho(\mathcal{X})^{-1}$.
On the contrary, the pseudo-inverse reconstructor $\Psi^{\dagger}(y)$ in
Equation (2.1) is unstable for any $\epsilon>0$ when $A$ is ill-conditioned.
Indeed,

$||\Psi^{\dagger}(Ax^{gt}+e)-x^{gt}||=||(A^{T}A)^{-1}(A^{T}A)x^{gt}+(A^{T}A)^{-1}A^{T}e-x^{gt}||=||(A^{T}A)^{-1}A^{T}e||.$

If $A=U\Sigma V^{T}$ is the SVD decomposition of $A$, then

$(A^{T}A)^{-1}A^{T}e=(V\Sigma^{2}V^{T})^{-1}V\Sigma
U^{T}e=V\Sigma^{-1}U^{T}e=\sum_{i=1}^{n}\frac{u_{i}^{T}e}{\sigma_{i}}v_{i}$

which implies that $||(A^{T}A)^{-1}A^{T}e||\gg||e||$ when $A$ is ill-
conditioned.

From these examples, it is clear that there is a trade-off between accuracy
and stability for a given reconstructor $\Psi$. In the next paragraphs, we
study such relationship.

###  2.1 Stability and Local Lipschitz constant

To simplify the computation, from now on we assume we are working with
(locally) Lipschitz continuous reconstructors, as it is common in the
literature. First of all, recover the definition of Local Lipschitz continuous
reconstructors.

######  Definition 2.6.

Given $\mathcal{Y}\subseteq\mathbb{R}^{m}$ and $\epsilon>0$, we define the
$\epsilon$-Lipschitz constant of $\Psi$ over $\mathcal{Y}$ as:

$\displaystyle
L^{\epsilon}(\Psi,\mathcal{Y})=\sup_{\begin{subarray}{c}y\in\mathcal{Y},z\in\mathbb{R}^{m}\\\
||z-y||\leq\epsilon\end{subarray}}\frac{||\Psi(z)-\Psi(y)||}{||z-y||}$ (2.4)

If $L^{\epsilon}(\Psi,\mathcal{Y})<\infty$ for some $\epsilon>0$, then $\Psi$
is said to be locally Lipschitz continuous.

Since we are working with problem of the form $y^{\delta}=Ax^{gt}+e$, we are
interested in the cases where $\mathcal{Y}=Rg(A,\mathcal{X})$, the range of
$A$ over $\mathcal{X}$. When this is the case, the definition of
$L^{\epsilon}(\Psi,\mathcal{Y})$ can be simplified. Indeed, since
$y\in\mathcal{Y}$, then $\exists x^{gt}\in\mathcal{X}$ such that $y=Ax^{gt}$.
Moreover, each $z\in\mathbb{R}^{m}$ with $||z-y||\leq\epsilon$ can be
characterized by $z=Ax^{gt}+e$ for some $e\in\mathbb{R}^{m}$ with
$||e||\leq\epsilon$. In this case, (2.4) can be rewritten as:

$\displaystyle
L^{\epsilon}(\Psi,\mathcal{Y})=\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
||e||\leq\epsilon\end{subarray}}\frac{||\Psi(Ax^{gt}+e)-\Psi(Ax^{gt})||}{||e||}$

From now on we will always assume $\mathcal{Y}=Rg(A,\mathcal{X})$.

The importance of the Local Lipschitz constant is that it is strongly related
to the stabiliy constant of $\Psi$. Indeed, if $y=Ax^{gt}\in\mathcal{Y}$ is
corrupted by additional noise $e$ with $||e||\leq\epsilon$, then
$L^{\epsilon}(\Psi,\mathcal{Y})$ represents the maximum variation of the
reconstruction via $\Psi$ around $y$, when corrupted with noise. This idea is
confirmed by the following Theorem.

######  Proposition 2.2.

If $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$ is a reconstructor with accuracy
$\eta^{-1}>0$, $\epsilon$-Lipschitz constant $L^{\epsilon}(\Psi,\mathcal{Y})$
with $\mathcal{Y}=Rg(A)$ and $||e||\leq\epsilon$, then

$\displaystyle||\Psi(Ax^{gt}+e)-x^{gt}||\leq\eta+L^{\epsilon}(\Psi,\mathcal{Y})||e||$
(2.5)

[1em]0em

###### Proof.

By the triangle inequality it follows that

$\displaystyle||\Psi(Ax^{gt}+e)-x^{gt}||\leq||\Psi(Ax^{gt}+e)-\Psi(Ax^{gt})||+||\Psi(Ax^{gt})-x^{gt}||$

since $||e||\leq\epsilon$, it follows from the definition of local Lipschitz
constant that

$\displaystyle||\Psi(Ax^{gt}+e)-\Psi(Ax^{gt})||\leq
L^{\epsilon}(\Psi,\mathcal{Y})||Ax^{gt}+e-Ax^{gt}||=L^{\epsilon}(\Psi,\mathcal{Y})||e||$

while the accuracy of $\Psi$ implies that

$\displaystyle||\Psi(Ax^{gt})-x^{gt}||\leq\eta$

Joining everything together gives

$\displaystyle||\Psi(Ax^{gt}+e)-x^{gt}||\leq
L^{\epsilon}(\Psi,\mathcal{Y})||e||+\eta$

$\blacksquare$

Consequently, we have

######  Corollary 2.1.

Under the assumptions of Proposition 2.2,

$C^{\epsilon}_{\Psi}\leq L^{\epsilon}(\Psi,\mathcal{Y})$

[1em]0em

###### Proof.

From inequality (2.5), we have

$||\Psi(Ax^{gt}+e)-x^{gt}||\leq\eta+L^{\epsilon}(\Psi,\mathcal{Y})||e||\iff
L^{\epsilon}(\Psi,\mathcal{Y})\geq\frac{||\Psi(Ax^{gt}+e)-x^{gt}||-\eta}{||e||}$

for any $x^{gt}\in\mathcal{X}$ and any $e\in\mathbb{R}^{m}$ with
$||e||\leq\epsilon$. Consequently, $L^{\epsilon}(\Psi,\mathcal{Y})$ is a
majorant of the set

$\Bigl{\\{}\frac{||\Psi(Ax^{gt}+e)-x^{gt}||-\eta}{||e||};x^{gt}\in\mathcal{X},||e||\leq\epsilon\Bigr{\\}}$

Since $C^{\epsilon}_{\Psi}$ is defined as the supremium of the same set, by
the minimality of the supremium we have $C^{\epsilon}_{\Psi}\leq
L^{\epsilon}(\Psi,\mathcal{Y})$. $\blacksquare$

We remark that we can state that $\Psi$ is $\epsilon$-stable if
$L^{\epsilon}(\Psi,\mathcal{Y})<1$. In addition, we prove in Appendix A an
important relation between the singular values of $A$,
$\\{\sigma_{1},\dots,\sigma_{n}\\}$ where
$\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{n}$, and the local Liptshitz
constant of an $\eta^{-1}$-accurate reconstructor $\Psi$:

$\displaystyle
L^{\epsilon}(\Psi,\mathcal{Y})\geq\frac{1}{\epsilon}\Bigl{(}\frac{\eta-2\sigma_{t}\eta}{\sigma_{t}}\Bigr{)}$
(2.6)

for any $t<rank(A)$. The inequality (2.6) shows a core issue appearing when we
apply an algorithm to solve an ill-posed linear problem $y^{\delta}=Ax^{gt}+e$
directly: there is no way to assess whether the the used solver is a stable
reconstructor for the noisy datum. On the other hand, as a consequence of
(2.6), a necessary condition for $L^{\epsilon}(\Psi,\mathcal{Y})<1$ is that
for any $t<rank(A)$,

$\displaystyle\frac{1}{\epsilon}\Bigl{(}\frac{\eta-2\sigma_{t}\eta}{\sigma_{t}}\Bigr{)}<1\iff\eta-2\sigma_{t}\eta<\epsilon\sigma_{t}\iff\sigma_{t}>\frac{\eta}{2\eta+\epsilon}.$
(2.7)

This necessary condition does not hold whenever the singular values of $A$
decay below the $\frac{\eta}{2\eta+\eta}$ threshold.

###  2.2 Accuracy vs. stability trade-off

In case of a full rank matrices $A$, we can also derive a relation between
accuracy and stability, becoming particularly interesting when $A$ is ill-
conditioned.

######  Lemma 2.1.

Let $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$ be an $\eta^{-1}$-accurate
reconstructor applied on the inverse problem $y^{\delta}=Ax^{gt}+e$ and let
$\epsilon>0$ be fixed. Then, if $rank(A)=n$, for any $x^{gt}\in\mathcal{X}$,
$\exists\tilde{e}\in\mathbb{R}^{m}$ with $||\tilde{e}||<\epsilon$ such that

$||\Psi(Ax^{gt}+\tilde{e})-x^{gt}||\geq||A^{\dagger}\tilde{e}||-\eta$ (2.8)

[1em]0em

###### Proof.

Since, for each $x^{gt}\in\mathcal{X}$, $Ax^{gt}\in\mathcal{Y}$ and since
$\mathcal{Y}$ is locally connected, then for any $\epsilon>0$ there is an
$\tilde{e}\in\mathbb{R}^{m}$ with $||\tilde{e}||\leq\epsilon$ such that
$Ax^{gt}+\tilde{e}\in\mathcal{Y}$. Thus, $\exists x^{\prime}\in\mathcal{X}$
such that $Ax^{gt}+\tilde{e}=Ax^{\prime}$. Consequently,

$\displaystyle||\Psi(Ax^{gt}+\tilde{e})-x^{gt}||=||\Psi(Ax^{\prime})-x^{gt}||\geq||x^{\prime}-x^{gt}||-||\Psi(Ax^{\prime})-x^{\prime}||\geq||x^{\prime}-x^{gt}||-\eta$

Since, by definition, $Ax^{gt}+\tilde{e}=Ax^{\prime}$, it holds that
$\tilde{e}=A(x^{gt}-x^{\prime})$, which implies that
$x^{gt}-x^{\prime}=A^{\dagger}\tilde{e}$. Consequently,

$\displaystyle||\Psi(Ax^{gt}+\tilde{e})-x^{gt}||\geq||x^{\prime}-x^{gt}||-\eta=||A^{\dagger}\tilde{e}||-\eta$

$\blacksquare$

######  Theorem 2.1 (Trade-off Theorem).

Under the same assumptions of Theorem 2.1, it holds that

$C^{\epsilon}_{\Psi}\geq\frac{||A^{\dagger}\tilde{e}||-2\eta}{||\tilde{e}||}=\frac{||A^{\dagger}\tilde{e}||}{||\tilde{e}||}-\frac{2\eta}{||\tilde{e}||}$
(2.9)

[1em]0em

###### Proof.

By the definition of $\epsilon$-stability constant, if
$||\tilde{e}||\leq\epsilon$ is the perturbation defined on Lemma 2.1, we have

$\displaystyle\begin{split}C^{\epsilon}_{\Psi}&=\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
e:||e||\leq\epsilon\end{subarray}}\frac{||\Psi(Ax^{gt}+e)-x^{gt}||-\eta}{||e||}\geq\sup_{x^{gt}\in\mathcal{X}}\frac{||\Psi(Ax^{gt}+\tilde{e})-x^{gt}||-\eta}{||\tilde{e}||}\geq\\\
&\geq\sup_{x^{gt}\in\mathcal{X}}\frac{||A^{\dagger}\tilde{e}||-2\eta}{||\tilde{e}||}=\frac{||A^{\dagger}\tilde{e}||-2\eta}{||\tilde{e}||}\end{split}$
(2.10)

$\blacksquare$

As a consequence of Theorem 2.1, there is a trade-off between $\eta$ and
$C^{\epsilon}_{\Psi}$, meaning that $C^{\epsilon}_{\Psi}$ grows if $\eta$
shrinks and, in particular, if $\eta\to 0$, then
$C^{\epsilon}_{\Psi}\to\frac{||A^{\dagger}e||}{||e||}$ if $A$ is ill-
conditioned.

##  3 Regularized reconstructors and Neural Networks

It is common to consider a reconstructor $\Psi$ with known accuracy and
stability properties and to approximate it with a sequence of approximators
$\\{\Psi_{\theta}\\}$. Examples of this practice are iterative algorithms and
neural networks approximating the solution of the inverse problem.
Interestingly, the stability of $\Psi_{\theta}$ is strongly related to the
stability of the approximated reconstructor $\Psi$. Indeed, the following
Theorem holds:

######  Theorem 3.1 (Approximation Theorem for Reconstructors).

Let $\Psi$ be an $\eta^{-1}$-accurate reconstructor and
$\\{\Psi_{\theta}\\}_{\theta\in\mathbb{R}^{s}}$ a set of reconstructors with
accuracy $\eta_{\theta}^{-1}$ for any $\theta$, we define:

  1. 1.

For any $\theta\in\mathbb{R}^{s}$

$\Delta(\theta):=\sup_{x^{gt}\in\mathcal{X}}||\Psi_{\theta}(Ax^{gt})-\Psi(Ax^{gt})||.$

  2. 2.

For any $\theta\in\mathbb{R}^{s}$

$\Delta_{\epsilon}(\theta):=\sup_{y^{\delta}\in\mathcal{Y}^{\epsilon}}||\Psi_{\theta}(y^{\delta})-\Psi(y^{\delta})||$

with $\epsilon>0$, $y^{\delta}=Ax^{gt}+e,e\sim\mathcal{N}(0,\delta^{2}I)$.

If $\Delta(\theta)\to 0$ when $\theta\to\theta^{*}$, then:

$\lim_{\theta\to\theta^{*}}\eta_{\theta}=\eta.$ (3.1)

Moreover, if $\Delta_{\epsilon}(\theta)\to 0$ when
$\theta\to\theta_{\epsilon}^{*}$, then:

$\lim_{\theta\to\theta_{\epsilon}^{*}}C^{\epsilon}_{\Psi_{\theta}}=C^{\epsilon}_{\Psi}.$
(3.2)

[1em]0em

###### Proof.

Consider $x^{gt}\in\mathcal{X}$. Since

$\displaystyle||\Psi_{\theta}(Ax^{gt})-x^{gt}||\leq||\Psi(Ax^{gt})-x^{gt}||+||\Psi_{\theta}(Ax^{gt})-\Psi(Ax^{gt})||$

and

$\displaystyle||\Psi_{\theta}(Ax^{gt})-x^{gt}||\geq||\Psi(Ax^{gt})-x^{gt}||-||\Psi_{\theta}(Ax^{gt})-\Psi(Ax^{gt})||$

it holds that

$|\>||\Psi_{\theta}(Ax^{gt})-x^{gt}||-||\Psi(Ax^{gt})-x^{gt}||\>|\leq||\Psi_{\theta}(Ax^{gt})-\Psi_{(}Ax^{gt})||\leq\Delta(\theta)$

which implies that
$||\Psi_{\theta}(Ax^{gt})-x^{gt}||\to||\Psi(Ax^{gt})-x^{gt}||$ as
$\Delta(\theta)\to 0$ and consequently, $\eta_{\theta}\to\eta$ as
$\Delta(\theta)\to 0$.

Now, consider $\epsilon>0$ and a corruption $e\in\mathbb{R}^{m}$ with
$||e||\leq\epsilon$. Following the computation above, it holds

$\displaystyle|\>||\Psi_{\theta}(Ax^{gt}+e)-x^{gt}||-||\Psi(Ax^{gt}+e)-x^{gt}||\>|\leq\Delta_{\epsilon}(\theta)$

which implies that
$||\Psi_{\theta}(Ax^{gt}+e)-x^{gt}||\to||\Psi(Ax^{gt}+e)-x^{gt}||$ for
$\Delta_{\epsilon}(\theta)\to 0$. Consequently, for
$\Delta_{\epsilon}(\theta)\to 0$,

$\displaystyle
C^{\epsilon}_{\Psi_{\theta}}=\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
e:||e||\leq\epsilon\end{subarray}}\frac{||\Psi_{\theta}(Ax^{gt}+e)-x^{gt}||-\eta_{\theta}}{||e||}\to\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
e:||e||\leq\epsilon\end{subarray}}\frac{||\Psi(Ax^{gt}+e)-x^{gt}||-\eta}{||e||}=C^{\epsilon}_{\Psi}$

which concludes the proof. $\blacksquare$

######  Corollary 3.1.

For any $\theta\in\mathbb{R}^{s}$, it holds:

$\eta_{\theta}\leq\eta+\Delta(\theta)$ (3.3)

[1em]0em

###### Proof.

Consider $x^{gt}\in\mathcal{X}$. Then

$\displaystyle||\Psi_{\theta}(Ax^{gt})-x^{gt}||\leq||\Psi_{\theta}(Ax^{gt})-\Psi(Ax^{gt})||+||\Psi(Ax^{gt})-x^{gt}||$

since $||\Psi_{\theta}(Ax^{gt})-\Psi(Ax^{gt})||\leq\Delta(\theta)$ by
hypothesis, while $||\Psi(Ax^{gt})-x^{gt}||\leq\eta$ since $\Psi$ is
$\eta^{-1}$-accurate, then

$||\Psi_{\theta}(Ax^{gt})-x^{gt}||\leq\eta+\Delta(\theta)$ (3.4)

which shows that $\eta^{\prime}\leq\eta+\Delta(\theta)$. $\blacksquare$

In the following, we will consider Neural Networks as the family
$\\{\Psi_{\theta}\\}_{\theta\in\mathbb{R}^{s}}$.

###  3.1 Neural Networks as reconstructors

Given $\mathcal{S}\subseteq\mathcal{X}$, consider the dataset
$\mathbb{D}=\\{(y^{\delta}_{i},x^{gt}_{i});x^{gt}_{i}\in\mathcal{S}\\}_{i=1}^{N_{\mathbb{D}}}$
of images such that, for any $i=1,\dots,N_{\mathbb{D}}$,
$y^{\delta}_{i}=Ax^{gt}_{i}+e_{i}$, $e_{i}\sim\mathcal{N}(0,\delta^{2}I)$.

######  Definition 3.1.

Given a neural network architecture $\mathcal{A}=(\nu,S)$ where
$\nu=(\nu_{0},\nu_{1},\dots,\nu_{L})\in\mathbb{N}^{L+1}$,
$\nu_{0}=m,\nu_{L}=n$, defines the width of each layer and
$S=(S_{1,1},\dots,S_{L,L}),S_{j,k}\in\mathbb{R}^{\nu_{j}\times\nu_{k}}$ is the
set of matrices representing the skip connections, we define the parametric
family of neural network reconstructors with architecture $\mathcal{A}$,
parameterized by $\theta\in\mathbb{R}^{s}$, as

$\mathcal{F}_{\theta}^{\mathcal{A}}=\\{\Psi_{\theta}:\mathbb{R}^{m}\to\mathbb{R}^{n};\theta\in\mathbb{R}^{s}\\}$

where

$\displaystyle\begin{cases}z^{0}=y\\\
z^{l+1}=\rho(W^{l}z^{l}+b^{l}+\sum_{k=1}^{l}S_{l,k}z^{k})\quad\forall
l=0,\dots,L-1\\\ \Psi_{\theta}(y^{\delta})=z^{L}\end{cases}$ (3.5)

and $W^{l}\in\mathbb{R}^{\nu_{l+1}\times\nu_{l}}$ is the weight matrix,
$b^{l}\in\mathbb{R}^{\nu_{l+1}}$ is the bias vector.

Training a neural network requires choosing a loss function
$\ell:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}_{+}$, and finding the
reconstructor $\Psi_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}$ solving the
minimization problem

$\min_{\Psi_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}\frac{1}{N_{\mathbb{D}}}\sum_{i=1}^{N_{\mathbb{D}}}\ell(\Psi_{\theta}(y^{\delta}_{i}),x^{gt}_{i})$
(3.6)

One of the most common choices for $\ell$ is simply the mean squared error
(MSE) function, defined as

$\ell(x,x^{\prime})=MSE(x,x^{\prime})=\frac{||x-x^{\prime}||_{2}^{2}}{n}$
(3.7)

Note that the minimization of the loss function
$MSE(\Psi_{\theta}(y^{\delta}),x^{gt})=MSE(\Psi_{\theta}(y^{\delta}),\Psi^{\dagger}(Ax^{gt}))$
is equivalent to minimizing the error $\Delta(\theta)$ introduced in Theorem
3.1. Since we have shown that $\Psi^{\dagger}$ is $\infty$-accurate, Theorem
3.1 implies that the accuracy of the neural network solving (3.6) tends to
infinity. However, as a consequence of Theorem 2.1:

$C^{\epsilon}_{\Psi_{\theta}}\approx\frac{||A^{\dagger}\tilde{e}||}{||\tilde{e}||}.$

We observe that if $A$ is ill-conditioned, $C^{\epsilon}_{\Psi_{\theta}}\gg
1$. It is thus clear that it is desirable to find a way to improve the
accuracy-stability trade-off in neural networks. Recently, more then one
approach has been proposed in this direction, not always with theoretical
guarantees.

A common technique to improve robustness in neural networks is to modify the
loss function $\ell$ to enforce consistency, as developed in [12, 16]. The
idea is to define $\Psi_{\theta}$ such as to

$\Psi_{\theta}(y^{\delta})\in\arg\min_{x\in\mathbb{R}^{n}}\frac{1}{2}||Ax-y^{\delta}||^{2}$
(3.8)

Unfortunately, it has been proved in [11], that this approach does not fully
solve instabilities.

Another common approach to reduce instabilities in neural network is
adversarial training [26], which can be implemented by data augmentation [21,
10]. The idea is to use equivariant transformations of $A$ to generate new
data. Unfortunately, for ill-condition matrices, enlarging the dataset
potentially increases the instability [11].

A third approach, called noise injection, is to add noise of different levels
to the input of the network during training. It has been proved in [6] that
this is equivalent to add Tikhonov regularization to the loss function of the
network. As explained in [3], this technique improves the stability of the
resulting network. Unfortunately, it is not clear to which extent the accuracy
of the resulting model is affected by noise injection and, moreover, how much
noise should be added to any input to maximize the trade-off between the
stability and accuracy.

In the following, we propose new approaches to improve the accuracy-stability
trade-off in neural networks as solvers of (1.2).

###  3.2 Tikhonov Regularization as Reconstructor

As an alternative to $\Psi^{\dagger}$ we introduce here a stable reconstructor
$\Psi$ of (1.2). A common choice is the Tikhonov regularized reconstructor
[24, 25], defined as:

$\Psi^{\lambda,L}(y^{\delta})=\arg\min_{x\in\mathbb{R}^{n}}\frac{1}{2}||Ax-y^{\delta}||^{2}+\frac{\lambda}{2}||Lx||^{2}$
(3.9)

where $\lambda>0$ is the regularization parameter, whereas
$L\in\mathbb{R}^{d\times n}$ is a matrix such that $\ker A\cap\ker L=\\{0\\}$.
$L$ is usually chosen to be the identity or the forward-difference operator.
The solution of (3.9) can be computed as the solution of the Normal Equations:

$\displaystyle\Psi^{\lambda,L}(y^{\delta})=(A^{T}A+\lambda
L^{T}L)^{-1}A^{T}y^{\delta}$ (3.10)

for any $y\in\mathcal{Y}^{\epsilon}$, where the kernel condition over $A$ and
$L$ is required to make $(A^{T}A+\lambda L^{T}L)$ invertible.

It is not hard to show that for any $L\in\mathbb{R}^{d\times n}$ and any
$\epsilon>0$, we can choose $\lambda>0$ such that $\Psi^{\lambda,L}$ is
$\epsilon$-stable. Indeed, it holds:

######  Proposition 3.1.

Let $\epsilon>0$ and $L\in\mathbb{R}^{d\times n}$. Then $\exists\lambda>0$
such that

$\displaystyle C^{\epsilon}_{\Psi^{\lambda,L}}\leq
L^{\epsilon}(\Psi^{\lambda,L},\mathcal{Y})<1$ (3.11)

[1em]0em

###### Proof.

For any $\lambda>0$ and any $y^{\delta}\in\mathcal{Y}^{\epsilon}$, similarly
to (3.10), it can be shown that

$\displaystyle\Psi^{\lambda,L}(y^{\delta})=\frac{1}{\lambda}\Bigl{(}\frac{1}{\lambda}A^{T}A+L^{T}L\Bigr{)}^{-1}A^{T}y^{\delta}$

Consequently, for any $y^{\delta}\in\mathcal{Y}^{\epsilon}$, it holds that

$\displaystyle\Psi^{\lambda,L}(y^{\delta})\to
0\quad\text{for}\quad\lambda\to\infty$

Then

$\displaystyle
L^{\epsilon}(\Psi^{\lambda,L},\mathcal{Y})=\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
||e||\leq\epsilon\end{subarray}}\frac{||\Psi^{\lambda,L}(Ax^{gt}+e)-\Psi^{\lambda,L}(Ax^{gt})||}{||e||}\to
0\quad\text{for}\quad\lambda\to\infty$

Which implies that $\forall\alpha>0$, $\exists\bar{\lambda}>0$ such that
$\forall\lambda>\bar{\lambda}$,
$L^{\epsilon}(\Psi^{\lambda,L},\mathcal{Y})<\alpha$. Choosing $\alpha=1$ gives
the required result. $\blacksquare$

###  3.3 Better conditioning implies better reconstructors: the advent of ReNN

Neural networks tend to be unstable since their training procedure
approximates the pseudo-inverse reconstructor $\Psi^{\dagger}$, as already
observed. As a consequence, we can improve neural networks stability by re-
defining the training problem, forcing the convergence of $\Psi_{\theta}$
towards Tikhonov regularized reconstructor $\Psi^{\lambda,L}$, which has been
shown to be $\epsilon$-stable for a suitable value of $\lambda$.

To this aim, we propose to modify the training procedure to learn
$\Psi_{\theta}$ solving

$\min_{\Psi_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}\sum_{i=1}^{N_{\mathbb{D}}}||\Psi_{\theta}(y^{\delta}_{i})-\Psi^{\lambda,L}(y^{\delta}_{i})||_{2}^{2},$
(3.12)

which corresponds to the minimization of $\Delta_{\epsilon}(\theta)$ in
Theorem 3.1. Hence, for Theorem 3.1 both the accuracy and the
$\epsilon$-stability constant of $\Psi_{\theta}$ approximately matches the
ones of $\Psi^{\lambda,L}$. In particular, for suitable values of $\lambda>0$,
$C^{\epsilon}_{\Psi_{\theta}}<1$. We will refer to such a network as the
Regularized Neural Network (ReNN) and we will indicate it as
$\Psi_{\theta}^{\lambda,L}$. The pseudocode of the algorithm to compute
$\Psi^{\lambda,L}_{\theta}$ is given in Algorithm 1.

Data: A collection
$\\{x^{gt}_{i}\\}_{i=1}^{N_{\mathbb{D}}}\subseteq\mathcal{X}$ of datapoints,
$A\in\mathbb{R}^{m\times n}$ and $\Psi^{\lambda,L}$.

Result: A trained ReNN $\Psi_{\theta}^{\lambda,L}$.

for  _$i\leftarrow 1:N_{\mathbb{D}}$_ do

Compute $y^{\delta}_{i}\leftarrow Ax^{gt}_{i}+e$;

Compute $\Psi^{\lambda,L}(y^{\delta}_{i})$;

Collect it into $\mathbb{D}^{\lambda,L}$;

end for

Given a neural network architecture $\mathcal{A}$;

Solve
$\arg\min_{\Psi_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}\sum_{i=1}^{N_{\mathbb{D}}}||\Psi_{\theta}(y^{\delta}_{i})-\Psi^{\lambda,L}(y^{\delta}_{i})||_{2}^{2}=\arg\min_{\Psi_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}\frac{1}{N_{\mathbb{D}}}\sum_{i=1}^{N_{\mathbb{D}}}MSE(\Psi_{\theta}(y^{\delta}_{i}),\Psi^{\lambda,L}(y^{\delta}_{i}))$

Algorithm 1 Regularized Neural Network (ReNN)

##  4 Stabilizers

We observe that $\Delta(\theta)$ strongly influences the assessment of the
$\Psi_{\theta}^{\lambda,L}$ stability which is in fact not guaranteed in case
of large values of $\Delta$. However, the Universal Approximation Theorem [22,
4] states that for any $\Delta(\theta)>0$ there always exists a neural network
$\Psi_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}$ (with suitable
$\mathcal{A}$), such that
$||\Psi^{\lambda,L}_{\theta}-\Psi^{\lambda,L}||\leq\Delta(\theta)$.
Interestingly, achieving the desired neural network may be undoable because of
the theory-practice gap of numerical computations, as observed in [7, 1]. It
becomes therefore necessary to design a new reconstructor, directly facing the
stability issue more than the approximation task. To this aim, we now
introduce a stabilization strategy.

######  Definition 4.1.

A continuous functions $\phi:\mathbb{R}^{m}\to\mathbb{R}^{t}$ is an
$\epsilon$-stabilizer of a reconstructor
$\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$ if:

  1. 1.

$\forall e\in\mathbb{R}^{m}$ with $||e||\leq\epsilon$, $\exists\
C^{\epsilon}_{\phi}\in[0,1)$ and $\exists\ e^{\prime}\in\mathbb{R}^{n}$ with
$||e^{\prime}||=C^{\epsilon}_{\phi}||e||$ such that

$\phi(Ax+e)=\phi(Ax)+e^{\prime}.$ (4.1)

  2. 2.

$\exists\ \gamma:\mathbb{R}^{t}\to\mathbb{R}^{n}$ such that
$\Psi=\gamma\circ\phi$.

In this case, the reconstructor $\Psi$ is said $\epsilon$-stabilized. The
smallest constant $C^{\epsilon}_{\phi}$ for which the definition holds is the
stability constant of the stabilizer $\phi$.

Note that, in the definition of $\epsilon$-stabilizer, we are requiring for
$\phi$ alone only a stability condition as in (4.1) and we do not ask anything
about its accuracy. We also remark that, by Definition 4.1, we disentangle the
problem of finding an accurate and a stable reconstructor.

Interestingly, given a $\epsilon$-stabilized reconstructor
$\Psi=\gamma\circ\phi$, we can estimate the $\epsilon$-stability constant
$C_{\Psi}^{\epsilon}$ of $\Psi$ by the constant $C_{\phi}^{\epsilon}$ defined
in Definition 4.1 and the local Lipschitz constant of $\gamma$ as proved in
the following proposition

######  Proposition 4.1.

Let $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$ be an $\epsilon$-stabilized
reconstructor, $\Psi=\gamma\circ\phi$. If $C_{\phi}^{\epsilon}$ is the
constant defined in (4.1), $L^{\epsilon}(\gamma,\mathcal{T})$ is the local
Lipschitz constant of $\gamma$ with $\mathcal{T}=\phi(\mathcal{Y})$, then

$C^{\epsilon}_{\Psi}\leq L^{\epsilon}(\gamma,\mathcal{T})C_{\phi}^{\epsilon}.$

[1em]0em

###### Proof.

Let $x^{gt}\in\mathcal{X}$ and $||e||\leq\epsilon$. Then,

$||\Psi(Ax^{gt}+e)-x^{gt}||=||\gamma(\phi(Ax^{gt}+e))-x^{gt}||$ (4.2)

Since $\phi$ is a stabilizer, $\phi(Ax^{gt}+e)=\phi(Ax^{gt})+e^{\prime}$ with
$||e^{\prime}||\leq C_{\phi}^{\epsilon}||e||$. Thus

$\displaystyle\begin{split}||(\gamma\circ\phi)(Ax^{gt}+e)-x^{gt}||&=||\gamma(\phi(Ax^{gt})+e^{\prime})-x^{gt}||\leq\\\
&\leq||\gamma(\phi(Ax^{gt})+e^{\prime})-\gamma(\phi(Ax^{gt}))||+||\gamma(\phi(Ax^{gt}))-x^{gt}||\leq\\\
&\leq\eta+L^{\epsilon}(\gamma,\mathcal{T})||e^{\prime}||\leq\eta+L^{\epsilon}(\gamma,\mathcal{T})C_{\phi}^{\epsilon}||e||\end{split}$
(4.3)

Which implies that $L^{\epsilon}(\gamma,\mathcal{T})C_{\phi}^{\epsilon}$ is a
majorant of the set

$\Bigl{\\{}\frac{||\Psi(Ax^{gt}+e)-x^{gt}||-\eta}{||e||};x^{gt}\in\mathcal{X},||e||\leq\epsilon\Bigr{\\}}$

Since $C^{\epsilon}_{\Psi}$ is defined as the supremium of the same set, by
the minimality of the supremium we have $C^{\epsilon}_{\Psi}\leq
L^{\epsilon}(\gamma,\mathcal{T})C_{\phi}^{\epsilon}$. $\blacksquare$

The Proposition above implies the following result:

######  Theorem 4.1 (Stabilization Condition).

For $\epsilon>0$, let $\phi_{1}:\mathbb{R}^{m}\to\mathbb{R}^{t}$ be an
$\epsilon$-stabilizer of $\Psi_{1}=\gamma_{1}\circ\phi_{1}$ and let $\Psi_{2}$
be another reconstructor such that $\Psi_{1}$, $\Psi_{2}$ have the same
accuracy $\eta^{-1}>0$. If

$\displaystyle
C^{\epsilon}_{\phi_{1}}\in\Bigl{[}0,\frac{C^{\epsilon}_{\Psi_{2}}}{L^{\epsilon}(\gamma_{1},\mathcal{T})}\Bigr{)}$
(4.4)

where $C_{\phi_{1}}^{\epsilon}$ is the constant defined in (4.1), then
$C_{\Psi_{1}}^{\epsilon}<C_{\Psi_{2}}^{\epsilon}$.

[1em]0em

###### Proof.

By hypothesis, both $\Psi_{1}=\gamma_{1}\circ\phi_{1}$ and $\Psi_{2}$ are
$\eta^{-1}$-accurate. Since

$C^{\epsilon}_{\phi_{1}}\in\Bigl{[}0,\frac{C^{\epsilon}_{\Psi_{2}}}{L^{\epsilon}(\gamma_{1},\mathcal{T})}\Bigr{)}$

and
$C^{\epsilon}_{\Psi_{1}}<L^{\epsilon}(\gamma_{1},\mathcal{T})C^{\epsilon}_{\phi_{1}}$
by Proposition 4.1, then

$C^{\epsilon}_{\Psi_{1}}\leq
L^{\epsilon}(\gamma_{1},\mathcal{T})C^{\epsilon}_{\phi_{1}}<L^{\epsilon}(\gamma_{1},\mathcal{T})\frac{C^{\epsilon}_{\Psi_{2}}}{L^{\epsilon}(\gamma_{1},\mathcal{T})}=C^{\epsilon}_{\Psi_{2}}$
(4.5)

which concludes the proof. $\blacksquare$

###  4.1 Disentangle accuracy and stability with stabilizers

The introduction of stabilizers open to the possibility of separately
designing a stable reconstructor and an accurate one. In particular, given
$\eta>0$, if one can define a function $\phi$ such that it exists a
reconstructor $\gamma$ with $\Psi=\gamma\circ\phi$ being $\eta^{-1}$-accurate,
then the lower $C^{\epsilon}_{\phi}$ is, the lower the stability constant of
$C^{\epsilon}_{\Psi}$ is.

Consequently, one can think that it would be great to consider
$\phi(Ax^{gt})=c$ for any $x^{gt}\in\mathcal{X}$. Since we already proven that
the constant reconstructor has stability constant equal to zero, than the
resulting $\epsilon$-stabilized reconstructor $\Psi=\gamma\circ\phi$ has to be
great. Unfortunately, as we will see, the accuracy of any
$\epsilon$-stabilized reconstructor, strongly depends on how far is $\phi$
from an injective function. Indeed, it holds

######  Proposition 4.2.

Let $\phi:\mathbb{R}^{m}\to\mathbb{R}^{t}$ be an $\epsilon$-stabilizer for an
$\eta^{-1}$-accurate reconstructor $\Psi=\gamma\circ\phi$. Let

$\displaystyle\sigma(\phi):=\sup\\{||x_{1}-x_{2}||;x_{1},x_{2}\in\mathcal{X},\phi(Ax_{1})=\phi(Ax_{2})\\}$
(4.6)

the maximum distance between elements in $\mathcal{X}$ that are reconstructed
as the same element by $\phi$. Then

$\eta^{-1}\leq\frac{2}{\sigma(\phi)}$ (4.7)

[1em]0em

###### Proof.

Let $x_{1},x_{2}\in\mathcal{X}$ such that $\phi(Ax_{1})=\phi(Ax_{2})$. Then,

$\displaystyle||x_{1}-x_{2}||$
$\displaystyle\leq||\phi(Ax_{1})-x_{1}||+||\phi(Ax_{1})-x_{2}||=$ (4.8)
$\displaystyle=||\phi(Ax_{1})-x_{1}||+||\phi(Ax_{2})-x_{2}||\leq 2\eta$ (4.9)

which implies that

$\eta\geq\frac{||x_{1}-x_{2}||}{2}.$ (4.10)

Since (4.10) holds for any $x_{1},x_{2}$ with $\phi(Ax_{1})=\phi(Ax_{2})$,
then it holds for $\sigma(\phi)$, concluding the proof. $\blacksquare$

Note that, if $\phi$ is the constant reconstructor, $\sigma(\phi)=\infty$,
which implies that for any $\gamma$, the accuracy of $\Psi=\gamma\circ\phi$
will be zero. Similarly, if $\sigma(\phi)<\infty$, then for any
$\eta^{-1}\in(0,\frac{2}{\sigma(\phi)}]$, if $\gamma$ is a suitable continuous
function, then $\Psi=\gamma\circ\phi$ is $\eta^{-1}$-accurate. A flexible way
to define such a $\gamma$, which is mathematically complex to describe, is to
use a neural network.

###  4.2 Iterative algorithm as stabilizers for neural networks

Consider the variational regularized reconstructor $\Psi^{\lambda,L}$
introduced in Section 3.2. Suppose to compute a sequence of approximators to
$\Psi^{\lambda,L}$, meaning a sequence $\\{\phi_{k}\\}_{k\in\mathbb{N}}$ such
that $\phi_{k}\to\Psi^{\lambda,L}$ for $k\to\infty$, i.e.

$\lim_{k\to\infty}||\phi_{k}-\Psi^{\lambda,L}||_{L^{2}}=0.$

This is not hard to do. A simple way is to consider a convergent iterative
algorithm:

$\begin{cases}x^{0}\in\mathbb{R}^{n}\\\
x^{k+1}=\mathcal{T}_{k}(x^{k}).\end{cases}$

For example, CGLS [14] algorithm converges to the Tikhonov regularized
reconstructor $\phi_{k}$. Then, for any $k\in\mathbb{N}$ we can define

$\phi_{k}=\circ_{i=1}^{k}\mathcal{T}_{i}$ (4.11)

Clearly, as a consequence of the convergence of $\mathcal{T}_{k}$,
$\\{\phi_{k}\\}_{k\in\mathbb{N}}$ is a sequence of approximators to
$\Psi^{\lambda,L}$. Now, given an algorithm
$\\{\mathcal{T}_{k}\\}_{k\in\mathbb{N}}$ and its associated sequence of
approximators $\\{\phi_{k}\\}_{k\in\mathbb{N}}$, we can define an
$\epsilon$-stabilizer by fixing $K\in\mathbb{N}$ and computing $\phi_{k}$.

######  Proposition 4.3.

Given a reconstructor $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$ with local
Lipschitz constant $L^{\epsilon}(\Psi,\mathcal{Y})<1$ and a sequence of
approximators $\\{\phi_{k}\\}_{k\in\mathbb{N}}$ to $\Psi$. Then there exists
$K\in\mathbb{N}$ such that for any $k\geq K$, $\phi_{k}$ is an
$\epsilon$-stabilizer.

[1em]0em

###### Proof.

Consider $x^{gt}\in X$ and $e\in\mathbb{R}^{m}$ with $||e||\leq\epsilon$. To
prove the Theorem, we need to show that for $k\geq K$,

$\phi_{k}(Ax^{gt}+e)=\phi_{k}(Ax^{gt})+e^{\prime}$

with $||e^{\prime}||=C^{\epsilon}_{\phi}||e||$ with
$C^{\epsilon}_{\phi}\in[0,1)$. Let
$e^{\prime}:=\phi_{k}(Ax^{gt}+e)-\phi_{k}(Ax^{gt})$. Then

$\displaystyle||e^{\prime}||=||\phi_{k}(Ax^{gt}+e)-\phi_{k}(Ax^{gt})||\leq
L^{\epsilon}(\phi_{k},\mathcal{Y})||Ax^{gt}+e-Ax^{gt}||=L^{\epsilon}(\phi_{k},\mathcal{Y})||e||$
(4.12)

which implies that $C^{\epsilon}_{\phi_{k}}\leq
L^{\epsilon}(\phi_{k},\mathcal{Y})$. Since $\\{\phi_{k}\\}_{k\in\mathbb{N}}$
is a sequence of approximators to $\Psi$, then for any $k\in\mathbb{N}$, there
is a constant $c_{k}$ such that $||\phi_{k}(y^{\delta})-\Psi(y^{\delta})||\leq
c_{k}$ with $c_{k}\to 0$ as $k\to\infty$. Consequently, it holds,

$\displaystyle\begin{split}&L^{\epsilon}(\phi_{k},\mathcal{Y})=\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
||e||\leq\epsilon\end{subarray}}\frac{||\phi_{k}(Ax^{gt}+e)-\phi_{k}(Ax^{gt})||}{||e||}\\\
&\leq\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
||e||\leq\epsilon\end{subarray}}\frac{||\phi_{k}(Ax^{gt}+e)-\Psi(Ax^{gt}+e)||+||\phi_{k}(Ax^{gt})-\Psi(Ax^{gt})||+||\Psi(Ax^{gt}+e)-\Psi(Ax^{gt})||}{||e||}\\\
&\leq\sup_{\begin{subarray}{c}x^{gt}\in\mathcal{X}\\\
||e||\leq\epsilon\end{subarray}}\frac{||\Psi(Ax^{gt}+e)-\Psi(Ax^{gt})||+2c_{k}}{||e||}\end{split}$
(4.13)

which implies that $L^{\epsilon}(\phi_{k},\mathcal{Y})\to
L^{\epsilon}(\Psi,\mathcal{Y})$ as $k\to\infty$. Since
$L^{\epsilon}(\Psi,\mathcal{Y})<1$, this implies that $\exists K\in\mathbb{N}$
such that for any $k\geq K$, $L^{\epsilon}(\phi_{k},\mathcal{Y})<1$. For those
values of $k$, $C^{\epsilon}_{\phi_{k}}\leq
L^{\epsilon}(\phi_{k},\mathcal{Y})<1$ which implies that $\phi_{k}$ is an
$\epsilon$-stabilizer. $\blacksquare$

Consider as a case of study the Tikhonov regularized reconstructor
$\Psi^{\lambda,L}$. We already proved in Proposition 3.1 that
$L^{\epsilon}(\Psi,\mathcal{Y})<1$ for some $\lambda>0$. Moreover, from (3.9)
we know that that $\Psi^{\lambda,L}(y^{\delta})$, is the solution of the
normal equation system

$(A^{T}A+\lambda L^{T}L)x=A^{T}y^{\delta}.$ (4.14)

A well-known and widely used iterative algorithm to solve (4.14) is the
Conjugate Gradient for Least Squares (CGLS) [14]. Moreover, by the convergence
proprieties of CLGS, $\\{\phi_{k}^{\lambda,L}\\}_{k\in\mathbb{N}}$ is a
sequence of approximators to $\Psi^{\lambda,L}$. Thus, by Proposition 4.3, it
defines a stabilizer for a neural network $\Psi_{\theta}$. Consequently, we
can join either NN or ReNN with this stabilizer, obtaining what we will call
StNN and StReNN, respectively. In Algorithms 2 and 3, an implementation of
StNN and StReNN can be found.

Data: A collection
$\mathbb{D}=\\{(y^{\delta}_{i},x^{gt}_{i}\\}_{i=1}^{N_{\mathbb{D}}}$ of
datapoints, a convergence algorithm $\mathcal{T}_{k}$ for $\Psi^{\lambda,L}$,
an integer $k\in\mathbb{N}$ such that $\phi_{k}$ is an $\epsilon$-stabilizer.

Result: A trained StNN $\Psi_{\theta}$.

for  _$i\leftarrow 1:N_{\mathbb{D}}$_ do

Compute $x^{k}_{i}=\phi_{k}(y^{\delta})$;

Collect $(x^{k}_{i},x_{i}^{gt})$ into $\mathbb{D}^{St}$;

end for

Given a neural network architecture $\mathcal{A}$;

Solve
$\gamma_{\theta}=\arg\min_{\gamma_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}\sum_{i=1}^{N_{\mathbb{D}}}||\gamma_{\theta}(x^{k}_{i})-x^{gt}_{i}||_{2}^{2}=\arg\min_{\gamma_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}=\frac{1}{N_{\mathbb{D}}}\sum_{i=1}^{N_{\mathbb{D}}}MSE(\gamma_{\theta}(x^{k}_{i}),x^{gt}_{i})\;$
Define $\Psi_{\theta}=\gamma_{\theta}\circ\phi_{k}$.

Algorithm 2 Stabilized Neural Network (StNN)

Data: A collection
$\mathbb{D}=\\{(y^{\delta}_{i},x^{gt}_{i})\\}_{i=1}^{N_{\mathbb{D}}}$ of
datapoints, a convergence algorithm $\mathcal{T}_{k}$ for $\Psi^{\lambda,L}$,
an integer $k\in\mathbb{N}$ such that $\phi_{k}$ is an $\epsilon$-stabilizer.

Result: A trained StReNN $\Psi_{\theta}$.

for  _$i\leftarrow 1:N_{\mathbb{D}}$_ do

Compute $x^{k}_{i}=\phi_{k}^{\lambda,L}(y^{\delta})$;

Compute $x^{\lambda,L}_{i}=\Psi^{\lambda,L}(y^{\delta}_{i})$;

Collect $(x^{k}_{i},x^{\lambda,L}_{i})$ them into $\mathbb{D}^{StRe}$;

end for

Given a neural network architecture $\mathcal{A}$;

Solve
$\gamma_{\theta}=\arg\min_{\gamma_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}\sum_{i=1}^{N_{\mathbb{D}}}||\gamma_{\theta}(x^{k}_{i})-x^{\lambda,L}_{i}||_{2}^{2}=\arg\min_{\gamma_{\theta}\in\mathcal{F}_{\theta}^{\mathcal{A}}}=\frac{1}{N_{\mathbb{D}}}\sum_{i=1}^{N_{\mathbb{D}}}MSE(\gamma_{\theta}(x^{k}_{i}),x^{\lambda,L}_{i})\;$
Define $\Psi_{\theta}=\gamma_{\theta}\circ\phi_{k}$.

Algorithm 3 Stabilized Regularized Neural Network (StReNN)

##  5 Experimental Setup

To test that the estimated results work in practice, we did three kind of
experiments, referred as A,B and C in the following, with the networks NN,
ReNN, StNN and StReNN, that we have defined in the previous sections. The
codes can be found in our GitHub repository for this work, at
https://github.com/loibo/stabilize_neural_networks.

We consider as a test case image deblurring [15], which is a fundamental
problem in imaging. In this particular case $Ax=K\ast x$, where $K$ is the
blurring operator. In our experiments, we use the $11\times 11$ Gaussian blur
kernel with variance $\sigma_{G}=1.3$, i.e.

$\displaystyle
K_{i,j}=\begin{cases}e^{-\frac{1}{2}\frac{i^{2}+j^{2}}{\sigma_{G}^{2}}}\quad&i,j\in\\{-10,\dots,10\\}^{2}\\\
0&\text{otherwise}\end{cases}$ (5.1)

Note that, by enforcing periodic boundary conditions over $x^{gt}$, both
$Ax^{gt}$ and $A^{T}y^{\delta}$ can be efficiently computed by Fast Fourier
Transform (FFT).

###  5.1 Dataset

Our results have been tested on the famous GoPro dataset
(https://seungjunnah.github.io/Datasets/gopro), introduced in [20].
Originally, the dataset was composed by high resolution RGB images, which we
will use as ground-truth. We also mention that in the dataset each image is
associated with a blurring kernel, that we didn’t use for our experiments. We
did the following preprocessing steps on the images, in order to create our
new dataset:

  * •

Each $x^{gt}$ is cropped into $256\times 256\times 3$ patches without
overlapping, and each patch is added to our dataset.

  * •

We transformed all the $256\times 256\times 3$ images into grey-scale and we
normalized them in $[0,1]$ by:

$x^{\prime}=\frac{x-\min(x)}{\max(x)-\min(x)}.$

  * •

We splitted the dataset into train and test sets, with $N_{train}=2503$ and
$N_{test}=1111$.

  * •

For $i=1,\dots,N_{train}$ we then computed the blurred version of $x^{gt}_{i}$
by sampling $e_{i}\sim\mathcal{N}(0,\delta^{2}\odot I)$, and computing
$y^{\hat{\delta}}_{i}=K\ast x^{gt}_{i}+e_{i}$, where $K$ is given by 5.1.

###  5.2 Iterative solver for Tikhonov reconstructor

We used, in our experiments, Tikhonov reconstructor as defined in 3.9 with $L$
is the discretization of the gradient, defined by $L=D\otimes D$ where

$\displaystyle\begin{split}&D=\begin{bmatrix}-2&1&0&\dots&\dots&0\\\
1&-2&1&0&\dots&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&\dots&0&1&-2&1\\\ 0&\dots&\dots&0&1&-2\end{bmatrix}\in\mathbb{R}^{n\times
n}\end{split}$ (5.2)

We choose $\lambda$ heuristically and we efficiently computed the minimum in
$\Psi^{\lambda,L}(y^{\delta})$ by using the CGLS algorithm [14]. We applied it
to any corrupted datum $y^{\delta}_{i}$, and we stopped at iteration $K$, when
one of the following conditions are satisfied:

$\displaystyle\begin{cases}&||y^{\delta}-Ax^{(k)}||_{2}^{2}<tol_{f}\\\
&||x^{(k)}-x^{(k-1)}||_{2}^{2}<tol_{x}\end{cases}$ (5.3)

where $tol_{f}=tol_{x}=10^{-6}$ in the experiments. We will refer to the
solution $x^{(K)}$ of CGLS as $\Psi_{K}^{\lambda,L}(y^{\delta})$.

###  5.3 Experiment A

In order to test the stability of our frameworks with respect to noise on the
data, we train the networks to solve the noiseless task $y=Ax^{gt}$ and we
compute the stability constant $C^{\epsilon}_{\Psi}$ as in (2.2), where $e$,
such that $||e||\leq\epsilon$, is computed by sampling
$e\sim\mathcal{N}(0,\delta I)$ for $\delta=0.01$. To do that, we created the
following data set:
$\mathbb{D}=\\{(Ax_{i}^{gt},x_{i}^{gt})\\}_{i=1}^{N_{train}}$ and
$\mathbb{D}^{\lambda,L}=\\{(Ax_{i}^{gt},\Psi_{K}^{\lambda,L}(Ax_{i}^{gt})\\}_{i=1}^{N_{train}}$.
We remark that by minimizing the MSE either on $\mathbb{D}$ or
$\mathbb{D}^{\lambda,L}$ we are minimizing $\Delta(\theta)$ as defined in
Theorem 3.1. The four models NN, StNN, ReNN and StReNN are based on the same
end-to-end U-net architecture (details on the neural network used and its
training setup can be found in Appendix B). In StNN and StReNN the input is
pre-processed by a stabilizer obtained with $k=3$ iterations of the CGLS
algorithm applied on (3.9). In particular, we underline that NN and StNN are
trained on the dataset $\mathbb{D}$, whereas ReNN and StReNN on
$\mathbb{D}^{\lambda,L}$. Finally, for a fair comparison we have computed the
accuracy and the stability constant of $\Psi^{\lambda,L}$.

###  5.4 Experiment B

In order to understand the behaviour of the proposed methods when a
stabilization technique is applied on them. Among the stabilizing solutions
listed in Section 3.1, we have chosen noise injection, which is considered the
most reliable. Hence we have modified experiment A only by adding gaussian
noise with variance $\delta^{I}=0.025$ to the input of each network. We
observe that in StNN and StReNN the noise has been added after stabilizing the
input. Thanks to the improved robustness introduced by noise injection, in
this experiment we have tested the reconstructors with $e$ computed by
sampling $\mathcal{N}(0,\delta I)$ for $\delta=0.075$ and $\delta=0.105$.
Clearly, in this experiment the comparison with Tikhonov reconstructor is
meaningless, hence it is not reported in the results.

###  5.5 Experiment C

Theorem 3.1 states in (3.2) that, by reducing $\Delta_{\epsilon}(\theta)$ the
stability constant $C_{\Psi_{\theta}}^{\epsilon}$ converges to
$C_{\Psi}^{\epsilon}$. In this experiment we show that, if we choose Tikhonov
as reconstructor $\Psi$ we can modify the training of ReNN and StReNN methods
to minimize $\Delta_{\epsilon}(\theta)$ in place of $\Delta(\theta)$, aiming
at dramatically improve their stability. To do that, we modified
$\mathbb{D}^{\lambda,L}$ by introducing noise in the target, thus obtaining
$\hat{\mathbb{D}}^{\lambda,L}=\\{(Ax_{i}^{gt}+\hat{e},\Psi^{\lambda,L}(Ax_{i}^{gt}+\hat{e})\\}_{i=1}^{N_{train}}$
with $e\sim\mathcal{N}(0,\hat{\delta}I)$ with $\hat{\delta}=0.025$. In this
experiment we have tested the reconstructors with $e$ computed by sampling
$\mathcal{N}(0,\delta I)$ for $\delta=0.075$ and $\delta=0.105$.

##  6 Numerical Results

In order to estimate the accuracy and the stability constant of a given
reconstrcutor $\Psi$, in our experiments, we compute the following parameters:
the empirical accuracy $\hat{\eta}^{-1}$ and the empirical stability constant
$\hat{C}^{\epsilon}_{\Psi}$, respectively defined as

$\hat{\eta}=\sup_{x^{gt}\in\mathcal{S}}||\Psi(Ax^{gt})-x^{gt}||$ (6.1)

and

$\hat{C^{\epsilon}_{\Psi}}=\sup_{x^{gt}\in\mathcal{S}}\frac{||\Psi(Ax^{gt}+e)-x^{gt}||-\hat{\eta}}{||e||}$
(6.2)

where $\mathcal{S}\subseteq\mathcal{X}$ is the test set and $e$ is a noise
realization from $\mathcal{N}(0,\delta^{2}I)$ with $||e||\leq\epsilon$ which
is considered different for any datum $x^{gt}\in\mathcal{S}$.

In order to care about the stochasticity of our tests, we repeated $T=20$
times the experiments on test sets, with different realization of noise. In
the following, we report the maximum value of the computed parameters
$\hat{\eta}$ and $\hat{C^{\epsilon}_{\Psi}}$ over the T experiments.

###  6.1 Results of experiment A

The results obtained in experiment A, described in Section 5.3, are reported
in Table 1 and graphically shown in Figure 1. In Figure 1(a) we plot the
values of the stability constant $\hat{C^{\epsilon}_{\Psi}}$ versus the
accuracy $\hat{\eta}^{-1}$ for the considered reconstructors NN, ReNN, StNN,
StReNN. As expected, NN has an excellent accuracy but it lacks in stability.
ReNN appears as a good trade-off between accuracy and stability: it shows a
slightly higher value of $\hat{\eta}^{-1}$ with respect to NN, whereas its
constant $\hat{C}^{\epsilon}_{\Psi}$ is drammatically lowered (from a value
greater than 32 in NN to about 10 in ReNN). Concerning the introduction of a
stabilizer, we observe that StNN and StReNN are $\epsilon$-stable, since
$\hat{C}^{\epsilon}_{\Psi}<1$. Finally, by comparing with the results from the
variational reconstructor $\Psi_{K}^{\lambda,L}$, denoted as IS in the
following, we remark that the proposed StNN and StReNN improve its accuracy
while preserving the stability.

In Figure 1(b) we plot the reconstruction error with respect to $x^{gt}$ when
the value of $\delta$ varies in the interval $[0,0.02]$ for all the considered
reconstructors. We observe that errors of NN and ReNN methods increase very
rapidly: for values of $\delta$ near zero they are the lowest but they soon
are overcome. Moreover, we underline that the stabilizer introduced in StNN
and StReNN is very effective, since the curves have a gentle slope with
increasing $\delta$.

To show the statistical behavior over all the T=20 tests, in Figure 2 we
report the box plots of $\hat{\eta}^{-1}$ and $\hat{C}^{\epsilon}_{\Psi}$
(Figures 2(a) and 2(b) respectively). They show that the values lies in a very
small range, except for very few exceptions, confirming the robustness of the
proposed methods.

$\hat{\delta}=0$ | NN | StNN | ReNN | StReNN | IS
---|---|---|---|---|---
$\hat{\eta}^{-1}$ | 0.1287 | 0.1029 | 0.1296 | 0.0980 | 0.0940
$\hat{C^{\epsilon}_{\Psi}}(\delta=0.01)$ | 32.6398 | 0.9269 | 10.9700 | 0.7046 | 0.5392
Table 1: Values of empirical accuracy and stability constant obtained in
experiment A.

(a)

(b)

Figure 1: Results obtained in Experiment A. Left. Scatterplot of the stability
constant versus the accuracy. Right. Plot of the reconstruction error when the
value of $\delta$ varies in the interval $[0,0.02]$.

(a)

(b)

Figure 2: Boxplot of the results obtained in experiment A for accuracy (Left)
and stability constant (Right) with $T=20$ executions.

###  6.2 Results of experiment B

The results obtained in experiment B, described in Section 5.4, are reported
in Table 2 and graphically shown in Figure 3. We remark that the difference
between this experiment and the previous one is noise injection in the input
of each network. Since this technique stabilizes the results, we have tested
the methods with greater values of the noise variance ($\delta=0.075$ and
$\delta=0.105$). In Figure 3(a), we observe that the stability constant
$C^{\epsilon}_{\Psi}$ is now less than one for all the reconstructors.
Moreover, from the values of $C^{\epsilon}_{\Psi}$ reported in Table 2, it is
evident that the stabilized reconstructors StNN and StReNN benefit a lot of
noise injection. In Figure 3, as for experiment A, we plot the reconstruction
error for increasing noise variance. We observe the curves for NN and ReNN to
remain constant for $\delta$ varying in the interval $[0,0.025]$, showing that
noise injection prevent instability up to the noise level introduced as an
input of the network. After $\delta=0.025$, the curves for NN and ReNN grows
rapidly, surpassing StNN and StReNN which grows consistently slow for the
complete considered interval, showing great stability agains high level of
noise.

$\delta^{I}=0.025$ | NN | StNN | ReNN | StReNN
---|---|---|---|---
$\hat{\eta}^{-1}$ | 0.0764 | 0.0710 | 0.0761 | 0.0705
$\hat{C^{\epsilon}_{\Psi}}(\delta=0.075)$ | 0.7719 | 0.2408 | 0.7674 | 0.2495
$\hat{C^{\epsilon}_{\Psi}}(\delta=0.105)$ | 0.9632 | 0.3608 | 0.9516 | 0.3846
Table 2: Values of empirical accuracy and stability constant obtained in
experiment B.

(a)

(b)

Figure 3: Results obtained in Experiment B. Left. Scatterplot of the stability
constant versus the accuracy. Right. Plot of the reconstruction error when the
value of $\delta$ varies in the interval $[0,0.1]$.

###  6.3 Results of experiment C

The results obtained in experiment C, described in Section 5.5, are reported
in Table 3 and graphically shown in Figure 4. We remark that we aim at
experimentally verify Theorem 3.1 by showing that, if we train ReNN and StReNN
with as described in section 5.5, their accuracy and stability constant are
close to those $\Psi^{\lambda,L}$. We test the reconstructors with
$\delta=0.075$ and $\delta=0.105$. Figure 4(a) shows that the markers relative
to ReNN and StReNN are indeed close to the purple ones, associated to IS.
Table 3 confirms that the accuracy of ReNN and StReNN, as well as the
stability constant of ReNN, well approximates those of IS. StReNN appears even
more stable than IS, since its stability constant is lower than the one of IS
by approximately 0.07. These results are confirmed by Figure 4.

$\hat{\delta}=0.025$ | ReNN | StReNN | IS
---|---|---|---
$\hat{\eta}^{-1}$ | 0.0633 | 0.0639 | 0.0605
$\hat{C^{\epsilon}_{\Psi}}(\delta=0.075)$ | 0.1500 | 0.1238 | 0.2034
$\hat{C^{\epsilon}_{\Psi}}(\delta=0.105)$ | 0.2536 | 0.1865 | 0.2566
Table 3: Values of empirical accuracy and stability constant obtained in
experiment C.

(a)

(b)

Figure 4: Results obtained in Experiment C. Left. Scatterplot of the stability
constant versus the accuracy. Right. Plot of the reconstruction error when the
value of $\delta$ varies in the interval $[0,0.1]$.

##  7 Conclusions

After formalizing the concept of accuracy and stability of a reconstructor for
the solution of discrete linear inverse problems, we have performed a
theoretical analysis on neural network reconstructors. In particular, we
proved in Theorem 2.1 that it is not possible to increase stability without
decreasing accuracy, depending on a constant characterizing the solution
process. Balancing the trade-off between stability and accuracy we have
proposed new reconstructors: the ReNN approach increases the stability by
inheriting, in the network training, regularization from a model-based scheme,
the StNN framework stabilizes the solution process by reducing the impact of
the noise on the input data performing few iterations of a model-based
algorithm and, combining the two previous approaches, the StReNN scheme
approximates the entire execution of a model-based iterative regularized
solver. In the numerical experiments we compared our proposal with a standard
neural network reconstructor, which exhibits high instability to noise
perturbation on the data. The results show that our reconstructors improve the
stability of neural networks by reducing its stability constant of around
$92\%$, especially when stabilizers are introduced in StNN and StReNN, with
minimal accuracy loss of about $10-20\%$ (in one case we even gain accuracy
with respect to neural networks). Our methods outperform end-to-end neural
networks even when stabilization techniques such as noise injection are
introduced. We believe that this new approach to solve noisy linear inverse
problems with stable deep learning based tools is relevant in this field and
it can be further theoretically extended to more general problems and applied
to the solution of real imaging applications.

## References

  * [1] Ben Adcock and Nick Dexter.  The gap between theory and practice in function approximation with deep neural networks.  SIAM Journal on Mathematics of Data Science, 3(2):624–655, 2021\.
  * [2] Jaweria Amjad, Jure Sokolić, and Miguel RD Rodrigues.  On deep learning for inverse problems.  In 2018 26th European Signal Processing Conference (EUSIPCO), pages 1895–1899. IEEE, 2018.
  * [3] Simon Arridge, Peter Maass, Ozan Öktem, and Carola-Bibiane Schönlieb.  Solving inverse problems using data-driven models.  Acta Numerica, 28:1–174, 2019.
  * [4] Julius Berner, Philipp Grohs, Gitta Kutyniok, and Philipp Petersen.  The modern mathematics of deep learning.  arXiv preprint arXiv:2105.04026, 2021.
  * [5] Mario Bertero, Patrizia Boccacci, and Christine De Mol.  Introduction to inverse problems in imaging.  CRC press, 2021.
  * [6] Chris M Bishop.  Training with noise is equivalent to tikhonov regularization.  Neural computation, 7(1):108–116, 1995.
  * [7] Matthew J Colbrook, Vegard Antun, and Anders C Hansen.  Can stable and accurate neural networks be computed?–on the barriers of deep learning and smale’s 18th problem.  arXiv preprint arXiv:2101.08286, 2021.
  * [8] Heinz Werner Engl, Martin Hanke, and Andreas Neubauer.  Regularization of inverse problems, volume 375.  Springer Science & Business Media, 1996.
  * [9] Davide Evangelista, Elena Morotti, and Elena Loli Piccolomini.  Rising a new framework for few-view tomographic image reconstruction with deep learning.  arXiv preprint arXiv:2201.09777, 2022.
  * [10] Zalan Fabian, Reinhard Heckel, and Mahdi Soltanolkotabi.  Data augmentation for deep learning based accelerated mri reconstruction with limited data.  In International Conference on Machine Learning, pages 3057–3067. PMLR, 2021.
  * [11] Nina M Gottschling, Vegard Antun, Ben Adcock, and Anders C Hansen.  The troublesome kernel: why deep learning for inverse problems is typically unstable.  arXiv preprint arXiv:2001.01258, 2020.
  * [12] Harshit Gupta, Kyong Hwan Jin, Ha Q Nguyen, Michael T McCann, and Michael Unser.  Cnn-based projected gradient descent for consistent ct image reconstruction.  IEEE transactions on medical imaging, 37(6):1440–1453, 2018.
  * [13] Per Christian Hansen.  Discrete inverse problems: insight and algorithms.  SIAM, 2010.
  * [14] Per Christian Hansen.  Discrete inverse problems: insight and algorithms.  SIAM, 2010.
  * [15] Per Christian Hansen, James G Nagy, and Dianne P O’leary.  Deblurring images: matrices, spectra, and filtering.  SIAM, 2006.
  * [16] Yixing Huang, Alexander Preuhs, Günter Lauritsch, Michael Manhart, Xiaolin Huang, and Andreas Maier.  Data consistent artifact reduction for limited angle tomography with deep learning prior.  In International workshop on machine learning for medical image reconstruction, pages 101–112. Springer, 2019.
  * [17] Chang Min Hyun, Seong Hyeon Baek, Mingyu Lee, Sung Min Lee, and Jin Keun Seo.  Deep learning-based solvability of underdetermined inverse problems in medical imaging.  Medical Image Analysis, 69:101967, 2021.
  * [18] Elena Morotti, Davide Evangelista, and Elena Loli Piccolomini.  A green prospective for learned post-processing in sparse-view tomographic reconstruction.  Journal of Imaging, 7(8):139, 2021.
  * [19] Jennifer L Mueller and Samuli Siltanen.  Linear and nonlinear inverse problems with practical applications.  SIAM, 2012.
  * [20] Seungjun Nah, Tae Hyun Kim, and Kyoung Mu Lee.  Deep multi-scale convolutional neural network for dynamic scene deblurring.  In CVPR, 07 2017.
  * [21] Daniel Obmann, Linh Nguyen, Johannes Schwab, and Markus Haltmeier.  Augmented nett regularization of inverse problems.  Journal of Physics Communications, 5(10):105002, 2021.
  * [22] Allan Pinkus.  Approximation theory of the mlp model in neural networks.  Acta numerica, 8:143–195, 1999.
  * [23] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [24] Otmar Scherzer, Markus Grasmair, Harald Grossauer, Markus Haltmeier, and Frank Lenzen.  Variational methods in imaging.  2009\.
  * [25] Andrei Nikolaevich Tikhonov, AV Goncharsky, VV Stepanov, and Anatoly G Yagola.  Numerical methods for the solution of ill-posed problems, volume 328.  Springer Science & Business Media, 1995.
  * [26] Zhengxia Zou, Tianyang Shi, Zhenwei Shi, and Jieping Ye.  Adversarial training for solving inverse problems in image processing.  IEEE Transactions on Image Processing, 30:2513–2525, 2021.

##  Appendix A Estimation of Local Lipschitz Constant

Consider a reconstructor $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$. We already
shown that its local Lipschitz constant $L^{\epsilon}(\Psi,\mathcal{Y})$ is
strongly related to its stability. In particular, we said that
$C_{\Psi}^{\epsilon}\leq L^{\epsilon}(\Psi,\mathcal{Y})$. Consequently, it is
mandatory to find a way to lower-bound $L^{\epsilon}(\Psi,\mathcal{Y})$, to
prove instability of reconstructors when $L^{\epsilon}(\Psi,\mathcal{Y})\geq
1$. To do that, assume $y\in\mathcal{Y}$. Then, we can define the local
Lipschitz constant of $\Psi$ over $y\in\mathcal{Y}$ as

$\displaystyle
L^{\epsilon}(\Psi,y):=L^{\epsilon}(\Psi,\\{y\\})=\sup_{z:||z-y||\leq\epsilon}\frac{||\Psi(z)-\Psi(y)||}{||z-y||}$

Note that, since
$L^{\epsilon}(\Psi,\mathcal{Y})=\sup_{y\in\mathcal{Y}}L^{\epsilon}(\Psi,y)$,
if we show that $L^{\epsilon}(\Psi,Ax^{gt})\geq 1$ for at least one
$x^{gt}\in\mathcal{X}$, then $L^{\epsilon}(\Psi,\mathcal{Y})\geq 1$.

As the following Theorem shows, the local Lipschitz constant of a
reconstructor $\Psi$, accurately solving an inverse problem
$y^{\delta}=Ax^{gt}+e$ is strongly related to the ill-conditioning of $A$.

######  Theorem A.1.

If $\Psi:\mathbb{R}^{m}\to\mathbb{R}^{n}$ is an $\eta^{-1}$-accurate
reconstructor, $x,x^{\prime}\in\mathcal{X}$ such that
$||A(x-x^{\prime})||<\eta$, then for any $\epsilon\geq\eta$, it holds

$\displaystyle
L^{\epsilon}(\Psi,Ax)\geq\frac{||x-x^{\prime}||-2\eta}{\epsilon}$ (A.1)

[1em]0em

###### Proof.

Let $y=Ax$. By the definition of $\epsilon$-Lipschitz constant we have

$\displaystyle L^{\epsilon}(\Psi,y)$
$\displaystyle=\sup_{z:||z-y||\leq\epsilon}\frac{||\Psi(z)-\Psi(y)||}{||z-y||}$

Since, by hypothesis, $||Ax-
Ax^{\prime}||=||A(x-x^{\prime})||<\eta\leq\epsilon$, we have:

$\displaystyle L^{\epsilon}(\Psi,y)$
$\displaystyle=\sup_{z:||z-y||\leq\epsilon}\frac{||\Psi(z)-\Psi(y)||}{||z-y||}\geq\frac{||\Psi(Ax^{\prime})-\Psi(Ax)||}{||A(x^{\prime}-x)||}=$
$\displaystyle=\frac{||\Psi(Ax^{\prime})-x^{\prime}-\Psi(Ax)+x+x^{\prime}-x||}{||A(x^{\prime}-x)||}\geq$
$\displaystyle\geq\frac{-||\Psi(Ax^{\prime})-x^{\prime}||-||\Psi(Ax)-x||+||x^{\prime}-x||}{||A(x^{\prime}-x)||}>$
$\displaystyle>\frac{||x^{\prime}-x||-2\eta}{\eta}\overset{(\eta\leq\epsilon)}{\geq}\frac{||x^{\prime}-x||-2\eta}{\epsilon}$

which implies that

$\displaystyle L^{\epsilon}(\Psi,y)>\frac{||x^{\prime}-x||-2\eta}{\epsilon}$

$\blacksquare$

Consequently, a sufficient condition that implies instability is that there
exists $x,x^{\prime}\in\mathcal{X}$ such that
$||A(x-x^{\prime})||<\eta\leq\epsilon$ and

$\displaystyle\begin{split}&\frac{||x-x^{\prime}||-2\eta}{\epsilon}\geq
1\iff\\\ &||x-x^{\prime}||\geq\epsilon+2\eta\end{split}$ (A.2)

This way, we moved the problem of proving instability of $\Psi$ to a problem
of finding points $x,x^{\prime}\in\mathcal{X}$ distant from each other such
that their difference lies close to the kernel of $A$. In the next Section, we
will prove that this condition is closely related to the conditioning of $A$
and, in particular, to how quick their singular values decay.

###  A.1 Existence of points far from each other lying close to the kernel

Let $x^{gt}\in\mathcal{X}$ and $y=Ax^{gt}$. Let
$(\sigma_{1},\dots,\sigma_{n})^{T}$ be the singular values of $A$ and
$(v_{1},\dots,v_{n})^{T}$ the left-singular vectors associated to
$(\sigma_{1},\dots,\sigma_{n})^{T}$. For any $r=1,\dots,n$, define
$V_{r}:=Span\\{v_{r},v_{r+1},\dots,v_{n}\\}$ the subspace defined by the left
singular vectors of $A$ smaller then $\sigma_{r}$. Consider the accuracy
$\eta^{-1}>0$ to be fixed.

######  Theorem A.2.

Let $t<n$ be a positive integer and $\eta>0$. Then there are points
$x,x^{\prime}\in\mathcal{X}$ such that
$||x-x^{\prime}||=\frac{\eta}{\sigma_{t}}$ but $||A(x-x^{\prime})||<\eta$. In
particular, for every $x\in\mathcal{X}$, the set of such points defines a
vector space of dimension $n-t$.

[1em]0em

###### Proof.

Let $t<n$ be fixed. Consider $x\in\mathcal{X}$, $\delta x\in V_{t+1}$ and
$x^{\prime}:=x+\delta x$.

Since $\delta x\in V_{t+1}$, then we can write it in coordinate with respect
to the base $\\{v_{t+1},\dots,v_{n}\\}$ of $V_{t+1}$ as

$\displaystyle\delta x=\sum_{i=t+1}^{n}a_{i}v_{i}$

Since the $V$ matrix of the SVD decomposition of $A$ is an orthogonal matrix,
its columns $\\{v_{1},\dots,v_{n}\\}$ forms an orthonormal base of
$\mathbb{R}^{n}$. Thus

$\displaystyle||x-x^{\prime}||=||\delta x||=\sqrt{\sum_{i=t+1}^{n}a_{i}^{2}}$

and

$\displaystyle||A(x^{\prime}-x)||=||A\delta x||=||U\Sigma V^{T}\delta x||$

but

$\displaystyle U\Sigma V^{T}\delta x=\sum_{i=1}^{n}\sigma_{i}(v_{i}^{T}\delta
x)u_{i}=\sum_{i=1}^{n}\sigma_{i}(v_{i}^{T}\sum_{j=t+1}^{n}a_{j}v_{j})u_{i}=\sum_{i=t+1}^{n}\sigma_{i}a_{i}u_{i}$

which implies that

$\displaystyle||A(x-x^{\prime})||=\sqrt{\sum_{i=t+1}^{n}\sigma_{i}^{2}a_{i}^{2}}\leq\sqrt{\sum_{i=t+1}^{n}\sigma_{t}^{2}a_{i}^{2}}=\sigma_{t}\sqrt{\sum_{i=t+1}^{n}a_{i}^{2}}=\sigma_{t}||\delta
x||$

As a consequence, if we choose $||\delta x||$ to be $||\delta
x||=\frac{\eta}{\sigma_{t}}$ we have that:

$\displaystyle||A(x-x^{\prime})||<\eta$
$\displaystyle||x-x^{\prime}||=||\delta x||=\frac{\eta}{\sigma_{t}}$

$\blacksquare$

In particular, given an $\eta^{-1}$-accurate reconstructor $\Psi$, we can find
$x,x^{\prime}\in\mathcal{X}$ such that $||A(x-x^{\prime})||<\eta$ and
$||x-x^{\prime}||=\frac{\eta}{\sigma_{t}}$ where $t$ is a positive integer
such that $t<n$. Consequently, joining Theorem A.1 and Theorem A.2, we have

$\displaystyle
L^{\epsilon}(\Psi,\tilde{y})\geq\frac{1}{\epsilon}\Bigl{(}||x-x^{\prime}||-2\eta\Bigr{)}=\frac{1}{\epsilon}\Bigl{(}\frac{\eta}{\sigma_{t}}-2\eta\Bigr{)}=\frac{1}{\epsilon}\Bigl{(}\frac{\eta-2\sigma_{t}\eta}{\sigma_{t}}\Bigr{)}$
(A.3)

which shows that any $\eta^{-1}$-accurate reconstructor $\Psi$ is unstable as
soon as

$\displaystyle\begin{split}&\frac{1}{\epsilon}\Bigl{(}\frac{\eta-2\sigma_{t}\eta}{\sigma_{t}}\Bigr{)}\geq
1\iff\\\ &\sigma_{t}\leq\frac{\eta}{\epsilon+2\eta}\end{split}$ (A.4)

Moreover, Theorem A.2 shows that given $x\in\mathcal{X}$ and $y=Ax$, the set
of points $x^{\prime}$ such that $||A(x-x^{\prime})||<\eta$ with distance
$\frac{\eta}{\sigma_{t}}$ from $x$ has dimension $n-t$ which implies that, if
the smallest $t$ such that (A.4) is low, the probability of getting a locally
unstable point $x^{\prime}$ is high.

##  Appendix B Neural Network Details

Figure 5: A diagram of the architecture used in the experiments.

As already mentioned, the neural network architecture was not a focus in this
work, since we was interested in understanding how to manipulate the stability
of the reconstructor by acting on the conditioning of the problem, with fixed
architecture. Moreover, in Appendix A we showed that the instability
phenomenon is unrelated with the architecture used to solve it, since it only
depends on the ill-conditioning of the matrix $A$ defining the forward problem
and the accuracy that the reconstructor achieve. Consequently, we decided to
define each network to be a slight modification of the original UNet [23]
architecture, where we added a residual connection between the input and the
output of the network, following what we already done in [18, 9].

###  B.1 Architecture design

As it is clear in the diagram in Figure 5, the UNet architecture is naturally
divided in resolution levels, defined by the successive application of Max
Pooling and Up Sampling layers. Let $L$ be the total number of levels.
Moreover, each level (except for the last one) is composed two symmetric
groups of convolutional layers, each containing $N_{c}$ layers with $N_{l}$,
$l=1,\dots,L$ filters. After each application of Max Pooling, the number of
channels is doubled, thus following the equation $N_{l+1}=2N_{l}$. Clearly, to
be consistent, the number of channel is halved after each Up Sampling. In each
level, the left-side group of convolution is connected to the right-side with
a skip connection, necessary to keep informations while we scale down over the
levels. Thus the triple $(L,N_{c},N_{1})$ uniquely defines a network of that
kind. In our experiments, we fixed $L=4$, $N_{c}=3$, $N_{1}=64$, as it is
represented in Figure 5.

###  B.2 Training details

Each network is trained for 50 epochs with Adam optimization algorithm and
fixed step size of $10^{-4}$. Following the theory developed in the previous
sections, the network is trained to minimize the Mean Squared Error (MSE)
between the output and the processed input, for each dataset, with no
additional modifications such as regularizers or the use of Dropout layers.
The batch size is set to the maximum value our GPU could handle, that is 8.
The training is executed on an NVIDIA RTX A4000 GPU card, requiring
approximately 1 to 2 hours to be completed.


# Objective Methods for Assessing Models for Wildfire Spread

Jeff Picka

###### Abstract

Models for wildfires must be stochastic if their ability to represent
wildfires is to be objectively assessed. The need for models to be stochastic
emerges naturally from the physics of the fire, and methods for assessing fit
are constructed to exploit information found in the time evolution of the burn
region.

##  1 Introduction

Management of wildfires requires the development of real-time forecasting
methods for the prediction of their evolution. Development of these models
requires acknowledging the difficulties associated with data acquisition, data
processing, and modeling of wildfire events. At present, many predictive
models are deterministic, and produce one forecast for every set of observed
fire conditions. When a single deterministic prediction is made, two major
questions arise: what does that single prediction represent, and how reliable
is it? There is no way to answer either of these questions unless the
deterministic model is replaced by a model which satisfies three conditions.

  1. 1.

The model predicts only a single aspect of the fire rather than predicting
many diverse aspects of the fire.

  2. 2.

The model must be stochastic, and make predictions in the form of estimated
probability distributions or probability maps.

  3. 3.

The fit of any model must be assessed by statistical methods after it has been
calibrated.

These conditions are imposed by our lack of ability to model the fire as a
physical process, by our inability to observe the details of the internal
structure of the fire, and by the difficulties associated with simulating
sequences of regions of random shape and size.

###  1.1 Terminology

A trajectory is a sequence of regions of irregular shape and extent which
describe the cumulative region burnt by a fire over time. Even if the
prediction problem does not directly involve the burned region, it will almost
certainly involve predicting a response statistic which is defined by the
shape and size of that region (e.g. cumulative value of timber burned,
expenditure required to extinguish the fire).

A stochastic model is a collection of stochastic processes whose outcomes may
represent fire trajectories. The parameters of the stochastic model are
collections of numbers which index the stochastic processes in the model. A
fitted (stochastic) model is the one stochastic process in the model which is
chosen to best represent the fire for which a forecast needs to be made. The
fitted stochastic model can be made to produce a sequence of realizations,
each of which may represent a fire trajectory.

Since fire trajectories are unpredictable, the burn regions at any fixed time
are unpredictable in shape, size, and orientation. At a fixed time, a sample
of these shapes is a sample from a random set process [1]. The stochastic
model for the fire spread is a time series of random set processes.

##  2 Stochastic Modeling of Fires

Useful models for fires must be stochastic because of the nature of wildfires
and our ability to model complex multiscale phenomena. The stochastic aspects
of wildfires are unavoidable, and cannot be ignored in model construction,
model fitting, or model assessment.

###  2.1 Models for Multiscale Physical Phenomena

The simplest version of a multiscale physical phenomenon consists of a
collection of fundamental elements at one scale (the microscale) which act
collectively to produce a distinctive pattern of behaviour at a much larger
scale of observation (the macroscale). At the microscale, the fundamental
elements may or may not be observable. Even if they are observable, it may be
difficult to model their behaviour. It is assumed that the macroscale
phenomena are always observable.

In the ideal case, the behaviour of the fundamental elements can be used to
derive the behaviour of the macroscale phenomena. The classical example of
this is the Gibbs canonical ensemble model of an ideal gas, from which the
basic results of classical thermodynamics for an ideal gas at equilibrium can
be derived [2]. If the state of an inert gas changes very slowly, then its
physical behaviour can be modeled well by results derived from the Gibbs
model. This model works well because the fundamental units are atoms of an
inert gas, which have simple structure and which interact rarely. The
simplicity of the microscale behaviour combined with the immense number of
atoms involved results in macroscale behaviour which can be comprehensively
described by a very small number of macroscopic state variables. Microscale
simplicity also allows the derivation of models for the relationships between
state variables, which do not need to include any consideration of macroscale
fluctuations induced by microscale behaviour.

Models for the steady-state laminar flow of Newtonian fluids and the steady-
state conduction of heat through solids are also successful multiscale models
[3]. In these models, the mechanisms of energy transfer through small volumes
is simple, and macroscale flow can be modeled using partial differential
equations. In the case of heat transfer by convection or turbulent fluid flow,
rigorous definition of macroscale properties from microscale behaviour is
impossible, since there are no mathematical models that can faithfully
represent the behaviour of the fundamental elements and their interactions. In
other cases such as the modeling of the flow of powders [4], the number of
fundamental elements involved is not sufficient to eliminate fluctuations in
properties at the macroscale.

In the modeling of mass, heat, and momentum transfer, it is rare to be able to
formulate a model which faithfully captures microscale phenomena and then
derives from the microscale model a useful description of the limiting
macroscale phenomena. Any model for macroscale behaviour depends on being able
to summarize that behaviour with appropriately defined state variables, and
then being able to construct a model which describes the relationships between
the state variables. To be of use, this model must not be seriously
inconsistent with the physics at the microscale, even if the microscale
behaviour is unobserved and unmodelable. This approach can be successful if
the model developer is modeling an inherently simple multiscale phenomenon,
manages to develop a useful set of state variables, manages to find a model to
relate them, and then finds that the resulting model produces predictions
which are close enough to reality to solve the problem of interest.

###  2.2 Modeling of Wildfires

Wildfires are not amenable to modeling by simple multiscale models which have
deterministic macroscale behaviour. The combustion takes place in a vapour
which forms above pyrolyzing fuel. The heat produced by the combustion drives
the formation of more vapour, as well as slowly changing the geometry of the
combustible materials from which the vapour is produced. The combustible
material occurs in many different forms on many different length scales, and
does not have homogenous structure. The heat evolved from combustion causes
turbulent flow of air and of combusting and non-combusting fuel vapour. These
turbulent flows are of sufficient strength to affect local weather. The
modeling of any one of these phenomena is very difficult on a small scale. The
violent and complex interactions of these processes cannot be observed except
on the largest scales.

The complexity of wildfires makes it impossible to imagine modeling them at
the microscale and then establishing the limiting macroscale behaviour. The
relatively small size of the macroscale relative to the microscale together
with the complexity of interactions between the fundamental elements ensures
that, no matter how much is known about the past history of that fire, there
will be many possible future fire trajectories which are consistent with that
history. Given what can be observed about a real fire, there is an ensemble of
fire trajectories which are consistent with these observations. Given our
inability to comprehensively model the fire, nature appears to choose at
random one of the fire trajectories from the ensemble. Any model for a fire
trajectory has to sample from that ensemble in the same way that the physical
process would, and must yield a probabilistic rather than a deterministic
prediction of future fire behaviour.

Stochastic modeling of wildfire behaviour is possible, since the large-scale
behaviour of the fires is more predictable than small-scale combustion
phenomena. Any model proposed will need to be stochastic, but there will not
be any rigorous theoretical argument suggesting exactly what form this model
should take. The trajectory of each real fire will be impossible to model, and
so any model must be based on a random sample of some form of approximate
trajectory. Construction of such models will require much creative programming
and physical intuition. When constructing models, it would be unreasonable to
assume that the first successful models will be useful in all prediction
problems, or that the selection among possible approximated fire trajectories
should be undertaken via a uniform distribution. Instead, it would be best to
focus on one prediction problem at a time and then try to unify approaches
once successful models have been found.

###  2.3 The Risks Associated With Not Using Stochastic Models

In many models now in use, the macroscale behaviour of a fire is simulated by
a set of partial differential equations whose numerical solution produces a
single approximate fire trajectory. This approach to modeling has many
inherent problems which can be avoided by using stochastic models.

The most difficult problem associated with a single predictive trajectory is
that it is not clear what that trajectory represents. If it is a single
trajectory from the ensemble, then is it a typical trajectory or one with very
unusual features when compared to other trajectories in the ensemble?
Alternatively, the deterministic model may produce some form of average of
many trajectories. This implies that at a fixed time, the burn region in the
trajectory is an average of many such regions from the ensemble. Any region
produced by a geometrical averaging process will be much simpler in shape than
any single burn region from a trajectory in the ensemble.

Even if the nature of the single prediction could be determined, in both cases
the single prediction overstates the accuracy of the prediction procedure. The
deterministic prediction can be used to produce graphics which can mislead a
naive observer as to how uncertain the prediction really is, and gives no
information whatsoever as to how typical that prediction is of others from
ensemble.

Given a single trajectory from any model, it is almost certain that the
prediction will differ from the observed trajectory of the fire. If this
difference is trivial to the eye, the model may be judged to be adequate in
the modeling of that fire. If major differences are observed, it may be
tempting to conclude that the model does not fit. To make such a judgment
objectively, especially in a case where the ensemble contains a great variety
of possible trajectories, it is necessary to make the judgment of model fit
based on comparison between the observed trajectory and many simulated
trajectories from the model. This process requires the use of a stochastic
model, and is the subject of the next section.

###  2.4 Constructing a Stochastic Model

If a deterministic model samples a single trajectory from the ensemble of
physically viable trajectories, then it can be used as the basis of a
stochastic model [5]. Stochastic elements can be added which represent lack of
ability to clearly observe initial and boundary conditions, but they can also
be added to compensate for aspects of fire propagation that are not fully
understood.

##  3 Assessment of Model Fit

Any stochastic model consists of a family of stochastic processes, only one of
which is chosen to represent the fire. There is no guarantee that any method
of calibration will produce a good choice, and so the fit of the chosen model
must be assessed by statistical means.

###  3.1 Fitting a stochastic model

Fitting a model involves finding a method to choose parameter values so that
the fitted model can provide useful predictions of the evolution of a
wildfire. This choice must be based on the past history of the fire.

The parameters can be divided into two classes. Physical parameters appear in
those aspects of the model derived from physics, and take the form of physical
constants such as the acceleration of gravity. These parameters are generally
set using values established by physical experiments, often in contexts very
different from wildfires. While these parameters are not generally allowed to
be set by a model-fitting procedure, allowing these values to be fit may be
useful in model assessment.

Calibration parameters are all parameters in the model which are not set by
physical arguments. They may be physical parameters which cannot be observed
or estimated, or they may be parameters which have no clear physical meaning
but which make the model flexible enough to fit the data. They can be set by
the guesses of experts who are attempting to make the model output look right,
or by statistical means.

Development of objective fitting procedures requires that the model predict
only one specific aspect of the fire. Once the specific aspect is chosen, it
is necessary to find a small set of response statistics which fully describes
that aspect. Seeking a model which works in all predictive problems would
require finding a small set of statistics which could describe or predict all
possible responses. It is doubtful that any such universal set of predictive
statistics exists.

Given a set of response statistics, a fitted stochastic model can be found by
comparing the observed response to simulated responses. A large number of
combinations of calibration parameters could be defined, and then a large
sample of realizations generated from each stochastic process. From each
sample, a mean response would be found. The fitted model would be the
stochastic process whose mean response was closest to that of the observed
fire. A response consisting of a time series of measurements would be best for
use in this type of fitting procedure.

It is also possible to fit a model by varying the calibration parameters until
the simulated responses look right. This method needs to be compared to
objective methods since, in the hands of an expert, it may work well in some
cases. This method may miss better parameter values, may be subject to
conscious or unconscious bias, or may produce conclusions which are impossible
to explain to others.

###  3.2 Assessing the fit of a fire spread model

Once a fire spread model is fit, the adequacy of the model must be assessed by
criteria appropriate to the prediction problem. This is necessary not only
because all models are descriptive and speculative, but also because there is
a severe risk that the fitting process will only match superficial
similarities between the fire and the model realizations.

Fit assessment for a fire spread model requires comparing a sample of
simulated trajectories with the trajectory of the fire that they are modeling.
If no evidence can be found of differences which could affect the response,
then the model can be accepted as being useful. If evidence of differences are
found, then these may suggest how the model can be improved.

To assess the fit, the first stage is to summarize aspects of the fire
trajectory using descriptive statistics. These statistics may be based on the
burn region at one time, or on a sequence of burn regions in time. They may be
based on the shape of the burn region or on observable and quantifiable
properties of the fire within it.

For a descriptive statistic to summarize essential elements of the prediction
problem, it is necessary for the response to depend on the statistic. This
dependency can only be investigated by conducting statistical tests on fitted
model output and then by attempting to identify associations. It would be
better to investigate associations using fire data, but this is impossible
unless the fire can be repeated many times under the same conditions. This
classification process can reduce the number of statistics which are used to
compare the model and the data. If a statistic is found not to affect the
response, then there is no need for its values to be consistent for the model
and the data. Since the classification is based on statistical inference,
there is always a risk that a descriptive statistic will be misclassified.

It is necessary to find a large number of descriptive statistics in any model
assessment problem. Many descriptive statistics will be spatial, and will
depend on the shape and size of regions. Statistics of this kind suffer from a
lack of power, and may not be able to distinguish between spatial processes
with visibly different realizations [6]. There is also no general theory to
suggest which statistics may be useful in finding differences between the
response in the fire and the response in the fitted model.

Finding descriptive statistics is a major challenge. Given expert concerns or
past model failings, it is possible to quantify these in a similar manner to
the development of response statistics. Using these statistics alone is not
good enough, since they may not identify differences that have not been
anticipated or observed in previous modeling efforts. Realizations must be
summarized by statistics that can identify aspects of model failure that are
difficult or impossible to see, even by experts. This requires developing a
large library of statistics, all of which describe different aspects of the
burn region. Analogues of second moment statistics, such as the pair
correlation function, the nearest neighbour function, the empty space
function, the $K-$function, and the 2-point correlation function [7] are
inadequate for this purpose, since the stochastic model is not Gaussian. These
statistics are very useful, but must be supplemented by other statistics
including (possibly) statistics based on triangulations, statistics based on
models for physical properties applied out of context [8], and statistics
which have yet to be invented.

Given a set of descriptive statistics which affect the response, it is
necessary to establish if they show any evidence that the model is wrong. This
requires comparing a sample of a collection of descriptive statistic values
from the observed fire with a sample which can be generated from many
realizations of the fitted model. This comparison could be made with one
descriptive statistic at a time, but better results are expected if small
collections of statistics are compared. These comparisons are difficult to
make reliably with the data from only a single fire, but this is unavoidable
since each fire takes place under different conditions. Comparisons could be
better made between data from controlled burns undertaken under fixed
experimental conditions, but there is no guarantee that these burns would
sample from ensembles similar to those associated with real fires.

The assessment of fit will be subject to unavoidable errors since it is based
on statistical tests. Investigations of association between descriptive
statistics and the response will misclassify statistics. There will be no
guarantee that the observed trajectory is typical of trajectories from
physical process, and so differences between the observed and model
trajectories may be found for a useful model. It is also possible that no
evidence of poor fit will be found in models which are seriously wrong. This
will happen if none of the statistics used to compare the model to the data
are capable of identifying the difference. These errors are also unavoidable
in fits by expert judgment.

##  4 Conclusions

Models which faithfully represent what is known and unknown about fire spread
must be stochastic. The goal of a stochastic model is not to reproduce any one
burn model exactly, but instead to produce simulated fire trajectories which
approximate fire trajectories from the ensemble of possible fires. They need
not capture every aspect of these fires, but they need to be able to capture
the distribution of the response over the ensemble of real fires. The
assessment of fit looks for evidence that the approximate trajectories are
inconsistent with the one trajectory that is observed. A model will be useful
if this general approach to fitting and assessment allows it to be predict
outcomes successfully from many different fires.

The objective approach to model fitting presents a number of challenges which
must be addressed before it can be usefully implemented. It requires the
development of new statistics for the description of fire trajectories. It
requires the development of a stochasticized version of a model like
Prometheus [9] which can usefully be fit to data. It will be necessary to use
historical data and simulations to determine what must be observed in order to
predict fire behaviour and to assess the fit of fire models. It may require
the acquisition or development of new equipment in order to gather the
information on which predictions are based.

Given the effort required, it may seem easier to continue fitting models by
eye. This approach implicitly assumes that experts will not be influenced by
any beliefs they may have about the fitness of model and that an expert will
be able to see all possible ways in which the model could be wrong. If these
assumptions are incorrect, then mistaken assumptions about fire behaviour
could live on in flawed models for far longer than they should. Without a
trajectory-based assessment of the quality of a fitted model, it will always
be possible that a fitted model will fail to capture the physics of fire
spread.

The use of stochastic methods is also essential to making useful predictions
about fire behaviour. If a model is purely deterministic, there is no way to
determine if its single prediction is typical of those from the ensemble of
possible fires or not. Also, there will no measure of uncertainty arising from
the model that can be assigned to the prediction. This could lead to a loss of
faith in modeling, due to failures in predictions for new fires. Deterministic
predictions could also mislead non-experts into overconfidence about what is
known and can be predicted about fire behaviour. Neither outcome would be
beneficial for the development of useful models for fire spread.

##  5 Acknowledgments

This work has been partly supported by NSERC but has mostly been supported by
GEOIDE.

## References

  * [1] I. Molchanov.  Theory of Random Sets.  Springer, 2005.
  * [2] R.A. Minlos.  Introduction to Mathematical Statistical Mechanics.  American Mathematical Society, 2000.
  * [3] R.B. Bird, W.E. Stewart, and E.N. Lightfoot.  Transport Phenomena.  Wiley, 1960.
  * [4] J. Kakalios.  Resource letter GP-1: Granular physics or non-linear dynamics in a sandbox.  Am. J. Phys., 73:8–22, 2003.
  * [5] T. Garcia, J. Braun, R. Bryce, and C. Tymstra.  Smoothing and bootstrapping the PROMETHEUS fire growth model.  Environmetrics, 19:836–848, 2008.
  * [6] A.J. Baddeley and B.W. Silverman.  A cautionary example on the use of 2nd-order methods for analyzing point patterns.  Biometrics, 40(4):1089–1093, 1984.
  * [7] D. Stoyan, W. S. Kendall, and J. Mecke.  Stochastic Geometry and Its Applications.  John Wiley & Sons, 2nd edition, 1995.
  * [8] J.D. Picka and T.D. Stewart.  A conduction-based statistic for description of orientation in anisotropic sphere packings.  In Powders and Grains 2005 Proceedings, pages 329–332. Balkema, 2005.
  * [9] C. Tymstra, R.W. Bryce, B.M. Wotton, and O.B. Armitage.  Development and structure of prometheus: the canadian wildland fire growth simulation model.  Technical Report NOR-X-417, Natural Resources Canada, Canadian Forestry Centre Edmonton, 2009.  forthcoming.