On the randomised stability constant for inverse problems

1.   Introduction

Inverse problems are the key to all experimental setups where the physical quantity of interest is not directly observable and must be recovered from indirect measurements. They appear in many different contexts including medical imaging, non-destructive testing, seismic imaging or signal processing. In mathematical terms, an inverse problem consists in the inversion of a linear or nonlinear map

which models how the quantity of interest $x$ belonging to a space $X$ is related to the measurements $y = T(x)$ in the space $Y$. The reader is referred to the many books on inverse problems for a comprehensive exposition (see, e.g., ^{[2,3,4,24,27,30,44,47]}).

Inverse problems can be ill-posed: the map $T$ may not be injective (i.e., two different $x_{1}$ and $x_{2}$ may correspond to the same measurement $T(x_{1}) = T(x_{2})$) or, when injective, $T^{-1}\colon{\rm{ran}} T\subseteq Y\to X$ may not be continuous (i.e., two different and not close $x_{1}$ and $x_{2}$ may correspond to almost identical measurements $T(x_{1})\approx T(x_{2})$). Various strategies have been introduced to tackle the issue of inversion in this setting, Tykhonov regularisation being the most famous method ^[17,29].

Our purpose is to investigate the role of randomisation in the resolution of inverse problems. By randomisation, we refer to the use of random measurements, or to the case of random unknowns. We do not mean the well-established statistical approaches in inverse problems, as in the Bayesian framework, where probability is used to assess the reliability of the reconstruction.

Even with this distinction, the wording "randomisation" may refer to many different concepts in the framework of inverse problems. A huge literature is devoted to the issue of randomisation in the measuring process: the unknown $x$ in $X$ is fixed and deterministic, and we choose the measurements randomly according to some suitable distribution. For example, compressed sensing ^[19] and passive imaging with ambient noise ^[20] belong to this class. Another very popular idea consists in randomising the unknown: in this scenario, we try to recover most unknowns $x$ in $X$, according to some distribution. For example, this is typically the situation when deep learning ^[21] is applied to inverse problems. We briefly review these instances in Appendix A.

This article deals with a more recent approach, introduced in the framework of observability or control theory: the so-called randomised observability constants ^[38,39,40]. We argue that, in contrast with the initial motivation for its introduction, the randomised observability constant is not necessarily indicative of the likelihood for randomised unknowns to be observable.

In the next section, we recall the notion of randomised observability constant and comment on its use in optimal design problems. In section 3, we reformulate the randomised observability constant in an abstract setting as randomised stability constant. In section 4 we show that when the classical (deterministic) stability constant is null, a positive randomised stability constant need not imply, as one could hope, that the inverse problem can be solved for most unknowns. In the course of our study, we make several observations on the properties of the randomised stability constant. Section 5 contains several concluding remarks and discusses possible future directions.

2.   Randomising initial data of PDEs

In this section we briefly review the randomised observability constant introduced in ^[38,39,40] and its main properties. This was motivated by the original idea of randomisation by Paley and Zygmund, which we now briefly discuss.

2.1.   On randomisation processes

In order to understand the principle of randomisation that will be at the core of this paper, it is useful to recall the historical result by Paley and Zygmund on Fourier series. Let $(c_{n})_{n\in\mathbb{Z}}$ be an element of $\ell^{2}(\mathbb{C})$ and $f$ be the Fourier series given by

where $\mathbb{T}$ denotes the torus $\mathbb{R}/(2\pi)$. According to the so-called Parseval identity, the function $f$ belongs to $L^{2}(\mathbb{T})$; furthermore, the coefficients $c_{n}$ can be chosen in such a way that $f$ does not belong to any $L^{q}(\mathbb{T})$ for $q > 2$. Some of the results obtained by Paley and Zygmund (see ^[35,36,37]) address the regularity of the Fourier series $f$. They show that if one changes randomly and independently the signs of the Fourier coefficients $c_{n}$, then the resulting random Fourier series belongs almost surely to any $L^{q}(\mathbb{T})$ for $q > 2$. More precisely, introducing a sequence $(\beta_{n}^{\nu})_{n\in\mathbb{Z}}$ of independent Bernoulli random variables on a probability space $(A, \mathcal{A}, \mathbb{P})$ such that

then, the Fourier series $f^{\nu}$ given by

belongs almost surely to $L^{q}(\mathbb{T})$ for all $q < +\infty$.

In ^[13], the effect of the randomisation on the initial data of solutions to dispersive equations was investigated. In that case, initial data and solutions of the PDE considered are expanded in a Hilbert basis of $L^{2}(\Omega)$ made of eigenfunctions $(\phi_{j})_{j\geq1}$ of the Laplace operator. The randomisation procedure consists hence in multiplying all the terms of the series decomposition by well-chosen independent random variables. In particular, regarding for instance the homogeneous wave equation with Dirichlet boundary conditions, one can show that for all initial data $(y^{0}, y^{1})\in H_{0}^{1}(\Omega)\times L^{2}(\Omega)$, the Bernoulli randomisation keeps the $H_{0}^{1}\times L^{2}$ norm constant. It is observed that many other choices of randomisation are possible. For instance a positive effect of the randomisation can be observed by considering independent centred Gaussian random variables with variance 1 (Gaussian randomisation). Indeed, it allows to generate a dense subset of the space of initial data $H_{0}^{1}(\Omega)\times L^{2}(\Omega)$ through the mapping

where $(y_{\nu}^{0}, y_{\nu}^{1})$ denotes the pair of randomised initial data, provided that all the coefficients in the series expansion of $(y^{0}, y^{1})$ are nonzero. Several other properties of these randomisation procedures are also established in ^[13].

2.2.   Randomised observability constant

