# Statistical inverse learning problems with random observations

Abhishake, Tapio Helin, Nicole Mücke 111corresponding author nicole.muecke@tu-
braunschweig.de

###### Abstract

We provide an overview of recent progress in statistical inverse problems with
random experimental design, covering both linear and nonlinear inverse
problems. Different regularization schemes have been studied to produce robust
and stable solutions. We discuss recent results in spectral regularization
methods and regularization by projection, exploring both approaches within the
context of Hilbert scales and presenting new insights particularly in
regularization by projection. Additionally, we overview recent advancements in
regularization using convex penalties. Convergence rates are analyzed in terms
of the sample size in a probabilistic sense, yielding minimax rates in both
expectation and probability. To achieve these results, the structure of
reproducing kernel Hilbert spaces is leveraged to establish minimax rates in
the statistical learning setting. We detail the assumptions underpinning these
key elements of our proofs. Finally, we demonstrate the application of these
concepts to nonlinear inverse problems in pharmacokinetic/pharmacodynamic
(PK/PD) models, where the task is to predict changes in drug concentrations in
patients.

###### Contents

  1. 1 Introduction
    1. 1.1 Literature overview
  2. 2 Statistical Inverse Learning Theory
    1. 2.1 Problem setting
    2. 2.2 Hilbert Scales
    3. 2.3 A Priori Assumptions
    4. 2.4 Minimax optimality
  3. 3 Linear problems
    1. 3.1 Discretization by Random Sampling
    2. 3.2 Optimality of spectral regularization methods
      1. 3.2.1 Spectral regularization methods
      2. 3.2.2 Classical Error Bounds
      3. 3.2.3 Error Bounds in Hilbert Scales
    3. 3.3 Regularization by projection
    4. 3.4 Convex regularization
    5. 3.5 Parameter choice rules and adaptivity
  4. 4 Non-linear problems
    1. 4.1 Classical Error Bounds
    2. 4.2 Error Bounds in Hilbert Scales
  5. 5 Applications
  6. 6 Conclusions and outlook for future directions

##  1 Introduction

Statistical inverse learning refers to the process of inferring an unknown
function $f^{\dagger}$ by using statistical methods on point evaluations of a
related function $g^{\dagger}$. These evaluations might be sparse or noisy,
and the connection between $f^{\dagger}$ and $g^{\dagger}$ is established
through an ill-posed mathematical model

$A(f^{\dagger})=g^{\dagger},\qquad\text{for}\quad
f^{\dagger}\in\mathcal{H}\quad\text{and}\quad
g^{\dagger}\in\mathcal{H}^{\prime},$ (1.1)

where $A$ can be a linear or nonlinear operator between a real separable
Hilbert spaces $\mathcal{H}$ and a reproducing kernel Hilber space
$\mathcal{H}^{\prime}$. In particular, $\mathcal{H}^{\prime}$ represents the
set of functions $g:X\to\mathbb{R}$, where $X$ is a standard Borel space (the
input space) and $\mathbb{R}$ is a space of real numbers (the output space).
The ill-posedness of (1.1) means the problem lacks a unique solution or is
sensitive to data perturbations such as measurement noise.

In this work, we focus on random design regression in inverse problems, where
the data
$\mathbf{z}=\mathopen{}\mathclose{{}\left((x_{i},y_{i})}\right)_{i=1}^{m}$
emerges due to random observations

$y_{i}=g^{\dagger}(x_{i})+\varepsilon_{i},\quad i=1,\dots,m,$ (1.2)

where $\varepsilon_{i}$ is the observational noise, and $m$ denotes the sample
size. Here, we assume that the random observations $(x_{i},y_{i})$ are i.i.d.
and follow the _unknown_ probability distribution $\rho$, defined on the
sample space $Z=X\times\mathbb{R}$. The key objective is to understand
interplay between uncertainty of the design measure and optimal reconstruction
rates.

Classical applications in inverse problems community often involve laboratory
settings, where deterministic design of observational points
$\mathopen{}\mathclose{{}\left(x_{i}}\right)_{i=1}^{m}$ is possible. However,
there is a growing literature considering inverse problems where observational
design contains inherent uncertainties related to areas such as machine
learning [13, 42, 67] and imaging [61, 62, 15]. Statistical inverse learning
provides a principled approach to analyse concentration in such context.

The starting point in the inverse learning paradigm with random design is that
we do not necessarily possess exact information regarding the distribution
$\rho$. In this regard, our work differs from other avenues of random design
regression in inverse problems such as [37, 62, 58], where a specific
(uniform) random design is assumed. Statistical inverse learning aims to
assess how assumptions on $\rho$ influence concentration of given statistical
estimators. That being said, it is known that under suitable regularity
assumptions on $\rho$, the observation scheme in (1.2) is asymptotically
equivalent (in a Le Cam sense of statistical experiments) to the full
observation of $g^{\dagger}$ corrupted by Gaussian noise with vanishing
variance [69].

The ill-posedness of inverse problems is usually mitigated by regularization
techniques [40] which stabilize the reconstruction and prevent overfitting.
These techniques introduce additional information or constraints to the
problem to guide the solution towards a more reasonable and stable outcome. In
this work, we overview recent concentration results for penalized least
squares estimators

$f_{\mathbf{z},\lambda}=\mathop{\rm
arg\,min}\limits_{f\in\mathcal{D}(A)\subset\mathcal{H}}J_{\lambda}(f),\quad
J_{\lambda}(f)=\frac{1}{m}\sum_{i=1}^{m}\mathopen{}\mathclose{{}\left(A(f)(x_{i})-y_{i}}\right)^{2}+\lambda
G(f),$ (1.3)

where the first term in $J_{\lambda}$ is called the empirical error, the
second term is the penalty term $G:\mathcal{H}\to\mathbb{R}_{+}$ and
$\lambda>0$ is the regularization parameter that balances both terms. From the
recent literature studying statistical inverse learning with random design, we
overview three different regularization strategies: classical Tikhonov-type
penalties, projection methods and general convex penalties. We separately
consider results appearing for linear and non-linear forward mapping $A$.

Concentration results in inverse learning require assumptions limiting the
complexity of the set of probability models that may produce the data
$\mathbf{z}$. Such assumptions are typically characterized by three
conditions: First, a _source condition_ limiting complexity of the source set
(set of possible solutions $f^{\dagger}$), second, mapping properties of the
forward operator, e.g. _smoothness condition_ for linear operators and, third,
_design measure_ $\nu$ generating points $(x_{i})_{i=1}^{m}\subset X^{m}$. As
we will see later, these conditions are often intertwined and the analysis of
different penalties has been coupled with varying assumptions on these
conditions. Here, with some cost of generality, we try to unify the exposition
for clarity.

###  1.1 Literature overview

The literature on statistical learning, reproducing kernel Hilbert spaces and
kernel regression, and regularization theory is vast and beyond the scope of
this treatise. We guide interested readers to [3, 27, 28, 36, 75, 82] for
classical sources on these subjects.

Connections of kernel regression methods to regularization theory were first
investigated in [32, 34, 53]. Early work on upper rates of convergence in a
reproducing kernel Hilbert space was carried out in [26] utilizing a covering
number technique. Subsequent work includes [72, 73, 4, 86, 18] culminating to
the minimax optimal rates demonstrated in [13]. In this paper, we focus on
results on random design regression in inverse problems introduced in [13, 68]
and extended to non-linear problems [67, 65, 66], and projection and convex
regularization strategies in [14, 43]. Here, let us note that convex penalties
of type
$G(f)=\frac{1}{p}\mathopen{}\mathclose{{}\left\|{f}}\right\|_{\mathcal{H}}^{p}$
were considered in statistical learning in [57, 76]. However, these results
are not directly comparable to the inverse problem literature as, among other
things, convergence takes place on a different domain
$\mathcal{H}\neq\mathcal{H}^{\prime}$.

The literature around statistical inverse problems has also expanded rapidly
during last two decades. For a general overview, we refer to [50, 77, 20, 61].
Early works involving spectral regularization methods include [8, 38, 4, 86,
33]. Projection methods are also well-studied and their minimax optimality is
understood in various settings with deterministic design [83, 63, 49, 35, 56,
55, 22, 24]. In the context of convex regularization strategies, Bregman
distances were first introduced in [17] and have been frequently used since as
an error measure for studying convergence rates in Banach spaces. Statistical
inverse problems with convex penalties have been studied in [16, 84]. For a
recent discussion on Bregman distances we refer to [6]. Furthermore,
adaptivity and parameter choice rules are studied in [19, 31, 54, 12, 9, 10,
74, 2].

For other related works involving non-Gaussian likelihood models, nonlinear
forward mapping or iterative methods, see e.g. [47, 46, 41, 48, 10, 25, 11].

This paper is organized as follows. In Section 2 we describe our problem
setting in detail and cover mathematical background to the field including
Hilbert scales and model class assumptions. In addition, we discuss minimax
optimality of statistical estimators. In Section 3 we overview the results on
spectral regularization (in the classical framework and in Hilbert scales),
projection methods and convex penalties. Moreover, literature on parameter
choice rules with random design is discussed. In Section 4 we describe non-
linear inverse problems isolating different scenarios on differentiability or
stability on the forward map in classical setup or in Hilbert scales. Finally,
we discuss an application in dynamic system monitoring in Section 5 and give a
brief outlook in Section 6.

##  2 Statistical Inverse Learning Theory

In this section, we present fundamental concepts in supervised inverse
statistical learning utilizing the least squares method [13, 67, 43].

###  2.1 Problem setting

Statistical Model and Noise Assumption. In the context of supervised inverse
learning, the objective is to learn a function from random samples of input-
output pairs $(x,y)\in X\times Y$. Here, we assume that the input space $X$ is
a Polish space, such as $\mathbb{R}^{d}$, and the output space is denoted as
$Y=\mathbb{R}$. The product space $Z=X\times\mathbb{R}$ is equipped with a
probability measure $\rho$, defined on the Borel $\sigma$-algebra of
$X\times\mathbb{R}$.

The probability distribution $\rho$ is only known through a training set
${\mathbf{z}}=\mathopen{}\mathclose{{}\left((x_{1},y_{1}),\ldots,(x_{m},y_{m})}\right)\in
Z^{m}$, consisting of pairs sampled independently and identically according to
$\rho$. Given this training set ${\mathbf{z}}$, the goal in supervised inverse
learning is to find an estimate $f_{{\bf z}}:X\to\mathbb{R}$ so that, for a
new unlabeled input $x\in X$, the value $A(f_{{\bf z}})(x)$ serves as a
reliable approximation of the true label $y$.

To formalize this objective, we introduce the notion of expected square loss
for any measurable function $f:X\to\mathbb{R}$:

${\mathcal{E}}(f)=\int_{X\times\mathbb{R}}(A(f)(x)-y)^{2}\,d\rho(x,y).$

This loss function captures the discrepancy between the predicted values and
the true labels, with pairs $(x,y)$ that are more likely to be sampled
exerting a greater contribution to the overall error. The expected loss can be
expressed in an alternative form, assuming that

$\int_{X\times\mathbb{R}}y^{2}\,d\rho(x,y)<\infty.$

We recall a useful integral decomposition. Let
$h:X\times\mathbb{R}\to\mathbb{R}$ be a measurable function. Then, the
integral of $h$ over $X\times\mathbb{R}$ can be decomposed as

$\int_{X\times\mathbb{R}}h(x,y)\,d\rho(x,y)=\int_{X}\int_{\mathbb{R}}h(x,y)\,d\rho(y|x)\,d\nu(x),$

where $\nu$ represents the marginal measure on $X$, and $\rho(\cdot|x)$ is the
conditional probability measure on $\mathbb{R}$ given $x\in X$. Recall that
the existence of $\rho(\cdot|x)$ is ensured by our assumption on $Z$, see e.g.
[71]. Denoting by

$g_{\rho}(x)=\int_{\mathbb{R}}y\,d\rho(y|x),$

the conditional mean function, the expected risk can be written as

$\mathcal{E}(f)=\mathopen{}\mathclose{{}\left\|{A(f)-g_{\rho}}}\right\|^{2}_{\nu}+{\mathcal{E}}(g_{\rho}),$

where we denote the $L^{2}$-norm on $X$ by
$\mathopen{}\mathclose{{}\left\|{g}}\right\|_{\nu}^{2}=\int_{X}|g(x)|^{2}\,d\nu(x)$.
Therefore, finding the minimizer of the expected risk is equivalent to
identifying the minimizer of the first term on the right side in the above
equation. In other words, the minimizers of expected risk typically coincide
with the classical minimum norm solutions of inverse problems with the
difference that data fidelity is weighted by the design measure. Moreover,
unlike in statistical learning we are interested in providing convergence
rates for the unknown $f$ assuming some ground truth.

######  Remark 2.1.

When there exists an element $f_{\rho}\in\mathcal{D}(A)~{}$ such that
$g_{\rho}=A(f_{\rho})$, then $f_{\rho}$ is a minimizer of the expected risk.
In the scenario of centered noise $\int_{\mathbb{R}}\varepsilon d\rho(y|x)=0$,
we obtain from the equations (1.1), (1.2) that
$\int_{\mathbb{R}}yd\rho(y|x)=g^{\dagger}(x)=A(f^{\dagger})(x)$. Consequently,
in this case for the injective operator $A$, the minimizer of the expected
risk, denoted as $f_{\rho}$, is equal to $f^{\dagger}$.

In what follows, we state our assumptions on the true solution of the inverse
problem (1.1) and the noise in (1.2), see also [13].

######  Assumption 1 (True solution).

The conditional expectation of $y$ given $x$, with respect to the probability
distribution $\rho$, is guaranteed to exist almost surely (a.s.). Furthermore,
there exists an element $f^{\dagger}\in\mathcal{D}(A)~{}$ such that:

$\int_{\mathbb{R}}y\,d\rho(y|x)=g^{\dagger}(x)=[A(f^{\dagger})](x),\text{ for
all }x\in X.$

The noise satisfies the following Bernstein condition:

######  Assumption 2 (Noise condition).

There exist positive constants $\Sigma>0$ and $M>0$ such that for any integer
$l\geq 2$

$\int_{\mathbb{R}}\mathopen{}\mathclose{{}\left\lvert
y-[A(f^{\dagger})](x)}\right\rvert^{l}d\rho(y|x)\leq\frac{1}{2}l!\Sigma^{2}M^{l-2},$

$\nu$-almost surely.

Assumption 2 characterizes the nature of the noise affecting the output $y$.
This equation remains valid under various scenarios, including cases where the
noise $\varepsilon$ is bounded and follows either a Gaussian distribution or a
sub-Gaussian distribution with zero mean and it is independent from $x$ [81].

RKHS Structure. We endow the space ${\mathcal{H}}^{\prime}$ with an additional
structure that allows us to exploit the rich theory of reproducing kernels.

######  Assumption 3 (RKHS Structure).

The space $\mathcal{H}^{\prime}$ is supposed to be a _reproducing kernel
Hilbert space_ (RKHS) of functions $g:X\to\mathbb{R}$ associated to a
measurable bounded kernel $K:X\times X\to\mathbb{R}$ with
$\kappa^{2}:=\sup\limits_{x\in X}K(x,x)<\infty.$ We additionally assume that
${\mathcal{H}}^{\prime}$ can be continuously embedded into
$\mathscr{L}^{2}(X,\nu)$. We denote the inclusion by
$I_{\nu}:{\mathcal{H}}^{\prime}\hookrightarrow\mathscr{L}^{2}(X,\nu)$.

We recall from [75], Chapter 4, that this assumption is equivalent to assuming
that for any $x\in X$, the evaluation functional

$\displaystyle S_{x}:{\mathcal{H}}^{\prime}$ $\displaystyle\to\mathbb{R}$
$\displaystyle g$ $\displaystyle\mapsto S_{x}(g):=g(x)$ (2.0)

is continuous. By Riesz’s representation theorem, there exists a unique
$F^{\prime}_{x}\in{\mathcal{H}}^{\prime}$ such that
$g(x)=\mathopen{}\mathclose{{}\left\langle
g,F^{\prime}_{x}}\right\rangle_{{\mathcal{H}}^{\prime}}$ and

$K(x,x^{\prime})=\mathopen{}\mathclose{{}\left\langle
F^{\prime}_{x},F^{\prime}_{x^{\prime}}}\right\rangle_{{\mathcal{H}}^{\prime}}\;,\quad
x,x^{\prime}\in X.$

In particular, for any $g\in\operatorname{Im}(A)$, we may write

$g(x)=[A(f)](x)=\mathopen{}\mathclose{{}\left\langle
A(f),F^{\prime}_{x}}\right\rangle_{{\mathcal{H}}^{\prime}},$

for some $f\in{\mathcal{H}}$ and for all $x\in X$. If $A$ is linear, we even
have

$[A(f)](x)=\mathopen{}\mathclose{{}\left\langle
f,F_{x}}\right\rangle_{{\mathcal{H}}},\quad F_{x}:=A^{*}F^{\prime}_{x}\;.$
(2.1)

From standard theory, we know that the feature map $x\mapsto F_{x}$ gives rise
to another reproducing kernel

$K_{A}(x,x^{\prime}):=\mathopen{}\mathclose{{}\left\langle
F_{x},F_{x^{\prime}}}\right\rangle_{{\mathcal{H}}},\quad x,x^{\prime}\in X,$

depending on the operator $A$. The associated RKHS ${\mathcal{H}}_{K_{A}}$
coincides with the image of ${\mathcal{H}}$ under $A$:
${\mathcal{H}}_{K_{A}}=\operatorname{Im}(A)$ and $A$ acts as a partial
isometry, see [13, Proposition 2.2], for an extended discussion. Hence, in the
linear case, we can forget about any structure on the space
${\mathcal{H}}^{\prime}$ and consider $A$ as an operator from ${\mathcal{H}}$
to ${\mathcal{H}}_{K_{A}}$ by assuming continuity of $S_{x}$ on
$\operatorname{Im}(A)$ only, for any $x\in X$. We formulate this as an
independent assumption and will analyze this case in Section 3 more deeply.

######  Assumption 4 (Linear Case).

The operator $A:{\mathcal{H}}\to{\mathcal{H}}^{\prime}:=\operatorname{Im}(A)$
is linear. Moreover, the evaluation functional
$S_{x}:{\mathcal{H}}^{\prime}\to\mathbb{R}$ defined by (2.1) is continuous for
all $x\in X$ and the map $x\mapsto S_{x}g$ is measurable for all
$g\in{\mathcal{H}}^{\prime}$.

We emphasize once more that in the general non-linear case in Assumption 3,
the kernel $K$ and thus the RKHS do not depend on the operator $A$, while in
the linear case, the operator $A$ is explicitly used to construct the kernel
$K_{A}$.

