Penalized wavelet estimation and robust denoising for irregular spaced data

Nonparametric univariate regression via wavelets is usually implemented under the assumptions of dyadic sample size, equally spaced fixed sample points, and i.i.d. normal errors. In this work, we propose, study and compare some wavelet based nonparametric estimation methods designed to recover a one-dimensional regression function for data that not necessary possess the above requirements. These methods use appropriate regularizations by penalizing the decomposition of the unknown regression function on a wavelet basis of functions evaluated on the sampling design. Exploiting the sparsity of wavelet decompositions for signals belonging to homogeneous Besov spaces, we use some efficient proximal gradient descent algorithms, available in recent literature, for computing the estimates with fast computation times. Our wavelet based procedures, in both the standard and the robust regression case have favorable theoretical properties, thanks in large part to the separability nature of the (non convex) regularization they are based on. We establish asymptotic global optimal rates of convergence under weak conditions. It is known that such rates are, in general, unattainable by smoothing splines or other linear nonparametric smoothers. Lastly, we present several experiments to examine the empirical performance of our procedures and their comparisons with other proposals available in the literature. An interesting regression analysis of some real data applications using these procedures unambiguously demonstrate their effectiveness.


Introduction
Wavelets based methods work quickly, efficiently, and attain fast rates of convergence in the problem of nonparametric univariate regression over equally spaced, fixed sample points on an interval, usually [0, 1] with a dyadic number of data points and independent noise: (1.1) In this model, the y i are the observed data, x i = (i − 1)/n are the sample points, and σ ε i is noise. The noise is assumed to be i.i.d. normal random variables in most cases. The function f is an unknown function of interest. We wish to estimate f and one can measure the performance of an estimatef by the expected global mean squared error, the goal being to construct estimates that have "small" risk R n . In order to have some meaningful estimate according to this criterion, one must assume certain regularity conditions on the unknown function f , such as f belongs to some Hölder classes, Sobolev classes, Besov classes and so forth. When the regression function f is sufficiently smooth, efficient smoothing methods such as kernel, splines and basis expansions have received considerable attention in the nonparametric literature (see, for example, Green and Silverman (1993), Eubank (1999) and Härdle (1990) and references given therein). In contrast, which is the case of this paper, when the unknown function f , mostly smooth, is suspected to have few discontinuities, sharp spikes and abrupt changes, wavelet methods are very popular. The application of wavelet theory to the field of statistical function estimation with dyadic, equispaced data, was pioneered in  and . The methodology includes a coherent set of procedures that are spatially adaptive and near optimal over a range of function spaces of inhomogeneous smoothness. Wavelet procedures achieve adaptivity through thresholding of the empirical wavelet coefficients. They enjoy excellent mean squared error properties when estimating functions that are only piecewise smooth and have near optimal convergence rates over large function classes. For example they attain optimal convergence rates of order log(n)n −2s/(2s+1) for the mean squared error risk when f is in a ball of a Besov space B s p,r with a smoothness parameter s, which can not be achieved by any linear estimator. For a thorough review of wavelet methods in statistics and an important list of references the reader is referred to Antoniadis (2007).
Despite their considerable advantages, however, standard wavelet procedures have limitations. It might be noticed that the vast majority of wavelet-based regression estimation have been conducted within the setting that the design points are fixed and equally spaced to enable the application of the discrete wavelet transform (DWT) to a compactly supported signal. Since this orthogonal transform is a matrix with dyadic dimensions, these methods have generally required a dyadic number of data points.
When the data does not meet these requirements, various modifications have been proposed. If n is not dyadic, it can be extended to a dyadic integer by reflection or periodizing. For data that is not equispaced, Cai and Brown (1998) investigated wavelet methods on samples with fixed designs via an approximation approach. They showed that applying the methods devised for equispaced data directly to nonequispaced data may lead to suboptimal estimators. They then proposed a method that was adaptive and near optimal, but used the assumptions of i.i.d. noise and dyadic number of sample points. Work has also been done with sample points that are not only nonequispaced, but random as well. Cai and Brown (1999) show that good rates can be obtained for certain random design schemes using the DWT.
We will deal with fixed, equispaced or nonequispaced, sampling designs of a size that is not supposed to be a power of 2. In the nonequispaced case, the design points x i are random design points drawn from a univariate distribution with a compactly supported density function f X , that satisfies certain conditions (essentially no flat parts in the density, or 'no holes' in the design). However, our wavelet-based regression procedures are defined using a standard least squares loss function which is unbounded and can be very sensitive to the presence of outlying observations. Moreover, it is challenging to detect outliers with standard wavelet-based regression procedures, especially when denoising functions with abrupt changes. Thus, robust estimators for such regression procedures are in need, for such cases. Therefore, a further relaxation of the constraints in (1.1) will be made, this time on the noise level assumptions by assuming that there are either sub-Gaussian or even, sometime, heavy-tailed, when a presence of aberrant observations (outliers) of the response variable is suspected. We will then also address robust estimation, meaning that our procedures remain valid even when there are aberrant observations of the response variable. The literature on this line of work, as far as we can tell, is not extensive. This paper is organized as follows. In Sect. 2 we describe the framework for the wavelet-based regression model with the basic concept of wavelets. A short description of the standard wavelet based regression procedures for the equispaced case via the DWT transform is also reviewed in this section. In Sect. 3 we present an aspect of wavelets described in Antoniadis and Fan (2001) crystallizing the penalized least squares approaches to wavelet nonparametric regression showing that they can be used to construct a set of basis functions over an arbitrary compact interval and that linear combinations of such basis functions are able to estimate particular, usually jagged, regression functions better than spline bases. A detailed description of the procedures, some completely new, and their asymptotic properties in the general fixed design case is presented in several subsections of this section. The problem of robust nonparametric wavelet smoothing is studied in Sect. 4. All procedures are investigated via various simulation settings and real data application examples in relevant subsections throughout the paper. Some concluding remarks are given in Sect. 5. Technical proofs are relegated to an "Appendix".
The R codes and their description implementing the methods in the paper together with the driver scripts for simulations, plots and the analysis of real examples are made available in a compressed archive as supplementary material.