We now review how randomisation appeared in the framework of inverse problems involving an observability inequality. The property of observability of a system is related to the following issue: How to recover the solutions of a PDE from the knowledge of partial measurements of the solutions. In what follows, we will concentrate on wave models, having in particular photoacoustic/thermoacoustic tomography imaging in mind.

Let $T > 0$ and $\Omega\subseteq\mathbb{R}^{d}$ be a bounded Lipschitz domain with outer unit normal $\nu$. We consider the homogeneous wave equation with Dirichlet boundary conditions

It is well known that, for all $(y^{0}, y^{1})\in H_{0}^{1}(\Omega)\times L^{2}(\Omega)$, there exists a unique solution $y\in C^{0}([0, T], H_{0}^{1}(\Omega))\cap C^{1}((0, T), L^{2}(\Omega))$ of (2.1) such that $y(0, x) = y^{0}(x)$ and $\partial_{t}y(0, x) = y^{1}(x)$ for almost every $x\in\Omega$. Let $\Gamma$ be a measurable subset of $\partial\Omega$, representing the domain occupied by some sensors, which take some measurements over a time horizon $[0, T]$.

The inverse problem under consideration reads as follows.

Inverse problem: Reconstruct the initial condition $(y^{0}, y^{1})$ from the knowledge of the partial boundary measurements

To solve this problem, we introduce the so-called observability constant: $C_{T}(\Gamma)$ is defined as the largest non-negative constant $C$ such that

for any solution $y$ of (2.1), where $H_{0}^{1}(\Omega)$ is equipped with the norm $\|u\|_{H_{0}^{1}(\Omega)} = \|\nabla u\|_{L^{2}(\Omega)}$. Then, the aforementioned inverse problem is well-posed if and only if $C_{T}(\Gamma) > 0$. In such a case, we will say that observability holds true in time $T$. Moreover, observability holds true within the class of $\mathcal{C}^{\infty}$ domains $\Omega$ if $(\Gamma, T)$ satisfies the Geometric Control Condition (GCC) (see ^[8]), and this sufficient condition is almost necessary.

Let us express the observability constant more explicitly. Fix an orthonormal basis (ONB) $(\phi_{j})_{j\geq1}$ of $L^{2}(\Omega)$ consisting of (real-valued) eigenfunctions of the Dirichlet-Laplacian operator on $\Omega$, associated with the negative eigenvalues $(-\lambda_{j}^{2})_{j\geq1}$. Then, any solution $y$ of (2.1) can be expanded as

where the coefficients $a_{j}$ and $b_{j}$ account for initial data. More precisely, we consider the ONB of $H_{0}^{1}(\Omega)\times L^{2}(\Omega)$ given by $\{\psi_{j}^{+}, \psi_{j}^{-}:j\ge1\}$, where

Expanding now the initial data with respect to this basis we can write $(y^{0}, y^{1}) = \sum_{j = 1}^{+\infty}a_{j}\psi_{j}^{+}+b_{j}\psi_{j}^{-}$, namely,

The corresponding solution to (2.1) is given by (2.3). In addition, Parseval's identity yields $\|(y^{0}, y^{1})\|_{H_{0}^{1}(\Omega)\times L^{2}(\Omega)}^{2} = \sum_{j = 1}^{+\infty}|a_{j}|^{2}+|b_{j}|^{2}$.

Then, the constant $C_{T}(\Gamma)$ rewrites

where $y(t, x)$ is given by (2.3).

The constant $C_{T}(\Gamma)$ is deterministic and takes into account any $(a_{j}), (b_{j})\in\ell^{2}(\mathbb{C})$, including the worst possible cases. Interpreting $C_{T}(\Gamma)$ as a quantitative measure of the well-posed character of the aforementioned inverse problem, one could expect that such worst cases do not occur too often; thus it would appear desirable to consider a notion of observation in average.

Motivated by the findings of Paley and Zygmund (see §2.1) and its recent use in another context ^[11,13], making a random selection of all possible initial data for the wave equation (2.1) consists in replacing $C_{T}(\Gamma)$ with the so-called randomised observability constant defined by

where

and $(\beta_{1, j}^{\nu})_{j\in\mathbb{N}}$ and $(\beta_{2, j}^{\nu})_{j\in\mathbb{N}}$ are two sequences of independent random variables of Bernoulli or Gaussian type, on a probability space $(A, \mathcal{A}, \mathbb{P})$ with mean 0 and variance 1. Here, $\mathbb{E}$ is the expectation in the probability space, and runs over all possible events $\nu$. In other words, we are randomising the Fourier coefficients $\{a_{j}, b_{j}\}_{j\ge1}$ of the initial data $(y^{0}, y^{1})$ with respect to the basis $\{\psi_{j}^{\pm}\}_{j\ge1}$.

The randomised observability constant was introduced in ^{[38,39,40,41,42]}. It can be expressed in terms of deterministic quantities (see [40,Theorem 2.2]).

Proposition 1. Let $\Gamma\subset\partial\Omega$ be measurable. We have

Proof. In view of Lemma 1 (see below), we have that

where $y_{j}^{\pm}$ is the solution to (2.1) with initial condition $\psi_{j}^{\pm}$. Thus

and in turn

which leads to our thesis.

We have $C_{T, \mathrm{rand}}(\Gamma)\geq C_{T}(\Gamma)$ (see Proposition 4 below). It has been noted in ^[41] that the observability inequality defining $C_{T, \mathrm{rand}}(\Gamma)$ is associated to a deterministic control problem for the wave equation (2.1), where the control has a particular form but acts in the whole domain $\Omega$.

