Optimal Recovery from Inaccurate Data in Hilbert Spaces:
Regularize, but what of the Parameter?  

This article is concerned with a central problem in Data Science, namely: a function $f$ is acquired through point evaluations

and these data should be used to learn $f$—or to recover it, with the terminology preferred in this article. Importantly, the evaluation points $x^{(1)},\ldots,x^{(m)}$ are considered fixed entities in our scenario: they cannot be chosen in a favorable way, as in Information-Based Complexity [Novak and Woźniakowski, 2008], nor do they occur as independent realizations of a random variable, as in Statistical Learning Theory [Hastie, Tibshirani, and Friedman, 2009]. In particular, without an underlying probability distribution, the performance of the recovery process cannot be assessed via generalization error. Instead, it is assessed via a notion of worst-case error, central to the theory of Optimal Recovery [Micchelli and Rivlin, 1977].

To outline this theory, we make the framework slightly more abstract. Precisely, given a normed space $F$, the unknown function is replaced by an element $f\in F$. This element is accessible only through a priori information expressing an educated belief about $f$ and a posteriori information akin to (1). In other words, our partial knowledge about $f$ is summed up via

• •

  the fact that $f\in\mathcal{K}$ for a subset $\mathcal{K}$ of $F$ called a model set;

• •

  the observational data $y_{i}=\lambda_{i}(f)$, $i=1,\ldots,m$, for some linear functionals $\lambda_{1},\ldots,\lambda_{m}\in F^{*}$ making up the observation map $\Lambda:g\in F\mapsto[\lambda_{1}(g);\ldots;\lambda_{m}(g)]\in\mathbb{R}^{m}$.

We wish to approximate $f$ by some $\widehat{f}\in F$ produced using this partial knowledge of $f$. Since the error $\|f-\widehat{f}\|$ involves the unknown $f$, which is only accessible via $f\in\mathcal{K}$ and $\Lambda(f)=y$, we take a worst-case perspective leading to the local worst-case error

Our objective consists in finding an element $\widehat{f}$ that minimizes ${\rm lwce}(y,\widehat{f})$. Such an $\widehat{f}$ can be described, almost tautologically, as a center of a smallest ball containing $\mathcal{K}\cap\Lambda^{-1}(\{y\})$. It is called a Chebyshev center of this set of model- and data-consistent elements. This remark, however, does not come with any practical construction of a Chebyshev center.

The term local was used above to make a distinction with the global worst-case error of a recovery map $\Delta(=\Delta_{\mathcal{K}}):\mathbb{R}^{m}\to F$, defined as

The minimal value of ${\rm gwce}(\Delta)$ is called the intrinsic error (of the observation map $\Lambda$ over the model set $\mathcal{K}$) and the maps $\Delta$ that achieve this minimal value are called globally optimal recovery maps. Our objective consists in constructing such maps—of course, the map that assigns to $y$ a Chebyshev center of $\mathcal{K}\cap\Lambda^{-1}(\{y\})$ is one of them, but it may be impractical. By contrast, for model sets that are convex and symmetric, the existence of linear maps among the set of globally optimal recovery maps is guaranteed by fundamental results from Optimal Recovery in at least two settings: when $F$ is a Hilbert space and when $F$ is an arbitrary normed space but the full recovery of $f$ gives way to the recovery of a quantity of interest $Q(f)$, $Q$ being a linear functional. We refer the readers to [Foucart, To appear, Chapter 9] for details.

The problem solved in this article is a quintessential Optimal Recovery problem—its specificity lies in the particular model set and in the incorporation of errors in the observation process. The underlying normed space $F$ is a Hilbert space and is therefore denoted by $H$ from now on. Reproducing kernel Hilbert spaces, whose usage is widespread in Data Science [Schölkopf and Smola, 2002], are of particular interest as point evaluations of type (1) make perfect sense there.

