Learning the optimal Tikhonov regularizer for inverse problems
Year: 2022

Abstract:
In this work, we consider the linear inverse problem y = Ax+ε, whereA : X → Y is a known linear operator between the separable Hilbert spaces X and Y , x is a random variable in X and ε is a zero-mean random process in Y . This setting covers several inverse problems in imaging including denoising, deblurring and Xray tomography. Within the classical framework of regularization, we focus on the case where the regularization functional is not given a priori, but learned from data. Our first result is a characterization of the optimal generalized Tikhonov regularizer, with respect to the mean squared error. We find that it is completely independent of the forward operator A and depends only on the mean and covariance of x. Then, we consider the problem of learning the regularizer from a finite training set in two different frameworks: one supervised, based on samples of both x and y, and one unsupervised, based only on samples of x. In both cases we prove generalization bounds, under some weak assumptions on the distribution of x and ε, including the case of sub-Gaussian variables. Our bounds hold in infinite-dimensional spaces, thereby showing that finer and finer discretizations do not make this learning problem harder. The results are validated through numerical simulations.

heading: 1 Introduction
text: The aim of an inverse problem is to recover information about a physical quantity from indirect measurements. Virtually all imaging problems and modalities fall within this framework, including denoising, deblurring [14], computed tomography [41] and magnetic resonance imaging [15]. Classical and general approaches to solve inverse problems consist in studying a variational (minimization) problem and can be divided into two classes.
The first is based on the so-called regularization theory [14]. The aim is to recover a single, deterministic, unknown x† from noisy data y = F (x†) + ε by solving a minimization problem
min x dY (F (x), y) + J(x) (1)
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
for a fidelity term dY : Y × Y → R and a regularization functional J : X → [0,+∞). The latter is chosen in order to mitigate the ill-posedness of the map F , and represent some a-priori knowledge on x. For instance, in the classical Tikhonov regularization we have J(x) = λ‖x‖2X for λ > 0. The second approach considers the unknown as a random variable and is based on statistical/Bayesian methods [23, 45]. In this case one can recover the unknown using point estimators such as the maximum a posteriori (MAP) estimator, or extract richer information on the probability distribution of the unknown. In practice, the MAP estimator is found by solving a minimization problem of the same form as (1). The main difference is that the fidelity term and the regularizer are tailored to the statistical properties of the unknown and the noise, which are usually assumed to be known.
In recent years, machine learning techniques, and especially deep learning, have shaken the field of inverse problems by providing the basis for data-driven methods that have outperformed the state-of-the-art in most imaging modalities [3, 43]. The most successful methods take inspiration from regularization theory [10, 1, 26, 2, 36, 21, 22, 27, 34, 39, 8, 17, 40]: the physical model given by the forward map F is assumed to be known while the regularizer (or the gradient updates related to it) is learned from a training set. While these approaches have shown impressive results in applications, a solid theory behind their successes is lacking. In view of the many sensitive applications where these methods are already being employed, e.g. in medical imaging, it is of utmost importance to fill this theoretical gap to better understand the strengths and limits of data-driven imaging modalities. Moreover, many inverse problems are naturally formulated in infinite-dimensional spaces [14, 41, 23, 15, 45, 37], and their discretization must be carefully treated due to their ill-posedness [25]. Hence, it is of main interest to provide a theoretical analysis in the infinite-dimensional setting.
In this work, we consider the problem of learning a regularizer for a linear inverse problem in the framework of statistical learning theory, which is the natural setting to derive precise theoretical guarantees. This is part of the growing research area of learning an operator between infinite dimensional spaces [46, 13, 32, 42, 31]. We study the case where the measurements are modeled by a linear, possibly ill-posed, forward map, and the penalty term is a generalized Tikhonov regularizer [47].
More precisely, let A : X → Y be a bounded linear operator between the separable real Hilbert spaces X and Y . We consider the inverse problem
y = Ax+ ε, (2)
which consists of the reconstruction of x from the knowledge of y, where ε represents noise. We assume that x is a random variable onX with mean µ and covariance Σx, and the noise ε, independent of the variable x, is a zero-mean random process on Y with covariance Σε, (see Section 2 for more details). The operator A is typically injective but its inverse may be unbounded: typical examples include denoising (A is the identity) and deblurring (A is a convolution operator).
We aim to recover the unknown via generalized Tikhonov regularization. For a quadratic fidelity term dY : Y × Y → R, the minimization problem
min x dY (Ax, y) + ‖B−1(x− h)‖2X , (3)
has a unique solution Rh,B(y), called the generalized Tikhonov reconstruction. Here the pair (h,B), where h ∈ X and B : X → X is a positive bounded operator, is considered as a free parameter that we call regularization pair. For example, the operator B can be a smoothing operator in L2, as a negative power of the Laplacian, so that B−1 is a differential operator, yielding a classical regularization in Sobolev spaces. We want to characterize and learn the optimal pair (h,B) with respect to the expected, or mean squared, error
L(h,B) = Ex,y‖Rh,B(y)− x‖2X .
Our first contribution is a complete characterization of the minimizers of L. In particular we find that (µ,Σ1/2x ) is a global minimizer (and is unique if A is injective), which shows that the best regularizer is completely independent of the forward operator A and depends only on the mean and covariance of x. This is consistent with the known linearized minimum mean squared error estimator in the finite-dimensional case [24], but it is usually not taken into account in the machine learning approaches to inverse problems mentioned above. The extension to the infinite-dimensional case is not straightforward due to the presence of unbounded operators in inverse problems.
Since the computation of the expected error requires the full distribution ρ of x and y, we study how this can be approximated from a finite training set. We suppose to have access to a sample of m pairs z = {(xj , yj)}mj=1 drawn independently from ρ. In view of the results on the expected error, we consider two alternative ways to learn a regularizer pair (ĥz, B̂z) from a training set z: either by minimizing the empirical risk [12] (supervised learning), or by using the empirical mean and covariance of {xj} (unsupervised learning). In both cases, we prove generalization bounds for the sample error |L(ĥz, B̂z)− L(µ,Σ1/2x )|. Under some natural compactness assumptions on the class of regularization pairs, we prove that the sample error has the asymptotic behavior
|L(ĥz, B̂z)− L(µ,Σ1/2x )| . 1√ m , (4)
with high probability, in both the supervised and the unsupervised approaches. We stress the point that these bounds hold in the infinite-dimensional setting, or, in other words, they do not depend on the discretization of the signal and of the measurements.
Finally, we complement our theoretical findings with some numerical experiments. For a 1D denoising problem (i.e. A is the identity operator) we replicate the asymptotic bound (4) at different discretization scales. Moreover, we find that the unsupervised approach, despite yielding the same rate (4), clearly outperforms the supervised one.
The paper is organized as follows. In Section 2 we introduce the main notation and technical assumptions that will be used throughout the paper, including several examples. Section 3 presents the main results for the minimization of the expected error, while Section 4 is devoted to the study of the sample error. Numerical experiments are the subject of Section 5. Concluding remarks and discussions are reserved for Section 6.