######  Remark 2.2 (Equivalence with Classical Kernel Learning Setting).

Under Assumption 4, the inverse problem (1.1) is naturally equivalent to the
classical kernel learning problem (the _direct_ problem) in non-parametric
regression in [4, 18, 34], that follows the model

$y_{i}={g^{\dagger}}(x_{i})+\varepsilon_{i}.$

The goal is to estimate the function $g^{\dagger}$. Kernel methods _posit_
that $g^{\dagger}$ belongs to some reproducing kernel Hilbert space and
construct an estimate $\hat{g}$ of $g^{\dagger}$ based on the observed data.
The reconstruction error $\hat{g}-g^{\dagger}$ is analyzed in
$L^{2}(\nu)$-norm or in RKHS-norm. In the setting of Assumption 4 and given
the model (1.1), the linear operator $A$ is a partial isometry and hence,
defining $\hat{f}:=A^{+}\hat{g}$, $f^{\dagger}:=A^{+}g^{\dagger}$ with
pseudoinverse $A^{+}$ gives

$\mathopen{}\mathclose{{}\left\|{f^{\dagger}-\hat{f}}}\right\|_{{\mathcal{H}}}=\mathopen{}\mathclose{{}\left\|{g^{\dagger}-\hat{g}}}\right\|_{{\mathcal{H}}^{\prime}},$

see [13]. Hence, a bound established for the direct learning setting in the
sense of the RKHS-norm reconstruction also applies to the inverse problem
reconstruction error. Furthermore, the eigenvalue decay conditions and the
source conditions involving the operator $T_{\nu}$ introduced in Section 2.3
are, via the same isometry, equivalent to similar conditions involving the
kernel integral operator in the direct learning setting, as considered, for
instance, in [4, 18, 34]. It follows that estimates in
${\mathcal{H}}^{\prime}$-norm available from those references are directly
applicable to the inverse learning setting.

######  Remark 2.3.

When employing kernel methods in statistical inverse learning problem (1.1),
one tackles the ensuing minimization problem:

$\inf\limits_{f\in\mathcal{H}}\int_{X\times\mathbb{R}}\mathopen{}\mathclose{{}\left((Af)(x)-y}\right)^{2}d\rho(x,y).$

Given Assumption 4 for $(Af)(x)=\mathopen{}\mathclose{{}\left\langle
f,F_{x}}\right\rangle_{\mathcal{H}}$, the previous expression can be
reformulated as:

$\inf\limits_{f\in\mathcal{H}}\int_{X\times\mathbb{R}}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\langle
f,F_{x}}\right\rangle_{\mathcal{H}}-y}\right)^{2}d\rho(x,y).$

By introducing an alternative probability measure
$\rho^{\prime}(F_{x},y)=\rho(x,y)$, the above simplifies to the form:

$\inf\limits_{f\in\mathcal{H}}\int_{X\times\mathbb{R}}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\langle
f,F_{x}}\right\rangle_{\mathcal{H}}-y}\right)^{2}d\rho^{\prime}(F_{x},y)=\inf\limits_{f\in\mathcal{H}}\int_{\mathcal{H}\times\mathbb{R}}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\langle
f,\xi}\right\rangle_{\mathcal{H}}-y}\right)^{2}d\rho^{\prime}(\xi,y).$

Therefore, minimizing the considered statistical inverse problem (1.2) in the
linear case is equivalent to solving the following problem:

$y_{i}=\mathopen{}\mathclose{{}\left\langle
f,\xi_{i}}\right\rangle_{\mathcal{H}}+\varepsilon_{i},\quad i=1,\ldots,m.$
(2.2)

###  2.2 Hilbert Scales

In the context of regularized statistical inverse problems, the utilization of
a Hilbert scale enables a delicate balance between the fidelity of the
estimate to the observed data and the smoothness of the estimated function
[60, 36]. Hilbert scales, comprised of a nested sequence of Hilbert spaces,
facilitate systematic control over the smoothness and complexity of functions.
Each space in the sequence is compactly embedded in the next, offering a
structured approach to managing these properties. Pioneering works by Tikhonov
[78, 79], extended the regularization process to a more general optimization
problem, minimizing a functional of the form

$\mathopen{}\mathclose{{}\left\|{A(f)-g^{\dagger}}}\right\|^{2}_{\nu}+\lambda\mathopen{}\mathclose{{}\left\|{Lf}}\right\|_{{\mathcal{H}}}^{2}\,,\quad
f\in{\mathcal{D}}(L)\;,$

where $L:{\mathcal{D}}(L)\to{\mathcal{H}}$ is an unbounded operator with
domain ${\mathcal{D}}(L)$. A notable instance occurs when $L$ is a
differential operator, effectively enforcing smoothness in the approximation,
given that the search space is a subset of ${\mathcal{D}}(L)$.

Remarkably, convergence of these regularized solutions can be achieved in the
original ${\mathcal{H}}$-norm even when no solution of (1.1) is an element of
${\mathcal{D}}(L)$. In this case, we are in the setting of _oversmoothing_. We
now give a rigorous definition:

######  Definition 2.4 (Hilbert Scale).

Let $a\in\mathbb{R}$ and let $L:{\mathcal{D}}(L)\to{\mathcal{H}}$ be a densely
defined unbounded self-adjoint positive operator, with
${\mathcal{D}}(L)\subset{\mathcal{H}}$. We set
$\tilde{\mathcal{H}}:=\bigcap_{a\in\mathbb{R}}{\mathcal{D}}(L^{a})$ and define
on $\tilde{\mathcal{H}}$ the inner product and norm

$\mathopen{}\mathclose{{}\left\langle
f_{1},f_{2}}\right\rangle_{a}:=\mathopen{}\mathclose{{}\left\langle
L^{a}f_{1},L^{a}f_{2}}\right\rangle_{{\mathcal{H}}}\;,\quad\mathopen{}\mathclose{{}\left\|{f}}\right\|_{a}:=\mathopen{}\mathclose{{}\left\|{L^{a}f}}\right\|_{\mathcal{H}}\;.$

The Hilbert space ${\mathcal{H}}_{a}$ is defined as the completion of
$\tilde{\mathcal{H}}$ with respect to the norm
$\mathopen{}\mathclose{{}\left\|{\cdot}}\right\|_{a}$. The scale
$({\mathcal{H}}_{a})_{a\in\mathbb{R}}$ is called the _Hilbert scale induced by
$L$_.

Note that by spectral theory, the powers $L^{a}$ are defined on
$\tilde{\mathcal{H}}$. Moreover, for any $-\infty<a<a^{\prime}<\infty$, we
have ${\mathcal{H}}_{a^{\prime}}\subset{\mathcal{H}}_{a}$, where the embedding
is dense and continuous, see [36, Section 8.4] for mathematical details.

In many instances, it is assumed that the scale of Hilbert spaces corresponds
to a scale of Sobolev spaces, and the smoothing properties of the underlying
operator $A$ are assessed in relation to this scale. Additionally, the
smoothness of the solution is characterized by assuming it belongs to a
specific space within this scale. Additional examples of Hilbert Scales
pertinent to learning in RKHSs can be found in the work of [59].

There has been growing interest in regularization within Hilbert scales for
statistical inverse problems, particularly in the context of the Bayesian
approach [39, 1]. Recently, this approach gained attention in statistical
learning and statistical inverse learning [65, 68, 59].

###  2.3 A Priori Assumptions

In classical statistical learning theory for non-parametric regression, the
convergence speed of an estimator to the target function can be exceedingly
slow as the sample size grows, a phenomenon known as the _No Free Lunch
Theorem_. To address this issue, the set of all data-generating distributions
is indirectly limited, for more efficient learning. A priori assumptions play
a critical role in this process by enabling the inclusion of prior knowledge
and information, ultimately leading to faster rates of convergence.

In what follows, we incorporate prior knowledge in two ways: First, we gauge
the _degree of ill-posedness_ of the inverse problem by assuming a certain
decay behavior of the eigenvalues of a second moment operator associated to
$A$. The faster the decay, the more ill-posed the problem is. Secondly, we
impose a smoothness/ sparsity constraint on the target function $f^{\dagger}$.
The smoother the target, the easier the reconstruction is. We now provide a
formalization of both concepts utilizing the following covariance operator
related to feature maps.

######  Definition 2.5 (Uncentered Covariance Operator).

Let $\Phi:X\to{\mathcal{H}}$ be a measurable feature map with

$\mathbb{E}_{\nu}[\mathopen{}\mathclose{{}\left\|{\Phi(x)}}\right\|_{\mathcal{H}}^{2}]<\infty.$
(2.3)

The _uncentered covariance operator_ associated to $\Phi$ is defined to be the
operator $T_{\nu}:{\mathcal{H}}\to{\mathcal{H}}$ with

$T_{\nu}:=\mathbb{E}_{\nu}[\Phi(x)\otimes\Phi(x)]=\int_{X}\mathopen{}\mathclose{{}\left\langle\cdot,\Phi(x)}\right\rangle\Phi(x)\;d\nu(x).$

The operator $T_{\nu}$ is positive, self-adjoint and compact, in particular
trace class by (2.3) with $\sum_{j\in\mathbb{N}}\sigma_{j}<\infty$, where
$(\sigma_{j})_{j\in\mathbb{N}}$ denotes the decreasing sequence of eigenvalues
of $T_{\nu}$.

Depending on the context, linear or non-linear problems, we will consider
different feature maps $\Phi$, giving rise to different operators $T_{\nu}$.
However, the non-linear approach contains the linear approach as a special
case, as exemplified below.

######  Example 2.6 (Linear Case).

If $A:{\mathcal{H}}\to{\mathcal{H}}^{\prime}$ is linear, we let
$\Phi(x):=A^{*}F_{x}^{\prime}=F_{x}$, see (2.1). Setting $B_{\nu}=I_{\nu}A$,
we also have the alternative representation of the covariance operator
$T_{\nu}=B_{\nu}^{*}B_{\nu}$, see e.g. [13, 18].

######  Example 2.7 (Non-Linear Case).

Assume that $A:{\mathcal{H}}\to{\mathcal{H}}^{\prime}$ is Fréchet
differentiable at $f^{\dagger}$ with derivative
$A^{\prime}(f^{\dagger}):{\mathcal{H}}\to{\mathcal{H}}^{\prime}$. We let
$\Phi(x):=A^{\prime}(f^{\dagger})^{*}F^{\prime}_{x}$. Setting
$B_{\nu}=I_{\nu}A^{\prime}(f^{\dagger})$, we again obtain
$T_{\nu}:=B_{\nu}^{*}B_{\nu}$. Note that we recover the linear case since
$A^{\prime}(f^{\dagger})=A$ if $A$ is linear.

In our given examples we observe that the second moment operator $T_{\nu}$
depends on both, the marginal distribution $\nu$ and the operator $A$ itself,
through the feature map $\Phi$. The properties of $T_{\nu}$ are crucial for
deriving fast rates.

######  Remark 2.8.

(Normal Equation) Let $A$ be a linear operator, and let $f^{\dagger}$ be the
expected risk minimizer, i.e., the least-squares solution of the model (1.1)
given by

$f^{\dagger}=\mathop{\rm
arg\,min}\limits_{f\in{\mathcal{H}}}\mathopen{}\mathclose{{}\left\|{B_{\nu}f-g^{\dagger}}}\right\|_{\nu}=\mathop{\rm
arg\,min}\limits_{f\in{\mathcal{H}}}{\mathcal{E}}(f)\;.$

This minimization is equivalent to the validity of the Gaussian normal
equation [36]

$T_{\nu}f^{\dagger}=B_{\nu}^{*}g^{\dagger}.$ (2.4)

Eigenvalue Decay. We impose a-priori assumptions on the set of marginal
distributions on $X$ by assuming a decay rate for the eigenvalues of the
covariance operator $T_{\nu}$, reflecting the severity of the ill-posed nature
of (1.1).

######  Assumption 5.

Suppose the eigenvalues
$\mathopen{}\mathclose{{}\left\\{\sigma_{i}}\right\\}_{i\in\mathbb{N}}$ of the
integral operator $T_{\nu}$ follow a polynomial decay: For fixed positive
constants $\alpha,\beta$ and $b<1$, we asume that

  1. (i)

$\sigma_{i}\leq\beta i^{-\frac{1}{b}}$, for all $i\in\mathbb{N}$,

  2. (ii)

$\alpha i^{-\frac{1}{b}}\leq\sigma_{i}$, for all $i\in\mathbb{N}$,

  3. (iii)

$\exists\gamma>0,~{}i_{0}\geq 1~{}~{}\text{s.t.}~{}~{}\sigma_{i}\leq
2^{\gamma}\sigma_{2i}$, for all $i\in\mathbb{N}$.

Since the covariance operator $T_{\nu}$ depends on the marginal $\nu$,
Assumption 5 imposes a restriction on the class of marginal distributions
${\mathcal{P}}(X)$ on $X$. We therefore introduce the following sub-classes:

${\mathcal{P}}_{\text{poly}}^{<}(\beta,b):=\\{\nu\in{\mathcal{P}}(X):\sigma_{i}\leq\beta
i^{-\frac{1}{b}}\;\;\forall\;\;i\in\mathbb{N}\\}\;,$ (2.5)

${\mathcal{P}}_{\text{poly}}^{>}(\alpha,b):=\\{\nu\in{\mathcal{P}}(X):\sigma_{i}\geq\alpha
i^{-\frac{1}{b}}\;\;\forall\;\;i\in\mathbb{N}\\}\;,$ (2.6)

${\mathcal{P}}_{\text{str}}^{>}(\gamma):=\\{\nu\in{\mathcal{P}}(X):\sigma_{i}\leq
2^{\gamma}\sigma_{2i}\;\;\forall\;\;i\geq
i_{0}\;,\;\;i_{0}\in\mathbb{N}\\}\;.$ (2.7)

######  Remark 2.9 (Effective Dimension).

In statistical learning theory, assumptions about marginal distributions are
expressed more broadly using the _effective dimension_ , also known as
_statistical dimension_ , representing the effective number of degrees of
freedom. In statistical modeling, particularly with high-dimensional data, the
effective dimension denotes the parameters or features significantly
influencing data variation. This is crucial when the actual data dimension is
high, yet numerous features may be redundant or carry little information. For
$\lambda>0$, it may be defined as

$\mathcal{N}(\lambda):=\operatorname{Tr}\mathopen{}\mathclose{{}\left((T_{\nu}+\lambda
I)^{-1}T_{\nu}}\right)\;.$

Since $T_{\nu}$ is trace class and bounded, we immediately get

$\mathcal{N}(\lambda)\leq\kappa^{2}\lambda^{-1}\;.$

However, an upper bound for the decay rate of the eigenvalues of $T_{\nu}$
implies an even sharper bound for ${\mathcal{N}}(\lambda)$. Specifically, for
any $\nu\in{\mathcal{P}}_{\text{poly}}^{<}(\beta,b)$ we find

$\mathcal{N}(\lambda)\leq\frac{\beta b}{b-1}\lambda^{-b},$

see [18, Proposition 3]. Note that in our context of inverse learning, bounds
on ${\mathcal{N}}(\lambda)$ also depend on $A$, incorporating properties of
the inverse problem into the statistical dimension naturally.

A-Priori Smoothness.  Another way to integrate prior knowledge into solving
inverse problems is by means of _(approximate) source conditions_. In
statistical inverse problems, a source condition describes the regularity or
smoothness of the true solution or unknown parameter being estimated. The
choice of a source condition is motivated by the specific characteristics of
the problem at hand and the available information about the true solution.

In the context of regularization of linear ill-posed inverse problems, the
a-priori-smoothness is expressed in terms of powers of the covariance operator
$T_{\nu}$, to be considered as a smoothing operator:

######  Assumption 6 (Source condition).

Let $r>0$, $R_{0}>0$. We assume that the solution $f^{\dagger}$ is a member of
the class $\Omega(r,R_{0},T_{\nu})$ with

$\Omega(r,R_{0},T_{\nu}):=\mathopen{}\mathclose{{}\left\\{f\in\mathcal{H}:f=T_{\nu}^{r}v\text{
for some }v\in{\mathcal{H}}\text{ with
}\mathopen{}\mathclose{{}\left\|{v}}\right\|_{\mathcal{H}}\leq
R_{0}}\right\\}.$

In the context of regularization schemes within Hilbert Scales, the source
condition is considered in terms of the known unbounded operator $L$, which is
involved in the regularization approach.

######  Assumption 7 (General Source Condition).

For $s>0$ and $R_{0}>0$, we let

$f^{\dagger}\in\Omega(s,R_{0},L^{-1}):=\mathopen{}\mathclose{{}\left\\{f\in\mathcal{H}:f=L^{-s}v\text{
for some }v\in{\mathcal{H}}\text{ with
}\mathopen{}\mathclose{{}\left\|{u}}\right\|_{\mathcal{H}}\leq
R_{0}}\right\\}.$

Source conditions are generally difficult to verify in practice or sometimes
may not be pre-specified. In this case, the concept of _distance functions_
can be employed to quantify the deviation from an assumed level of smoothness,
serving as a benchmark. This benchmark can be chosen by the user. For a given
_benchmark smoothness_ $p\geq 1$, we define the corresponding distance
function as specified below.

######  Definition 2.10 (Distance Function).

Let $R>0$, $p\geq 1$ and $L$ be a generator of a Hilbert scale. We define the
_distance function_ $d_{p}:[0,\infty)\to[0,\infty)$ by

$\displaystyle
d_{p}(R)=\inf\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left\|{L(f-f^{\dagger})}}\right\|_{\mathcal{H}}:f=L^{-p}u\text{
for some }u\in{\mathcal{H}}\text{ with
}\mathopen{}\mathclose{{}\left\|{u}}\right\|_{\mathcal{H}}\leq R}\right\\}.$
(2.8)

The distance function is positive, decreasing, convex, and continuous for all
values of $0\leq R<1$. Hence, the minimizer exists. As $R$ approaches
infinity, the distance function approaches $0$, as detailed in [44]. Note that
the distance function can be bounded explicitly in terms of $R_{0}>0$ provided
$f^{\dagger}$ satisfies a more general source condition Assumption 7. More
precisely, for $s<p$ we have

$d_{p}(R)\leq(R_{0})^{\frac{p-1}{p-s}}R^{\frac{1-s}{p-s}},$ (2.9)

see [45, Theorem 5.9]. If $f^{\dagger}\in{\mathcal{H}}_{p}$, i.e.
$f^{\dagger}=L^{-p}u$ for some $u\in{\mathcal{H}}$, then obviously
$d_{p}(R)=0$.