Regarding $C_{T, \mathrm{rand}}(\Gamma)$, we refer to [42,Section 4] for a discussion on the positivity of this constant. The authors show that if $\Omega$ is either a hypercube or a disk, then $C_{T, \mathrm{rand}}(\Gamma) > 0$ for every relatively non-empty open subset $\Gamma$ of $\partial\Omega$. In particular, in some cases $C_{T}(\Gamma) = 0$ while $C_{T, \mathrm{rand}}(\Gamma) > 0$. This raised hopes that, even if recovering all unknowns is an unstable process, recovering most unknowns could be feasible, since apparently most unknowns are observable. This heuristic argument, mentioned amongst the motivations for the study of the optimisation of the randomised observability constant, was not investigated further in the aforementioned papers. This matter will be studied in the following sections, dedicated more generally on the possible use of such a constant for investigating the well-posed character of general inverse problems.

2.2.1.   Applications to optimal design problems

A larger observability constant $C_{T}(\Gamma)$ in (2.2) leads to a smaller Lipschitz norm bound of the inverse map. Therefore $C_{T}(\Gamma)$ can be used as the quantity to maximise when searching for optimal sensors' positions. However, this turns out to be somewhat impractical. When implementing a reconstruction process, one has to carry out in general a very large number of measures; likewise, when implementing a control procedure, the control strategy is expected to be efficient in general, but possibly not for all cases. Thus, one aims at exhibiting an observation domain designed to be the best possible in average, that is, over a large number of experiments. Adopting this point of view, it appeared relevant to consider an average over random initial data. In ^[38,39,40], the best observation is modelled in terms of maximising a randomised observability constant, which coincides with $C_{T, \mathrm{rand}}(\Gamma)$ when dealing with the boundary observation of the wave equation.

When dealing with internal observation of the wave equation on a closed manifold, it has been shown in ^[26] that the related observability constant reads as the minimum of two quantities: The infimum of the randomised observability constants over every orthonormal eigenbasis and a purely geometric criterion standing for the minimal average time spent by a geodesic in the observation set.

However, one should keep in mind that a large randomised constant may not be associated with a reconstruction method (see section 4).

3.   The randomised stability constant for abstract inverse problems

It is convenient to generalise the construction of the previous section to an abstract setting. In what follows, none of the arguments require the precise form of the forward operator related to the wave equation, as they rely solely on the structure of the randomised constant.

For the remainder of this paper, we let $X$ and $Y$ be separable infinite-dimensional Hilbert Spaces, and $P\colon X\to Y$ be an injective bounded linear operator.

If $P^{-1}\colon\operatorname{ran}P\to X$ is a bounded operator, the inverse problem of finding $x$ from $P(x)$ can be solved in a stable way for all $x\in X$, without the need for randomisation. This can be measured quantitatively by the constant

The smaller $C_{\mathrm{det}}$ is, the more ill-conditioned the inversion becomes.

On the other hand, when $P^{-1}$ is unbounded, the situation is different, and the inverse problem is ill-posed ^[24,30]. In this case, although the kernel of $P$ reduces to $\{0\}$, we have

Examples of such maps abound. The example that motivated our study was that introduced in section 2.2, with $X = H_{0}^{1}(\Omega)\times L^{2}(\Omega)$, $Y = L^{2}([0, T]\times\partial\Omega)$ and

where $\Gamma\subseteq\partial\Omega$ and $y$ is the solution to (2.1) with initial condition $(y^{0}, y^{1})$: If $\Gamma$ is not large enough, $P$ is still injective but $P^{-1}$ is unbounded ^[31]. Any injective compact linear operator satisfies (3.1).

Let us now introduce the randomised stability constant, which generalises the randomised observability constant to this general setting. We consider the class of random variables introduced in the last section. Choose an ONB $\mathrm{e} = \{\mathrm{e}_{k}\}_{k\in\mathbb{N}}$ of $X$ and write $x = \sum_{k = 1}^{\infty}x_{k}\mathrm{e}_{k}\in X$. We consider random variables of the form

where $\beta_{k}^{\nu}$ are i.i.d. complex-valued random variables on a probability space $(A, \mathcal{A}, \mathbb{P})$ with vanishing mean and variance $1$, so that $\mathbb{E}(|\beta_{k}^{\nu}|^{2}) = 1$ for every $k$. These include the Bernoulli and Gaussian random variables considered in the previous section. It is worth observing that, in the case of Bernoulli random variables, we have $|\beta_{k}^{\nu}|^{2} = 1$ for every $k$, so that $\|x^{\nu}\|_{X} = \|x\|_{X}$.

Definition 1. The randomised stability constant is defined as

As in the previous section, this constant should represent stability of the inverse problem for most unknowns. By definition, we have $C_{\mathrm{rand}}(\mathrm{e})\ge C_{\mathrm{det}}$. In general, this is a strict inequality: We will provide examples in section 4. This can be heuristically seen also by the following deterministic expression for the randomised stability constant.

Lemma 1. There holds

Proof. Since $\left(Y, \left\Vert \cdot\right\Vert _{Y}\right)$ is also a Hilbert space, we find

Since $\beta_{k}^{\nu}$ are i.i.d. with vanishing mean and such that $\mathbb{E}(|\beta_{k}^{\nu}|^{2}) = 1$, we obtain

which means $C_{\mathrm{rand}}(\mathrm{e})\ge\inf_{k}\left\Vert P\left(\mathrm{e}_{k}\right)\right\Vert _{Y}^{2}$. Choosing $x = \mathrm{e}_{k}$, we obtain (3.2).

4.   Exploiting $C_{\mathrm{rand}}$ in inverse problems

The aim of this section is to discuss the impact of the randomised observability constant in inverse problems. We would like to address the following question: How can the positivity of $C_{ \rm{rand}}$ be exploited in the solution of an inverse ill-posed problem? We will not fully address this issue, but rather provide a few positive and negative partial results.