Concerning the model set, we concentrate on an approximation-based choice that is increasingly scrutinized, see e.g. [Maday, Patera, Penn, and Yano, 2015], [DeVore, Petrova, and Wojtaszczyk, 2017] and [Cohen, Dahmen, Mula, and Nichols, 2020]. Depending on a linear subspace $\mathcal{V}$ of $H$ and on a parameter $\varepsilon>0$, it takes the form

Binev, Cohen, Dahmen, DeVore, Petrova, and Wojtaszczyk [2017] completely solved the Optimal Recovery problem with exact data in this situation (locally and globally). Precisely, they showed that the solution $\widehat{f}$ to

which clearly belongs to the model- and data-consistent set $\mathcal{K}\cap\Lambda^{-1}(\{y\})$, turns out to be its Chebyshev center. Moreover, with $P_{\mathcal{V}}$ and $P_{\mathcal{V}^{\perp}}$ denoting the orthogonal projectors onto $\mathcal{V}$ and onto the orthogonal complement $\mathcal{V}^{\perp}$ of $\mathcal{V}$, the fact that ${\rm dist}(f,\mathcal{V})=\|f-P_{\mathcal{V}}f\|=\|P_{\mathcal{V}^{\perp}}f\|$ makes the optimization program (4) tractable. It can actually be seen that $\Delta:y\mapsto\widehat{f}$ is a linear map. This is a significant advantage because $\Delta$ can then be precomputed in an offline stage knowing only $\mathcal{V}$ and $\Lambda$ and the program (4) need not be solved afresh for each new data $y\in\mathbb{R}^{m}$ arriving in an online stage.

Concerning the observation process, instead of exact data $y=\Lambda(f)\in\mathbb{R}^{m}$, it is now assumed that

for some unknown error vector $e\in\mathbb{R}^{m}$. This error vector is not modeled as random noise but through the deterministic $\ell_{2}$-bound $\|e\|_{2}\leq\eta$. Although other $\ell_{p}$-norms can the considered for the optimal recovery of $Q(f_{0})$ when $Q$ is a linear functional on an arbitrary normed space $F$ (see [Ettehad and Foucart, 2021]), here the arguments rely critically on $\mathbb{R}^{m}$ being endowed with the $\ell_{2}$-norm. It will be written simply as $\|\cdot\|$ below, hoping that it does not create confusion with the Hilbert norm on $H$.

For our specific problem, the worst-case recovery errors (2) and (3) need to be adjusted. The local worst-case recovery error at $y$ for $\widehat{f}$ becomes

As for the global worst-case error of $\Delta:\mathbb{R}^{m}\to H$, it reads

Note that both worst-case errors are infinite if one can find a nonzero $h$ in $\mathcal{V}\cap\ker(\Lambda)$. Indeed, the element $f_{t}:=f+th$, $t\in\mathbb{R}$, obeys $\|P_{\mathcal{V}^{\perp}}f_{t}\|=\|P_{\mathcal{V}^{\perp}}f\|\leq\varepsilon$ and $\|y-\Lambda(f_{t})\|=\|y-\Lambda(f)\|\leq\eta$, so for instance ${\rm lwce}(y,\widehat{f})\geq\sup_{t\in\mathbb{R}}\|f_{t}-\widehat{f}\|=+\infty$. Thus, we always make the assumption that

We keep in mind that the latter forces $n:=\dim(\mathcal{V})\leq m$, as can be seen by dimension arguments. With $\Lambda^{*}$ denoting the Hermitian adjoint of $\Lambda$, another assumption that we sometimes make reads