In the Hilbert scale scenario when the source condition is considered in terms
of $L^{-1}$, the following Link condition is essential to transfer the
smoothness in terms of the covariance operator $L^{-1}T_{\nu}L^{-1}$ to the
smoothness with respect to the operator $L^{-1}$ and vice versa.

######  Assumption 8 (Link condition).

For some $p\geq 1$, $0<a<\frac{1}{2}$ and $1\leq c<\infty$ we have

$\mathopen{}\mathclose{{}\left\|{L^{-p}u}}\right\|_{\mathcal{H}}\leq\mathopen{}\mathclose{{}\left\|{(L^{-1}T_{\nu}L^{-1})^{ap}u}}\right\|_{\mathcal{H}}\leq
c^{p}\mathopen{}\mathclose{{}\left\|{L^{-p}u}}\right\|_{\mathcal{H}},\quad
u\in\mathcal{H}.$ (2.10)

######  Remark 2.11.

The Weyl Monotonicity Theorem [7, Cor. III.2.3] establishes a relation between
the eigenvalues $\sigma_{j}(L^{-1})$ and
$\sigma_{j}((L^{-1}T_{\nu}L^{-1})^{a})$, namely,

$\sigma_{j}(L^{-1})\leq\sigma_{j}((L^{-1}T_{\nu}L^{-1})^{a})\leq\sigma_{j}(\beta
L^{-1})\;,\quad j\in\mathbb{N}.$

Hence, we observe that when the eigenvalues of $L^{-1}$ exhibit a decay rate,
say, ${\mathcal{O}}(j^{-w})$ for some $w>0$, the eigenvalues of
$L^{-1}T_{\nu}L^{-1}$ will follow a decay rate of
${\mathcal{O}}(j^{-\frac{w}{a}})$. What is more, one can show that $T_{\nu}$
has a spectral decay ${\mathcal{O}}(j^{(2-\frac{1}{a})w})$.

Class of Models. Let $\Omega\subseteq{\mathcal{H}}$ and denote by
${\mathcal{K}}(\Omega)$ the set of regular conditional distributions
$\rho(\cdot|\cdot)$ on ${\mathcal{B}}(\mathbb{R})\times X$ such that (1) and
(2) hold for some $f^{\dagger}\in\Omega$. Given $r>0$, $R_{0}>0$, the class of
models considered in most cases below is defined by

${\mathcal{M}}_{\theta}:=\\{\rho(dx,dy)=\rho(dy|x)\nu(dx)\;:\;\rho(\cdot|\cdot)\in{\mathcal{K}}(\Omega)\;,\nu\in{\mathcal{P}}^{\prime}\\}\;,$
(2.11)

where ${\mathcal{P}}^{\prime}$ specifies the admissible design measures. This
varies depending on the regularization method, see e.g. definitions in
(2.5)-(2.7) utilized in spectral regularization. The variable $\theta$
represents parametrization for data-generating distributions and our study
focuses on finding parameter choice rules that yield uniform concentration
rates over the domain of $\theta$. For the rigorous definition of this idea,
see next section for minimax optimality in learning. In practise, we mostly
consider $\theta$ as the pair
$(\Gamma,R_{0})\in\mathbb{R}_{+}\times\mathbb{R}_{+}$, i.e. the noise variance
and a priori norm bound to $f^{\dagger}$, but other parametrizations could be
considered as well.

We note that while ${\mathcal{M}}_{\theta}$ is specified by two main elements,
the source set $\Omega$ and the admissible design measures
${\mathcal{P}}^{\prime}$, their definitions are often intertwined and are both
often related to the mapping properties of the forward mapping. In particular,
for the convex penalties we will not specify $\Omega$ and
${\mathcal{P}}^{\prime}$ but instead consider ${\mathcal{M}}_{\theta}$
generated by a set of pairs $(f^{\dagger},\nu)$ satisfying suitable
conditions.

###  2.4 Minimax optimality

In statistical (inverse) learning theory, the speed at which a learning
algorithm $(\mathbf{z},\lambda)\mapsto f_{\mathbf{z},\lambda}$ converges to
the true underlying model as the sample size increases is referred to as the
_rate of convergence_. It quantifies how quickly the learned model approaches
the correct solution as more data becomes available. The rate of convergence
is a crucial concept in assessing the efficiency and performance of
statistical learning methods. A faster rate of convergence implies that the
estimator requires fewer observations to provide a more accurate estimate or
model. Typically, rates are discussed by means of probabilistic tail
inequalities, involving the confidence level $\eta\in(0,1]$ and the sample
size $m\in\mathbb{N}$. More specifically, they take the form

$\rho^{m}\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\leq
Ca_{m}\log\mathopen{}\mathclose{{}\left(\frac{1}{\eta}}\right)}\right\\}\geq
1-\eta.$ (2.12)

Here, $\rho^{m}$ is the joint distribution, the function $m\mapsto a_{m}$ is a
positive decreasing function and describes the specific rate as $m\to\infty$.
Recall from the previous section that we restrict the class of admissible
distributions to a subclass ${\mathcal{M}}_{\theta}$, where $\theta$ is a
parameter belonging to a parameter space $\Theta$. Then, (2.12) allows to
write

$\limsup\limits_{m\to\infty}\sup\limits_{\rho\in{\mathcal{M}}_{\theta}}\rho^{m}\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\geq\tau
a_{m}}\right\\}\leq e^{-\frac{\tau}{C}}\;,$

with $\tau=C\log\mathopen{}\mathclose{{}\left(\frac{1}{\eta}}\right)$. Hence,

$\lim\limits_{\tau\to\infty}\limsup\limits_{m\to\infty}\sup\limits_{\rho\in{\mathcal{M}}_{\theta}}\rho^{m}\mathopen{}\mathclose{{}\left\\{\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\geq\tau
a_{m}}\right\\}=0.$

Integration immediately leads to a bound in expectation:

$\displaystyle\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}^{p}}\right]^{\frac{1}{p}}=$
$\displaystyle\mathopen{}\mathclose{{}\left(\int\limits_{0}^{\infty}\rho^{m}\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}^{p}>t}\right)dt}\right)^{\frac{1}{p}}\leq
C_{p}a_{m},$

for some $C_{p}<\infty$, which implies

$\limsup\limits_{m\to\infty}\sup\limits_{\rho\in{\mathcal{M}}_{\theta}}\frac{\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}^{p}}\right]^{\frac{1}{p}}}{a_{m}}\leq
C_{p}<\infty.$

This motivates the following definitions. We refer to [4, 18, 13] for more
details.

######  Definition 2.12 (Upper Rate of Convergence).

A family of sequences $(a_{m,\theta})_{(m,\theta)\in\mathbb{N}\times\Theta}$
of positive numbers is called an _upper rate of convergence_ in
$\mathscr{L}^{p}(X,\nu)$ for the reconstruction error, over the family of
models $({\mathcal{M}}_{\theta})_{\theta\in\Theta}$, for the sequence of
estimated solutions
$(f_{\mathbf{z},\lambda_{m,\theta}})_{(m,\theta)\in\mathbb{N}\times\Theta}$,
using regularization parameters
$(\lambda_{m,\theta})_{(m,\theta)\in\mathbb{N}\times\Theta}$, if

$\sup\limits_{\theta\in\Theta}\limsup\limits_{m\to\infty}\sup\limits_{\rho\in{\mathcal{M}}_{\theta}}\frac{\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda_{m,\theta}}-f^{\dagger}}}\right\|_{\mathcal{H}}^{p}}\right]^{\frac{1}{p}}}{a_{m,\theta}}<\infty.$

######  Definition 2.13 (Weak and Strong Minimax Lower Rate of Convergence).

A family of sequences $(a_{m,\theta})_{(m,\theta)\in\mathbb{N}\times\Theta}$
of positive numbers is called a _weak minimax lower rate of convergence_ in
$\mathscr{L}^{p}(X,\nu)$ for the reconstruction error, over the family of
models $({\mathcal{M}}_{\theta})_{\theta\in\Theta}$, if

$\sup\limits_{\theta\in\Theta}\limsup\limits_{m\to\infty}\inf\limits_{f_{\cdot}}\sup\limits_{\rho\in{\mathcal{M}}_{\theta}}\frac{\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z}}-f^{\dagger}}}\right\|_{\mathcal{H}}^{p}}\right]^{\frac{1}{p}}}{a_{m,\theta}}>0.$

where the infimum is taken over all estimators, i.e., measurable mappings
$f_{\cdot}:(X\times\mathbb{R})^{m}\to\mathcal{H}$. It is called a _strong
minimax lower rate of convergence_ in $L^{p}$ if

$\sup\limits_{\theta\in\Theta}\liminf\limits_{m\to\infty}\inf\limits_{f_{\cdot}}\sup\limits_{\rho\in{\mathcal{M}}_{\theta}}\frac{\mathbb{E}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z}}-f^{\dagger}}}\right\|_{\mathcal{H}}^{p}}\right]^{\frac{1}{p}}}{a_{m,\theta}}>0.$

######  Definition 2.14 (Minimax Optimal Rate of Convergence).

The sequence of estimated solutions $(f_{\mathbf{z},\lambda_{m,\theta}})_{m}$
using the regularization parameters
$(\lambda_{m},\theta)_{(m,\theta)\in\mathbb{N}\times\Theta}$ is called
_weak/strong minimax optimal_ in $\mathscr{L}^{p}(X,\nu)$ for the
interpolation norm of parameter $s\in[0,\frac{1}{2}]$, over the model family
$({\mathcal{M}}_{\theta})_{\theta\in\Theta}$, with rate of convergence given
by the sequence $(a_{m,\theta})_{(m,\theta)\in\mathbb{N}\times\Theta}$, if the
latter is a weak/strong minimax lower rate as well as an upper rate for
$(f_{\mathbf{z},\lambda_{m,\theta}})_{(m,\theta)\in\mathbb{N}\times\Theta}$.

##  3 Linear problems

In this Section we present the theory for regularization of linear inverse
learning problems. We begin with introducing a suitable discretization
strategy based on random sampling that is able to provide practical
implementations of our regularized solutions. We move forward with presenting
three major classes of regularization methods: _spectral regularization_ ,
_regularization by projection_ , and more generally _convex regularization_.
For each of these methods we state convergence results and discuss their
minimax optimality.

###  3.1 Discretization by Random Sampling

Discretization is a crucial step in solving inverse problems because it
transforms complex, continuous mathematical models into a computationally
tractable format that can be approached using numerical methods. In our
framework, we employ the random sample
$\mathopen{}\mathclose{{}\left((x_{i},y_{i})}\right)_{i=1}^{m}$ to construct a
discretization of the operator $A$ in the following way: We define the
_sampling operator_ $S_{\mathbf{x}}:{\mathcal{H}}^{\prime}\to\mathbb{R}^{m}$
by

$S_{\mathbf{x}}g:=(g(x_{1}),...,g(x_{m})).$

We further define
$B_{\mathbf{x}}:=S_{\mathbf{x}}A:\mathcal{H}\to\mathbb{R}^{m}$, yielding an
empirical approximation to $B_{\nu}$ and
$T_{\mathbf{x}}:=B_{\mathbf{x}}^{*}B_{\mathbf{x}}$, yielding an approximation
to the covariance $T_{\nu}$.

The empirical expected risk based on the sample is given by

$\mathcal{E}_{\mathbf{z}}(f)=\frac{1}{m}\sum_{i=1}^{m}\mathopen{}\mathclose{{}\left([A(f)](x_{i})-y_{i}}\right)^{2}=\frac{1}{m}\mathopen{}\mathclose{{}\left\|{B_{\mathbf{x}}f-\mathbf{y}}}\right\|_{2}^{2}.$
(3.1)

Minimizing ${\mathcal{E}}_{\mathbf{z}}$ amounts to finding $f_{{\bf
z}}\in{\mathcal{H}}$ such that $B_{\mathbf{x}}f_{{\bf z}}=\mathbf{y}$, which
can be seen as a noisy, randomly discretized version of (1.1) and one would
expect that with increasing sample size, model (1.1) will be better
approximated by its empirical counterpart. However, it is not straightforward
to measure the accuracy of the random approximation as the ranges of the
operators defining the inverse problems are different. Considering the normal
equations for both problems

$T_{\nu}f^{\dagger}=B_{\nu}^{*}g^{\dagger},\quad T_{\mathbf{x}}f_{{\bf
z}}=B_{\mathbf{x}}^{*}\mathbf{y}$ (3.2)

suggest to compare both inverse problems by bounding the norms

$\mathopen{}\mathclose{{}\left\|{T_{\nu}-T_{\mathbf{x}}}}\right\|,\quad\mathopen{}\mathclose{{}\left\|{B_{\nu}^{*}g^{\dagger}-B_{\mathbf{x}}^{*}\mathbf{y}}}\right\|_{{\mathcal{H}}}.$

Since both norms are random quantities, estimating the perturbation due to
random noise and random sampling requires probabilistic tools, such as
concentration of measure results, see [18, 34]. We provide a result that shows
that $T_{\nu}$ and $B_{\nu}^{*}g^{\dagger}$ can be well approximated by their
respective empirical counterparts.

######  Proposition 3.1 (Lemma 9 in [4]).

Let $\eta\in(0,1]$ and suppose that Assumption 2 is satisfied. With
probability at least $1-\eta$ we have

$\displaystyle\mathopen{}\mathclose{{}\left\|{T_{\nu}-T_{\mathbf{x}}}}\right\|$
$\displaystyle\leq\frac{2\sqrt{2}\kappa^{2}}{\sqrt{n}}\log(4/\eta),$
$\displaystyle\mathopen{}\mathclose{{}\left\|{B_{\nu}^{*}g^{\dagger}-B_{\mathbf{x}}^{*}\mathbf{y}}}\right\|_{{\mathcal{H}}}$
$\displaystyle\leq
2\kappa\mathopen{}\mathclose{{}\left(\frac{M}{n}+\frac{\sqrt{\Sigma}}{\sqrt{n}}}\right)\log(4/\eta)\;.$

######  Remark 3.2 (Measure theoretic point of view).

If we let $\hat{\nu}=\frac{1}{m}\sum_{j=1}^{m}\delta_{x_{j}}$ denote the
empirical measure of the input data, we can identify
$\mathscr{L}^{2}(X,\hat{\nu})$ with $\mathbb{R}^{m}$, equipped with the inner
product
$\langle{\mathbf{y}},{\mathbf{y}^{\prime}}\rangle_{m}:=\frac{1}{m}\sum_{i=1}^{m}y_{i}y^{\prime}_{i}$.
It is then straightforward to see that the operators $B_{\nu}$, $T_{\nu}$
reduce to their respective empirical counterparts when replacing $\nu$ with
$\hat{\nu}$, i.e. $B_{\hat{\nu}}=B_{\mathbf{x}}$ and
$T_{\hat{\nu}}=T_{\mathbf{x}}$.

###  3.2 Optimality of spectral regularization methods

Regularization strategies are essential in solving inverse problems to address
issues related to ill-posedness, instability, and noise sensitivity in the
input data. Regularization also helps to prevent overfitting by penalizing
overly complex solutions, encouraging smoother or simpler solutions that are
more likely to generalize well.

Common regularization strategies in the context of inverse problems include
Tikhonov regularization (also known as ridge regression), total variation
regularization, and sparse regularization. These techniques add mathematical
constraints or penalties to the optimization problem, promoting solutions that
are both faithful to the data and have desirable properties such as smoothness
or sparsity. The choice of regularization method depends on the specific
characteristics of the problem and the desired properties of the solution. In
what follows, we will focus on spectral methods.

####  3.2.1 Spectral regularization methods

We give now the general definition of spectral methods, followed by classical
examples.

######  Definition 3.3 (Spectral regularization).

We say that a family of functions $\\{g_{\lambda}\\}_{\lambda>0}$ with
$g_{\lambda}:[0,\kappa^{2}]\to\mathbb{R}$ is a _spectral regularization
scheme_ or _spectral filter_ if there exist constants $D<\infty$, $E<\infty$,
and $\gamma_{0}<\infty$ such that for all $\lambda>0$

  1. (i)

$\sup\limits_{t\in[0,\kappa^{2}]}\mathopen{}\mathclose{{}\left\lvert
tg_{\lambda}(t)}\right\rvert\leq D$,

  2. (ii)

$\sup\limits_{t\in[0,\kappa^{2}]}\mathopen{}\mathclose{{}\left\lvert
g_{\lambda}(t)}\right\rvert\leq\frac{E}{\lambda}$ and

  3. (iii)

$\sup\limits_{t\in[0,\kappa^{2}]}\mathopen{}\mathclose{{}\left\lvert
r_{\lambda}(t)}\right\rvert\leq\gamma_{0}$, where
$r_{\lambda}(t)=1-g_{\lambda}(t)t$.

Moreover, the maximal $q$ satisfying the condition

$\sup\limits_{t\in[0,\kappa^{2}]}\mathopen{}\mathclose{{}\left\lvert
r_{\lambda}(t)}\right\rvert t^{q}\leq\gamma_{q}\lambda^{q},$

for some constant $\gamma_{q}$ and any $\lambda>0$ is said to be the
_qualification_ of the regularization scheme $g_{\lambda}$.

######  Example 3.4 (Spectral cut-off).

Consider the spectral cut-off or truncated singular value decomposition (TSVD)
defined by