We remind the reader that the randomisation introduced in the last section when (3.1) holds, was based on the point of view that the ratio in (3.1) is not "usually" small, and that, hopefully, in most cases inversion "should" be possible. It is worthwhile to observe that the subset of the $x\in X$ such that $\left\Vert Px\right\Vert _{Y}\geq c\left\Vert x\right\Vert _{X}$ for some fixed $c > 0$ is never generic in $X$, since $\left\{ x\in X:\left\Vert Px\right\Vert _{Y} < c\left\Vert x\right\Vert _{X}\right\}$ is open and non empty by (3.1). This caveat in mind, we nevertheless wish to test if some evidence can be given to support our optimistic approach that in most cases, inversion should be possible.

Proposition 2. For every $\epsilon > 0$ and $x\in X$, there exists $c > 0$ such that

Proof. Take $x\in X$. Define the real-valued map $g(c) = \mathbb{P}\left(\left\Vert Px^{\nu}\right\Vert _{Y}\geq c\left\Vert x\right\Vert _{X}\right).$ Take a sequence $c_{n}\searrow0$. It is enough to show that

We write

Note that $f_{n}$ is monotone increasing and, since $\ker P = \left\{ 0\right\}$, $\lim_{n\to\infty}f_{n}(\nu) = 1$ for every $\nu$. Thus by the Monotone Convergence Theorem,

This thus shows that, for a fixed $x$, our intuition is vindicated: in the vast majority of cases, the inequality $\left\Vert Px^{\nu}\right\Vert _{Y}\geq c\left\Vert x\right\Vert _{X}$ holds true. This is true independently of $C_{\mathrm{rand}}(\mathrm{e})$; we now investigate whether the positivity of $C_{\mathrm{rand}}(\mathrm{e})$ may yield a stronger estimate.

4.1.   Large deviation inequalities

The next step is to estimate the probability that, for a given $x\in X\setminus\{0\}$, the square of the ratio $\frac{\left\Vert P\left(x^{\nu}\right)\right\Vert _{Y}}{\left\Vert x\right\Vert _{X}}$ used in Definition 1 is close to its mean value. The large deviation result we could derive describes the deviation from an upper bound to $C_{\mathrm{rand}}(\mathrm{e})$, namely the constant $K_{\mathrm{rand}}(\mathrm{e})$ defined by

Theorem 1 (large deviation estimate). Assume that $Y = L^{2}(\Sigma, \mu)$, where $(\Sigma, S, \mu)$ is a measure space. Let $(\beta_{k}^{\nu})_{k\in\mathbb{N}}$ be a sequence of independent random variables of Bernoulli type, on a probability space $(A, \mathcal{A}, \mathbb{P})$ with mean 0 and variance 1. Let $x\in X\setminus\{0\}$ and $x^{\nu} = \sum_{k = 1}^{\infty}\beta_{k}^{\nu}\mathrm{e}_{k}x_{k}$. Then, for every $\delta > 0$ we have

The proof of Theorem 1 is postponed to Appendix B. The argument follows the same lines as the one of [11,Theorem 2.1] and the general method introduced in ^[13].

Remark 1 (Application to a wave system). Considering the wave equation (2.1) and adopting the framework of section 2.2 leads to choose $\{\psi_{j}^{\pm}\}_{j\geq1}$ defined by (2.4) as the orthonormal basis $\mathrm{e}$. In that case, $X = H_{0}^{1}(\Omega)\times L^{2}(\Omega)$, $\Sigma = [0, T]\times\partial\Omega$, $d\mu = dt\, d\mathcal{H}^{d-1}$ and $Y = L^{2}(\Sigma)$. Following the discussion in section 2.2, the map $P$ is given by $P(y^{0}, y^{1}) = { \mathbb{1}}_{\Gamma}\frac{\partial y}{\partial\nu}$, where $y$ is the unique solution of (2.1). Further, we have

where the first equality is given in Proposition 1, and the second one follows by applying the same argument.

Note that according to the so-called Rellich identity^*, we have $0 < K_{\text{rand}}\left(\mathrm{e}\right)\leq\frac{T}{2}\operatorname{diam}(\Omega)$, under additional mild assumptions on the domain $\Omega$.

This identity, discovered by Rellich in 1940 ^[43], reads

for every eigenpair $(\lambda, \phi)$ of the Laplacian-Dirichlet operator, $\Omega$ being a bounded connected domain of $\mathbb{R}^{n}$ either convex or with a $C^{1, 1}$ boundary.

The estimate given in Theorem 1 is on the "wrong side", since we show that the ratio related to the inversion is much bigger than $K_{\mathrm{rand}}(\mathrm{e})$ with low probability. The issue is not the boundedness of $P$, which is given a priori, but of its inverse. This would correspond to a result of the type

namely, a quantification of the estimate given in Proposition 2, uniform in $x$. If such a bound held, it would show that $C_{\mathrm{rand}}(\mathrm{e})$ is a reliable estimator of the behaviour of the ratio $\frac{\|P(x)\|_{Y}^{2}}{\|x\|_{X}^{2}}$ in general. Notice that, in the favourable case when $P^{-1}$ is bounded, there exists $\delta_0\in [0, C_{\mathrm{rand}}(\mathrm{e}))$ such that

for all $\delta\in [\delta_0, C_{\mathrm{rand}}(\mathrm{e}))$.

In this general framework, estimate (4.3) does not hold, see Example 2. Using a concentration inequality, a weaker bound can be derived.

Proposition 3. Assume that $Y = L^{2}(\Sigma, \mu)$, where $(\Sigma, S, \mu)$ is a measure space. Let $(\beta_{k}^{\nu})_{k\in\mathbb{N}}$ be a sequence of independent random variables of Bernoulli type, on a probability space $(A, \mathcal{A}, \mathbb{P})$ with mean 0 and variance 1. Let $x\in X\setminus\{0\}$ and $x^{\nu} = \sum_{k = 1}^{\infty}\beta_{k}^{\nu}\mathrm{e}_{k}x_{k}$. Then, for every $\delta > 0$ we have