Wavelet regression for equispaced data
Let us first consider a particular univariate regression problem stated as with a non-stochastic equidistant design t i = i/n, i = 1, . . . , n of size n = 2 J for some positive integer J , noise variables i that are i.i.d. Gaussian N (0, σ 2 ) and with a potentially non-smooth function f , defined on [0, 1], that may present a wide range of irregular effects. There is not any loss of generality to assume that μ = 1 0 f (x)dx = 0, since this parameter can always be estimated with a root-n rate by the empirical mean of the observations. We will therefore assume that the data has been centered.

A short background on wavelets
Wavelets denoising procedures have the ability to estimate inhomogeneous functions showing abrupt changes and discontinuities. Models for f in (2.1), that allow a wide range of irregular effects, are through the sequence space representation of (inhomogeneous) Besov spaces B s d,q ([0, 1]), 1 ≤ d ≤ ∞, 1 ≤ q ≤ ∞ (see e.g. ). To capture key characteristics of variations in f and to exploit its sparse wavelet coefficients representation, we will assume that f belongs to a "ball" of B s d,q ([0, 1]) of radius M with 1 ≤ q ≤ ∞ and s + 1/d − 1/2 > 0. The last condition ensures in particular that evaluation of f at a given point makes sense.
To deal with the estimation of f in model (2.1) several wavelet based estimation procedures have been developed in the literature (see, e.g. Antoniadis (2007) for a review) and we will summarize in this subsection the main ingredients most relevant to our work.
Let us adopt a vector-matrix form for model (2.1), After applying a linear and orthogonal wavelet transform, the discretized model becomes where Z = W n×n Y, γ = W n×n f and˜ = W n×n . The orthogonality of the DWT matrix W n×n ensures that the transformed noise vector˜ is still distributed as a Gaussian white noise with variance σ 2 I n . Hence, the representation of the model in the wavelet domain allows to exploit in an efficient way the sparsity of the wavelet coefficients γ in the representation of the nonparametric component including the case of hard sparsity, meaning that γ has exactly K * non-zero entries, or soft sparsity, based on the decay rate on the ordered entries of γ of wavelet coefficients is bounded above by say R. By modeling the unknown signal with its DWT, we over-parameterized the linear model (2.2) derived from the observations. The prior knowledge on the sparsity of the wavelet coefficients of f , is used to regularize efficiently the over-parameterized linear model by minimizing with respect to γ a separable penalized least squares expression: where, for λ > 0, p λ (| · |), is an appropriate penalty that preserves the sparsity (see e.g. Antoniadis and Fan (2001)). For example, using a SCAD penalty and an appropriate sequence of regularization parameters λ n one can show, thatf has an optimal asymptotic mean-squared error rate Thus, for regularly spaced data with n = 2 J , wavelet regularizations provide consistent nonparametric estimators with nearly optimal rates.