heading: 2 Setting the stage

heading: 2.1 The random objects x and ε
text: As mentioned in the introduction, we formulate (2) as a statistical inverse problem, where x and ε are not deterministic but random. Let us start with the description of the prior on x. Assumption 2.1. Let x be a random variable on a probability space (Ω,P) taking values in X . More precisely, x is square-integrable, so that its expectation µ ∈ X and its covariance Σx : X → X is a trace-class operator. We assume that Σx is injective.
Without loss of generality we can always assume that Σx is injective (i.e., x is non-degenerate), since otherwise it would be enough to consider the inverse problem only in (ker Σx)⊥ ( X . Let us consider some common examples of priors arising in inverse problems. Example 2.2 (Gaussian random variables). A general class of priors arises when considering Gaussian random variables. We recall that x is a Gaussian random variable if for all v ∈ X , 〈x, v〉 is a real Gaussian random variable and, by Fernique’s theorem, x is square-integrable [7]. Since Σx : X → X is self-adjoint, positive and trace-class, we can write its singular value decomposition (SVD) as
Σxv = ∑ k σ2k〈v, ek〉Xek, v ∈ X,
where {ek}k is an orthonormal basis of X , ∑ k σ 2 k < +∞ and 〈x, ek〉 ∼ N (µk, σk), where
µ = ∑ k µkek is the mean of x. In other words,
x = µ+ ∑ k σkakek, (5)
where ak are i.i.d. standard Gaussian variables. This shows that, in infinite dimension, since σk → 0, the variations of x along the direction ek become smaller and smaller as k → +∞.
This abstract construction reduces to a smoothness prior by suitably choosing the covariance operator. Example 2.3 (Smoothing priors). Let X = L2(Td), where Td = Rd/Zd is the d-dimensional torus and d ≥ 1. Let ∆ denote the Laplace-Beltrami operator on Td, which is simply the classical Laplace
operator on [0, 1]d with periodic boundary conditions. For s > d2 , the operator (I −∆)−s : L2(Td)→ L2(Td) is trace class, and can be used to define the Gaussian distribution N (0, (I −∆)−s). In the notation of Example 2.2, the SVD of (I −∆)−s is given by σ2k = (1 + 4π|k|2)−s and ek(t) = e2πik·t with k ∈ Zd. This enforces a smoothness prior on x, depending on the parameter s, which controls the decay of the Fourier coefficients of x (see [33, Appendix B] and [7]).
Let us now discuss the model for the noise ε. Assumption 2.4. Let ε = (ε)v∈Y be a (linear) random process on Y with zero mean and such that its covariance Σε : Y → Y defined by E[εvεw] = 〈Σεv, w〉Y is bounded and injective.
Notice that the injectivity of Σε implies that noise is present in all directions of Y , whereas the boundedness of Σε allows us to regard the random process ε as a bounded (linear) operator from Y into L2(Ω,P). The reader is referred to [16] and to Appendix A.1 for additional details on random processes in Hilbert spaces. We mention here some basic properties that will be needed for the following discussion. Remark 2.5. Even if ε may not belong to Y almost surely (see Example 2.7 below), it is always possible to view it as an element of a larger space, as we now discuss. Let K be a separable Hilbert space and ι : K → Y be an injective linear map such that ι(K) is dense in Y and
ι∗ ◦ Σε ◦ ι : K → K∗ is trace-class,1 (6) where we identify Y ∗ = Y , but we do not identify K∗ with K, and we regard ι∗ as the canonical embedding Y → K∗. The restriction of ε to K is a Hilbert-Schmidt operator from K into L2(Ω,P), hence there exists a unique square-integrable random vector ε taking values in K∗ such that εv = 〈ε, v〉K∗×K for v ∈ K. It is easy to show that the random vector ε has zero mean and its covariance operator is ι∗ ◦ Σε ◦ ι : K → K∗, since
E[〈ε, v〉K∗×K〈ε, w〉K∗×K ] = E[ει(v)ει(w)] = 〈Σειv, ιw〉Y , v, w ∈ K. (7)
A random variable is always a random process, as we now describe. Example 2.6. A simple example of this abstract construction consists of considering a random variable ε. In this case Σε is trace-class itself, so that we can choose K = Y and ι = I and ε ∈ Y almost surely. As discussed in Example 2.2 for Gaussian variables, this means that, since σk → 0, the expected amplitude of the noise in the direction ek goes to 0 as k →∞. For instance, the choice of Σε = (I −∆)−s as in Example 2.3 corresponds to smaller noise levels for higher Fourier modes.
A random process allows for considering noise that is uniformly distributed in all directions. Example 2.7 (White noise). The Gaussian white noise ε is a random process on Y such that for any v ∈ Y it holds that εv is a standard Gaussian variable (mean 0 and variance 1), so that Σε = I . Heuristically, in the notation of Example 2.2, this corresponds to σk = 1 for every k, and so by (5)
ε = ∑ k akek, ak ∼ N (0, 1),
so that ε /∈ Y with probability 1 whenever Y is infinite dimensional (see, e.g., [16]). In view of Remark 2.5, it is possible to consider a larger space K∗ so that ε ∈ K∗ almost surely. For concreteness of explanation, we focus on the case when Y = L2(Td), a typical framework in imaging. A possible choice for the space K is the Sobolev space Hs(Td) with s > d/2 (see [25]), so that the canonical embedding ι : Hs(Td)→ L2(Td) is a Hilbert-Schmidt operator, hence (6) is satisfied and ε can naturally be seen as an element of H−s(Td) = Hs(Td)∗.

heading: 2.2 The new formulation of the inverse problem and of the regularization
text: As a consequence of Assumptions 2.4, since ε may not belong to Y , the inverse problem (2) must be interpreted from a different perspective, namely considering y as the stochastic process yv = 〈Ax, v〉Y + εv on Y or by formulating the problem as an equation in K∗:
y = ι∗Ax+ ε,
1It is worth observing that K and ι always exist: it is enough to choose them, independently of ε, so that the embedding ι : K → Y is Hilbert-Schmidt, which implies that ι∗ ◦ Σε ◦ ι is trace class, since Σε is bounded.
i.e. 〈y, v〉K∗×K = 〈Ax, ι(v)〉Y + 〈ε, v〉K∗×K for v ∈ K, where ι∗ : Y → K∗ is the natural embedding. We denote the joint probability distribution of (x, y) on X ×K∗ by ρ. We now provide a consistent formulation of the quadratic functional appearing in (3). The goal is to replicate what would be the natural choice in a finite dimensional context, i.e.
min x ‖Σ−1/2ε (Ax− y)‖2Y + ‖B−1(x− h)‖2X . (8)
Unfortunately, if Y is infinite dimensional, the first factor is in general not well-defined, since for example for Gaussian processes ε ∈ Im Σ1/2ε with probability 0 [7]. Thus, we need to write this minimization problem in a different formulation. We start by stating the assumptions we make on B. Assumption 2.8. Let us assume that B : X → X is a bounded positive operator such that
Im(AB) ⊆ Im(Σει). (9)
It is worth observing that, whenever Y is infinite dimensional, then Im(Σει) ( Y since Σει is compact. Furthermore, (9) requires, in some sense, the operator AB : X → Y to be at least as “smoothing” as the operator Σει : K → Y . For instance, in the case when AB and Σει have the same left-singular vectors, this condition means that the singular values of AB should go to 0 at least as fast as the singular values of Σει.
We are now ready to rewrite the functional in (8). The penalty term involving B−1(x− h) suggests the change of variables x = h+Bx′. The corresponding minimization problem for x′ reads
min x′∈X
‖Σ−1/2ε (A(h+Bx′)− y)‖2Y + ‖x′‖2X . (10)
This expression does not require the injectivity of B. By formally expanding the first factor we obtain
‖Σ−1/2ε ABx′‖2Y − 2〈Σ−1/2ε (y −Ah),Σ−1/2ε ABx′〉Y + ‖Σ−1/2ε (y −Ah)‖2Y .
Let us analyze these three terms separately:
1. Since Σε is self-adjoint we have ‖Σ−1/2ε ABx′‖2Y = 〈Σ−1ε ABx′, AB′x〉Y , which is well-defined because Im(AB) ⊆ Im(Σε) thanks to (9). 2. The second factor is formally equivalent to −2〈y − Ah,Σ−1ε ABx′〉Y , which is not welldefined as scalar product in Y since y may not belong to Y . However, since y ∈ K∗ and Σ−1ε ABx
′ ∈ ι(K) by (9), this scalar product can be interpreted as the duality pairing −2〈y − ι∗Ah, ι−1Σ−1ε ABx′〉K∗×K .
3. The third factor ‖Σ−1/2ε (y−Ah)‖2Y is independent of x, and so it irrelevant for the minimization task: thus, we remove it. This is a key step in infinite dimension, since, as mentioned above, ‖Σ−1/2ε y‖2Y = +∞ almost surely. See [45, Remark 3.8] for additional details on this aspect.
This discussion motivates the introduction of the following functional, formally equivalent to (10). For y ∈ K∗, we define the regularized solution of the inverse problem as x̂ = h+Bx̂′, where
x̂′ = arg min x′∈X ‖Σ−1/2ε ABx′‖2Y − 2〈y − ι∗Ah, (Σει)−1ABx′〉K∗×K + ‖x′‖2X . (11)
The minimizer exists and is unique, and gives the following expression for x̂:
Rh,B(y) := h+Bx̂ ′ = h+B(BA∗Σ−1ε AB + I) −1((Σει) −1AB)∗(y − ι∗Ah), (12)
where Rh,B : K∗ → X is a bounded affine map. See Proposition A.2 for all the details. Note that ((Σει) −1AB)∗ : K∗ → X is well-defined thanks to (9).

heading: 3 The optimal regularizer
text: The regularization approach described so far is based on the choice of h and B. In classical regularization theory, these are chosen depending on the prior knowledge of the problem under consideration. In the data-driven approach we consider in this work, h and B are learned from training data. In this section, we let learning come into play and consider the problem of determining
the optimal h and B, under the assumptions that the distributions of x and ε are fully known. More precisely, this allows for the explicit computation of the expected error
L(h,B) = E(x,y)∼ρ‖Rh,B(y)− x‖2X = Ex,ε‖Rh,B(ι∗Ax+ ε)− x‖2X , which quantifies the mean square error that our regularization functional (11) yields. Optimal choices of h? and B? are those that minimize L(h,B), and are characterized in the following result. Theorem 3.1. Let X and Y and be separable real Hilbert spaces, A : X → Y be a bounded linear operator, x and ε satisfy Assumptions 2.1 and 2.4 and be independent, and K and ι be as in Remark 2.5. Suppose B = Σ1/2x satisfies Assumption 2.8.
Consider the minimization problem
min h,B {E(x,y)∼ρ
[ ‖Rh,B(y)− x‖2X ] }, (13)
where the minimum is taken over all B satisfying Assumption 2.8 and over all h ∈ X . Then (B?, h?) is a global minimizer of (13) if and only if
h? = µ and B2|(kerA)⊥ = Σx|(kerA)⊥ .
In particular, B? = Σ1/2x is always a global minimizer, and is unique if A is injective. Furthermore, for every minimizer (h?, B?), the corresponding reconstruction map is independent of B? and, for all y ∈ K∗, is given by
R?(y) = µ+ Σ1/2x (Σ 1/2 x A ∗Σ−1ε AΣ 1/2 x + IX) −1((Σει) −1AΣ1/2x ) ∗(y − ι∗Aµ) (14) = µ+ ΣxA ∗(ι∗(AΣxA ∗ + Σε)) −1(y − ι∗Aµ). (15)
The proof is in Appendix A.3. Some comments on this result are in order.
• By assumption ι∗(AΣxA∗ + Σε) is an injective compact operator from Y to K∗, so that its inverse is not bounded, however it is possible to prove that ΣxA∗(ι∗(AΣxA∗ + Σε))−1 extends to a bounded operator from K∗ into X . With a slight abuse of notation, we denote this extension in the same way, so that (15) makes sense for all y ∈ K∗.
• To prove this result, we first consider the minimization in (13) over all possibile affine maps, which yields the so-called Linearized Minimum Mean Square Error (LMMSE) estimator of x. Then, it is possible to show that such optimal affine functional is of the form Rh?,B? , for suitable B? and h?. In a finite-dimensional context, such a result is a direct consequence of the expression of the LMMSE estimator (see, e.g., [24, Theorem 12.1]). Theorem 3.1 generalizes this result to the infinite-dimensional case.
• In the case of Gaussian random variables, the expression of the optimal regularizer R? coincides with the maximum a posteriori (MAP) estimator. Nevertheless, our result does not require any assumptions on x and ε being Gaussian (see the discussion in [19, 20]).
• The minimum expected loss can be computed as L(h?, B?) = tr ( Σ1/2x (Σ 1/2 x A ∗Σ−1ε AΣ 1/2 x + IX) −1Σ1/2x ) , (16)
as it is reported in Appendix A.3.
• It is worth observing that the optimal regularization parameters B? = Σ1/2x and h? = µ are independent of A and ε, and depend only on the mean and the covariance of x.

heading: 4 Finding the optimal regularizer: the sample error
text: The computation of the optimal regularizer proposed in the previous section through the minimization of the expected loss L requires the knowledge of the joint probability distribution ρ of x and y. In this section, we suppose that ρ is unknown, 2 but we have access to a sample z = {(xj , yj)}mj=1 of m pairs (xj , yj) ∈ Z = X ×K∗ drawn independently from the joint probability distribution ρ, and we study how to learn an estimator (ĥz, B̂z) of the optimal parameters (h?, B?). We propose two alternative ways to learn an estimator based on a training sample z. For the ease of notation, from now on we omit the dependence on z.
2More precisely, we only assume that Σε is known.
1. Supervised learning: (ĥS , B̂S) is determined by minimizing the empirical risk L̂, namely
(ĥS , B̂S) = argmin (h,B)∈Θ
L̂(h,B), L̂(h,B) = 1
m m∑ j=1 ‖Rh,B(yj)− xj‖2X , (17)
where Θ is a suitable subset of X × L(X,X).
2. Unsupervised learning: since the best parameters are h? = µ and B? = Σ1/2x , a natural estimator is provided by means of the sample {xj} alone as follows:
ĥU = µ̂ = 1
m m∑ j=1 xj , B̂U = Σ̂x 1/2 , Σ̂x = 1 m m∑ j=1 (xj − µ̂)⊗ (xj − µ̂). (18)
In both cases, we evaluate the quality of (ĥ, B̂) in terms of its excess error L(ĥ, B̂)− L(h?, B?).

heading: 4.1 Supervised learning: empirical risk minimization
text: There exist several techniques to show the convergence of the empirical risk minimizer to the optimal parameter, involving tools such as the VC dimension and the Rademacher complexity (see, e.g., [44]), which require some compactness assumption on Θ. Here, we fix a Hilbert space H with a compact embedding j : H → X having dense range. For %1 > 0, set
Θ1 = {j(h̄) : h̄ ∈ H, ‖h̄‖H ≤ %1}, Θ2 = {jB̄j∗ : B̄ ∈ HS(H∗, H), ‖B̄‖HS(H∗,H) ≤ %1}, (19)
and define Θ as the set of pairs {(h,B) ∈ Θ1 ×Θ2 : B is positive}. Here, HS(H∗, H) denotes the space of Hilbert-Schmidt operators from H∗ to H . We further assume that
a) the map j can be decomposed as j = j2 ◦ j1, where j1 : H → X and j2 : X → X are compact and satisfy
sk(j1) . k −s, s > 0, being sk(j1) the singular values of j1; (20)
whereas j2 is such that Im(Aj2) ⊆ Im(Σει). (21)
b) The optimal parameter (h?, B?) = (µ,Σ1/2x ) belongs to Θ.
Assumption a), and in particular (20), allows us to explicitly compute the covering numbers, whereas (21) ensures that Assumption 2.8 holds uniformly for each positive operator B ∈ Θ2. For example, when H = Hσ1(T1) and X = Hσ2(T1) are Sobolev spaces on the one-dimensional torus, assumption a) is satisfied if s = σ1 − σ2 > 0. As a consequence, assumption b) can be interpreted as an a priori regularity assumption on the problem. Such hypothesis can be relaxed by introducing the approximation error, namely, the rate at which the space Θ approximates X ×L(X,X) as the radius %1 grows to∞. Such an analysis, which easily follows from the range-density property of j, is not treated here.
Finally, we assume that both the inputs and the outputs are bounded, i.e.
supp(ρ) ⊂ BZ(%2), a ball of Z = X ×K∗ of radius %2. (22)
Theorem 4.1. Under the above conditions, let (ĥS , B̂S) be defined by (17) and take τ > 0. We have |L(ĥS , B̂S)− L(h?, B?)| ≤ ( c1 + c2 √ τ√
m
)1− 1 2s′+1
, m ≥ m0, (23)
with probability exceeding 1− e−τ , being 0 < s′ < s, where c1, c2,m0 are independent of m and τ .
The proof of Theorem 4.1 is reported in the Appendices A.4 and A.5. The approach is inspired by [12, Proposition 4] and is suited for a much broader class of learning problems: by adapting Lemma A.6, it is possible to extend the current approach to non-quadratic regularization functionals.
A prominent example of H satisfying (20) comes from Sobolev spaces. Consider, e.g., X = L2(Td), where Td is the d-dimensional torus, and H = Hσ(Td). If σ > 0, the embedding of H in X is compact, and its singular values show a polynomial decay (20) with s = σ/d.