This result is based on an appropriate Hoeffding inequality; its proof is postponed to Appendix B. Note that (4.4) is not a large deviation result: the quantity under consideration is bounded between 0 and 1, and the upper bound obtained is not small. This is unavoidable, see Example 2.

4.2.   Can you reconstruct two numbers from their sum?

We collect here several observations that suggest that the positivity of the randomised stability constant may not be helpful for solving the inverse problems, not even for most unknowns.

4.2.1.   Instability arises for every $x\in X$

We remind the reader why (3.1) renders inversion unstable. Hypothesis (3.1) implies that there exists a sequence $\left(x_{n}\right)_{n\in\mathbb{N}}$ such that

Suppose that our measurements are not perfect, and are affected by a low level of noise $\delta$, $\left\Vert \delta\right\Vert _{Y}\leq\epsilon$, with $\epsilon > 0$. Then, for every $n$ such that $n\epsilon > 1$, we have

hence $x$ and $x+x_{n}$ correspond to the same measured data, even if $\left\Vert \left(x+x_{n}\right)-x\right\Vert _{X} = 1$. This is an unavoidable consequence of the unboundedness of $P^{-1}$, and is true for every $x\in X$, even if the randomised stability constant were positive (and possibly large).

4.2.2.   The dependence of $C_{\mathrm{rand}}$ on the basis

Lemma 1 shows that the ONB used to randomise our input plays a role, as it appears explicitly in the formula (3.2). The following proposition underscores that point. Namely, if we consider all possible randomisations with respect to all ONB of $X$ we recover the deterministic stability constant $C_{\mathrm{det}}$.

Proposition 4. We have

where the infimum is taken over all ONB of $X$. In particular, if $P^{-1}$ is unbounded then $\inf_{\mathrm{e}} C_{\mathrm{rand}}(\mathrm{e}) = 0$.

Proof. By definition of $C_{\mathrm{rand}}(\mathrm{e})$, we have that $C_{\mathrm{rand}}(\mathrm{e})\ge C_{\mathrm{det}}$ for every ONB $\mathrm{e}$, and so it remains to prove that

By definition of $C_{\mathrm{det}}$, we can find a sequence $x_{n}\in X$ such that $\|x_{n}\|_{X} = 1$ for every $n$ and $\|Px_{n}\|_{Y}^{2}\to C_{\mathrm{det}}.$ For every $n$, complete $x_{n}$ to an ONB of X, which we call $\mathrm{e}^{(n)}$. By Lemma 1 we have $C_{\mathrm{rand}}(\mathrm{e}^{(n)})\le\|Px_{n}\|_{Y}^{2}$, and so

This result shows that, in general, the randomised stability constant strongly depends on the choice of the basis. There will always be bases for which it becomes arbitrarily small when $P^{-1}$ is unbounded.

It is also worth observing that for compact operators, which arise frequently in inverse problems, the randomised stability constant is always zero.

Lemma 2. If $P$ is compact then $C_{\mathrm{rand}}(\mathrm{e}) = 0$ for every ONB $\mathrm{e}$ of $X$.

Proof. Since $\mathrm{e}_{k}$ tends to zero weakly in $X$, by the compactness of $P$ we deduce that $P(\mathrm{e}_{k})$ tends to zero strongly in $Y$. Thus, by Lemma 1 we have

as desired.

4.2.3.   Examples

Let us now consider some examples. The first example is finite dimensional and the kernel of the operator is not trivial. Note that the definition of $C_{\mathrm{rand}}(\mathrm{e})$ and all results above, except Proposition 2 and Lemma 2, are valid also in this case, with suitable changes due to the finiteness of the ONB.

Example 1. Choose $X = \mathbb{R}^{2}$ and $Y = \mathbb{R}$, and consider the map

The associated inverse problem can be phrased: Find the two numbers whose sum is given. This problem is ill-posed and impossible to solve. The deterministic stability constant vanishes

and $S^{-1}$ does not exist. However, the randomised constant obtained using the canonical basis is positive. Indeed, $\left|P\left(1, 0\right)\right| = \left|P\left(0, 1\right)\right| = 1$, therefore

The positivity of this constant does not imply the existence of any useful method to perform the reconstruction of $x$ and $y$ from $x+y$, even for most $(x, y)\in\mathbb{R}^{2}$.

Had we chosen as orthonormal vectors $\frac{1}{\sqrt{2}}\left(1, 1\right)$ and $\frac{1}{\sqrt{2}}\left(1, -1\right)$, since $\left|P\left(1, -1\right)\right| = 0$, we would have found

One may wonder whether the features highlighted above are due to the fact that the kernel is not trivial. That is not the case, as the following infinite-dimensional generalisation with trivial kernel shows.

Example 2. Consider the case when $X = Y = \ell^{2}$, equipped with the canonical euclidean norm. Let $\mathrm{e} = \{\mathrm{e}_{k}\}_{k = 0}^{+\infty}$ denote the canonical ONB of $\ell^{2}$. Take a sequence $\left(\eta_{n}\right)_{n\in\mathbb{N}_{0}}$ such that $\eta_{n} > 0$ for all $n$, and $\lim_{n\to\infty}\eta_{n} = 0$. We consider the operator $P$ defined by

The operator $P$ may be represented with respect to the canonical basis $\mathrm{e}$ by the block-diagonal matrix

In other words, $P$ may be expressed as

We note that $\ker T = \left\{ 0\right\}$ and that its inverse is given by

which is an unbounded operator since $\eta_{n}^{-1}\to+\infty$. Given the block diagonal structure of this map, the inversion consists of solving countably many inverse problems (i.e., linear systems) of the form