General fixed design case
Our task again is to estimate an unknown function f defined on [0, 1], given noisy data y i at sampling points x 1 , . . . , x n : We will now assume that the sampling points x 1 , . . . , x n are i.i.d realizations from an input distribution with a continuous density on [0, 1] which is bounded below by a constant b 0 > 0. However, the resulting x 1 , . . . , x n are regarded as deterministic points in [0, 1], so the only randomness involved is the sampling of the y i 's. When the underlying unknown signal f is smooth enough, several nonparametric approaches (local polynomials, kernel based procedures, spline smoothing or penalized splines) exist and their properties are well studied. However these methods are suboptimal when f belongs to some spaces of inhomogeneous smoothness such as the ball of the inhomogeneous Besov space B s d,q ([0, 1]) we used before. We shortly reviewed in the previous subsection some efficient wavelet based procedures for estimating f . However, the task is not so easy when n is not necessarily a power of 2 and especially when we are also dealing with irregularly spaced data which is the realistic case when sampling points occur from a general input density.
Using the same notation as before, and once an r -regular multi-resolution analysis of L 2 [0, 1] has been chosen, the unknown univariate function f has a derived wavelet basis expansion. For a given n, since f belongs to B s d,q (M), we can approximate the function f by its truncated expansion on the wavelet basis functions, say {W } : where K is an appropriate truncation index (depending on s) that is allowed to increase to infinity with n. Under our assumptions on f , the values of f on the design points are then well approximated by the above expansion and their estimation is therefore equivalent in estimating the (sparse) wavelet coefficient vector γ = (γ 1 , . . . , γ K ) T . As alluded to in Antoniadis and Fan (2001) and implemented in Wand and Ormerod (2011) we can also define the design matrices containing wavelet basis functions evaluated at the design points. Using this construction, we denote by W the corresponding wavelet regression n×K matrix of the wavelet basis functions evaluated at the observed sampling points, i.e.
Adopt again a vector-matrix form of the regression model, given by equation (2.1) to get: The function f is therefore characterised by the K -dimensional sparse coefficient vector γ . If the response Y is not centered, the intercept may be efficiently estimated by the empirical median or mean of the observations and therefore without any loss of generality we may assume that the columns of W are centered, which, in particular, means that the basis generating W doesn't involve scaling functions. Moreover, there is no loss of generality in assuming hereafter that W has its columns normalized, .i.e. W i 2 = 1, for i = 1, . . . , K .

Remark 1
To get an appropriate approximation of the function f and achieve optimal asymptotic rates for the resulting estimator, the index K will be allowed to be of the order O(n q ) with q ∈]0, 1]. Therefore in the linear regression model (3.1) the matrix W may have less rows than columns. In this case, where n K , problem (3.1) is severely underdetermined (and ill-posed), and hence it is challenging to obtain a meaningful solution. The sparsity of γ , looking for an approximative solution with many zero entries is then essential for estimating the unknown regression function.
Under the assumption that γ is sparse in the sense that the number K * of important nonzero elements of γ is small relative to n, we can estimate γ with several approaches proposed in the literature. Among them a popular method for realizing sparsity constraints is the Lasso (Tibshirani 1996). It leads to a convex but nonsmooth optimization problem: where · 1 denotes the 1 -norm of a vector, and λ > 0 is a regularization parameter. Since its introduction (Chen et al. 1998;Tibshirani 1996), problem (3.2) has gained immense popularity in many diverse disciplines, which can largely be attributed to the fact that it admits efficient numerical solution. Under certain regularity conditions (e.g., level wise restricted isometry property and level wise restricted eigenvalue condition) on the design matrix W and the sparsity level of the true vector γ , it can produce models with good estimation and prediction accuracy [see e.g. Zhao and Yu (2006) and Wainwright (2009)]. However, Lasso tends to over-shrink large coefficients, which leads to biased estimates (Antoniadis and Fan 2001;Fan and Li 2001;Fan and Peng 2004) and this implies that there does not exist a sequence λ n of the tuning parameter λ which can lead to both variable selection consistency and asymptotic normality [see Fan and Li (2001) and Zhao and Yu (2006)]. In other words this means that the Lasso estimates do not have the oracle property even when the minimum of the nonzero coefficients of γ is bounded below. A way to improve variable selection accuracy and gain oracle properties is by reducing the bias of Lasso via the adaptive Lasso procedure (Zou 2006) which solves the following weighted l 1 regularization problem for some α ∈ (0, 1) :γ whereŵ is an estimator of γ (for example, the solution of the standard unweighted Lasso with regularisation parameter λ n or a ridge estimate of γ ). A low-dimensional analysis in Zou (2006) shows that the adaptive Lasso solution can achieve the oracle property asymptotically. A high dimensional analysis of this procedure was given in Huang et al. (2008). For variable selection consistency and oracle properties to hold, the adaptive lasso requires strong non optimal conditions in terms of the minimal strength of nonzero components of the true wavelet coefficients γ and some extra conditions on the design matrix.

