Nonparametric denoising Signals of Unknown Local Structure, II: Nonparametric Regression Estimation

Let ${\mathbf{F}}=(\Omega,\Sigma,P)$ be a probability space. We consider the problem of recovering unknown complex-valued random field $(s_{\tau}=s_{\tau}(\xi))_{{\tau\in\mathbb{Z}^{d}\atop\xi\in\Omega}}$ over $\mathbb{Z}^{d}$ from noisy observations

We assume that the field $(e_{\tau})$ of observation noises is independent of $(s_{\tau})$ and is of the form $e_{\tau}=\sigma\epsilon_{\tau}$, where $(\epsilon_{\tau})$ are independent of each other standard Gaussian complex-valued variables; the adjective “standard” means that $\Re(\epsilon_{\tau})$, $\Im(\epsilon_{\tau})$ are independent of each other ${\mathbf{N}}(0,1)$ random variables.

We suppose that the observations (1) come from a function (“signal”) $f$ of continuous argument (which we assume to vary in the $d$-dimensional unit cube $[0,1]^{d}$); this function is observed in noise along an $n$-point equidistant grid in $[0,1]^{d}$, and the problem is to recover $f$ via these observations. This problem fits the framework of nonparametric regression estimation with a “traditional setting” as follows:

A. The objective is to recover an unknown smooth function $f:\,[0,1]^{d}\to{\mathbb{R}}$, which is sampled on the observation grid $\Gamma_{n}=\{x_{\tau}=m^{-1}\tau:\,0\leq\tau_{1},...,\tau_{d}\leq m\}$ with $(m+1)^{d}=n$, so that $s_{\tau}=f(x_{\tau})$. The error of recovery is measured with some functional norm (or a semi-norm) $\|\cdot\|$ on $[0,1]^{d}$, and the risk of recovery $\widehat{f}$ of $f$ is the expectation $E_{f}\|\widehat{f}-f\|^{2}$;

B. The estimation routines are aimed at recovering smooth signals, and their quality is measured by their maximal risks, the maximum being taken over $f$ running through natural families of smooth signals, e.g., Hölder or Sobolev balls;

C. The focus is on the asymptotic, as the volume of observations $n$ goes to infinity, behavior of the estimation routines, with emphasis on asymptotically minimax (nearly) optimal estimates – those which attain (nearly) best possible rates of convergence of the risks to 0 as the observation sample size $n\to\infty$.

Initially, the research was focused on recovering smooth signals with a priori known smoothness parameters and the estimation routines were tuned to these parameters (see, e.g., [23, 34, 38, 24, 2, 31, 39, 22, 36, 21, 27]). Later on, there was a significant research on adaptive estimation. Adaptive estimation methods are free of a priori assumptions on the smoothness parameters of the signal to be recovered, and the primary goal is to develop the routines which exhibit asymptotically (nearly) optimal behavior on a wide variety of families of smooth functions (cf. [35, 28, 29, 30, 6, 8, 9, 25, 3, 7, 19]). For a more compete overview of results on smooth nonparametric regression estimation see, for instance, [33].³³3Our “brief outline” of adaptive approach to nonparametric regression would be severely incomplete without mentioning a novel approach aimed at recovering nonsmooth signals possessing sparse representations in properly constructed functional systems [5, 10, 4, 11, 12, 13, 14, 15, 16, 17, 37, 18]. This promising approach is completely beyond the scope of our paper.
The traditional focus on recovering smooth signals ultimately comes from the fact that such a signal locally can be well-approximated by a polynomial of a fixed order $r$, and such a polynomial is an “easy to estimate” entity. Specifically, for every integer $T\geq 0$, the value of a polynomial $p$ at an observation point $x_{t}$ can be recovered via $(2T+1)^{d}$ neighboring observations $\{x_{\tau}:|\tau_{j}-t_{j}|\leq T,1\leq j\leq d\}$ “at a parametric rate” – with the expected squared error $C\sigma^{2}(2T+1)^{-d}$ which is inverse proportional to the amount $(2T+1)^{d}$ of the observations used by the estimate. The coefficient $C$ depends solely on the order $r$ and the dimensionality $d$ of the polynomial. The corresponding estimate $\widehat{p}(x_{t})$ of $p(x_{t})$ is pretty simple: it is given by a “time-invariant filter”, that is, by convolution of observations with an appropriate discrete kernel $q^{(T)}=(q^{(T)}_{\tau})_{\tau\in\mathbb{Z}^{d}}$ vanishing outside the box ${\mathbf{O}}_{T}=\{\tau\in\mathbb{Z}^{d}:|\tau_{j}|\leq T,1\leq j\leq d\}$:

then the estimation $\widehat{f}$ of $f(x_{t})$ is taken as $\widehat{f}=\widehat{p}(x_{t})$.

Note that the kernel $q^{(T)}$ is readily given by the degree $r$ of the approximating polynomial, $T$ and dimension $d$. The “classical” adaptation routines takes care of choosing “good” values of the approximation parameters (namely, $T$ and $r$). On the other hand, the polynomial approximation “mechanism” is supposed to be fixed once for ever. Thus, in those procedures the “form” of the kernel is considered as given in advance.

In the companion paper [26] (referred hereafter as Part I) we have introduced the notion of a well-filtered signal. In brief, the signal $(s_{\tau})_{\tau\in\mathbb{Z}^{d}}$ is $T$-well-filtered for some $T\in\mathbb{N}_{+}$ if there is a filter (kernel) $q=q^{(T)}\in\mathbf{O}_{T}$ which recovers $(s_{\tau})$ in the box $\{u:|u-t|\leq 3T\}$ with the mean square error comparable with $\sigma T^{-d/2}$:

The universe of these signals is much wider than the one of smooth signals.As we have seen in Part I that it contains, in particular, “modulated smooth signals” – sums of a fixed number of products of smooth functions and multivariate harmonic oscillations of unknown (and arbitrarily high) frequencies. We have shown in Part I that whenever a discrete time signal (that is, a signal defined on a regular discrete grid) is well-filtered, we can recover this signal at a “nearly parametric” rate without a priori knowledge of the associated filter. In other words, a well-filtered signal can be recovered on the observation grid basically as well as if it were an algebraic polynomial of a given order.

We are about to demonstrate that the results of Part I on recovering well-filtered signals of unknown structure can be applied to recovering nonparametric signals which admit well-filtered local approximations. Such an extension has an unavoidable price – now we cannot hope to recover the signal well outside of the observation grid (a highly oscillating signal can merely vanish on the observation grid and be arbitrarily large outside it). As a result, in what follows we are interested in recovering the signals along the observation grid only and, consequently, replace the error measures based on the functional norms on $[0,1]^{d}$ by their grid analogies.

The estimates to be developed will be “double adaptive”, that is, adaptive with respect to both the unknown in advance structures of well-filtered approximations of our signals and to the unknown in advance “approximation rate” – the dependence between the size of a neighborhood of a point where the signal in question is approximated and the quality of approximation in this neighborhood. Note that in the case of smooth signals, this approximation rate is exactly what is given by the smoothness parameters. The results to follow can be seen as extensions of the results of [32, 20] (see also [33]) dealing with the particular case of univariate signals satisfying differential inequalities with unknown differential operators.

We start with the formal description of the components of the nonparametric regression problem.

Let for $\tau\in\mathbb{Z}^{d}$, $|\tau|=\max\{|\tau_{1}|,...,|\tau_{d}|\}$, and let $\tau\leq m$ for some $a\in\mathbb{N}$ denote $\tau_{i}\leq m,\;i=1,...,d$. Let $m$ be a positive integer, $n=(m+1)^{d}$, and let $\Gamma_{n}=\{x=m^{-1}\alpha:\alpha\in\mathbb{Z}^{d},0\leq\alpha,|\alpha|\leq m\}$.

Let $C([0,1]^{d})$ be the linear space of complex-valued fields over $[0,1]^{d}$. We associate with a signal $f\in C([0,1]^{d})$ its observations along $\Gamma_{n}$:

where $\{\epsilon_{\tau}\}_{\tau\in\mathbb{Z}^{d}}$ are independent standard Gaussian complex-valued random noises. Our goal is to recover $f\big{|}_{\Gamma_{n}}$ from observations (2). In what follows, we write