As soon as $\tilde{n}$ is such that $\eta_{\tilde{n}}$ becomes smaller than the noise level, all the following inverse problems for $n\ge\tilde{n}$ are impossible to be solved, since they reduce to the "sum of two numbers" discussed in Example 1.

Note that

therefore

If we choose instead the rotated orthonormal basis,

then $P\left(v_{2n+1}\right) = -\eta_{n}\frac{1}{\sqrt{2}}\mathrm{e}_{2n+1},$ and so

We now turn to (4.3) and (4.4). For some $k\geq0$, consider $x = \mathrm{e}_{2k}$. Then $P(x^{\nu}) = {\beta_{2k}^{\nu}}(\mathrm{e}_{2k}+\mathrm{e}_{2k+1})$ and therefore $\|P(x^{\nu})\|_{Y} = \|P(x)\|_{Y}$: there is no deviation as

Thus, the probabilities in (4.3) and (4.4) are equal to $0$, and the inequalities are trivial.

However, alternatively, consider $x = \mathrm{e}_{2k}+\mathrm{e}_{2k+1}$. Then $\|P(x^{\nu})\|_{Y}^{2} = 4+(2+\eta_{2k})^{2}$ if $\beta_{1}^{\nu}\beta_{2}^{\nu} = 1$ and $\|P(x^{\nu})\|_{Y}^{2} = \eta_{2k}^{2}$ if $\beta_{1}^{\nu}\beta_{2}^{\nu} = -1$. Therefore

As a consequence, (4.3) cannot be true in general for every $x$. Similarly, we have

and the left-hand side of (4.4) can indeed be large for some $x$.

It is worth observing that a very similar example is considered in ^[32] to show that particular complications may arise when using neural networks for solving some inverse problems, even naive and small scale (cfr. §A.3).

These examples show that considering the observability constant for a particular basis sheds little light on a potential stable inversion of the problem in average, and that considering all possible randomisations leads to the same conclusion as the deterministic case (confirming Proposition 4).

4.3.   Linear versus nonlinear problems

The pitfalls we encountered when we tried to make use of the randomised stability constant all stem from the linearity of the problems we are considering. The seminal work of Burq and Tzvetkov ^[12], which showed existence of solutions in super-critical regimes for a semilinear problem did not involve tinkering with associated linear operator (the wave equation); it is the nonlinearity that controlled the critical threshold. In both compressed sensing and passive imaging with random noise sources, nonlinearity plays a key role; further, deep networks are nonlinear maps (cfr. Appendix A).

The naive intuition we discussed earlier, namely, that extreme situations do not occur often, is more plausible for nonlinear maps where pathological behaviour is local.

Example 3. As a toy finite-dimensional example, consider the map

for some small $\epsilon > 0$. Then $T$ can be stably inverted outside of a region with size of order $\epsilon$, since there the inverse is continuous. Thus, a random initial data has little chance of falling precisely in the problematic region. Such a case cannot happen with linear maps.

Example 4. Let $A\colon H\to H$ be an unbounded linear operator on a Hilbert space with compact resolvent, so that the spectrum of $A$ is discrete. Define the nonlinear map

Note that $A+\lambda I$ is invertible with probability $1$ if $\lambda$ is chosen uniformly in $[0, 1]$. Thus, if $H\times[0, 1]$ is equipped with a product probability measure whose factor on $[0, 1]$ is the uniform probability, then $x$ may be reconstructed from $T(x)$ with probability $1$.

5.   Concluding remarks

In this paper we focused on the randomised stability constant for linear inverse problems, which we introduced as a generalisation of the randomised observability constant.

We argue that, despite its intuitive and simple definition, the randomised stability constant has no implications in the practical solution of inverse problems, even for most unknowns. As the examples provided show, this may be due to the linearity of the problem. With nonlinear problems, the situation is expected to be completely different. It could be that the randomised stability constant is meaningful in the context of a nonlinear inversion process, involving for example a hierarchical decomposition ^[34,46], but we do not know of results in that direction: this is left for future research.

Conflict of interest

The authors declare no conflict of interest.

A.   Examples of techniques based on randomisation

In this appendix we briefly review three different techniques for solving inverse problems where randomisation plays a crucial role. We do not aim to provide an exhaustive overview, or to report on the most recent advances, or to discuss the many variants that have been studied. The examples we present are used to contrast possible different approaches, and the level of mathematical understanding associated with them.

A.1.   Compressed sensing

Since the seminal works ^[15,16], compressed sensing (CS) has provided a theoretical and numerical framework to overcome Nyqvist criterion in sampling theory for the reconstruction of sparse signals. In other words, sparse signals in $\mathbb{C}^{n}$ may be reconstructed from $k$ discrete Fourier measurements, with $k$ smaller than $n$ and directly proportional to the sparsity of the signal (up to log factors, see eq. (A.2) below). Let us give a quick overview of the main aspects of CS, in order to show how it fits in the general framework of section 1. For additional details, the reader is referred to the book ^[19], and to the references therein.

Given $s\in\mathbb{N} = \{1, 2, \dots\}$, let $X$ be the set of $s$-sparse signals in $\mathbb{C}^{n}$, namely

Let $F\colon\mathbb{C}^{n}\to\mathbb{C}^{n}$ denote the discrete Fourier transform. In fact, any unitary map may be considered, by means of the notion of incoherence ^[14]. In any case, the Fourier transform is a key example for the applications to Magnetic Resonance Imaging and Computerised Tomography (via the Fourier Slice Theorem). In order to subsample the Fourier measurements, we consider subsets $S_{a}$ of cardinality $k$ of $\{1, 2, \dots, n\}$ and parametrise them with $a\in\{1, 2, \dots, {\binom{n}{k}}\}$. Let $Y = \mathbb{C}^{k}$ and $P_{a}\colon\mathbb{C}^{n}\to\mathbb{C}^{k}$ be the projection selecting the entries corresponding to $S_{a}$. We then define the measurement map