Remark 2
The restricted isometry condition (RIP) has been introduced first in the compressive sensing literature to analyze 1 -regularized recovery of a sparse β from its random projection X β with i.i.d. N (0, 1) entries in X . While conditions for RIP are specialized to hold for random designs with specific covariance matrices, related conditions using sparse eigenvalues can be defined (for example, the sparse Riesz condition (SRC)) and lead to 2 -norm estimation error bounds of the optimal order for gives the necessary and sufficient conditions on a penalty function for the solution of the penalised least-squares problem to be thresholding, continuous and approximately unbiased for large values of their argument. In the subsections that follow we will therefore adopt such penalties.

A non-convex wavelet based approach
Our nonconvex approach for wavelet based denoising leads to the following optimization problem: where ρ λ,τ is a nonconvex penalty, λ > 0 is a regularization parameter, and τ controls the degree of concavity of the penalty. We already mentioned important examples of such penalties in this section. The nonconvexity and nonsmoothess of the penalties ρ λ,τ pose some challenge for the optimization problem (3.3), but their attractive theoretical properties and empirical successes justify our choice. We focus on the most commonly used nonconvex penalties, i.e., SCAD, MCP, and capped 1 which is a linear approximation of the SCAD penalty [but similar results can be derived for some of the nonquadratic penalties discussed in Antoniadis et al. (2011)], for recovering the unknown function f ; see Table 1 for explicit formulas (and the associated thresholding operators, to be defined later).
It is worth mentioning here that the penalties ρ λ,τ (t) in the above table are all (μ, ξ )-amenable (see Loh and Wainwright (2015)) as soon as τ > 2. Recall in particular, that if ρ λ,τ (t) is (μ, ξ )-amenable, then q λ,τ (t) := λ|t| − ρ λ,τ (t) is everywhere differentiable. Defining the vector version q λ,τ : IR K → IR accordingly, it is easy to see that μ 2 γ 2 2 − q λ,τ (γ ) is convex. A consequence of the above facts is that for any of the nonconvex penalties in Table 1, there exists at least one global minimizer to problem (3.3). This essentially follows from the boundedness of L from below by zero, and by the lower semi-continuity of each of the penalties.
We now recall the expressions of thresholding operators for the penalties in Table 1 which form the basis of many existing algorithms, e.g., coordinate descent, iterative thresholding and first-order proximal mappings, which may be found in the literature. We recall also here an unified derivation of the thresholding operators associated to these penalties. For any of the penalties ρ(t) in Table 1 (the subscripts λ and τ are omitted for simplicity), define a function g(t) It is easy to show that T * = inf t>0 g(t) is attained at some point t * ≥ 0 and that one has (t * , T * ) = (0, λ) for the capped-1 , SCAD and MCP.  Antoniadis and Fan (2001) The thresholding operators S ρ are defined by and can potentially be set-valued. Explicit expressions of the thresholding operators S ρ associated with the three nonconvex penalties (capped-1 , SCAD and MCP) are as given in Table 1. Note that they are singled-valued, except at v = λ(τ + 1 2 ) for the capped-1 penalty.
In order to obtain stationary points of the objective function (3.3), we can use the composite gradient descent algorithm (see Nesterov (2007)) and the fact that by our assumptions on f , the vector γ has a bounded 1 norm by M. Denotinḡ Then the composite gradient iterates are given by where η is a step-size parameter. Defining the soft-thresholding operator S λ/η (β) component-wise according to a simple calculation shows that the iterates (3.5) take the form Theorem 3 of Loh (2017) guarantees that the composite gradient descent algorithm will converge at a linear rate to the global minimizer as long as the initial point γ 0 is chosen close enough to the true γ , since L has an almost-everywhere bounded second derivative and since q λ,τ is convex, as is the case for the penalties in Table 1.
We may now examine the finite sample properties and the empirical quadratic risk of the univariate estimatorf := Wγ , whereγ is a stationary point of the objective function (3.3), close to the true wavelet coefficient vector γ . Our first result below establishes near minimax convergence rates for the quadratic risk of the estimator. Its proof is given in the "Appendix". for the n × K -wavelet design matrix W associated to an r -regular orthogonal wavelet basis with r > max{1, s}. Further, . Assume that f 0 ∈ B s d,q (M) and that the mother wavelet ψ, has r null moments and r continuous derivatives with r > max{1, s}. Suppose λ ≥ c 1 η 2 + 2 log K for some η > 0 and some constant c 1 > 0 that may depend on M and the distribution of i . Then, for sufficiently large K (specifically K ≥ c 1 n 1/(2s+1) ), with probability at least 1 − 2 exp(−η 2 /2), we have