This is not extremely stringent: assuming the surjectivity of $\Lambda$ is quite natural, otherwise certain observations need not be collected; then the map $\Lambda$ can be preprocessed into another map $\widetilde{\Lambda}$ satisfying $\widetilde{\Lambda}\widetilde{\Lambda}^{*}=\mathrm{Id}_{\mathbb{R}^{m}}$ by setting $\widetilde{\Lambda}=(\Lambda\Lambda^{*})^{-1/2}\Lambda$. Incidentally, if $u_{1},\ldots,u_{m}\in H$ represent the Riesz representers of the observation functionals $\lambda_{1},\ldots,\lambda_{m}\in H^{*}$, characterized by $\langle u_{i},f\rangle=\lambda_{i}(f)$ for all $f\in H$, then the assumption (6) is equivalent to the orthonormality of the system $(u_{1},\ldots,u_{m})$. In a reproducing kernel Hilbert space with kernel $K$, if the $\lambda_{i}$’s are point evaluations at some $x^{(i)}$’s, so that $u_{i}=K(\cdot,x^{(i)})$, then (6) is equivalent to $K(x^{(i)},x^{(j)})=\delta_{i,j}$ for all $i,j=1,\ldots,m$. This occurs e.g. for the Paley–Wiener space of functions with Fourier transform supported on $[-\pi,\pi]$ when the evaluations points come from an integer grid, since the kernel is given by $K(x,x^{\prime})={\rm sinc}(\pi(x-x^{\prime}))$, $x,x^{\prime}\in\mathbb{R}$.

There are previous works on Optimal Recovery in Hilbert spaces in the presence of observation error bounded in $\ell_{2}$. Notably, [Beck and Eldar, 2007] dealt with the local setting, while [Melkman and Micchelli, 1979] and [Micchelli, 1993] dealt with the global setting. These works underline the importance of regularization, which is prominent in many other settings [Chen and Haykin, 2002]. They establish that the optimal recovery maps are obtained by solving the unconstrained program

for some $\tau\in[0,1]$. It is the precise choice of this regularization parameter $\tau$ which is the purpose of this article. Assuming from now on that $H$ is finite dimensional²²2It is likely that the results are still valid in the infinite-dimensional case. But then it is unclear how to solve (8) and (9) numerically, so the infinite-dimensional case is not given proper scrutiny in the rest of the article., we provide a complete (almost) picture of the local and global Optimal Recovery solutions, as summarized in the four points below, three of them being new:

1. L1.

   With $H$ restricted here to be a complex Hilbert space, the Chebyshev center of the set $\{f\in H:\|P_{\mathcal{V}^{\perp}}f\|\leq\varepsilon,\|\Lambda f-y\|\leq\eta\}$ is the minimizer of (7) for the choice $\tau=d_{\sharp}/(c_{\sharp}+d_{\sharp})$, where $c_{\sharp},d_{\sharp}$ are solutions to the semidefinite program³³3In the statement of this semidefinite program and elsewhere, the notation $T\succeq 0$ means that an operator $T$ is positive semidefinite on $H$, i.e., that $\langle Tf,f\rangle\geq 0$ for all $f\in H$.

   	$\displaystyle\underset{c,d,t\geq 0}{\rm minimize}\,\;\varepsilon^{2}c+(\eta^{2}-\|y\|^{2})d+t$	s.to	$\displaystyle cP_{\mathcal{V}^{\perp}}+d\Lambda^{*}\Lambda\succeq\mathrm{Id},$
   		and	$\displaystyle\begin{bmatrix}cP_{\mathcal{V}^{\perp}}+d\Lambda^{*}\Lambda&|&-d\Lambda^{*}y\\ \hline\cr-d(\Lambda^{*}y)^{*}&|&t\end{bmatrix}\succeq 0.$