In other words, $T_{a}$ is the partial Fourier transform, since only the frequencies in $S_{a}$ are measured, and $\#S_{a} = k\le n$.

Given an unknown signal $x_{0}\in X$, we need to reconstruct it from the partial knowledge of its Fourier measurements represented by $y: = T_{a}(x_{0})$. The sparsity of $x_{0}$ has to play a crucial role in the reconstruction, since as soon as $k < n$ the map $P_{a}\circ F\colon\mathbb{C}^{n}\to\mathbb{C}^{k}$ necessarily has a non-trivial kernel. It is worth observing that sparsity is a nonlinear condition: If $X$ were a linear subspace of $\mathbb{C}^{n}$, the problem would be either trivial or impossible, depending on $\ker(P_{a}\circ F)\cap X$. Thus nonlinearity plays a crucial role here.

The simplest reconstruction algorithm is to look for the sparsest solution to $T_{a}x = y$, namely to solve the minimisation problem

where $\|x\|_{0} = \#{\rm{supp}} x$. However, this problem is NP complex, and its direct resolution impractical. Considering the convex relaxation $\|\cdot\|_{1}$ of $\|\cdot\|_{0}$ leads to a well-defined minimisation problem

whose solution may be easily found by convex optimisation (in fact, by linear programming).

The theory of CS guarantees exact reconstruction. More precisely, if $\tilde{x}$ is a minimiser of (A.1), then $\tilde{x} = x_{0}$ with high probability, provided that

and that $a$ is chosen uniformly at random in $\{1, 2, \dots, {\binom{n}{k}}\}$ (namely, the subset $S_{a}$ is chosen uniformly at random among all the subsets of cardinality $k$ of $\{1, 2, \dots, n\}$) ^[15]. In addition, in the noisy case, with measurements of the form $y = T_{a}(x_{0})+\eta$ where $\|\eta\|_{2}\le\epsilon$, by relaxing the equality "$T_{a}x = y$" to the inequality "$\|T_{a}x-y\|_{2}\le\epsilon$" in (A.1), one obtains the linear convergence rate $\|x_{0}-\tilde{x}\|_{2}\le C\epsilon$, namely, the solution is stable.

In summary, CS allows for the stable reconstruction of all sparse signals from partial Fourier measurements, for most choices of the measured frequencies. The corresponding forward map $T_{a}\colon X\to\mathbb{C}^{k}$ is nonlinear, simply because $X$ is not a vector space.

A.2.   Passive imaging with random noise sources

The material presented in this part is taken from ^[20], to which the reader is referred for more detailed discussion on this topic.

A typical multistatic imaging problem is the recovery of some properties of a medium with velocity of propagation $c(x) > 0$ from some measurements at locations $x_{j}\in\mathbb{R}^{3}$ of the solution $u(t, x)$ of the wave equation

where $f(t)$ is the source pulse located at $y$. One of the major applications of this setup is geophysical imaging, where one wants to recover properties of the structure of the earth from measurements taken on the surface. Generating sources in this context is expensive and disruptive. Earthquakes are often used as sources, but they are rare and isolated events. Yet, noisy signals, as they may be recorded by seismographers, even if low in amplitude, may be relevant and useful even in absence of important events.

The key idea is to consider the data generated by random sources (e.g., in seismology, those related to the waves of the sea). The equation becomes

where the source term $n(t, x)$ is a zero-mean stationary random process that models the ambient noise sources. We assume that its autocorrelation function is

where $F$ is the time correlation function (normalised so that $F(0) = 1$) and $K$ characterises the spatial support of the sources. The presence of $\delta(y_{1}-y_{2})$ makes the process $n$ delta-correlated in space.

The reconstruction is based on the calculation of the empirical cross correlation of the signals recorded at $x_{1}$ and $x_{2}$ up to time $T$:

Its expectation is the statistical cross correlation

which is given by

where $\hat{\cdot}$ denotes the Fourier transform in time and $G(t, x, y)$ is the time-dependent Green's function. Moreover, $C_{T}$ is a self-averaging quantity, namely

in probability.

The role of randomised sources is now clear: from the measured empirical cross correlation $C_{T}$ with large values of $T$ it is possible to estimate, with high probability, the statistical cross correlation $C^{(1)}$. Using (A.3), from $C^{(1)}(\tau, x_{1}, x_{2})$ it is possible to recover (some properties of) the Green function $G$, which yield useful information about the medium, such as travel times.

A.3.   Deep Learning in inverse problems for PDE

Convolutional Neural Networks have recently been used for a variety of imaging and parameter reconstruction problems ^[22], including Electrical Impedance Tomography (EIT) ^[23,33,48], optical tomography ^[18], inverse problems with internal data ^[9], diffusion problems in imaging ^[7], computerised tomography ^[10,28], photoacoustic tomography ^[5,25] and magnetic resonance imaging ^[49,50]. In the following brief discussion, we decided to focus on inverse problems for partial differential equations (PDE), and in particular on EIT, but similar considerations are valid for most methods cited above.

Significant improvement has been observed in EIT with deep learning compared to previous imaging approaches. Let $\Omega\subseteq\mathbb{R}^{d}$, $d\ge2$, be a bounded Lipschitz domain with outer unit normal $\nu$. The data in EIT is (a part of) of the Dirichlet-to-Neumann map

where $u(x)$ denotes the unique solution of the elliptic problem

and $\sigma(x) > 0$ is the unknown conductivity. The experimental data is usually part of the inverse map, namely the Neumann-to-Dirichlet map $\Lambda_{\sigma}^{-1}$. In two dimensions, provided that the electrodes are equally separated on the unit disk, the data may be modelled by