An OEM perspective
We consider here a similar setting as the one in Antoniadis and Fan (2001). Since n is not necessarily a power of 2, and the design is not regular, we start by embedding the observed data points {x 1 , . . . , x n } into a fine equispaced grid by taking N = 2 J dyadic grid points with J is a chosen resolution, which we denote t , = 1, . . . , N (since we are working on the unit interval [0, 1] this would mean taking t = /N ). The number N should be much larger than the sample size n (i.e N n). Embedding here means that we are approximating the set of sampling points {x 1 , . . . , x n } with a subset of dyadic points at the scale J , that is we assume that x i ≈ n i /N for i = 1, . . . , N for some 1 ≤ n i ≤ N (i.e. n i /N is closest design point to the point x i ). It is easy to see that when N is much larger than n, and given that the sampling points x i are i.i.d. realizations of a probability distribution on [0, 1] with a density that is continuous and bounded below by a strict positive constant, the approximation errors by moving non dyadic points to dyadic points is negligible. Letf N be the underlying signal f collected at all dyadic points { /N , = 1, . . . , N }, that isf N = ( f (t 1 ), . . . , f (t N )) T and let this time W be a given N × N orthogonal discrete wavelet transform (DWT). The wavelet transform off N is the vector β of discrete wavelet coefficients given by β = Wf N . By orthonormality of the discrete wavelet transform we havef N = W β.
Let S be the n × N sampling matrix which extracts the [n 1 , . . . , n n ]th entry from f N , then f = Sf N is the discrete signal at the sampling points. Denoting by A = SW the observed data, up to a negligible approximation error, obeys the linear model where is the noise vector which again is assumed to be independent of the entries of A. The linear model derived from the observations is over-parameterized, but we have already noticed that for functions belonging to the Besov bodies adopted in this work, their wavelet coefficients are sparse. We can therefore estimate β by regularization.
We have already concisely reviewed in our previous subsection a coordinate descent Newton-type method that solves problem (3.6). In a similar spirit, for penalty functions such as the capped-1 , the MCP and the SCAD, Xiong et al. (2017) proposed a novel orthogonalizing expectation maximization (OEM) algorithm to solve the least squares problems in linear regression when the regularizers are separable over the variables.
Inspired by Daubechies et al. (2004) where the idea below was mentioned for the first time and by a recent work of Liu et al. (2020), we suggest a new perspective of OEM algorithm with emphasis on the capability of separating variables to avoid the computation of matrix inverse. With this, an algorithm is proposed to solve problem (3.6) and to apply it for wavelet denoising. The advantage again is that the minimization of (3.6) will be carried out by decomposing it as a sequence of univariate minimizations that then can be efficiently computed via the univariate thresholding operators reviewed above.
For a problem in (3.6) with a separable penalty P λ,τ (β) = 2 K i=1 ρ λ,τ (β i ) denote byβ the final solution and introduce a function where F = d I p − A A and d is a positive integer such that F is positive semidefinite. Quite clearly, RSS(β) can also be minimized atβ. On the other hand, by some simple algebraic computation, it is equally clear that, up to some additive constant whereβ i is the ith element of the dualβ = A T Y/d + Fβ/d. This gives rise to the following relationship: which involves only a single variable and can be done by applying component wise in an iterative manner the thresholding operator (3.4). From the discussion from Xiong et al. (2017), any choice of d ≥ σ max A A will make the algorithm work in the sense that, under our assumptions on f and the regularity of the wavelet basis used in constructing W, the thresholding iterations converge again to a stationary point of problem (3.6). We may, again, state the following proposition (very similar to Proposition 1 stated before), whose proof is omitted since it is a direct application of Theorem 4 of Antoniadis and Fan (2001). We may conclude that both versions of the univariate wavelet based estimators studied above, under appropriate assumptions, present optimal prediction performance. The next Sect. 3.3 reports an numerical comparison of their empirical performances, compared also to another recently developed wavelet based estimator, called Wavemesh (see Haris et al. (2018)).