2. L2.

   Under the orthonormal observations assumption (6) but without the above restriction on $H$, the Chebyshev center of the set $\{f\in H:\|P_{\mathcal{V}^{\perp}}f\|\leq\varepsilon,\|\Lambda f-y\|\leq\eta\}$ is the minimizer of (7) for the choice $\tau$ that satisfies

   (8)

   $$\lambda_{\min}((1-\tau)P_{\mathcal{V}^{\perp}}+\tau\Lambda^{*}\Lambda)=\frac{(1-\tau)^{2}\varepsilon^{2}-\tau^{2}\eta^{2}}{(1-\tau)\varepsilon^{2}-\tau\eta^{2}+(1-\tau)\tau(1-2\tau)\delta^{2}},$$

   where $\delta$ is precomputed as $\delta=\min\{\|P_{\mathcal{V}^{\perp}}f\|:\Lambda f=y\}=\min\{\|\Lambda f-y\|:f\in\mathcal{V}\}$. For the distinct case $\mathcal{V}=\{0\}$, the best choice of parameter is more simply $\tau=\max\{1-\eta/\|y\|,0\}$.

3. G1.

   A globally optimal recovery map is provided by the linear map sending $y\in\mathbb{R}^{m}$ to the minimizer of (7) with parameter $\tau=d_{\flat}/(c_{\flat}+d_{\flat})$, where $c_{\flat},d_{\flat}$ are solutions to the semidefinite program

   (9)

   $$\underset{c,d}{\rm minimize}\,\;\varepsilon^{2}c+\eta^{2}d\qquad\mbox{s.to}\quad cP_{\mathcal{V}^{\perp}}+d\Lambda^{*}\Lambda\succeq\mathrm{Id}.$$

4. G2.

   Under the orthonormal observations assumption (6), the linear map sending $y\in\mathbb{R}^{m}$ to the minimizer of (7) is a globally optimal recovery map for any choice of parameter $\tau\in[0,1]$.

Before entering the technicalities, a few comments are in order to put these results in context. Item L1 is the result of [Beck and Eldar, 2007] (see Corollary 3.2 there) adapted to our situation. It relies on an extension of the S-lemma involving two quadratic constraints. This extension is valid in the complex finite-dimensional setting, but not necessarily in the real setting, hence the restriction on $H$ (this does not preclude the validity of the result in the real setting, though). It is worth pointing out the nonlinearity of the map that sends $y\in\mathbb{R}^{m}$ to the above Chebyshev center. Incidentally, we can safely talk about the Chebyshev center, because it is known [Garkavi, 1962] that a bounded set in a uniformly convex Banach space has exactly one Chebyshev center. A sketch of the argument adapted to our situation is presented in the appendix.

For item L2, working with an observation map $\Lambda$ satisfying $\Lambda\Lambda^{*}=\mathrm{Id}_{\mathbb{R}^{m}}$ allows us to construct the Chebyshev center even in the setting of a real Hilbert space. This is possible because our argument does not rely on the extension of the S-lemma—it just uses the obvious implication. As for equation (8), it is easily solved using the bisection method or the Newton/secant method. Moreover, it gives some insight on the value of the optimal parameter $\tau$. For instance, the proof reveals that $\tau$ is always between $1/2$ and $\varepsilon/(\varepsilon+\eta)$. When $\varepsilon\geq\eta$, say, the optimal parameter should then satisfy $\tau\geq 1/2$, which is somewhat intuitive: $\varepsilon\geq\eta$ means that there is more model mismatch than data mismatch, so the regularization should penalize model fidelity less than data fidelity by taking $1-\tau\leq\tau$, i.e., $\tau\geq 1/2$. As an aside, we point out that, here too, the map that sends $y\in\mathbb{R}^{m}$ to the Chebyshev center is not a linear map—if it was, then the optimal parameter should be independent of $y$.

In contrast, the globally optimal recovery map of item G1 is linear. It is one of several globally optimal recovery maps, since the locally optimal one (which is nonlinear) is also globally optimal. However, as revealed in the reproducible⁴⁴4matlab and Python files illustrating the findings of this article are located at https://github.com/foucart/COR. accompanying this article, it is in general the only regularization map that turns out to be globally optimal. The fact that regularization produces globally optimal recovery maps was recognized by Micchelli, who wrote in the abstract of [Micchelli, 1993] that “the regularization parameter must be chosen with care”. However, a recipe for selecting the parameter was not given there, except on a specific example. The closest to a nonexhaustive search is found in [Plaskota, 1996, Lemma 2.6.2] for the case $\mathcal{V}=\{0\}$, but even this result does not translate into a numerically tractable recipe. The selection stemming from (9) does, at least when $H$ is finite-dimensional, which is assumed here. Semidefinite programs can indeed be solved in matlab using CVX [Grant and Boyd, 2014] and in Python using CVXPY [Diamond and Boyd, 2016].