heading: 4.2 Unsupervised learning: empirical mean and covariance
text: As pointed out in (18), it is possible to recover an approximation of the optimal parameter (h?, B?) only by taking advantage of a sample of the output variable {xj}mj=1. Since this technique does not require matched couples of inputs and outputs, we refer to it as an unsupervised learning approach. In order to prove a bound in probability for the sample error, in this section we assume that x is a κ-sub-Gaussian random vector, i.e.,
∀v ∈ X, 〈x, v〉X is a real sub-Gaussian r.v., i.e. ‖〈x, v〉X‖p ≤ κ √ p‖〈x, v〉X‖2, ∀p > 1, (24)
where ‖〈x, v〉‖pp = E[|〈x, v〉|p]. It its known [48] that E[‖x‖p] is finite for all p > 0, so that x has finite mean and its covariance operator Σx is trace-class. Gaussian random variables are a particular instance of sub-Gaussian random variables by Fernique’s theorem [7]. Note that, in infinite-dimensional spaces, bounded random vectors in general are not sub-Gaussian.
We further assume that the injective operator Σε has a bounded inverse, thus Σε + AΣxA∗ is invertible. This is satisfied for example if ε is the white-noise, since Σε = I. We also require that A∗(ι∗(Σε +AΣxA ∗))−1, defined on ι∗(Y ) ⊂ K∗, extends to a bounded operator from K∗ into X .
Theorem 4.2. Under the above conditions, let (ĥU , B̂U ) be defined by (18) and take τ > 0. Then,
|L(ĥU , B̂U )− L(h?, B?)| ≤ c3 + c4 √ τ√
m , m ≥ m0, (25)
with probability exceeding 1− e−τ , where m0, c3 and c4 depend only on Σx, Σε and A.
The proof of Theorem 4.2 is based on several concentration estimates reported in Appendix A.6. The rates we obtain can be meaningfully compared with recent results in supervised learning: see [6, 35].