Experimental comparison setup
The simulation setup to illustrate the above univariate wavelet based denoising procedures is as follows. Three testing functions, namely, Heavisine, Blip and Corner are used that represent a variety of function characteristics (see Fig. 1). A relatively low level of signal-to-noise ratio (SNR) equal to 4 is chosen, where SNR is defined as SNR=std( f )/σ . Two different sample sizes n = 100 and n = 300 are used. Nonequispaced sample points x i , i = 1, . . . , n, are generated according to, either a uniform distribution on [0, 1] (case (I)) or a Gaussian distribution N (2, 1) whose values are where i 's are i.i.d. standard normal distributed and σ is chosen to achieve a signalto-noise ratio S N R = 4. For each simulated data set Y h , we applied each of the three methods of this subsection to estimate the test signal with optimal tuning parameters chosen from pre-specified grids. The root mean square error (RMSE) of the solutionf and the true signal f at the nonequispaced design points x i is used to assess the quality of the estimates.
Boxplots of the RMSEs of the estimates obtained by the three methods for the design case (I) are presented in Figs. 2 and 3 for Heavisine, Blip, and Corner, respectively, for each sample size used.
When the sample size is small (n = 100) our proximal method is almost equivalent to the OEM based method and slightly worse than wavemesh. This is due to the fact that with such a sample size both proximal and OEM tend to over smooth the cusps and discontinuities present in the signals especially for the Heavisine function. The behaviour of the OEM based method is also expected since it is very close to the ROSE one-step estimator of Antoniadis and Fan (2001), whose optimal properties are mainly asymptotic. When n = 300, all three methods have an excellent behavior. Contrary to the small sample case, our proximal estimator performs better than the two other methods for the signal Heavisine and Blip and is slightly worse than wavemesh for the signal Corner.
When the design grid is generated according to case (II), each of three methods of this subsection to estimate the test signal with optimal tuning parameters chosen from pre-specified grids is less sensitive to boundary effects since the tails of data within the interval [0, 1] are lighter by design.
Boxplots of the RMSEs of the estimates obtained by the three methods for the design case (II) are presented in Figs. 4 and 5 for Heavisine, Blip, and Corner, respectively, for each sample size used. When the sample size is small (n = 100) our proximal method is better than the the two other methods.
When n = 300, contrary to the small sample case, our proximal estimator performs always much better than the two other methods and wavemesh is better than OEM. It is worth to outlining that, in the average, the fastest procedure in CPU time is OEM followed by wavemesh and then proximal.

Robust wavelet smoothing
The nonparametric univariate regression procedures developed and studied in our previous subsections are well adapted to the cases of not uniformly spaced design points, sample sizes that are not a power of 2 and a standard Gaussian or Gaussian alike noise term, and mainly derived as penalized regression estimations using a penalized least squares loss function.The attraction of such procedures is that the unknown "inhomogeneous" regression function can then be estimated with optimal convergence rates. In practice however, some extreme observations may occur, and estimation using an unbounded loss function suffers from a lack of robustness, meaning that the estimated functions can be distorted by the outliers. Both the nonparametric function estimates themselves and the choice of the penalization parameters associated to them are affected. Indeed, when noise has a non-normal distribution (or a sub-Gaussian one as in Proposition 1), for instance, heavy-tailed distribution, even for data observed on a regular grid, the classical wavelet shrinkage is not efficient for estimating the true function. We may mention here some robust wavelet procedures that have been previously studied in the equidistant dyadic sample case. Bruce et al. (1994) proposed a robust wavelet transform to estimate wavelet coefficients at each resolution level which is robust to outliers. Kovac and Silverman (2000) developed a direct approach. Their approach can be summarized as identifying outliers with a classical robustness test, removing the outliers, and applying their thresholding procedure for unequally spaced data. However, their approach has a drawback in that information is lost by removing data. Sardy et al. (2001) proposed a robust wavelet de-noising estimator based on a robust M-type loss function. Their procedure is computationally intensive as it involves solving a non trivial nonlinear optimization problem. Averkamp and Houdré (2003) extend the minimax theory for wavelet thresholding to some broader classes of symmetric heavy-tailed distributions, but is only applicable if the distribution is known. Oh et al. (2007) studied a robust wavelet shrinkage algorithm that is based on the concept of pseudo data. However the selection of a robust threshold and the asymptotic properties of the resulting estimators were not addressed in their paper.
To address the lack of robustness, we introduce a new class of wavelet based robust estimators for performing estimation in sparse "inhomogeneous" univariate regression models. Our emphasis is on M-type penalized smoothing in a spirit much similar to the one adopted in our previous subsections. The class of estimators derived here are minimizers of a criterion that balances fidelity of a wavelet basis expansion to the data with sparseness/smoothness in the estimate. We will use for such a task, the methods that have been developed (Amato et al. 2021), which under high-dimensional sparse linear models, and for appropriate robust loss functions, are able in handling optimal, both computationally and asymptotically, predictions.
We will focus again at the linear wavelet regression model, stated, with a negligible approximation error, in (3.1): is the vector of responses, W ∈ IR n×K is the wavelet basis regression design matrix, and γ ∈ IR K is the vector of wavelet regression coefficients. The difference here is that we assume that the noise vector = ( 1 , ..., n ) T ∈ IR n is having i.i.d components with a symmetric around zero distribution F = F with a finite second moment which deviates from the Gaussian distribution, in the sense that it is a more or less heavy-tailed distribution producing outliers in the responses. As before the design points may be either deterministic or random. Therefore the parameter of interest Wγ in (3.1) is the mean of Y in the deterministic case or the conditional mean of Y on X in the random design case.
We will assume again that the unknown regression f to be recovered belongs to a Besov body B s d,q (M) which implies in particular that the weak sparsity K = K n , while allowed to tend to infinity with n is such that n ≥ c 0 K n log n, for a sufficiently large constant c 0 . By known information-theoretic results of Wainwright (2009), this type of lower bound is required for any method to recover the support of a K -sparse signal, hence is not a limiting restriction.
Since the error distribution F deviates from the normal distribution, the standard unpenalized squared error loss function used in the previous subsections is typically changed to the log-likelihood − log f , where f is the density of F , leading to penalised log-likelihood estimators. Unfortunately, in real life situations the error distribution F is unknown and methods that adapt to many different distributions are needed. A large class of penalized M-estimators has been recently proposed by Amato et al. (2021). The methods in this class penalise an empirical loss function L n in the following mannerγ for a chosen scaled loss and for a suitable penalty function p λ . The side condition γ 1 ≤ R for some positive number R guarantees the existence of local/global optima in cases when the loss and the penalty may be nonconvex. Such a condition is satisfied here by the fact that f ∈ B s p,q (M) so the true regression vector γ is feasible for the optimisation problem (4.1).
The loss functions that we will adopt in this subsection include Tukey's biweight and Welsh's loss (see Amato et al. (2021), Appendix 1). One could also use some other popular losses such as the least absolute loss or a Huber loss, but their comparison with Tukey's biweight and Welsh's loss made in Amato et al. (2021), favor the latter. An unifying property of these loss functions is that the derivative is bounded and odd in each case which turns out to be an important property for robustness of the resulting estimator because it is lessening the contamination effect of gross outliers.
Since for the class R of robust loss functions that we consider also have (u) > −κ for all u ∈ IR and some constant κ > 0, and the density of the design is continuous and bounded below by a strict positive constant, the local RSC condition, already discussed in our Remark 2, for these loss functions holds with probability at least 1 − c exp(−c 2 log K n ) for some constant c > 0 as soon as K n tends to infinity and n > c 0 K n log n. Since the penalties used are μ-amenable, by Theorem 4.1 of Amato et al. (2021) the resulting wavelet regression estimator is asymptotically consistent. Moreover, the estimate can be again computed with a variant of the numerical algorithms described in Amato et al. (2021).