$g_{\lambda}(t)=\mathopen{}\mathclose{{}\left\\{\begin{array}[]{ll}t^{-1}&,\quad\text{if}\quad
t\geq\lambda;\\\ 0&,\quad\text{if}\quad t<\lambda\end{array}}\right.$

The qualification $q$ could be any positive number and
$D=E=\gamma_{0}=\gamma_{q}=1$.

######  Example 3.5 (Explicit regularization: (Iterated) ridge regression).

For $l\in\mathbb{N}$, the iterated Tikhonov regularization is given by the
function

$g_{\lambda}(t)=\sum\limits_{i=1}^{l}\lambda^{i-1}(\lambda+t)^{-i}=\frac{1}{t}\mathopen{}\mathclose{{}\left(1-\frac{\lambda^{l}}{(\lambda+t)^{l}}}\right)$

In this case we have $q=l$, $D=\gamma_{0}=\gamma_{q}=1$, and $E=l$. For $l=1$,
this scheme corresponds to the Tikhonov regularization or ridge regression.

######  Example 3.6 (Implicit regularization: Gradient decent).

Gradient descent with constant stepsize $\eta>0$, a.k.a. Landweber iteration,
to the function $g_{\lambda}(t)=\sum\limits_{k=1}^{u}\eta(1-\eta t)^{u-k}$
with $\eta\in(0,\kappa^{2}]$ with the choice $\lambda=(\eta u)^{-1}$. Here,
$D=E=\gamma_{0}=1$ and the qualification $q$ can be any positive number with
$\gamma_{q}=1$ if $0<q\leq 1$ and $\gamma_{q}=q^{q}$ if $q>1$.

####  3.2.2 Classical Error Bounds

Equation (3.2) suggests to define our regularized empirical solution as given
below.

######  Definition 3.7.

Let $\\{g_{\lambda}\\}_{\lambda>0}$ be a family of filter functions defined on
the spectrum of $T_{\nu}$ and let $((x_{j},y_{j}))_{j=1}^{m}$ be an i.i.d.
sample. The regularized empirical estimator for the solution to the linear
inverse problem (1.1) is defined by

$f_{\mathbf{z},\lambda}=g_{\lambda}(T_{\mathbf{x}})B_{\mathbf{x}}^{*}\mathbf{y}.$
(3.3)

In what follows we present the rates of convergence for the reconstruction
error for the estimator defined above.

######  Theorem 3.8 (Blanchard et al. [13]).

Suppose Assumptions 1, 2, 4, 5 (i), 6 are satisfied. Suppose further that
$\lambda\geq m^{-\frac{1}{b+1}}$ and the qualification is $q\geq r$. For all
$\eta\in(0,1]$, the regularized solution $f_{\mathbf{z},\lambda}$ (3.3)
satisfies with probability at least $1-\eta$ the upper bound

$\displaystyle\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\;\leq$
$\displaystyle\;C\mathopen{}\mathclose{{}\left(R_{0}\lambda^{r}+\frac{\Sigma}{\sqrt{m}\lambda^{\frac{b+1}{2}}}+\frac{R_{0}}{\sqrt{m}}+\frac{M}{m\lambda}}\right)\log^{3}\mathopen{}\mathclose{{}\left(\frac{8}{\eta}}\right),$
(3.4)

where $C$ depends on all model parameters
$r,~{}\kappa,~{}\gamma_{0},\gamma_{q},~{}D,~{}E,~{}b,~{}\beta$.

Note that this upper bound is given in high probability. Deriving the final
rate of convergence requires to carefully select the regularization parameter
$\lambda$, depending on the sample size and the a-priori assumptions. Setting
$\theta=(\Sigma,R_{0})$, a straightforward calculation shows that the choice

$\lambda_{m,\theta}=\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{2r+b+1}}$
(3.5)

equalizes the first and second term of (3.4). Moreover, the remaining terms
can be shown to be of lower order. We summarize these observations below.

######  Corollary 3.9.

Suppose the assumptions of Theorem 3.8 are satisfied. Choose $\lambda$
according to (3.5). With probability at least $1-\eta$, we have

$\displaystyle\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\;\leq$
$\displaystyle\;CR_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2r}{2r+b+1}}\log^{3}\mathopen{}\mathclose{{}\left(\frac{8}{\eta}}\right),$

provided that $m$ is sufficiently large and where $C$ depends on all model
parameters $r$, $\kappa$, $\gamma_{0}$, $\gamma_{q}$, $D$, $E$, $b$, $\beta$.

Following the lines of Section 2.4, we immediately get the upper rate of
convergence through integration.

######  Theorem 3.10.

Suppose the assumptions of Corollary 3.9 hold true. Then, the sequence

$\displaystyle
a_{m,\theta}=R_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2r}{2r+b+1}}$
(3.6)

is an upper rate of convergence in $\mathscr{L}^{p}(X,\nu)$ for all $p>0$, for
the reconstruction error, over the family of models
${\mathcal{M}}(r,R_{0},{\mathcal{P}}_{\text{poly}}^{<}(\beta,b))$ for the
sequence of solutions of the regularization scheme (3.3) with regularization
parameter (3.5).

Notably, the sequence given in (3.6) can also shown to be a lower bound for
the reconstruction error:

######  Theorem 3.11.

Let Assumptions 1, 2, 4, and 6 hold true. Then the sequence
$(a_{m,\theta})_{\theta}$ defined in (3.6) is a weak lower rate of convergence
in $\mathscr{L}^{p}(X,\nu)$ for all $p>0$ for the reconstruction error, over
the family of models
${\mathcal{M}}(r,R_{0},{\mathcal{P}}_{\text{poly}}^{>}(\alpha,b))$. Moreover,
the sequence $(a_{m,\theta})_{\theta}$ is a strong lower rate of convergence
in $\mathscr{L}^{p}(X,\nu)$ for all $p>0$ over the family of models
${\mathcal{M}}(r,R_{0},{\mathcal{P}}_{\text{str}}^{>}(\gamma))$.

Finally, by Definition 2.14, we may conclude:

######  Theorem 3.12.

Suppose the Assumptions of Theorem 3.10 and Theorem 3.11 are satisfied. Then
the sequence of solutions $(f_{\mathbf{z},\lambda})_{m}$ using the
regularization parameters given in (3.5) is weak and strong minimax optimal in
$\mathscr{L}^{p}(X,\nu)$, for all $p>0$, for the reconstruction error, over
the model family ${\mathcal{M}}(r,R_{0},{\mathcal{P}}^{\prime})$, with
${\mathcal{P}}^{\prime}={\mathcal{P}}_{\text{poly}}^{<}(\beta,b)\cap{\mathcal{P}}_{\text{poly}}^{>}(\alpha,b))$,
${\mathcal{P}}^{\prime}={\mathcal{P}}_{\text{poly}}^{<}(\beta,b)\cap{\mathcal{P}}_{\text{str}}^{>}(\gamma)$,
respectively.

####  3.2.3 Error Bounds in Hilbert Scales

Now, we analyze the regularization schemes in Hilbert Scales for the
considered statistical inverse problem (1.1). The background of these schemes
is discussed in Section 2.2. To this end, let $L$ be an unbounded operator
with domain ${\mathcal{D}}(L)\subset{\mathcal{H}}$, generating a Hilbert
scale. Tikhonov regularization in Hilbert Scales is defined as:

$f_{\mathbf{z},\lambda}=\mathop{\rm
arg\,min}\limits\limits_{f\in\mathcal{D}(L)}\mathopen{}\mathclose{{}\left\\{\frac{1}{m}\sum_{i=1}^{m}\mathopen{}\mathclose{{}\left([Af](x_{i})-y_{i}}\right)^{2}+\lambda\mathopen{}\mathclose{{}\left\|{Lf}}\right\|_{\mathcal{H}}^{2}}\right\\},$
(3.7)

with explicit solution

$f_{\mathbf{z},\lambda}=L^{-1}(L^{-1}T_{\mathbf{x}}L^{-1}+\lambda
I)^{-1}L^{-1}B_{\mathbf{x}}^{*}\mathbf{y}.$

Clearly, the minimizer $f_{\mathbf{z},\lambda}$ belongs to $\mathcal{D}(L)$,
and formally we may introduce
$u_{\mathbf{z},\lambda}:=Lf_{\mathbf{z},\lambda}\in\mathcal{H}$. When
$f^{\dagger}\in\mathcal{D}(L)$, which we call the _regular case_ , we may set
$u^{\dagger}:=Lf^{\dagger}\in\mathcal{H}$. With this notation, we can rewrite
(1.1) as

$g^{\dagger}=Af^{\dagger}=AL^{-1}u^{\dagger},\quad
u^{\dagger}\in\mathcal{D}(L).$

The Tikhonov minimization problem (3.7) would reduce to the standard one

$\min_{u\in{\mathcal{H}}}\;\;\frac{1}{m}\sum\limits_{i=1}^{m}\mathopen{}\mathclose{{}\left([AL^{-1}u](x_{i})-y_{i}}\right)^{2}+\lambda\mathopen{}\mathclose{{}\left\|{u}}\right\|_{\mathcal{H}}^{2},$

albeit for a different operator $AL^{-1}$. The reconstruction error is given
by

$\mathopen{}\mathclose{{}\left\|{f^{\dagger}-f_{\mathbf{z},\lambda}}}\right\|_{\mathcal{H}}=\mathopen{}\mathclose{{}\left\|{L^{-1}(u^{\dagger}-u_{\mathbf{z},\lambda})}}\right\|_{\mathcal{H}}.$

Therefore, error bounds for $u^{\dagger}-u_{\mathbf{z},\lambda}$ in the weak
norm, i.e. in $\mathcal{H}_{-1}$, yield bounds for
$f^{\dagger}-f_{\mathbf{z},\lambda}$.

The above observations motivate us to define more general spectral
regularization schemes from Section 3.2.1 as given in the next definition.

######  Definition 3.13.

Let $\\{g_{\lambda}\\}_{\lambda>0}$ be a family of filter functions defined on
the spectrum of $L^{-1}T_{\nu}L^{-1}$ and let $((x_{j},y_{j}))_{j=1}^{m}$ be
an i.i.d. sample. We define the estimated solutions to the linear inverse
problem (1.1) in the framework of Hilbert scales by

$f_{\mathbf{z},\lambda}=L^{-1}g_{\lambda}(L^{-1}T_{\mathbf{x}}L^{-1})L^{-1}B_{\mathbf{x}}^{*}\mathbf{y}\;.$

As discussed in [68], convergence of these regularized solutions can be
achieved in the original ${\mathcal{H}}$-norm even when no solution of (1.1)
is an element of ${\mathcal{D}}(L)$. In this case, we are in the setting of
_oversmoothing_.

In the following theorem, the error bounds are discussed using the approximate
source condition in the absence of a source condition for the regular case. As
discussed above, the covariance operator in Hilbert scales is
$L^{-1}T_{\nu}L^{-1}$. Therefore, we consider the eigenvalue decay condition 5
(i) and the link condition 8 in terms of $L^{-1}T_{\nu}L^{-1}$. We assume that
the eigenvalues $\sigma_{i}$ of the operator $L^{-1}T_{\nu}L^{-1}$ decay
polynomially, i.e., $\sigma_{i}\leq\beta i^{-\frac{1}{b}}$. The smoothness
between the operators $L^{-1}$ $L^{-1}T_{\nu}L^{-1}$ is linked by $a$ from
Assumption 8, i.e.,
$\mathcal{R}(L^{-1})=\mathcal{R}((L^{-1}T_{\nu}L^{-1})^{a})$.

######  Theorem 3.14 (Rastogi et al. [68]).

Suppose Assumptions 1–3, 5 (i), 8 are satisfied and let $\eta\in(0,1]$. Let
further the qualification $q\geq a(p-1)$. Define the regularized empirical
solutions according to Definition 3.13. Then, with probability at least
$1-\eta$ we have

$\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\;\leq\;C\lambda^{a}\mathopen{}\mathclose{{}\left(d_{p}(R)+R\lambda^{a(p-1)}+\frac{\Sigma}{\sqrt{m}\lambda^{\frac{b+1}{2}}}+\frac{R}{\sqrt{m}}+\frac{M}{m\lambda}}\right)\log^{4}\mathopen{}\mathclose{{}\left(\frac{4}{\eta}}\right),$

where $C$ depends on the parameters $\kappa$, $\gamma_{0}$, $\gamma_{q}$, $D$,
$E$, $b$, $\beta$ and $d_{p}(R)$ is defined in (2.8).

The distance function is bounded explicitly in (2.9) under the source
condition 7. Then, we get the explicit rate of convergence in terms of the
sample size by balancing the error terms for the regularization parameter
choice. These rates are presented in the last two rows of Table 1 under the
different condition on the benchmark smoothness $p$. Note that if
$f^{\dagger}\in{\mathcal{D}}(L)$, we may write $f^{\dagger}=L^{-1}u$ for some
$u\in\mathcal{H}$. Thus, if $s<1$, we are in the oversmoothing case while for
$s\geq 1$, we are in the regular case.

Now, we summarize the results in Table 1. We present a comprehensive picture
of the convergence rates with the corresponding parameter choices and the
conditions in both regular and oversmoothing scenarios. The table shows
convergence rates $a_{m}$ described in (2.12), corresponding orders of the
regularization parameter choice $\lambda_{m}$, smoothness of the true solution
$s$, and additional constraints on the benchmark smoothness $p$. Here, $n$ is
some natural number which reflects in the multiplicative constant in the error
bounds. The first row represents the oversmoothing case, and the last two rows
represent the regular case. In the regular case ($p>1$), the convergence rates
depend on the benchmark smoothness $p$.

In the regular case, the authors [68] showed that the convergence rates are
optimal by relating the assumptions considered in standard regularization and
regularization in Hilbert scales. The optimal rates are achieved in the last
row of Table 1, provided the benchmark smoothness $p$ is set appropriately,
i.e., $s\leq p\leq s+\frac{b+1}{2a}$. In contrast, the rates are sub-optimal
in the oversmoothing case.

Table 1: Convergence rates of the regularized solution $f_{\mathbf{z},\lambda}$ for $a<\frac{1}{2}$, $a(p-1)\leq{q}$. Case | Convergence rates | Reg. param. | True | Conditions
---|---|---|---|---
| $a_{m}$ | $\lambda_{m}$ | Smooth. |
Oversmooth. | $R_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2as}{b+1}}$ | $\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{b+1}}$ | $s<1$ | $a\in[\frac{1}{n+1},\frac{1}{2}),~{}n\in\mathbb{N}$
Regular | $R_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{s}{p-1}}$ | $\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{1}{a(p-1)}}$ | $s\geq 1$ | $p\geq s+\frac{b+1}{2a}$
$R_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2as}{2as+b+1-2a}}$ | $\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{2as+b+1-2a}}$ | $s\leq p\leq s+\frac{b+1}{2a}$

######  Remark 3.15.

The regularization schemes have a saturation effect, meaning that the rates
cannot improved corresponding to the smoothness of the solution beyond the
qualification of the regularization scheme. Here, in the regular case, the
optimal rates are established for the range $as\leq q$, provided that the
scheme has qualification $q$. For standard regularization schemes, this would
correspond to the range $\frac{as}{1-2a}\leq q$, only.

###  3.3 Regularization by projection

Regularization by projection is a classical method that has been studied
widely [36]. In statistical inverse learning, it was studied recently in [43].
Notice that in [43] an additive Gaussian noise model is assumed. Therefore, in
this section, $\epsilon_{i}$ is zero-centered Gaussian with variance $\Gamma$.
Here, we modify this recent result to align with the Hilbert scale structure
discussed above. To set the stage, we assume that there is a sequence of
nested subspaces of $\mathcal{H}$ that provide a consistent approximation
scheme.

######  Definition 3.16.

Consider a series of finite-dimensional subspaces of $\mathcal{H}$ denoted as
$V_{n}$, $n\geq 1$ with the condition $\text{dim}V_{n}=n$. Additionally, let
$P_{n}:\mathcal{H}\to V_{n}\subset\mathcal{H}$ represent an orthogonal
projection. The sequence
$\mathopen{}\mathclose{{}\left\\{V_{n}}\right\\}_{n=1}^{\infty}$ is termed
‘admissible subspaces’ if it meets the criteria that $V_{n}\subset V_{n+1}$
for all $n\in\mathbb{N}$ and $\cup_{n=1}^{\infty}V_{n}=\mathcal{H}$.

Given training data with $m$ samples, we aim to find an approximation space
$V_{n}$ and an estimator $f_{n,m}\in V_{n}$ that concentrates to the true
solution $\hat{f}$ with a proper parameter choice rule $n=n(m)$. Let us define
the maximum likelihood (ML) estimator within $V_{n}$ as follows:

$f_{n,m}=\mathop{\rm arg\,min}\limits\limits_{f\in
V_{n}}\mathopen{}\mathclose{{}\left\|{S_{\mathbf{x}}[A(f)]-\mathbf{y}}}\right\|_{m}^{2}\quad\text{a.s.}$
(3.8)

In order to ensure that the estimator $f_{n,m}$ is almost surely well-defined,
i.e. the minimization problem has a unique solution, we need to further assume
that the design measure is assumed to be atomless, $A$ is injective and the
kernel $k$ is strictly positive-definite. For what follows, recall that
$B^{+}$ refers to the pseudoinverse of an operator
$B:\mathcal{H}\to\mathcal{H}$.

Here, we assume that the source set is specified by the approximation scheme
$\\{V_{n}\\}_{n=1}^{\infty}$, namely, we set

$\Theta(s,R_{0})=\mathopen{}\mathclose{{}\left\\{f\in\mathcal{H}:\mathopen{}\mathclose{{}\left\|{(I-P_{n})f}}\right\|_{\mathcal{H}}\leq
R_{0}(n+1)^{-s}\quad\forall~{}~{}n\geq 0}\right\\}\subset\mathcal{H}.$ (3.9)

where $s,R_{0}>0$ and $V_{n}$, $n\geq 1$ are admissible subspaces defined in
3.16. With this groundwork, the probabilistic reconstruction error is defined
through the subsequent theorem:

######  Theorem 3.17.

Suppose $L$ specifies a Hilbert scale according to definition 2.4 such that
the spectrum of $L$ satisfies $\sigma_{j}(L^{-1})\geq
D_{1}j^{-\frac{t}{\gamma}}$, $t,\gamma>0$, and $V_{n}\subset\mathcal{D}(L)$
for any $n\geq 0$. We assume that $\nu$ satisfies Assumption 8 with $p=1$ and
$a<\frac{1}{2}$ such that $\gamma=1/a-2$ and suppose $T_{\nu}$ is a Hilbert-
Schmidt operator satisfying (3.10) with
${\mathcal{R}}(T_{\nu})\subset\mathcal{D}(L^{\gamma})$. Given a sequence of
admissible subspaces
$\mathopen{}\mathclose{{}\left\\{V_{n}}\right\\}_{n=1}^{\infty}$, we assume
that for some $D_{2}>0$ we have

$\max\mathopen{}\mathclose{{}\left\\{\sup_{n\in{\mathbb{N}}}\mathopen{}\mathclose{{}\left\|{(P_{n}T_{\nu}P_{n})^{+}P_{n}L^{-\gamma}P_{n}}}\right\|,\mathopen{}\mathclose{{}\left\|{L^{\gamma}T_{\nu}}}\right\|}\right\\}<D_{2}$
(3.10)

and

$\sup_{n\in{\mathbb{N}}}\mathopen{}\mathclose{{}\left\|{(P_{n}L^{-\gamma}P_{n})^{+}L^{-\gamma}}}\right\|\leq
D_{2}$ (3.11)

If $f^{\dagger}\in\Theta(s,R)$, there exists a constant $C$ depending on
$D_{j},~{}j=1,2$ such that for $m\geq n$ the ML estimator $f_{n,m}$ defined by
identity (3.8) satisfies the following bound with probability greater than
$1-\eta$ for $0<\eta<1$,

$\mathopen{}\mathclose{{}\left\|{f_{n,m}-f^{\dagger}}}\right\|_{\mathcal{H}}\leq
C\mathopen{}\mathclose{{}\left[R_{0}n^{-s}+\log\mathopen{}\mathclose{{}\left(\frac{8}{\eta}}\right)\delta\mathopen{}\mathclose{{}\left(\frac{n^{t}}{m}+\sqrt{\frac{n^{t+1}}{m}}}\right)}\right]$
(3.12)