heading: 5 Numerical simulations
text: We report some numerical results obtained from the supervised and unsupervised strategies for a denoising problem, using synthetic data. The goal of these experiments is twofold: on one hand, we want to study the asymptotic properties of the regularizers learned with the techniques proposed in Section 4 as the sample size m grows, verifying Theorems 4.1 and 4.2. On the other hand, we want to assess that those properties, obtained in an infinite-dimensional setting, do not suffer from the curse of dimensionality. We do so by introducing finer and finer discretizations, and showing that the theoretical bounds do not degrade as the dimension of the problem increases.

heading: 5.1 Problem formulation and discretization
text: We consider a denoising problem on X = Y = L2(T1), being T1 = R/Z the one-dimensional torus, which consists in determining a signal x from the noisy measurement y = x+ ε and thus corresponds to the case A = I . We define a statistical model both for ε and for x, which we use for the generation of the training data. In the learning process, though, we do not take advantage of the knowledge of the introduced probability distributions, apart from the covariance operator of the noise Σε. In accordance with Assumption 2.4, we assume that ε is a random process on Y , and in particular we consider a white noise process, i.e. with zero mean and Σε = σ2I . We consider a noise level of 5%, namely, the standard deviation σ is set to the 5% of the peak value of the average signal. In different tests, we employ different white noise processes with different distributions, including the Gaussian (cfr. Example 2.7) and the uniform distributions. Regarding x, we assume a Gaussian distribution with fixed mean µ and covariance Σx, where µ = 1− |2x− 1| and Σ1/2x is a convolution operator. In order to discretize the described problem, we fix N > 0 and approximate the space X by means of the N−dimensional space generated by a 1D-pixel basis. As a consequence, the functions in X and Y are approximated by vectors in RN , and the linear operators by matrices in RN×N . More details on the discretizations and on the random process generation are reported in Appendix A.7.