Below we use the following notations. For a set $B\subset[0,1]^{d}$, we denote by $\mathbb{Z}(B)$ the set of all $t\in\mathbb{Z}^{d}$ such that $m^{-1}t\in B$. We denote $\|\cdot\|_{q,B}$ the standard $L_{p}$-norm on $B$:

and $|g|_{q,B}$ its discrete analogy, so that

We set

Let $x=m^{-1}t\in\Gamma_{n}^{o}$. We say that a nonempty open cube

centered at $x$ is admissible for $x$, if $B_{h}(x)\subset[0,1]^{n}$. For such a cube, $T_{h}(x)$ denotes the largest nonnegative integer $T$ such that

For a cube

$D(B)=h$ stands for the edge of $B$. For $\gamma\in(0,1)$ we denote

the $\gamma$-shrinkage of $B$ to the center of $B$.

The recovering routine we are about to build is aimed at estimating functions from classes ${\mathbf{F}}^{k,\rho,p}(B,R)$ with unknown in advance parameters $k,\rho,p,B,R$. The only design parameters of the routine is an a priori upper bound $\mu$ on the parameter $\rho$ and a $\gamma\in(0,1)$.

Our main result is as follows:

Note that the rates of convergence to 0, as $n\to\infty$, of the risks $\widehat{{\mathbf{R}}}_{q}\left(\widehat{f}_{n};{\mathbf{F}}\right)$ of our adaptive estimate on the families ${\mathbf{F}}={\mathbf{F}}^{k,\rho,p}(B;R)$ are exactly the same as those stated by Theorem 3 from [31] (see also [30, 9, 19, 33]) in the case of recovering non-parametric smooth regression functions from Sobolev balls. It is well-known that in the smooth case the latter rates are optimal in order, up to logarithmic in $n$ factors. Since the families of locally well-filtered signals are much wider than local Sobolev balls (smooth signals are trivial examples of modulated smooth signals!), it follows that the rates of convergence stated by Theorem 3 also are nearly optimal.

In this section we present the results of a small simulation study of the adaptive filtering algorithm applied to the 2-dimensional de-noising problem. The simulation setting is as follows: we consider real-valued signals

$e_{(1,1)},...,\,e_{(m,m)}$ being independent standard Gaussian random variables. The problem is to estimate, given observations $(y_{\tau})$, the values of the signal $(f_{x_{\tau}})$ on the grid $\Gamma_{m}=\{m^{-1}\tau,\;1\leq\tau_{1},\tau_{2}\leq m\}$, and $f({x_{\tau}})=s_{\tau}$. The value $m=128$ is common to all experiments.

We consider signals which are sums of three harmonic components:

the frequencies $\omega_{i}$ and the phase shifts $\theta_{i}$, $i=1,...,3$ are drawn randomly from the uniform distribution over, respectively, $[0,\omega_{\max}]^{3}$ and $[0,1]^{3}$ and the coefficient $\alpha$ is chosen to obtain the signal-to-noise ratio equal to one.

In the simulations we present here we compared the result of adaptive recovery with $T=10$ to that of a “standard nonparametric recovery”, i.e. the recovery by the locally linear estimator with square window. We have done $k=100$ independent runs for each of eight values of $\omega_{\max}$,

In Table 1 we summarize the results for the mean integrated squared error (MISE) of the estimation,

The observed phenomenon is rather expectable: for slowly oscillating signals the quality of the adaptive recovery is slightly worse than that of “standard recovery”, which are tuned for estimation of regular signals. When we rise the frequency of the signal components, the adaptive recovery stably outperforms the standard recovery. Finally, standard recovery is clearly unable to recover highly oscillating signals (cf Figures 1-4).

We denote $C(\mathbb{Z}^{d})$ the linear space of complex-valued fields over $\mathbb{Z}^{d}$. A field $r\in C(\mathbb{Z}^{d})$ with finitely many nonzero entries $r_{\tau}$ is called a filter. We use the commun notation $\Delta_{j}$, $j=1,...,d$, for the “basic shift operators” on $C(\mathbb{Z}^{d})$:

and denote $r(\Delta)x$ the output of a filter $r$, the input to the filter being a field $x\in C(\mathbb{Z}^{d})$, so that $(r(\Delta)x)_{t}=\sum\limits_{\tau}r_{\tau}x_{t-\tau}.$