provided that

$\log\mathopen{}\mathclose{{}\left(\frac{8}{\eta}}\right)\leq\frac{\sqrt{m}}{12}n^{-t}.$

###### Proof.

Below, we describe the result [43, Thm. 1.2] and prove that conditions therein
are satisfied by the our assumptions. In the following, we employ the
convention $P_{0}=0$. In [43] the author considered following sub-classes of
design measures:

${\mathcal{P}}_{\text{proj}}^{>}(t,C)=\mathopen{}\mathclose{{}\left\\{\nu\in\mathbb{P}(X):\sigma_{\min}(P_{n}T_{\nu}P_{n})\geq
Cn^{-t}\;\text{for all}\;n\in\mathbb{N}}\right\\}$

and

${\mathcal{P}}_{\text{proj}}^{\times}(C)=\mathopen{}\mathclose{{}\left\\{\nu\in\mathbb{P}(X):\mathopen{}\mathclose{{}\left\|{(P_{n}T_{\nu}P_{n})^{+}T_{\nu}(I-P_{n})}}\right\|_{\mathcal{H}}\leq
C\;\text{for all}\;n\in\mathbb{N}}\right\\}$

where $\sigma_{\min}(P_{n}T_{\nu}P_{n})$ refers to the smallest eigenvalue of
the finite-dimensional operator $P_{n}T_{\nu}P_{n}:V_{n}\to V_{n}$. Now [43,
Thm. 1.2] states that if $f^{\dagger}\in\Theta(s,R)$ and
$\nu\in{\mathcal{P}}_{\text{proj}}^{>}(t,C_{1})\cap{\mathcal{P}}_{\text{proj}}^{\times}(C_{2})$,
then inequality (3.12) holds with probability

$\log\mathopen{}\mathclose{{}\left(\frac{8}{\eta}}\right)\leq\frac{\sqrt{m}}{12}\sigma_{\min}(P_{n}T_{\nu}P_{n})$

for a constant $C>0$ depending on $C_{1}$ and $C_{2}$.

Let us now deduce that our assumptions imply
$\nu\in{\mathcal{P}}_{\text{proj}}^{>}(t,C)\cap{\mathcal{P}}_{\text{proj}}^{\times}(C)$
for some $C>0$ depending on $D_{1}$ and $D_{2}$. Applying the operator concave
function $t\mapsto t^{1/2a}$ to inequality (2.10) in assumption 8 respects the
partial ordering, and we obtain that

$\langle{L^{-1}T_{\nu}L^{-1}u},{u}\rangle_{\mathcal{H}}\geq\langle{L^{-1/a}u},{u}\rangle_{\mathcal{H}}.$
(3.13)

Suppose $u:=LP_{n}v\in\mathcal{H}$ and notice that by assumption
$P_{n}v\in\mathcal{D}(L)$. Then we observe that

$\langle{T_{\nu}P_{n}v},{P_{n}v}\rangle_{\mathcal{H}}\geq\langle{L^{-\gamma}P_{n}v},{P_{n}v}\rangle_{\mathcal{H}}.$

Since $\gamma>0$, the operator $L^{-\gamma}$ is bounded, and by the Weyl
Monotonicity Theorem [7, Cor. III.2.3] we see that

$\sigma_{n}\mathopen{}\mathclose{{}\left(P_{n}T_{\nu}P_{n}}\right)\geq\sigma_{n}\mathopen{}\mathclose{{}\left(P_{n}L^{-\gamma}P_{n}}\right).$
(3.14)

In particular, by our assumption on the spectral decay of $L^{-1}$ we have
that

$\sigma_{\min}\mathopen{}\mathclose{{}\left(P_{n}T_{\nu}P_{n}}\right)\geq\sigma_{\min}\mathopen{}\mathclose{{}\left(P_{n}L^{-\gamma}P_{n}}\right)\geq
D_{1}n^{-t}.$ (3.15)

Next, we observe that

$\displaystyle\mathopen{}\mathclose{{}\left\|{(P_{n}T_{\nu}P_{n})^{+}T_{\nu}(I-P_{n})}}\right\|_{\mathcal{H}}$
$\displaystyle=$
$\displaystyle\mathopen{}\mathclose{{}\left\|{(P_{n}T_{\nu}P_{n})^{+}P_{n}L^{-\gamma}P_{n}(P_{n}L^{-\gamma}P_{n})^{+}L^{-\gamma}L^{\gamma}T_{\nu}(I-P_{n})}}\right\|_{\mathcal{H}}$
(3.16) $\displaystyle\leq$
$\displaystyle\mathopen{}\mathclose{{}\left\|{(P_{n}T_{\nu}P_{n})^{+}P_{n}L^{-\gamma}P_{n}}}\right\|_{\mathcal{H}}\mathopen{}\mathclose{{}\left\|{(P_{n}L^{-\gamma}P_{n})^{+}L^{-\gamma}}}\right\|_{\mathcal{H}}\mathopen{}\mathclose{{}\left\|{L^{\gamma}T_{\nu}}}\right\|_{\mathcal{H}}.$

Each of the three terms on the right-hand side of (3.16) is bounded due to
(3.10) and (3.11). In consequence, we obtain the result applying [43, Thm.
1.2]. ∎

We introduce a non-linear truncated estimator
$T_{\widehat{R}}:\mathcal{H}\to\mathcal{H}$ such that

$T_{\widehat{R}}(f)=\begin{cases}f,&\text{if
}\mathopen{}\mathclose{{}\left\|{f}}\right\|_{\mathcal{H}}\leq{\widehat{R}}\\\
0,&\text{otherwise}\end{cases}$

where a fixed non-negative value ${\widehat{R}}\geq 0$ is involved. We
formulate a non-linear estimator $g_{n,m}^{\widehat{R}}$ as follows:

$g_{n,m}^{\widehat{R}}=T_{\widehat{R}}(f_{n,m})$ (3.17)

where the selection of $\widehat{R}$ depends upon the parameters outlined
below. Now, we establish the upper convergence rate in $L_{p}$ for the
estimator $g_{n,m}^{\widehat{R}}$.

######  Theorem 3.18.

Suppose $L$ specifies a Hilbert scale according to definition 2.4 such that
the spectrum of $L$ satisfies $\sigma_{j}(L^{-1})\geq
D_{1}j^{-\frac{t}{\gamma}}$, $t,\gamma>0$, and $V_{n}\subset\mathcal{D}(L)$
for any $n\geq 0$. Let ${\mathcal{P}}^{\prime}$ be the set of probability
measures $\nu$ that

  1. 1.

satisfy Assumption 8 with $p=1$ and $a<\frac{1}{2}$ with $\gamma=1/a-2$,

  2. 2.

for which $T_{\nu}$ is a Hilbert-Schmidt operator satisfying (3.10) and (3.11)
with some $D_{2}>0$, and

  3. 3.

${\mathcal{R}}(T_{\nu})\subset\mathcal{D}(L^{-\gamma})$.

Let $2s-t+1>0$ and set $g_{n,m}^{\widehat{R}}$ as the ML estimator defined by
identity (3.17). For the sequence of solutions corresponding to the parameter
choice

$n=\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{-\frac{2}{2s+t+1}}$
(3.18)

and $\widehat{R}=\widehat{R}(n,\Sigma)$ given by

$\widehat{R}=C\mathopen{}\mathclose{{}\left(\frac{\Sigma}{\sigma_{\min}(P_{n}T_{\nu}P_{n})}+1}\right)$
(3.19)

with suitably large constant $C$ depending on $D_{2}$ and $R$, the sequence

$a_{m,R,\Sigma}=R_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2s}{2s+t+1}}$

is an upper rate of convergence in $L^{p}$ for all $p>0$ for the
reconstruction error over the family of
$\mathcal{M}(\Theta(s,R_{0}),{\mathcal{P}}^{\prime})$.

Let us also note that the paper [43] includes a proof of minimax optimality
regarding the above rates obtained using the class of models specified by
${\mathcal{P}}_{\text{proj}}^{>}$, ${\mathcal{P}}_{\text{proj}}^{\times}$ and
a set ${\mathcal{P}}_{\text{proj}}^{<}$ corresponding to an upper bound to the
spectral decay. However, extending the proof of lower bounds to the Hilbert
scale setting is beyond the scope of this paper.

###  3.4 Convex regularization

Convex variational regularization techniques are of interest in various
inverse problems, where certain features of the unknown, such as edges in
imaging or sparsity in signal processing, needs to be emphasized. In this
section, we deviate briefly from the earlier setup and assume that the forward
linear operator $A$ maps a separable Banach space ${\mathcal{B}}$ to a RKHS
${\mathcal{H}}^{\prime}$.

Below, we discuss the characteristics of regulated solutions
$f_{\mathbf{z},\lambda}$ obtained through the variational problem

$f_{\mathbf{z},\lambda}=\mathop{\rm
arg\,min}\limits\limits_{f\in\mathcal{H}}\mathopen{}\mathclose{{}\left\\{\frac{1}{2m}\sum\limits_{i=1}^{m}\mathopen{}\mathclose{{}\left(A(f)(x_{i})-y_{i}}\right)^{2}+\lambda
G(f)}\right\\}$ (3.20)

Here,
$G:{\mathcal{B}}\to\mathbb{R}\cup\mathopen{}\mathclose{{}\left\\{\infty}\right\\}$
represents a convex functional. Let us point out that penalties such as

$G(f)=\frac{1}{p}\mathopen{}\mathclose{{}\left\|{f}}\right\|_{{\mathcal{B}}}^{p}$
(3.21)

have been studied in statistical learning (see e.g. [57, 76]). However, in the
inverse problems context much less is known and to the best of our knowledge
their study with random design is so far limited to [14, 5].

The first article [14] leverages ideas from deterministic convex
regularization theory, where Bregman distances [6] have become popular during
last decade. While those techniques extend to more general penalties, here, we
restrict to the $p$-homogeneous functional (3.21). Due to strict convexity for
$p\geq 1$, the regularized solution $f_{\mathbf{z},\lambda}$ is unique. In
addition, the mappings

$(f,{\bf
z})\mapsto\frac{1}{2m}\sum\limits_{i=1}^{m}\mathopen{}\mathclose{{}\left(A(f)(x_{i})-y_{i}}\right)^{2}+\lambda
G(f)\quad\text{and}\quad{\mathbf{z}}\mapsto f_{\mathbf{z},\lambda}$ (3.22)

are continuous and measurable, respectively (see comment in [14, Sec. 4]).
Therefore, the learning problem in (3.20) becomes well-posed. Notice that [14]
assumes more general range (Banach space) of operator $A$ with similar
continuous pointwise evaluation functionals. As we limit to RKHS-valued range,
this simplifies the assumptions needed in [14].

The optimality criterion associated with (3.20) is given by

$B_{\mathbf{x}}^{*}(B_{\mathbf{x}}f_{\mathbf{z},\lambda}-\mathbf{y})+\lambda
r_{{\bf z},\lambda}=0$ (3.23)

for $r_{{\bf z},\lambda}\in\partial G(f_{\mathbf{z},\lambda})$, where the dual
operator $B_{\mathbf{x}}^{*}$ maps from the RKHS ${\mathcal{H}}^{\prime}$ to
the Banach dual ${\mathcal{B}}^{*}$ and $\partial G$ denotes the
subdifferential set

$\partial G(f)=\\{r\in{\mathcal{B}}^{*}\;|\;G(f)-G(\tilde{f})\leq\langle
r,f-\tilde{f}\rangle_{{\mathcal{B}}^{*}\times{\mathcal{B}}}\;\text{for
all}\;\tilde{f}\in{\mathcal{B}}\\}.$

For $r_{f}\in\partial G(f)$ and $r_{\tilde{f}}\in\partial G(f)$ we define
symmetric Bregman distance between $f$ and $\tilde{f}$ as

$D_{G}^{r_{f},r_{\tilde{f}}}(f,\tilde{f})=\langle
r_{f}-r_{\tilde{f}},f-\tilde{f}\rangle_{{\mathcal{B}}^{*}\times{\mathcal{B}}}$

Following [14] we restrict to the case $1<p<2$. In consequence, the
subdifferential elements $r_{f}\in\partial G(f)$ are unique and in what
follows we drop the dependence of Bregman distance on the subgradients and
simplify notation writing $D_{G}(f,\tilde{f})$.

The crux of analysing the non-quadratic objective functional in (3.20) is that
it does not give rise to similar explicit formula for the minimizer
$f_{\mathbf{z},\lambda}$ as e.g. in the Tikhonov case in (3.3). However, as is
shown in [14], the Bregman-based reconstruction error can be separated into
bias-variance type of decomposition. The bias is represented by the term

${\mathcal{R}}(\beta,{\mathbf{z}};f^{\dagger})=\inf_{{\bf
w}\in\mathbb{R}^{m}}\mathopen{}\mathclose{{}\left[G^{\ast}(R_{0}-B_{{\mathbf{z}}}^{*}{\bf
w})+\frac{\beta}{2m}\mathopen{}\mathclose{{}\left\|{{\bf
w}}}\right\|^{2}}\right],$ (3.24)

where $G^{\ast}$ is the convex conjugate of $G$. The functional
${\mathcal{R}}$ is closely related to the distance function in Definition
2.10. However, notice carefully that unlike the distance function, the
functional ${\mathcal{R}}$ is dependent on the random sampling. Therefore,
conditions for the source set $\Theta$ and set of design measures
${\mathcal{P}}^{\prime}$ become intertwined without further assumptions on the
structure of $G$. Moreover, the error variance is represented by
$G^{\ast}(B_{{\mathbf{z}}}^{*}{\bf\epsilon})$. The general concentration
result is given as follows:

######  Theorem 3.19 ([14, Thm. 4.11]).

Let $G$ be defined in (3.21) and let ${\mathcal{M}}_{\theta}$,
$\theta=(\Sigma,R_{0})$, consist of all data-generating distributions induced
by conditions

$G(f^{\dagger})\leq R_{0}$

and

$\mathbb{E}{\mathcal{R}}(\beta,{\mathbf{z}};f^{\dagger})\leq
D_{1}\beta+D_{2}m^{-Q}$

for some $R_{0},Q,D_{1},D_{2}>0$. Moreover, let us assume that

$\mathbb{E}G^{\ast}(B_{{\mathbf{z}}}^{*}{\bf\epsilon})\leq
D_{3}m^{-\frac{q}{2}}$

for some $D_{3}>0$, where $q$ is the Hölder conjugate of $1<p<2$. Then it
follows that for the parameter choice rule

$\lambda=\mathopen{}\mathclose{{}\left(\frac{D_{3}^{\frac{2}{q}}R_{0}^{-\frac{q+2}{q}}}{D_{1}}}\right)^{\frac{1}{3}}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{3}}$
(3.25)

we have

$\sup_{(\Sigma,R_{0})\in\mathbb{R}^{2}}\limsup_{m\to\infty}\sup_{\rho\in{\mathcal{M}}}\frac{\mathbb{E}D_{G}(f_{\mathbf{z},\lambda},f^{\dagger})}{a_{m,R_{0},\Sigma,p}}\leq
C,$ (3.26)

where the concentration rate is given by

$a_{m,R_{0},\Sigma,p}=\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{3}}$
(3.27)

and the constant $C$ depends on $D_{j}$, $j=1,2,3$, and $p$.

For a more concrete result, let us consider wavelet decompositions in
$\mathbb{R}^{d}$. Let
${\mathcal{B}}=B_{pp}^{s}(\mathbb{R}^{d}):=B_{p}^{s}(\mathbb{R}^{d})$ be a
Besov space [29] and

$G(f)=\frac{1}{p}\mathopen{}\mathclose{{}\left\|{f}}\right\|_{B_{p}^{s}}^{p}:=\frac{1}{p}\sum_{\alpha=1}^{\infty}c_{\alpha,p,s,d}|\langle
f,\psi_{\alpha}\rangle|^{p}$ (3.28)

for some $1<p<2$, where

$c_{\alpha,p,s,d}=2^{|\alpha|d\mathopen{}\mathclose{{}\left(p\mathopen{}\mathclose{{}\left(\frac{s}{d}+\frac{1}{2}}\right)-1}\right)}.$
(3.29)

Here, $\psi_{\alpha}:\mathbb{R}^{d}\to\mathbb{R}$, with $\alpha=1,...,\infty$,
are suitably regular functions that form an orthonormal wavelet basis for
$\mathscr{L}^{2}(\mathbb{R}^{d})$ with global indexing $\alpha$. The notation
$|\alpha|$ is used to denote the scale of the wavelet basis associated with
the index $\alpha$. Rephrasing [14, Cor. 5.6] in terms of the upper rate of
convergence yields the following statement.

######  Theorem 3.20.

Let $G$ be defined by (3.21) and let the model set ${\mathcal{M}}$ be induced
by conditions $G(f)\leq R_{0}$ and

$\partial
G(f)=\\{B^{*}_{\nu}w\\}\;\text{for}\;\mathopen{}\mathclose{{}\left\|{w}}\right\|\leq
R_{0}.$

Moreover, assume

$\sum_{\alpha=1}^{\infty}c_{\alpha,q,-s,d}\mathopen{}\mathclose{{}\left\|{A\psi_{\alpha}}}\right\|_{{\mathcal{H}}^{\prime}}^{q}<\infty,$

where coefficients $c_{\alpha,q,-s,d}$ are defined by (3.29). Then for
parameter choice rule

$\lambda=R_{0}^{\frac{2}{3p}-\frac{1}{3}}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{3}}$
(3.30)

we have

$\sup_{(\Sigma,R_{0})\in\mathbb{R}^{2}}\limsup_{m\to\infty}\sup_{\rho\in{\mathcal{M}}}\frac{\mathbb{E}D_{G}(f_{\mathbf{z},\lambda},f^{\dagger})}{a_{m,R_{0},\Sigma,p}}\leq
C$ (3.31)

with the concentration rate

$a_{m,R_{0},\Sigma,p}=R_{0}^{\frac{2+5p}{3p}}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{3}}.$
(3.32)

######  Remark 3.21.

As mentioned above, convex penalties in inverse problems, such as Lasso
regularization, are often motivated by sparsity-type assumption regarding the
unknown. While Theorems 3.19 and 3.20 provide first attempt to formalize
theory of convex penalties in statistical inverse learning with random design,
the highly interesting case of $p=1$ remains open and is not easily answered
by the Bregman distance based techniques. Moreover, the Bregman distance for
$1$-homogeneous functional (consider e.g. the $\ell^{1}$-norm) is typically
weak as the subgradient set is not a singleton. We point out that variational
source conditions have been successful in proving concentration rates in the
norm for deterministic design in statistical inversion [84] and it would be
tempting to use similar ideas for the random design problem.