heading: 5.2 Implementation and results
text: We denote θ = (h,B). The workflow of the numerical experiments is described as follows:
1. Fix the discretization size N , define the optimal regularizer Rθ? . Compute L(θ?).
2. For each selected value of the sample size m,
• generate the samples {xj}mj=1, {εj}mj=1; • minimize the empirical risk L̂(θ) to find θ̂S ; • compute the empirical mean and covariance µ̂, Σ̂x to find θ̂U ; • compute the excess risks |L(θ̂S)− L(θ?)| and |L(θ̂U )− L(θ?)|.
3. Show the decay of both computed quantities as m increases.
We compute the mean squared errors L(θ?), L(θ̂S) and L(θ̂U ) according to the definition of L, thus avoiding the use of a test set. Moreover, we perform the minimization of the empirical risk analytically, thanks to the explicit expression of the regularization functional provided by (12) (see Appendix A.7). As a final remark, the generalization bounds in Theorems 4.1 and 4.2 can be reformulated in expectation. Thus, to verify the expected decay, we repeat the same experiment 30 times, with different training samples for each size m and taking the average in each repetition.3
In Figure 1, we present the outcome of the numerical experiments, conducted under different statistical models for x and ε. The sample size ranges between 3 · 103 and 3 · 105. In all the presented scenarios, the decay of the excess risk both in the supervised and in the unsupervised cases agrees with the theoretical estimates, showing a decay of the order 1/ √ m. Finally, in Figure 2, we show that the theoretical results are equivalently matched by numerics when the discretization size is increased.
Additional details regarding the results of the numerical experiments are reported in Appendix A.7. Moreover, in Appendix A.8 we replicate the presented numerical study for a different example,
3All computations were implemented with Matlab R2019a, running on a laptop with 16GB of RAM and 2.2 GHz Intel Core i7 CPU. All the codes are available at https://github.com/LearnTikhonov/Code
namely, a deconvolution problem for 1D signals. In this case, A is a convolution operator with respect to a continuous kernel, whose inverse is in general unbounded.