Some supplementary simulations
This subsection reports the results from a small simulation study designed to assess the practical performance of the two robust methods proposed above and compare them to the robust wavelet shrinkage method of Oh et al. (2007). Since the Oh procedure is only designed for a dyadic equispaced grid the simulations below are designed with an equispaced deterministic design with a sample size equal to 256. Note however, that our robust methods can be applied to any type of design which is not true for Oh et al's method.
The experimental set-up is as follows. The three test functions 'Heavisine', 'Blip' and 'Corner' that we have also used before were selected for this illustration. A sample size, n = 256, a dimension K = n 2 and a contaminated Gaussian mixture, 0.9N(0, σ )+ 0 .1N (0, 4σ )), as a noise distribution, were considered with a value of σ such that, when no outliers are present, a SNR of 4 is realized. Note here that such a distribution for the noise covers the assumptions of Theorem 4.1 of Amato et al. (2021). For each combination of test function and noise, 50 samples were generated. Then, for each generated sample, three regression estimates were obtained by applying, respectively, the robust estimator based on Tukey's Biweight loss function and penalized with an MCP penalty, the robust estimator based on Welsh's loss function and penalized with an MCP penalty and the robust method by Oh et al.. To scale the loss functions an estimate of the noise scaleσ 2 adopting a good robust estimate of the noise variance (see Amato et al. (2021)) is used. The optimal penalty is chosen over a fine grid of values. For each value of the penalty parameter on the grid, a robust model selection , was calculated for each regression estimate. Figures 6 and 7 display some typical fits obtained by each of these methods for each of the test functions. The boxplots of the resulting RMSE for each scenario are displayed in Fig. 8.
The following major observations can be made. First, the proposed procedures, while less variable, are slightly worse than the Oh method for the Heavisine example, but remain very comparable. However, they both outperform the Oh method for the other two test functions. A likely explanation is that the Heavisine function has a less sparse representation in the wavelet domain compared to Blip and Corner. Having in mind that the other robust estimates can handle non dyadic and non regular designs they can be therefore considered superior.