###  3.5 Parameter choice rules and adaptivity

We recall that while tuning the regularization parameter is essential for
regularization to work well, an a-priori choice of the regularization
parameter is in general not feasible in statistical problems since the choice
necessarily depends on unknown structural properties such as smoothness of
$f^{\dagger}$ or eigenvalue decay, see e.g. equations (3.5), (3.18), Table 1
or equations (3.25), (3.30). Parameter choice rules, also known as _model
selection_ or _hyperparameter selection_ , are crucial in the context of
inverse problems and play a vital role in finding stable, reliable, and
interpretable solutions, that depend on the data only and are free from
additional unknown a-priori assumptions.

The aim is to construct estimators, i.e. to find a sequence of regularization
parameters $(\lambda_{\mathbf{z}})_{m}$, without knowledge of a-priori
assumptions, but depending on the data ${\mathbf{z}}$, such that the sequence
of regularized solutions $(f_{{\mathbf{z}},\lambda_{\mathbf{z}}})_{m}$ is
minimax optimal. Such a sequence of estimators is called _minimax optimal
adaptive_.

We give a very brief overview about some classical model selection criteria
used in inverse problems.

Unbiased risk estimation (URE). A very natural idea is to estimate the unknown
risk by a functional $U({\mathbf{z}},\lambda)$ of the observations, and then
to minimize this estimator of the risk. Ideally, $U$ is an _unbiased_
estimator of the risk, meaning that

${\mathcal{E}}(f_{{\mathbf{z}},\lambda})=\mathbb{E}[U({\mathbf{z}},\lambda)]\;,\quad\forall\;\;\lambda>0\;.$

The principle of URE suggests to minimize over $\lambda\in\Lambda$, where
$\Lambda$ is a well-chosen finite grid of possible values for the
regularization parameter:

$\lambda_{\mathbf{z}}:=\mathop{\rm
arg\,min}\limits_{\lambda\in\Lambda}U({\mathbf{z}},\lambda)\;.$

For inverse problems, this method was studied in [21], where exact oracle
inequalities for the mean square risk were obtained.

Hold-out and (generalized) cross-validation (GCV). Cross-validation is a
classical approach to select hyperparameters that is based on _data-
splitting_. The idea behind GCV is to use the data to both train and validate
the model, but with a focus on finding the optimal hyperparameters. The
process involves iteratively fitting the model with different hyperparameter
values and evaluating the model’s performance using a specific criterion. The
simplist version is the _hold-out method_ , where the available data are split
into a training set ${\mathbf{z}}^{t}$ of size $m_{1}$ and validation set
${\mathbf{z}}^{v}$ of size $m_{2}$, with $m=m_{1}+m_{2}$. The estimator
$f_{{\mathbf{z}}^{t},\lambda}$ is defined by means of the training data and is
evaluated using the empirical risk with respect to the validation data:

$\lambda_{{\mathbf{z}}}:=\mathop{\rm
arg\,min}\limits_{\lambda\in\Lambda}\frac{1}{m_{2}}\sum_{j=1}^{m_{2}}(f_{{\mathbf{z}}^{t},\lambda}(x_{j}^{v})-y_{j}^{v})^{2}\;.$

This method works well in practice and [19] show minimax optimal rates for the
$\mathscr{L}^{2}(X,\nu)$ error under polynomial eigendecay and a Hölder source
condition.

Lepski’s balancing principle. The _Lepski principle_ , a.k.a. the
_Goldenshluger-Lepski method_ was introduced in the seminal paper [52] and
does not use any data splitting. It is merely based on comparing estimators
for different levels of regularization. Surprisingly, the monotonicity of the
approximation error suffices to provide the knowledge to choose the
regularization parameter adaptively, provided an empirical approximation of
the sample error ${\mathcal{S}}_{\mathbf{z}}$ is available. More precisely,
let $\Lambda$ be a finite grid of possible values for $\lambda$. Set

$J_{{\mathbf{z}}}(\Lambda)=\mathopen{}\mathclose{{}\left\\{\lambda\in\Lambda:\mathopen{}\mathclose{{}\left\|(B_{\mathbf{x}}+\lambda^{\prime})^{\frac{1}{2}}(f_{{\mathbf{z}},\lambda}-f_{{\mathbf{z}},\lambda^{\prime}})}\right\|_{{\mathcal{H}}}\leq
c\;\sqrt{\lambda^{\prime}}\;{\mathcal{S}}_{\mathbf{z}}(\lambda^{\prime}),\quad\forall\lambda^{\prime}\in\Lambda,\,\lambda^{\prime}\leq\lambda}\right\\}\;.$

The parameter choice based on Lepski’s principle is given by

${\lambda}_{{\mathbf{z}}}:=\max_{\lambda\in\Lambda}J_{{\mathbf{z}}}(\Lambda)\;.$

In the context of learning theory, this method was analyzed in [31], [54] and
was finally improved by [12]. This parameter choice rule is minimax optimal
adaptive under various types of source conditions and eigenvalue decay.

Discrepancy principle and early stopping rules. The above mentioned criteria
require the computation of all estimators to be compared against each other.
For problems with high dimensions, this can incur a computationally
prohibitive cost. An alternative approach is the use of early stopping rules
based on the discrepancy principle. These rules terminate the procedure at an
iteration $t_{\mathbf{z}}$ based solely on iterates with index $t\leq
t_{\mathbf{z}}$ and potentially other quantities computed up to that point.
Because they require the calculation of significantly fewer iterates, early
stopping rules offer the potential to achieve both computational and
statistical efficiency simultaneously. The stopping criterion is defined based
on a validation performance: If the performance does not improve for a certain
number of consecutive evaluations, the training process is halted. A popular
validation criterion, for instance, is based on the smoothed empirical
residuals

$\hat{R}_{t,s}:=||(B_{\mathbf{x}}B_{\mathbf{x}}^{*})^{s}({\mathbf{y}}-B_{\mathbf{x}}f_{{\mathbf{z}},t})||_{2}\;,$

where $s\geq 0$ and $f_{{\mathbf{z}},t}$ is an iteratively defined estimator,
as such as the last iterate of a gradient based method or truncated SVD, see
e.g. [86, 74, 9, 10]. The iterative process is stopped at the smallest
$t_{\mathbf{z}}$, such that

$\hat{R}_{t_{\mathbf{z}},s}\leq\tau\;,$

for some pre-specified $\tau>0$. Provided that the inverse problem is
moderately ill-posed, it is known that this criterion leads to optimal
adaptive bounds for $f_{{\mathbf{z}},t_{{\mathbf{z}}}}$ for moderate
smoothing, i.e. $s<b/2$, over Sobolev-type ellipsoids [74].

We finally mention some recent results related to establishing data-driven
approaches for finding the optimal regularization parameter, highlighting
particularly the current need of a better understanding of such methods in
modern machine learning: [2] focus on the case where the regularization
functional is not given a priori, but learned from data, [23] study a
statistical learning approach, based on empirical risk minimization, and [30]
propose a novel method based on supervised learning to approximate the high-
dimensional function, mapping noisy data into a good approximation to the
optimal Tikhonov regularization parameter.

##  4 Non-linear problems

In this section, we discuss the statistical learning problem (1.1) when $A$ is
a nonlinear operator. However, in contrast to Section 3, we only focus on
Tikhonov regularization. We present classical error bounds for the
reconstruction error in Section 4.1 as well as bounds for Tikhonov
regularization in Hilbert scales in Section 4.2.

For analyzing the non-linear inverse problem (1.1) we require the forward
operator to be Fréchet differentiable (an assumption that will even be further
relaxed below).

######  Assumption 9 (Non-linearity of the operator I).

We impose the following assumptions:

  1. (i)

The domain $\mathcal{D}(A)$ is convex with nonempty interior.

  2. (ii)

The operator $A$ is Fréchet differentiable.

  3. (iii)

The Fréchet derivative $A^{\prime}(f)$ of $A$ at $f$ is bounded in a ball
$\mathcal{B}_{d}(f^{\dagger})$ of radius
$d:=4\mathopen{}\mathclose{{}\left\|{f^{\dagger}-\bar{f}}}\right\|_{\mathcal{H}}$,
i.e., there exists $\mu_{1}<\infty$ such that

$\mathopen{}\mathclose{{}\left\|{A^{\prime}(f)}}\right\|_{\mathcal{H}\to\mathcal{H}^{\prime}}\leq\mu_{1}\qquad\forall~{}f\in\mathcal{B}_{d}(f^{\dagger})\cap\mathcal{D}(A)\subset\mathcal{H}.~{}$

Recall from Example 2.7, the definition of the covariance operator

$T_{\nu}:=B_{\nu}^{*}B_{\nu},\quad B_{\nu}=I_{\nu}A^{\prime}(f^{\dagger}).$

###  4.1 Classical Error Bounds

The Tikhonov regularization scheme for this problem is again given by an
optimization problem

$f_{\mathbf{z},\lambda}=\mathop{\rm
arg\,min}\limits_{f\in\mathcal{D}(A)\subset\mathcal{H}}\mathopen{}\mathclose{{}\left\\{\frac{1}{m}\sum\limits_{i=1}^{m}\mathopen{}\mathclose{{}\left([A(f)](x_{i})-y_{i}}\right)^{2}+\lambda\mathopen{}\mathclose{{}\left\|{f-\bar{f}}}\right\|_{\mathcal{H}}^{2}}\right\\}.$
(4.1)

Here, $\bar{f}\in{\mathcal{H}}$ denotes some initial guess of the solution
$f^{\dagger}$, which offers the possibility to incorporate a-priori
information. Note, however, due to the non-linearity of $A$, the solution to
the minimization problem is not explicitly given and does not need to be
unique.

To establish our error bounds we need additional assumptions on the non-
linearity of $A$.

######  Assumption 10 (Non-linearity of the operator II).

We assume:

  1. (i)

The operator
$A:\mathcal{D}(A)\subset\mathcal{H}\to\mathcal{H}^{\prime}\hookrightarrow\mathscr{L}^{2}(X,\nu)$
is weakly sequentially closed and one-to-one.

  2. (ii)

There exists $\mu_{2}\geq 0$ such that for all
$f\in\mathcal{B}_{d}(f^{\dagger})\cap\mathcal{D}(A)\subset\mathcal{H}$ we have

$\mathopen{}\mathclose{{}\left\|{A(f)-A(f^{\dagger})-A^{\prime}(f^{\dagger})(f-f^{\dagger})}}\right\|_{\nu}\leq\frac{\mu_{2}}{2}\mathopen{}\mathclose{{}\left\|{f-f^{\dagger}}}\right\|_{\mathcal{H}}^{2}.$

The first assumption ensures that there exists a global minimizer of the
Tikhonov functional (4.1). If $\mathcal{D}(A)$ is weakly closed (closed and
convex) and $A$ is weakly continuous, then weak sequential closedness is
guaranteed. On the other hand, the second assumption helps to control the
remainder term of the linearized problem. The second condition also holds true
when the Fréchet derivative is Lipschitz in the operator norm (see [36, Chapt.
10]).

As in the linear case, we impose a source condition for $f^{\dagger}$ to
facilitate convergence.

######  Assumption 11 (Source condition).

Let $\bar{f}\in{\mathcal{H}}$. We assume that
$f^{\dagger}-\bar{f}\in\Omega(r,R_{0},T_{\nu})$, for $r\in[\frac{1}{2},1]$ and
where the source set is defined in Assumption 6.

Under the source condition $f^{\dagger}-\bar{f}=T_{\nu}^{r}u$ for
$r\in[\frac{1}{2},1]$, we have that
$f^{\dagger}-\bar{f}=T_{\nu}^{\frac{1}{2}}w$ for $w=T_{\nu}^{r-\frac{1}{2}}u$
satisfying

$\mathopen{}\mathclose{{}\left\|{w}}\right\|_{\mathcal{H}_{1}}\leq\mathopen{}\mathclose{{}\left\|{T_{\nu}}}\right\|^{r-\frac{1}{2}}\cdot\mathopen{}\mathclose{{}\left\|{u}}\right\|_{\mathcal{H}}\leq\kappa^{2r-1}R_{0}\leq\bar{\kappa}R_{0}=:\mu_{3}\;,$

with $\bar{\kappa}=\max\\{1,\kappa\\}$. We additionally assume that

$R_{0}<\frac{1}{2\mu_{2}\bar{\kappa}},$ (4.2)

implying that $2\mu_{2}\mu_{3}<1$. This assumption is a smallness condition
that imposes a constraint between $\mu_{3}$ and the non-linearity as measured
by the parameter $\mu_{2}$ in Assumption 9 (iii). It ensures that the initial
guess is close enough to the true solution and the residual error on
linearizing the nonlinear operator $A$ at the solution is sufficiently small.

We need an additional restriction to achieve the lower bounds:

######  Assumption 12.

Let us define $T=A^{\prime}(f)^{*}I_{\nu}^{*}I_{\nu}A^{\prime}(f)$ and
$\overline{T}=A^{\prime}(\bar{f})^{*}I_{\nu}^{*}I_{\nu}A^{\prime}(\bar{f})$.
There exists a family of bounded linear operators
$R_{f}:\mathcal{H}\to\mathcal{H}$ and $\zeta>0$ such that

$T^{r}=R_{f}\overline{T}^{r}\text{ and
}\mathopen{}\mathclose{{}\left\|{R_{f}-I}}\right\|\leq\zeta\mathopen{}\mathclose{{}\left\|{f-\bar{f}}}\right\|_{\mathcal{H}},$

where $f$ belongs to the sufficiently large ball and $\bar{f}$ is the initial
guess.

######  Theorem 4.1 (Rastogi et al. [67]).

Suppose Assumptions 1–3, 5, 9–11 are satisfied. Assume further that

$8\kappa^{2}\max\mathopen{}\mathclose{{}\left\\{1,\frac{\mu_{1}(M+\Sigma)}{\kappa
d}}\right\\}\log\mathopen{}\mathclose{{}\left(\frac{4}{\eta}}\right)\leq\sqrt{m}\lambda.$
(4.3)

Let (4.2) hold true with $r\in[\frac{1}{2},1]$. Then, for all $\eta\in(0,1]$,
for a regularized least-squares estimator $f_{\mathbf{z},\lambda}$ in (4.1),
the following upper bound holds with probability at least $1-\eta$:

$\|f_{\mathbf{z},\lambda}-f^{\dagger}\|_{\mathcal{H}}\;\leq\;C\mathopen{}\mathclose{{}\left(R_{0}\lambda^{r}+\frac{\Sigma}{\sqrt{m}\lambda^{\frac{b+1}{2}}}+\frac{M}{m\lambda}}\right)\log\mathopen{}\mathclose{{}\left(\frac{4}{\eta}}\right),$

where $C$ depends on the parameters
$\kappa,~{}\mu_{1},~{}\mu_{2},~{}\mu_{3},~{}b,~{}\beta$.

Again, as in the linear case, the choice of the regularization parameter is
crucial. By balancing the first and second term, we obtain the parameter
choice

$\lambda_{m,\theta}=\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2}{2r+b+1}}\;,$
(4.4)

with $\theta=(R_{0},\Sigma)$. A short calculation then shows that the third
term is of lower order and that (4.3) is satisfied for $m$ sufficiently large.
We summarize our findings in the next corollary.

######  Corollary 4.2.

Under the same assumptions of Theorem 4.1 with the parameter choice (4.4) we
have with probability at least $1-\eta$,

$\|f_{\mathbf{z},\lambda}-f^{\dagger}\|_{\mathcal{H}}\leq
CR_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2r}{2r+b+1}}\log\mathopen{}\mathclose{{}\left(\frac{4}{\eta}}\right)$

where $C$ depends on the parameters
$\kappa,~{}\mu_{1},~{}\mu_{2},~{}\mu_{3},~{}b,~{}\beta$ and provided that $m$
is sufficiently large.

Finally, proper integration using Corollary C.2 in [13] gives the upper rate
of convergence.

######  Theorem 4.3.

Suppose the assumptions of Corollary 4.2 hold true. Then, the sequence

$\displaystyle
a_{m,\theta}=R_{0}\mathopen{}\mathclose{{}\left(\frac{\Sigma}{R_{0}\sqrt{m}}}\right)^{\frac{2r}{2r+b+1}}$
(4.5)

is an upper rate of convergence in $L^{p}$ for all $p>0$, for the
reconstruction error, over the family of models
${\mathcal{M}}(r,R_{0},{\mathcal{P}}_{\text{poly}}^{<}(\beta,b))$ for the
sequence of solutions (4.1) with regularization parameter (4.4).

######  Theorem 4.4.

Let Assumptions 1–3, 5, 9–12 hold true. Then the sequence
$(a_{m,\theta})_{\theta}$ defined in (4.5) is a weak lower rate of convergence
in $L^{p}$ for all $p>0$ for the reconstruction error, over the family of
models ${\mathcal{M}}(r,R_{0},{\mathcal{P}}_{\text{poly}}^{>}(\alpha,b))$.

Finally, by Definition 2.14, we may conclude:

######  Theorem 4.5.

Suppose the Assumptions of Theorem 4.3 and Theorem 4.4 are satisfied. Then the
sequence of solutions $(f_{\mathbf{z},\lambda})_{m}$ using the regularization
parameters given in (4.4) is weak minimax optimal in $L^{p}$, for all $p>0$,
for the reconstruction error, over the model family
${\mathcal{M}}(r,R_{0},{\mathcal{P}}^{\prime})$, with
${\mathcal{P}}^{\prime}={\mathcal{P}}_{\text{poly}}^{<}(\beta,b)\cap{\mathcal{P}}_{\text{poly}}^{>}(\alpha,b))$.

###  4.2 Error Bounds in Hilbert Scales

Let $L:{\mathcal{D}}(L)\subset{\mathcal{H}}\to{\mathcal{H}}$ be an unbounded
operator, generating a Hilbert scale. The Tikhonov regularization scheme in
Hilbert Scales amounts to minimize the regularized empirical risk

$f_{\mathbf{z},\lambda}=\mathop{\rm
arg\,min}\limits_{f\in\mathcal{D}(A)\cap\mathcal{D}(L)}\mathopen{}\mathclose{{}\left\\{\frac{1}{m}\sum\limits_{i=1}^{m}\mathopen{}\mathclose{{}\left([A(f)](x_{i})-y_{i}}\right)^{2}+\lambda\mathopen{}\mathclose{{}\left\|{L(f-\bar{f})}}\right\|_{\mathcal{H}}^{2}}\right\\},$
(4.6)