heading: 6 Conclusions and limitations
text: We studied the problem of learning a regularization functional for an inverse problem between infinite dimensional spaces. This problem has received huge interest in recent years due to the successes in several imaging modalities. Our work provides theoretical support to machine learning approaches in sensitive applications such as medical imaging.
We have considered the case of a linear inverse problem that is solved via generalized Tikhonov regularization. This involves an unknown operator B and a signal h, both to be learned from data. We proposed two learning strategies, one supervised and one unsupervised. Surprisingly, we found that the regularizer learned with the unsupervised strategy has the same (or slightly better) generalization bounds than the supervised one. Furthermore, the unsupervised approach does not need the knowledge of the forward operator A nor that of the distribution of the noise ε. This motivates the development of more advanced unsupervised approaches to the problem, e.g. with deep learning methods (see [26, 27, 34, 36, 39]).
The analysis presented here was possible thanks to the simple form of the regularizer. Our work does not cover, for instance, the case of sparsity promoting regularization functionals [18] or more general convex or non-convex penalty terms arising from deep learning methods. Some results regarding optimal (non-quadratic) regularizers associated with different priors can be found e.g. in [19, 20, 9], which nevertheless deal with a finite-dimensional setting. Extensions to more general regularizers will be the subject of future studies.
Acknowledgments and Disclosure of Funding
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA8655-20-1-7027. GSA, EDV, LR and MS are members of GNAMPA, INdAM. GSA is supported by a UniGe starting grant “curiosity driven”. ML is supported by Academy of Finland, grants 273979 and 284715.