Finally, a surprise arises in item G2. Working with an observation map $\Lambda$ satisfying $\Lambda\Lambda^{*}=\mathrm{Id}_{\mathbb{R}^{m}}$, the latter indeed reveals that the regularization parameter does not need to be chosen with care after all, since regularization maps are globally optimal no matter how the parameter $\tau\in[0,1]$ is chosen. The precise interpretation of the choices $\tau=0$ and $\tau=1$ will be elucidated later.

The rest of this article is organized as follows. Section 2 gathers some auxiliary results that are used in the proofs of the main results. Section 3 elucidates item L1 and establishes item L2—in other words, it is concerned with local optimality. Section 4, which is concerned with global optimality, is the place where items G1 and G2 are proved. Lastly, a short appendix containing some side information is included after the bibliography.

Loosely speaking, the S-procedure is a relaxation technique expressing the fact that a quadratic inequality is a consequence of some quadratic constraints. In case of a single quadratic constraint, the relaxation turns out to be exact. This result, known as the S-lemma, can be stated as follows: given quadratic functions $q_{0}$ and $q_{1}$ defined on $\mathbb{K}^{N}$, with $\mathbb{K}=\mathbb{R}$ or $\mathbb{K}=\mathbb{C}$,

provided $q_{1}(\widetilde{x})<0$ for some $\widetilde{x}\in\mathbb{K}^{N}$. With more than one quadratic constraint, $q_{1},\ldots,q_{k}$, say, $q_{0}(x)\leq 0$ whenever $q_{1}(x)\leq 0,\ldots,q_{k}(x)\leq 0$ is still a consequence of $q_{0}\leq a_{1}q_{1}+\cdots+a_{k}q_{k}$ for some $a_{1},\ldots,a_{k}\geq 0$, but the reverse implication does not hold anymore. There is a subtlety when $k=2$, as the reverse implication holds for $\mathbb{K}=\mathbb{C}$ but not for $\mathbb{K}=\mathbb{R}$, see [Pólik and Terlaky, 2007, Section 3]. However, if the quadratic constraints do not feature linear terms, then the reverse implication holds for $k=2$ also when $\mathbb{K}=\mathbb{R}$. Since this result of [Polyak, 1998, Theorem 4.1] is to be invoked later, we state it formally below.

In this subsection, we take a closer look at the regularization program (7). The result below shows that its solution depends linearly on $y\in\mathbb{R}^{m}$. In fact, the result covers a slightly more general program and the linearity claim follows by taking $R=P_{\mathcal{V}^{\perp}}$, $r=0$, $S=\Lambda$, and $s=y$.

The expression (12) is not always the most convenient one. Under extra conditions on $R$ and $S$, we shall see that $f_{\tau}$, $\tau\in[0,1]$, can in fact be expressed as the convex combination $f_{\tau}=(1-\tau)f_{0}+\tau f_{1}$. The elements $f_{0}$ and $f_{1}$ should be interpreted⁵⁵5Intuitively, the solution to the program (10) written as the minimization of $\|Rf-r\|^{2}+(\tau/(1-\tau))\|Sf-s\|^{2}$ becomes, as $\tau\to 1$, the mininizer of $\|Rf-r\|^{2}$ subject to $\|Sf-s\|^{2}=0$. This explains the interpretation of $f_{1}$. A similar argument explains the interpretation of $f_{0}$. as