enforcing smoothness and with initial guess
$\bar{f}\in{\mathcal{D}}(A)\cap{\mathcal{D}}(L)$. As above, due to the non-
linearity of $A$, the solution to the minimization problem is not explicitly
given and does not need to be unique. To establish error bounds also in this
case, the nonlinear structure of the operator needs to be further specified.

######  Assumption 13 (Non-linearity of the operator III).

We suppose that:

  1. (i)

The operator $A:\mathcal{D}(A)\cap\mathcal{D}(L)\to\mathcal{H}^{\prime}$ is
weakly sequentially closed.

  2. (ii)

Link condition: There exist constants $p\geq 1$ and $0<c_{1}\leq c_{2}$ such
that for all $g\in\mathcal{H}$,

$c_{1}\|g\|_{{-(p-1)}}\leq\|A^{\prime}(f^{\dagger})g\|_{\nu}\leq
c_{2}\|g\|_{{-(p-1)}}.$

  3. (iii)

The Fréchet derivative $A^{\prime}(f)$ is Lipschitz continuous for all
$f\in\mathcal{D}(A)\cap\mathcal{D}(L)~{}$, i.e., there exists a constant
$\mu_{4}$ such that

$\|A^{\prime}(f^{\dagger})-A^{\prime}(f)\|_{-(p-1),\nu}\leq\mu_{4}\|f^{\dagger}-f\|_{\mathcal{H}}\leq\frac{c_{1}^{2}}{2c_{2}},$

where the norm $\mathopen{}\mathclose{{}\left\|{\cdot}}\right\|_{-(p-1),\nu}$
is defined for the operators $\mathcal{H}_{-(p-1)}$ to
$\mathscr{L}^{2}(X,\nu)$.

The link condition, described in Assumption 13 (ii) involves the interplay
between the operator $L^{-1}$ and the Fréchet derivative of $A$.

We incorporate a-priori smoothness in terms of a more general source
condition.

######  Assumption 14 (Source condition).

Let $\bar{f}\in{\mathcal{D}}(L)$, $s\in[1,p+1]$ for some benchmark smoothness
$p\geq 1$. We assume that $f^{\dagger}-\bar{f}\in\Omega(s,R_{0},L^{-1})$,
where the source set is defined in Assumption 7.

Having now established all assumptions, we are ready to state an upper bound
for the reconstruction error.

######  Theorem 4.6 (Rastogi et al. [65]).

Suppose Assumptions 1–3, 5, 9, 13 and 14 hold true. Then, for the Tikhonov
estimator $f_{\mathbf{z},\lambda}~{}$ in (4.6) with the a-priori choice of the
regularization parameter
$\lambda=\mathopen{}\mathclose{{}\left(\frac{1}{\sqrt{m}}}\right)^{\frac{2p}{p+s-1+bp}}$,
for all $0<\eta<1$, the following error bound holds with the confidence
$1-\eta$:

$\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\leq
C\mathopen{}\mathclose{{}\left(\frac{1}{\sqrt{m}}}\right)^{\frac{s}{p+s-1+bp}}\log^{2}\mathopen{}\mathclose{{}\left(\frac{4}{\eta}}\right).$

We will discuss the optimality of convergence rates later in this section.

Conditional stability estimates. The above results require the
differentiability of the operator $A$. This condition can be relaxed by
imposing _conditional stability_. A conditional stability estimate provides a
tool to characterize the degree of ill-posedness in inverse problems and
measures the sensitivity of a solution with respect to changes in the input or
observational data. This assumption implicitly addresses both the nonlinearity
conditions and solution smoothness within the underlying nonlinear inverse
problem. We will also establish the connection between link condition 13 (ii)
and the conditional stability estimates 15 (iv).

######  Assumption 15.

  1. (i)

The domain $\mathcal{D}(A)$ of $A$ is a convex and closed subset of
$\mathcal{H}$.

  2. (ii)

The operator $A:\mathcal{D}(A)\to\mathcal{H}^{\prime}$ is weak-to-weak
sequentially continuous222i.e., $f_{n}\rightharpoonup\hat{f}\in\mathcal{H}$
with $f_{n}\in\mathcal{D}(A)$, $n\in\mathbb{N}$, and
$\hat{f}\in\mathcal{D}(A)$ implies $A(f_{n})\rightharpoonup
A(\hat{f})\in\mathcal{H}^{\prime}$..

  3. (iii)

The operator $A$ is Lipschitz continuous with a Lipschitz constant
$\ell_{A}<\infty$ in a sufficiently large ball $\mathcal{B}_{d}(f^{\dagger})$,

$\mathopen{}\mathclose{{}\left\|{A(f)-A(\tilde{f})}}\right\|_{\mathcal{H}^{\prime}}\leq\ell_{A}\mathopen{}\mathclose{{}\left\|{f-\tilde{f}}}\right\|_{\mathcal{H}}\qquad\forall
f,\tilde{f}\in\mathcal{B}_{d}(f^{\dagger})\cap\mathcal{D}(A)\subset\mathcal{H},$

  4. (iv)

Conditional stability estimate: There exist constants $p\geq 1$, $\alpha>0$,
$d>0$, $\theta\geq 0$ and
$Q\subset\mathcal{D}_{d}^{\theta}(f^{\dagger})\cap\mathcal{D}(A)$ such that

$\mathopen{}\mathclose{{}\left\|{f-f^{\dagger}}}\right\|_{-(p-1)}\leq\alpha\mathopen{}\mathclose{{}\left\|{A(f)-A(f^{\dagger})}}\right\|_{\nu}$

for all $f\in Q$, where the constant $\alpha$ may depend on $p$, and $Q$.

Conditional stability estimates require additional regularization for
obtaining stable approximate solutions if the validity area $Q$ of such
estimates is not completely known. Conditional stability estimates can (in the
absence of differentiability) be verified by means of global inequalities of
the forward operator $A$. For examples, we refer the reader to [85].

######  Theorem 4.7 (Rastogi et al. [66]).

Suppose Assumptions 1–3, 5, 14 and 15 hold true and
$f_{\mathbf{z},\lambda},f^{\dagger}\in Q$ (for a sufficiently large sample
size $m$) for $p$, $Q$, and $s$ defined in Assumption 15 (iv). Then, for the
Tikhonov estimator $f_{\mathbf{z},\lambda}~{}$ in (4.6) with the a-priori
choice of the regularization parameter
$\lambda=\mathopen{}\mathclose{{}\left(\frac{1}{\sqrt{m}}}\right)^{\frac{2p}{p+s-1+bp}}$,
for all $0<\eta<1$, the following error bound holds with the confidence
$1-\eta$:

$\mathopen{}\mathclose{{}\left\|{f_{\mathbf{z},\lambda}-f^{\dagger}}}\right\|_{\mathcal{H}}\leq
C\mathopen{}\mathclose{{}\left(\frac{1}{\sqrt{m}}}\right)^{\frac{s}{p+s-1+bp}}\log^{2}\mathopen{}\mathclose{{}\left(\frac{4}{\eta}}\right).$

Optimality. Now, we discuss the optimality of the convergence rates
established in this section. To do so, we rely on the lower bounds for the
linear inverse problem presented in [68], which pertain to a similar setting,
albeit with a link condition.

First, we establish a connection between Assumption 15 (iv) and the link
conditions 13 (ii), 8. The conjunction of Assumptions 15 (iv) and

$\mathopen{}\mathclose{{}\left\|{A(f)-A(f^{\dagger})}}\right\|_{\nu}\leq\beta\mathopen{}\mathclose{{}\left\|{f-f^{\dagger}}}\right\|_{\mathcal{H}_{-(p-1)}},\qquad\beta>0$

leads to the following implication:

$\mathopen{}\mathclose{{}\left\|{f-f^{\dagger}}}\right\|_{\mathcal{H}_{-(p-1)}}\asymp\mathopen{}\mathclose{{}\left\|{A(f)-A(f^{\dagger})}}\right\|_{\nu}.$
(4.7)

Here, we write
$\mathopen{}\mathclose{{}\left\|{B_{1}(f)}}\right\|\asymp\mathopen{}\mathclose{{}\left\|{B_{2}(f)}}\right\|$
if there exist $c_{1}\leq c_{2}<\infty$, such that
$c_{1}\mathopen{}\mathclose{{}\left\|{B_{1}(f)}}\right\|\leq\mathopen{}\mathclose{{}\left\|{B_{2}(f)}}\right\|\leq
c_{2}\mathopen{}\mathclose{{}\left\|{B_{1}(f)}}\right\|$.

Let us next introduce an additional assumption that is utilized in both [40]
and [51]:

######  Assumption 16.

For $f,f^{\dagger}\in\mathcal{D}(A)$,

$\mathopen{}\mathclose{{}\left\|{A(f)-A(f^{\dagger})}}\right\|_{\nu}\asymp\mathopen{}\mathclose{{}\left\|{A^{\prime}(f^{\dagger})\mathopen{}\mathclose{{}\left(f-f^{\dagger}}\right)}}\right\|_{\nu}.$

The combination of Assumption 16 and equation (4.7) leads to link condition 13
(ii):

$\mathopen{}\mathclose{{}\left\|{f-f^{\dagger}}}\right\|_{\mathcal{H}_{-(p-1)}}\asymp\mathopen{}\mathclose{{}\left\|{A^{\prime}(f^{\dagger})\mathopen{}\mathclose{{}\left(f-f^{\dagger}}\right)}}\right\|_{\nu}.$

Next, by substituting $f-f^{\dagger}=L^{-1}h$, we derive:

$\mathopen{}\mathclose{{}\left\|{L^{-1}h}}\right\|_{\mathcal{H}_{-(p-1)}}\asymp\mathopen{}\mathclose{{}\left\|{A^{\prime}(f^{\dagger})L^{-1}h}}\right\|_{\nu}$

which further implies:

$\mathopen{}\mathclose{{}\left\|{h}}\right\|_{\mathcal{H}_{-p}}\asymp\mathopen{}\mathclose{{}\left\|{(L^{-1}T_{\nu}L^{-1})^{\frac{1}{2}}h}}\right\|_{\nu},\quad\text{for}\quad
T_{\nu}=(I_{\nu}A^{\prime}(f^{\dagger}))^{*}(I_{\nu}A^{\prime}(f^{\dagger})).$

The link condition 8 considered in the linear case coincides with the above
inequality for $p=\frac{1}{2a}$.

Hence, through a comparison of the aforementioned link relationships, we
deduce that we need to assign $\frac{1}{2a}\leftarrow p$.

As reported in [68] for the statistical inverse problem, assuming smoothness
$s$ in the source condition 7, and considering eigenvalue decay condition 5
with parameter $b$, the optimal rate is of the order
$\mathopen{}\mathclose{{}\left(\frac{1}{\sqrt{m}}}\right)^{\frac{s}{s-1+\frac{1}{2a}+\frac{b}{2a}}}$.
In this context, we have to assign $\frac{1}{2a}\leftarrow p$. This yields a
lower bound of the order

$\mathopen{}\mathclose{{}\left(\frac{1}{\sqrt{m}}}\right)^{\frac{s}{p+s-1+bp}}.$

This corresponds to the upper bound discussed in Theorems 4.6, 4.7,
demonstrating that the rates are of optimal order.

##  5 Applications

Pharmacokinetic/pharmacodynamic (PK/PD) models are employed to predict the
evolution of drug concentrations and effects in patients undergoing treatment.
Covariate models tell us how patient attributes (covariates) influence model
parameters and play a crucial role in PK/PD model construction. Common
covariates are body weight, age, and levels of crucial biomarkers (e.g.,
plasma creatinine). The problem of interest is the determination of the
functional relationship between covariates and parameters.

Mechanistic model $G$.

In the two-compartment Pharmacokinetic model [64], we consider the human body
into two compartments: the central and peripheral compartments. It can be
described by the system of ODEs:

$\displaystyle V_{1}\frac{dC_{1}}{dt}=$ $\displaystyle
Q(C_{2}-C_{1})-\text{CL}\cdot C_{1},$ $\displaystyle V_{2}\frac{dC_{2}}{dt}=$
$\displaystyle Q(C_{1}-C_{2}).$

The central and peripheral compartments contain drug concentrations
represented by $C_{1}$ and $C_{2}$ in mg/L, respectively. These compartments
have volumes $V_{1}$ and $V_{2}$ in L, and are connected by intercompartmental
flow $Q$ in L/day, with clearance CL in L/day. Initial conditions are set as
$C_{1}(0)=D_{\text{rel}}\cdot\frac{w}{V_{1}}$ and $C_{2}(0)=0$, where
$D_{\text{rel}}$ is an intravenous bolus dose of 15 mg/kg body weight. The
covariate-to-parameter relationship includes a maturation component
($\text{mat}(a)$) dependent on age $a$ and an allometric component dependent
on weight $w$.

$\displaystyle
V_{1}=V_{1}^{*}\mathopen{}\mathclose{{}\left(\frac{w}{w_{\text{ref}}}}\right);$
$\displaystyle
V_{2}=V_{2}^{*}\mathopen{}\mathclose{{}\left(\frac{w}{w_{\text{ref}}}}\right);$
$\displaystyle
Q=Q^{*}\mathopen{}\mathclose{{}\left(\frac{w}{w_{\text{ref}}}}\right)^{\frac{3}{4}};$
$\displaystyle\text{CL}=\text{CL}^{*}\cdot\text{mat}(a)\mathopen{}\mathclose{{}\left(\frac{w}{w_{\text{ref}}}}\right)^{\frac{3}{4}};$

with reference body weight $w_{\text{ref}}=70$ kg. The weight-normalized
parameters
$\mathopen{}\mathclose{{}\left(V_{1}^{*},V_{2}^{*},Q^{*},\text{CL}^{*}}\right)$
are unknown.

The model’s parameters are denoted by
$\theta=\mathopen{}\mathclose{{}\left(V_{1},V_{2},Q,\text{CL}}\right)$. The
observed quantity

$G(x,\theta)=\mathopen{}\mathclose{{}\left(\text{ln}~{}C_{1}(t_{1}),\ldots,\text{ln}~{}C_{1}(t_{q})}\right)$

was measured at fixed time points $t_{1},\ldots,t_{q}$.

Hence, this problem can be described as a statistical inverse problem in [42]:

$y_{i}=G(x_{i},\theta_{i})+\varepsilon_{i},\qquad\theta_{i}=f(x_{i}),~{}~{}~{}~{}i\in\mathopen{}\mathclose{{}\left\\{1,\ldots,m}\right\\},$

where the observations consist of covariates (age $a$, weight $w$) $x_{i}\in
X\subset\mathbb{R}^{2}$ and (noisy) drug concentration measurements
$y_{i}\in\mathbb{R}^{q}$. There are unobserved parameters
$\theta_{i}\in\Theta\subset\mathbb{R}^{4}$, along with a mechanistic model
$G:X\times\Theta\to\mathbb{R}^{q}$, which is a nonlinear function. The mapping
$f$ takes covariates to parameters. Moreover, the noise $\varepsilon_{i}$ is
assumed to be normally distributed centered noise.

The mechanistic model $G$ represents a solution to a system of ordinary
differential equations (ODEs) governed by a two-compartment model observed at
different time points $t_{1},\ldots,t_{q}$. The direct dependence of $G$ on
$x_{i}$ allows us to model personalized doses and focus on specific aspects of
a covariate model. The operator $G$ exhibits a nonlinear dependence on the
parameters. Here, the aim is to estimate the unknown covariate model $f$.

Hartung et al. [42] used Tikhonov regularization to estimate this covariate
model in RKHS. The kernel is chosen empirically. Frequently, a specific
parametric category of covariate-to-parameter relationships (e.g., linear,
exponential, etc.) is chosen on an ad hoc basis. Therefore, the authors
conducted hypothesis tests to check the structure of the covariate model $f$,
i.e., whether it belongs to this parametric class or nonparametric class as a
kernelized Tikhonov regularized. Then, the authors assessed the effect of age
on the clearance of a monoclonal antibody.

Now, we check the validity of the assumptions considered for Theorem 4.3 for
the upper rates of convergence. First, we focus on the nonlinear structure
assumptions 9 and 10 for the non-linear operator $A(f)(x)=G(x,f(x))$. Let us
record the following theorem from [80] regarding its properties.

######  Theorem 5.1 ([80, Thm. 3.1]).

Let us assume that $\Theta$ is a nonempty, open, bounded subset of
$\mathbb{R}^{p}$ with suitably smooth boundary and $X\subset\mathbb{R}$. For
any $f\in\mathcal{H}=W^{k,2}(\mathbb{R})$ (see [70, Sec. 1.3.5]) we have
$G\in\mathcal{C}^{k+1}(\mathbb{R}\times\overline{\Theta})$. Moreover, we
observe that the operator $A:\mathcal{H}\to\mathcal{H}$ is continuous and the
Fréchet derivative of $A$ at $f$ is given by

$[A^{\prime}(f)g](x)=g(x)D_{\theta}G(x,f(x)),$

where $D_{\theta}$ denotes the derivative of $G$ with respect to the second
coordinate. In addition, the Fréchet derivative is Lipschitz continuous.

Since the derivative of $G$ is Lipschitz continuous in the first coordinate
and $\Theta$ is bounded. Therefore,

$\mathopen{}\mathclose{{}\left\lvert[A^{\prime}(f)g](x)}\right\rvert\leq\mathopen{}\mathclose{{}\left\lvert
D_{\theta}G(x,f(x))}\right\rvert\mathopen{}\mathclose{{}\left\lvert
g(x)}\right\rvert\leq
c_{G}\mathopen{}\mathclose{{}\left\|{x}}\right\|\mathopen{}\mathclose{{}\left\lvert\mathopen{}\mathclose{{}\left\langle
K_{x},g}\right\rangle}\right\rvert\leq\kappa
c_{G}\mathopen{}\mathclose{{}\left\|{g}}\right\|\mathopen{}\mathclose{{}\left\|{x}}\right\|$

implies that

$\mathopen{}\mathclose{{}\left\|{A^{\prime}(f)}}\right\|\leq\kappa c_{G}.$

Hence, the Fréchet derivative is bounded. The Lipschitz continuity of the
Fréchet derivative of $A$ implies the condition 10 (ii) (see [36, Chapt. 10]).

The normally distributed centered noise satisfies the assumptions 1, 2. The
Sobolev space $W^{k,2}(\mathbb{R})$ satisfies the assumption 3. In case the
design measure satisfies the source condition 11 and the eigenvalue decay
condition 5 (i), we obtain the minimax rates of convergence for the parameter
choice provided in Theorem 4.3.

##  6 Conclusions and outlook for future directions