where $T_{N}$ is the $L^{2}$ projection on span$\{\theta\mapsto\cos(n\theta):1\leq n\leq N\}$, where $N$ is the number of electrodes: It is the partial Fourier transform limited to the first $N$ coefficients.

Direct neural network inversion approaches suffer from drawbacks alike direct non-regularised inversion attempts: the output is very sensitive to measurement errors and small variations. Successful approaches to Deep Learning EIT ^[23,33,48], and to other parameter identification problems in PDE, often involve two steps.

The first step consists in the derivation of an approximate conductivity $\sigma$ by a stable, albeit blurry, regularised inversion method. For instance, in ^[23] the "D-bar" equation is used, while in ^[33] a one-step Gauss-Newton method is used. In both cases, the output of this step is a representation of the conductivity coefficient, which depends on the inversion method used. This first step is deterministic and its analysis is well understood. The forward problem, relating the conductivity to the Dirichlet-to-Neumann map, is nonlinear, independently of the inversion algorithm used. Indeed, the map $\Lambda_{\sigma}$ is a linear operator, but $\sigma\mapsto\Lambda_{\sigma}$ is nonlinear.

The second, post-processing, step uses a neural network to "deblur" the image, and in fact restores details that were not identifiable after the first step.

The second step is not unlike other successful usage of deep-learning approaches for image classification; in general they are known to be successful only with very high probability (and in turn for random unknowns). More precisely, since the findings of ^[45], deep networks are known to be vulnerable to so-called "adversarial perturbations" (see the review article ^[1] and the references therein). Given an image $x$ that is correctly classified by the network with high confidence, an adversarial perturbation is a small perturbation $p$ such that the images $x$ and $y = x+p$ are visually indistinguishable but the perturbed image $y$ is misclassified by the network, possibly with high confidence too. State-of-the-art classification networks are successful for the vast majority of natural images, but are very often vulnerable to such perturbations.

These instabilities are not specific to image classification problems; they appear in the same way in image reconstruction ^[6]. In this case, given an image that is well-reconstructed by the network, it is possible to create another image that is visually indistinguishable from the original one, but that is not well-reconstructed by the network.

A full mathematical understanding of deep networks is still lacking, and the reasons of this phenomenon are not fully known. However, the large Lipschitz constant of the network certainly plays a role, since it is a sign of potential instability: in order for the network to be effective, the weights of its linear steps need to be chosen large enough, and the composition of several layers yields an exponentially large constant.

B.   Proofs of the large deviation estimates

The proofs of Theorem 1 and Proposition 3 rest upon a classical large deviation estimate, the so-called Hoeffding inequality, see e.g. [11,Prop. 2.2], whose proof is recalled for the convenience of the reader.

Proposition 5. Let $(\alpha_{n}^{\nu})_{n\geq1}$ be a sequence of sequence of independent random variables of Bernoulli type, on a probability space $(A, \mathcal{A}, \mathbb{P})$ with mean 0 and variance 1. Then, for any $t > 0$ and any sequence $(u_{n})_{n\geq1}\in\ell^{2}(\mathbb{C})$, we have

Proof. There holds $\mathbb{E}\left(\exp\left(\alpha_{n}^{\nu}v_{n}\right)\right) = \mathbb{E}\left(\sum_{k = 0}^{\infty}\frac{\left(\alpha_{n}^{\nu}v_{n}\right)^{k}}{k!}\right) = \sum_{k = 0}^{\infty}\frac{1}{k!}\mathbb{E}\left(\left(\alpha_{n}^{\nu}v_{n}\right)^{k}\right)$. All odd powers of $k$ vanish as $\alpha_{n}^{\nu}$ has zero mean and is symmetrical. Therefore for any $\lambda > 0$

Applying Chernoff's inequality, we obtain for any $\lambda > 0$

Choosing $\lambda = \frac{t}{\sum v_{n}^{2}}$ this yields

as desired.

We are now ready to prove Theorem 1.

Proof of Theorem 1. Fix $r\geq2$ and set $\mathcal{Y}^{\nu} = \frac{P\left(x^{\nu}\right)}{\left\Vert x^{\nu}\right\Vert _{X}}$. Markov's inequality yields

Let us denote by $\mathcal{L}_{\nu}^{r}$ the standard Lebesgue space with respect to the probability measure $d \mathbb{P}$. Recall that $Y = L^{2}(\Sigma, {\mu})$. To provide an estimate of the right-hand side, notice that

by using Jensen's inequality.

Furthermore, for a.e. $s\in\Sigma$, we have

Here, we use that if $X$ denotes a non-negative random variable and $\varphi\colon\mathbb{R}_{+}\to\mathbb{R}_{+}$, then

and by using Proposition 5 and the fact that $\mathcal{Y}^{\nu}(s)$ reads

one gets

As a consequence, by using the change of variable

we get

with

An elementary computation yields $C(r) < r^{r}$. Therefore,

According to (B.2), we infer that

From (4.2) and estimate (B.1), we get

using the triangle inequality. Minimising the upper bound with respect to $r$, that is, choosing

in the inequality above, one finally obtains

Note that we have assumed $r\geq2$, thus implicitly posited that ${e^{-2}\delta^{2}}\geq4{ K_{\mathrm{rand}}(\mathrm{e})}$; we multiply the bound by $\exp(2)$ to cover the other case.

We conclude by proving Proposition 3.

Proof of Proposition 3. As in (3.3), we have

Observing that the family $\alpha_{k, l}^{\nu} = \{\beta_{k}^{\nu}\beta_{l}^{\nu}\}_{\substack{k, l\\ k < l} }$ is made of independent Bernoulli variables, equal almost surely to $-1$ or 1, with $1/2$ as probability of success, we apply Proposition 5 with

and obtain for all $\delta > 0$

with