The requirements that $r\in{\rm ran}(R)$ and $s\in{\rm ran}(S)$ need to be imposed for $f_{0}$ and $f_{1}$ to even exist and the condition $\ker(R)\cap\ker(S)=\{0\}$ easily guarantees that $f_{0}$ and $f_{1}$ are unique. They obey

For instance, the identity $Rf_{0}=r$ reflects the constraint in the optimization program defining $f_{0}$, while $S^{*}(Sf_{0}-s)\in\ker(R)^{\perp}$ is obtained by expanding $\|S(f_{0}+tu)-s\|^{2}\geq\|S(f_{0})-s\|^{2}$ around $t=0$ for any $u\in\ker(R)$. At this point, we are ready to establish our claim under extra conditions on $R$ and $S$, namely that they are orthogonal projections. These conditions will be in place when the observation map satisfies $\Lambda\Lambda^{*}=\mathrm{Id}_{\mathbb{R}^{m}}$. Indeed, in view of $\|w\|^{2}=\langle w,\Lambda\Lambda^{*}w\rangle=\|\Lambda^{*}w\|^{2}$ for any $w\in\mathbb{R}^{m}$, the regularization program (7) also reads

where both $P_{\mathcal{V}^{\perp}}$ and $\Lambda^{*}\Lambda$ are orthogonal projections. The result below will then be applied with $R=P_{\mathcal{V}^{\perp}}$, $r=0$, $S=\Lambda^{*}\Lambda$, and $s=\Lambda^{*}y$.

We complement Proposition 3 with a few additional pieces of information.

In this subsection, we reproduce a result from [Beck and Eldar, 2007], albeit with different notation, and explain how it implies the statement of item L1. The result in question, namely Corollary 3.2, relies on the S-procedure with two constraints, and as such cannot be claimed in the real setting.

The statement made in item L1 is of course derived by taking $R=P_{\mathcal{V}^{\perp}}$, $r=0$, $S=\Lambda$, and $s=y$. Theorem 4 is indeed applicable, as $\widetilde{f}=(f_{0}+f_{1})/2$ satisfies the strict feasibility conditions, while the positive definiteness condition is not only fulfilled for some $\tau\in[0,1]$, but for all $\tau\in(0,1)$, since $\langle((1-\tau)P_{\mathcal{V}^{\perp}}+\tau\Lambda^{*}\Lambda)f,f\rangle=(1-\tau)\|P_{\mathcal{V}^{\perp}}f\|^{2}+\tau\|\Lambda f\|^{2}\geq 0$, with equality only possible if $f\in\mathcal{V}\cap\ker(\Lambda)$, i.e., if $f=0$ thanks to the assumption (5). We also note that, by virtue of (12), the element $f_{\sharp}$ defined above is nothing else than the regularized solution with parameter $\tau=d_{\sharp}/(c_{\sharp}+d_{\sharp})$.

In this subsection, we place ourselves in the situation of an observation map satisfying $\Lambda\Lambda^{*}=\mathrm{Id}_{\mathbb{R}^{m}}$ and we provide a proof of the statements made item L2. In fact, we prove some slightly more general results and L2 follows by taking $R=P_{\mathcal{V}^{\perp}}$, $r=0$, $S=\Lambda^{*}\Lambda$, and $s=\Lambda^{*}y$. Note that we must separate the cases where $R=\mathrm{Id}$ (corresponding to $\mathcal{V}=\{0\}$) and where $R$ is a proper orthogonal projection (corresponding to $\mathcal{V}\not=\{0\}$). We emphasize that, in each of these two cases, the optimal parameter $\tau_{\sharp}$ is not independent of $y$. Therefore, in view of (15) and of the linear dependence of $f_{0}$ and $f_{1}$ on $y$, the regularized solution $f_{\tau_{\sharp}}$ does not depend linearly on $y$. In other words, the locally optimal recovery map is not a linear map. The following two simple lemmas will be used to deal with both cases.