Motivated by the growing interest in supervised inverse learning with random
design, we hope to have presented a thorough introduction to this modern
research field. To analyze both linear and non-linear inverse problems, we
have defined a general mathematical framework within reproducing kernel
Hilbert spaces, allowing us to appropriately consider the stochastic nature of
the data.

For tackling linear inverse problems, we adopt a least squares approach with
various types of regularization: spectral regularization, regularization by
projection, and convex regularization. We discuss optimal rates of convergence
for our estimators, initially dependent on unknown structural model
assumptions such as the smoothness of the ground truth and the eigenvalue
decay of the uncentered covariance operator (both describing the degree of
ill-posedness of the inverse problem), ultimately leading to fast rates of
convergence. To prevent overfitting, a proper choice of the regularization
parameter is crucial. Since a-priori choices are not feasible in practical
applications, we gave a brief overview of common data-driven approaches.

While inverse learning with a known linear forward operator is well developed
for a fair amount of regularization methods, much less is known for inverse
learning with a known forward non-linear operator. We discussed error bounds
only for Tikhonov regularization in both the classical setting and within
Hilbert scales using the least squares loss. A Standard approach for handling
non-linear problems is to consider linear approximations of the operator:
Under the assumption of Fréchet differentiability of $A$ and some Lipschitz
conditions, the non-linearity is mild enough to obtain optimal bounds,
matching those of the linear case. Notably, in the case where $A$ is not
Fréchet differentiable, error bounds can still be derived by means of
conditional stability estimates.

Unfortunaltely, in the non-linear case, the Tikhonov regularized solution is
not directly computable, and hence, other regularization approaches that take
the non-linear structure into account and are efficiently implementable need
to be analyzed further. Due to their practicability, the hope is to
successfully apply kernel methods also to obtain approximate solutions that
maintain statistical optimality. A major focus should be on the development of
data-driven approaches that are also efficient, such as early stopping rules.
We hope that our contribution can provide an impetus to further research in
this area.

List of symbols and notation

Samples and sample space:

  * •

$X=\mathbb{R}^{d}$ input space, $Y=\mathbb{R}$ output space, $Z=X\times Y$
sample space

  * •

$\mathbf{x}=(x_{1},\ldots,x_{m})$ input data,
$\mathbf{y}=(y_{1},\ldots,y_{m})$ output data and
$\mathbf{z}=((x_{1},y_{1}),\ldots,(x_{m},y_{m}))$ full training data

  * •

$m$ sample size

Spaces and norms:

  * •

$\mathcal{H}$, $\mathcal{H}^{\prime}$ Hilbert spaces with inner products
$\mathopen{}\mathclose{{}\left\langle\cdot,\cdot}\right\rangle_{\mathcal{H}}$
and
$\mathopen{}\mathclose{{}\left\langle\cdot,\cdot}\right\rangle_{\mathcal{H}^{\prime}}$,
respectively

  * •

$(\mathcal{H}_{a})_{a\in\mathbb{R}}$, Hilbert scale with the inner product
$\mathopen{}\mathclose{{}\left\langle\cdot,\cdot}\right\rangle_{a}$ and norm
$||\cdot||_{a}$

  * •

$||\cdot||$ operator norm

  * •

$||\cdot||_{2}$ euclidean norm

  * •

$\mathscr{L}^{p}(X,\nu)$ space of $p$-integrable functions on $X$, the
$\mathscr{L}^{2}(X,\nu)$-norm is denoted by $||\cdot||_{\nu}$

Operators:

  * •

$\mathcal{D}(A)$ domain of an operator $A$, $\mathcal{R}(A)$ range of an
operator $A$

  * •

$A:\mathcal{D}(A)\subset{\mathcal{H}}\to{\mathcal{H}}^{\prime}$
linear/nonlinear forward operator

  * •

$A^{\prime}(f)$ Fréchet derivative of $A$ at $f$

  * •

$L$ unbounded operator generating a Hilbert scale

  * •

$I_{\nu}:\mathcal{H}^{\prime}\hookrightarrow\mathscr{L}^{2}(X,\nu)$ canonical
injection operator

  * •

$B_{\nu}=I_{\nu}\circ A^{\prime}(f)$, $T_{\nu}=B_{\nu}^{*}B_{\nu}$ uncentered
covariance operator

  * •

$S_{\mathbf{x}}$ sampling operator, $B_{\mathbf{x}}=S_{\mathbf{x}}\circ
A^{\prime}(f)$, $T_{\mathbf{x}}=B_{\mathbf{x}}^{*}B_{\mathbf{x}}$

  * •

$P_{n}$ orthogonal projection onto $V_{n}$, an $n$-dimensional subspace of
$\mathcal{H}$

  * •

$B^{\dagger}$ pseudoinverse of an operator $B$

Measures:

  * •

$\rho$ is the distribution of $(x,y)$

  * •

$\nu$ marginal distribution on $X$

  * •

${\mathcal{P}}(X)$ class of marginal distribution on Borel $\sigma$-field
${\mathcal{B}}(X)$

  * •

${\mathcal{P}}_{\text{poly}}^{<}(\beta,b)$,
${\mathcal{P}}_{\text{poly}}^{>}(\alpha,b)$,
${\mathcal{P}}_{\text{str}}^{>}(\gamma)$ sub-classes of ${\mathcal{P}}(X)$

  * •

$\mathcal{M}_{\theta}$ model classes

Functions:

  * •

$f^{\dagger}$, $g^{\dagger}$: ground truths with $\dagger$ symbol

  * •

$\bar{f}$ initial guess of the ground truth

  * •

$g_{\lambda}$ spectral regularization

  * •

$f_{{\bf z}}$ estimated solution based on samples

  * •

$f_{\mathbf{z},\lambda}$ regularized solution

  * •

$d_{p}$ distance function

Parameter:

  * •

$\varepsilon$ observation noise

  * •

$\kappa^{2}$ upper bound for the reproducing kernel $K$ on $X$

  * •

$\sigma_{i}$, $\sigma_{\min}$ eigenvalues/spectral decay etc

  * •

$t$ spectral decay of $T_{\nu}$

  * •

$\lambda$ the regularization parameter

  * •

$b$ eigenvalue decay parameter

  * •

$a$ and $p$ for link condition

  * •

$a$, $p$ parameters in link condition and $p$ is also the benchmark smoothness
in distance function

  * •

$q$ qualification of the regularization

  * •

$\Sigma$, $M$ parameters in the noise condition 2

  * •

$R_{0}$ norm bound

  * •

$s,p,r$ all relevant to smoothness in source conditions

## References

  * [1] Sergios Agapiou and Peter Mathé.  Designing truncated priors for direct and inverse Bayesian problems.  Electronic Journal of Statistics, 16(1):158–200, 2022.
  * [2] Giovanni S. Alberti, Ernesto De Vito, Matti Lassas, Luca Ratti, and Matteo Santacesaria.  Learning the optimal Tikhonov regularizer for inverse problems.  Advances in Neural Information Processing Systems, 34:25205–25216, 2021.
  * [3] Nachman Aronszajn.  Theory of reproducing kernels.  Transactions of the American Mathematical Society, 68:337–404, 1950\.
  * [4] Frank Bauer, Sergei Pereverzev, and Lorenzo Rosasco.  On regularization algorithms in learning theory.  Journal of Complexity, 23(1):52–72, 2007.
  * [5] Martin Benning, Tatiana A. Bubba, Luca Ratti, and Danilo Riccio.  Trust your source: quantifying source condition elements for variational regularisation methods.  arXiv preprint arXiv:2303.00696, 2023.
  * [6] Martin Benning and Martin Burger.  Modern regularization methods for inverse problems.  Acta numerica, 27:1–111, 2018.
  * [7] Rajendra Bhatia.  Matrix analysis.  In Graduate Texts in Mathematics, volume 169. Springer-Verlag, New York, 1997.
  * [8] Nicolai Bissantz, Thorsten Hohage, Axel Munk, and Frits Ruymgaart.  Convergence rates of general regularization methods for statistical inverse problems and applications.  SIAM Journal on Numerical Analysis, 45(6):2610–2636, 2007.
  * [9] Gilles Blanchard, Marc Hoffmann, and Markus Reiß.  Early stopping for statistical inverse problems via truncated SVD estimation.  Electronic Journal of Statistics, 12(2), 2018.
  * [10] Gilles Blanchard, Marc Hoffmann, and Markus Reiß.  Optimal adaptation for early stopping in statistical inverse problems.  SIAM/ASA Journal on Uncertainty Quantification, 6(3):1043–1075, 2018.
  * [11] Gilles Blanchard and Peter Mathé.  Discrepancy principle for statistical inverse problems with application to conjugate gradient iteration.  Inverse problems, 28(11):115011, 2012.
  * [12] Gilles Blanchard, Peter Mathé, and Nicole Mücke.  Lepskii principle in supervised learning.  arXiv preprint arXiv:1905.10764, 2019.
  * [13] Gilles Blanchard and Nicole Mücke.  Optimal rates for regularization of statistical inverse learning problems.  Foundations of Computational Mathematics, 18(4):971–1013, 2018\.
  * [14] Tatiana A. Bubba, Martin Burger, Tapio Helin, and Luca Ratti.  Convex regularization in statistical inverse learning problems.  Inverse Problems and Imaging, 17(6):1193–1225, 2023.
  * [15] Martin Burger, Hendrik Dirks, Lena Frerking, Andreas Hauptmann, Tapio Helin, and Samuli Siltanen.  A variational reconstruction method for undersampled dynamic x-ray tomography based on physical motion models.  Inverse Problems, 33(12):124008, 2017.
  * [16] Martin Burger, Tapio Helin, and Hanne Kekkonen.  Large noise in variational regularization.  Transactions of Mathematics and its Applications, 2(1):1–45, 2018\.
  * [17] Martin Burger and Stanley Osher.  Convergence rates of convex variational regularization.  Inverse problems, 20(5):1411, 2004.
  * [18] Andrea Caponnetto and Ernesto De Vito.  Optimal rates for the regularized least-squares algorithm.  Foundations of Computational Mathematics, 7(3):331–368, 2007.
  * [19] Andrea Caponnetto and Yuan Yao.  Cross-validation based adaptation for regularization operators in learning theory.  Analysis and Applications, 8(2):161–183, 2010.
  * [20] Laurent Cavalier.  Inverse problems in statistics.  In Inverse Problems and High-Dimensional Estimation: Stats in the Château Summer School, August 31-September 4, 2009, pages 3–96. Springer, 2011.
  * [21] Laurent Cavalier, Georgi K. Golubev, Dominique Picard, and Alexandre B. Tsybakov.  Oracle inequalities for inverse problems.  The Annals of Statistics, 30(3):843–874, 2002.
  * [22] Laurent Cavalier and Alexandre Tsybakov.  Sharp adaptation for inverse problems with random noise.  Probability Theory and Related Fields, 123(3):323–354, 2002.
  * [23] Jonathan Chirinos Rodriguez, Ernesto De Vito, Cesare Molinari, Lorenzo Rosasco, and Silvia Villa.  On learning the optimal regularization parameter in inverse problems.  arXiv preprint arXiv:2311.15845, 2023.
  * [24] Pao-Liu Chow, Ildar A. Ibragimov, and Rafail Z. Khasminskii.  Statistical approach to some ill-posed problems for linear partial differential equations.  Probability theory and related fields, 113(3):421–441, 1999.
  * [25] Christian Clason and Bangti Jin.  A semismooth newton method for nonlinear parameter identification problems with impulsive noise.  SIAM journal on imaging sciences, 5(2):505–536, 2012.
  * [26] Felipe Cucker and Steve Smale.  Best choices for regularization parameters in learning theory: on the bias-variance problem.  Foundations of Computational Mathematics, 2(4):413–428, 2002.
  * [27] Felipe Cucker and Steve Smale.  On the mathematical foundations of learning.  Bull. Am. Math. Soc., 39(1):1–49, 2002.
  * [28] Felipe Cucker and Ding-Xuan Zhou.  Learning theory: An approximation theory viewpoint, volume 24.  Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, UK, 2007.
  * [29] Ingrid Daubechies.  Ten lectures on wavelets.  Philadelphia, PA, USA: SIAM, 1992.
  * [30] Ernesto De Vito, Massimo Fornasier, and Valeriya Naumova.  A machine learning approach to optimal Tikhonov regularization I: affine manifolds.  Analysis and Applications, 20(02):353–400, 2022.
  * [31] Ernesto De Vito, Sergei Pereverzyev, and Lorenzo Rosasco.  Adaptive kernel methods using the balancing principle.  Foundations of Computational Mathematics, 10:455–479, 2010.
  * [32] Ernesto De Vito, Lorenzo Rosasco, and Andrea Caponnetto.  Discretization error analysis for Tikhonov regularization.  Analysis and Applications, 4(01):81–99, 2006.
  * [33] Ernesto De Vito, Lorenzo Rosasco, Andrea Caponnetto, Umberto De Giovannini, Francesca Odone, and Peter Bartlett.  Learning from examples as an inverse problem.  Journal of Machine Learning Research, 6(5), 2005.
  * [34] Ernesto De Vito, Lorenzo Rosasco, Andrea Caponnetto, Umberto De Giovannini, and Francesca Odone.  Learning from examples as an inverse problem.  J. Mach. Learn. Res., 6:883–904, 2005.
  * [35] David L. Donoho.  Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition.  Applied and computational harmonic analysis, 2(2):101–126, 1995\.
  * [36] Heinz W. Engl, Martin Hanke, and Andreas Neubauer.  Regularization of inverse problems, volume 375.  Mathematics and its applications, Kluwer Academic Publishers Group, Dordrecht, The Netherlands, 1996.
  * [37] Matteo Giordano and Richard Nickl.  Consistency of Bayesian inference with Gaussian process priors in an elliptic inverse problem.  Inverse Problems, 36(8):085001, 2020.
  * [38] Sabrina Guastavino and Federico Benvenuto.  Convergence rates of spectral regularization methods: A comparison between ill-posed inverse problems and statistical kernel learning.  SIAM Journal on Numerical Analysis, 58(6):3504–3529, 2020.
  * [39] Shota Gugushvili, Aad Van Der Vaart, and Dong Yan.  Bayesian linear inverse problems in regularity scales.  Annales de l’Institut Henri Poincaré Probabilités et Statistiques, 56(3):2081–2107, 2020.
  * [40] Martin Hanke, Andreas Neubauer, and Otmar Scherzer.  A convergence analysis of the Landweber iteration for nonlinear ill-posed problems.  Numerische Mathematik, 72(1):21–37, 1995.
  * [41] Bastian Harrach, Tim Jahn, and Roland Potthast.  Beyond the Bakushinkii veto: regularising linear inverse problems without knowing the noise distribution.  Numerische Mathematik, 145:581–603, 2020.
  * [42] Niklas Hartung, Martin Wahl, Abhishake Rastogi, and Wilhelm Huisinga.  Nonparametric goodness-of-fit testing for parametric covariate models in pharmacometric analyses.  CPT: Pharmacometrics & Systems Pharmacology, 10(6):564–576, 2021\.
  * [43] Tapio Helin.  Least squares approximations in linear statistical inverse learning problems.  arXiv preprint arXiv:2211.12121, 2022.
  * [44] Bernd Hofmann.  Approximate source conditions in Tikhonov-Phillips regularization and consequences for inverse problems with multiplication operators.  Mathematical Methods in the Applied Sciences, 29(3):351–371, 2006\.
  * [45] Bernd Hofmann and Peter Mathé.  Analysis of profile functions for general linear regularization methods.  SIAM Journal on Numerical Analysis, 45(3):1122–1141, 2007.
  * [46] Thorsten Hohage and Mihaela Pricop.  Nonlinear Tikhonov regularization in Hilbert scales for inverse boundary value problems with random noise.  Inverse Problems and Imaging, 2(271-290):96, 2008.
  * [47] Thorsten Hohage and Frank Werner.  Inverse problems with poisson data: statistical regularization theory, applications and algorithms.  Inverse Problems, 32(9):093001, 2016.
  * [48] Bangti Jin and Xiliang Lu.  On the regularizing property of stochastic gradient descent.  Inverse Problems, 35(1):015004, 2018.
  * [49] Iain M. Johnstone and Bernard W. Silverman.  Discretization effects in statistical inverse problems.  Journal of complexity, 7(1):1–34, 1991.
  * [50] Jari Kaipio and Erkki Somersalo.  Statistical inverse problems: discretization, model reduction and inverse crimes.  Journal of Computational and Applied Mathematics, 198(2):493–504, 2007.
  * [51] Barbara Kaltenbacher, Andreas Neubauer, and Otmar Scherzer.  Iterative regularization methods for nonlinear ill-posed problems.  Walter de Gruyter, 2008.
  * [52] O. V. Lepskii.  On a problem of adaptive estimation in Gaussian white noise.  Theory of Probability & Its Applications, 35(3):454–466, 1991\.
  * [53] L. Lo Gerfo, Lorenzo Rosasco, Francesca Odone, Ernesto De Vito, and Alessandro Verri.  Spectral algorithms for supervised learning.  Neural Computation, 20(7):1873–1897, 2008.
  * [54] Shuai Lu, Peter Mathé, and Sergei V. Pereverzev.  Balancing principle in supervised learning for a general regularization scheme.  Applied and Computational Harmonic Analysis, 48(1):123–148, 2020\.
  * [55] Mark A. Lukas.  Comparisons of parameter choice methods for regularization with discrete noisy data.  Inverse Problems, 14(1):161, 1998.
  * [56] Bernard A. Mair and Frits H. Ruymgaart.  Statistical inverse estimation in Hilbert scales.  SIAM Journal on Applied Mathematics, 56(5):1424–1444, 1996.
  * [57] Shahar Mendelson and Joseph Neeman.  Regularization in kernel learning.  The Annals of Statistics, 38(1):526 – 565, 2010.
  * [58] François Monard, Richard Nickl, and Gabriel P. Paternain.  Consistent inversion of noisy non-abelian x-ray transforms.  Communications on Pure and Applied Mathematics, 74(5):1045–1099, 2021.
  * [59] Nicole Mücke and Enrico Reiss.  Stochastic gradient descent in Hilbert scales: smoothness, preconditioning and earlier stopping.  arXiv preprint arXiv:2006.10840, 2020.
  * [60] Frank Natterer.  Error bounds for Tikhonov regularization in Hilbert scales.  Applicable Analysis, 18(1-2):29–37, 1984.
  * [61] Richard Nickl.  Bayesian non-linear statistical inverse problems.  Zurich Lectures in Advanced Mathematics. EMS Press, 2023.
  * [62] Richard Nickl, Sara van de Geer, and Sven Wang.  Convergence rates for penalized least squares estimators in PDE constrained regression problems.  SIAM/ASA Journal on Uncertainty Quantification, 8(1):374–413, 2020\.