Review summaries:
This paper studies the learning of the optimal generalized Tikhonov regularizer for linear inverse problems. Optimality is investigated with regards to the mean-squared error.  The authors show that the optimal parameters are independent of the forward operator and additive noise. Next they look into how the optimal parameters can be learned from finite data for the supervised and unsupervised setting. They bound the excess risk for both settings. Both these bounds turn out to be similar. Lastly, the authors verify the theoretical bounds through numerical simulations. 

This paper considers solving the inverse problem with a quadratic regularizer. It shows that the optimal parameters for the regularizer, in terms of achieving the smallest reconstruction error, are given by the mean and covariance of the true signal distribution. The paper also provides upper bounds for excess risk (expected error) for the estimated parameters learned from a training data set in either a supervised or unsupervised manner.

The paper considers solving linear inverse problems of the form $y=Ax+\epsilon$, where $A:X\to X$ denotes a bounded linear operator. The setup is a Bayesian setup with known distribution of source and noise. To solve this problem the paper focuses on using generalized Tikhonov regularization and solving the following optimization: $\min_x  d_Y(Ax,y)+\|\|B^{-1}(x-h)\|\|^2$. Let  $R_{h,B}(x)$  denote the unique solution of the mentioned optimization and define $L(h,B)=E_{x,y}\|\|R_{h,B}(x)-x\|\|^2$. For a given problem, the question is what is the optimal choice of $B$ and $h$ in terms of minimizing $L(h,B)$. The main result of the paper is to provide a full characterization the minimizers of $L(h,B)$ and show that $(\mu,\Sigma_x^{1/2})$ is a global minimizer of $L(h,B)$. (Here, $\mu$ and $\Sigma_x$ denote the mean and covariance of the input, respectively.) The paper also studies the case where instead of the actual source and noise distributions we have access to training data. In that case it studies two different methods for learning $(h,B)$ from data. 

The paper studies the linear MMSE in infinite-dimensional systems. Specifically, given an observation $Ax +\varepsilon$ where $A$ is a known linear operator and $\varepsilon$ is some form of noise, one wants to estimate $x$, which minimizes the mean squared error.  Focusing on regularization, the paper asks what is the optimal regularizer with respect to the mean squared error. It is shown that, when constrained to Tikhonov regularization, the optimal regularizer has a rather simple form which depends only on the statistical properties of the system and, in particular, is independent from the operator $A$.   The authors also address the problem of approximating the regularizer, based on a finite sample. They show two approximations schemes: A "supervised-learning" scheme, based on minimizing the empirical risk and an "unsupervised-learning" scheme, which is based on the particular form the optimal regularizer was found to have.  For both schemes, the authors show that as the sample size $m$ tend to infinity, the approximation error decays at the expected rate $\frac{1}{\sqrt{m}}$.