Real data examples
The proposed univariate wavelet procedures were first applied to the glint dataset analyzed by Sardy et al. (2001). The data are radar glint observations from a target captured at N = 512 angles of an object in degrees. The signal contains a number of glint spikes, causing the apparent signal to behave erratically. From physical considerations, a good model for the true signal is a frequency oscillation about 0 • . The data together with the regression estimates are displayed in Fig. 9 We have also applied the proposed methods to the balloon data used by Kovac and Silverman (2000). The data are radiation measurements from the sun, performed from a flight of a weather balloon. Since the measurement device was occasionally cut off from the sun, individual outliers and large outlier patches are introduced. We select every second or third observation from the 4984 observations to reduce the sample size to 2048. Although the data are very closely located to a regular grid, they are not exactly equally spaced. Therefore our wavelet procedures on an irregular design are appropriate to fit the data. Figure 10 shows the balloon data with a sample size of 2048 points and the estimates by a moving median filter with window size 14, the non robust MCP penalized wavelet procedure and the robust fitting based on MCP BiWeight, respectively. Note that the window size of the median filter is selected by eye to be "good". Increasing the window size does not help capturing the valley around 0.77 and decreasing the size makes the fit quite wiggly. The fit by the penalized MCP is still affected by the outliers. On the other hand, the robust fitting based MCP BiWeight reacts to the major pattern of the observations. Finally we show two use cases arising in the Horizon 2020 MADEin4 project. The project aims at developing next generation metrology tools, machine learning methods and applications in support of Industry 4.0 high volume manufacturing in the semiconductor manufacturing industry. It is headed by STMicroelectronics, leader company in semiconductor industry.
The first use case concerns monitoring some moving holders that block silicon wafers into an equipment during the growth of semiconductors in a particular process step for semiconductor device fabrication, with the objective of predicting their malfunctioning. Monitoring occurs by collecting images of the holders and processing them. As result of this analysis, an indicator has been devised that measures how much an image is different from a reference one. Presence of artifacts on the holders generated during the process unduly increase the error indicator. Therefore it is needed to track behaviour of the indicator as a function of time, cleaned by these artifacts. Since measurements are not taken every day, a nonequispaced time series is generated. Figure 11 shows the sequence of data points and the fit obtained by a Median filter and All fits show to be unaffected by the large outlier with indicator above 9, with Welsh loss function yielding back a smoother function.
In the second use case on data provided by STMicroelectronics, physical parameters were measured and collected over time, during another particular process step of semiconductor device fabrication. Data have to be further processed in a Virtual Metrology setting. The generated time sequences are not equispaced in time. It is interesting to fit such sequences, removing noise eventually, and approximating on an equispaced grid, introducing least correlation than interpolation, so to have more flexibility in processing them and extracting specific Features. Figure 12 shows data for a typical sequence of one specific parameter. The length of the sequence is 1714 data. The figure also shows the fits obtained by Penalized wavelet estimation (MCP) and robust Wavelet fitting (with BiWeight and Welsh loss function based on MCP).
Indeed the data do not show outliers, therefore a nonrobust approach seems more suited to represent the data with a greater accuracy. The results from all these real examples suggest that our methods possess promising empirical properties.

Conclusions and discussion
This paper introduced some wavelet-based method for nonparametric estimation of general univariate nonlinear regression models. The nonlinear regression functions in such models present a wide range of irregular effects in terms of smoothness and are approximated by appropriate truncated expansions on wavelet basis functions evaluated at the design sampling points. Unlike traditional methods the methods developed here do not require that the design points are uniformly spaced on the unit interval, nor do they require that the sample size is a power of 2. Adopting a matrix notation, we provide a unifying framework for achieving the estimation of the unknown regression functions by adapting general algorithms for penalized linear models and provide a solid theoretical framework for establishing mean-squared estimation consistency. Although in the paper we restrict to a univariate regression setting, the methodology presented can be extended to more flexible settings, which could be needed when dealing with more complex data, such as, for example, additive nonlinear regression models with irregular components. We would like to shortly discuss here a possible way to achieve it. To fix the ideas, assume that we now have at our disposal observations (y i , x i ), i = 1, . . . , n, from an additive regression model of the form Y i = μ + p j=1 f j (X j i )+ε i , i = 1, . . . , n, where μ ∈ IR is an overall mean parameter, each f j is a univariate function, a mean-zero error vector ε = (ε 1 , . . . , ε n ) T with independent components which is also independent of the covariates. Each of the p univariate functions f j could then be approximated by their approximation in a wavelet basis as in this paper. By stacking all involved wavelet basis matrices in a big design matrix, the problem to be studied can be written as for the linear setting used in our paper and estimation could follow by adapting our methods using grouped penalization methods and extensions of the definition of proximal operators to vector-valued arguments. We postpone a detailed study of such an extension for future research.

Supplementary Information
The online version contains supplementary material available at https://doi. org/10.1007/s00180-021-01174-4. coefficients f 0 . A nice feature of the above optimization problem is that it is exactly the lasso problem with design matrix W and loss function the convex functionalL.
The proof is completed by noting that for functions within the assumed Besov body the approximation error 1 n f 0 − Wγ 2 2 is of order O(2 −K s ).