Oracle Inequalities

and Nonparametric Function Estimation

Iain M. Johnstone

Abstract. In non-parametric function estimation, partial prior information about the unknown function is often expressed by a family of models or estimators, among which a choice must be made. Oracle inequalities bound the mean squared error of a given estimator in terms of the (unknowable) best possible choice of model for the unknown function. This survey concentrates on three examples: the James Stein estimator, soft thresholding, and complexity penalized least squares and as illustrations, we describe some consequences for adaptive estimation.

1991 Mathematics Subject Classification: 62G07, 62G20, 62C20

Keywords and Phrases: adaptive estimation, complexity penalty, James-Stein estimator, minimax estimation, thresholding, unconditional basis, wavelet shrinkage

1 Introduction

Statistical theory aims in part to articulate when and why certain applied methods of data analysis succeed. With emergence of large, often instrumentally acquired datasets, recent decades have seen a focus on “nonparametric” models in which the number of model parameters grows with the size of available data. Here we focus on the estimation (or “recovery” or “denoising”) of functions observed in additive noise and describe some relatively simple inequalities that encode information on the effect of sparse representation on the quality of estimation.

A common caricature is to posit observed data y ∈ R^{n} with structure y = μ + ϵz. Here μ is an unknown function which one desires to “estimate” or “recover”, and z ∈ R^{n} is a vector of standard Gaussian noise of known scale ϵ. When expressed in terms of coefficients in an orthonormal basis {ψ_{i}, I ∈ I_{n}}, the model becomes

y_{I} = μ_{I} + ϵz_{I}  I ∈ I_{n}.  (1)

Here z_{I} are independent Gaussian noises of mean zero and variance one. This sequence form of the “Gaussian white noise model”, whether finite as here, or infinite, as in Section 1.1 below, is the conceptually and technically simplest model of nonparametric estimation. Extensions to correlated noise and indirect data y = Kμ + ϵz are possible, but not covered here.

Examples 1. (a) The equispaced, fixed design in which y_{I} = f(t_{I}) + σz_{I}, with f an unknown function defined on [0, 1] and t_{I} = I/n, I = 1, … n.

(b) An initial segment of the continuous Gaussian white noise model. Suppose that W_{t} is a standard Brownian motion (or sheet) and that one observes dY_{t} = f(t)dt + ϵdW_{t} for t ∈ D, a compact set in ℝ^{d}. If d = 1 and D = [0, 1], this may be interpreted as Y_{t} = ∫_{0}^{t} f(s)ds + ϵW_{t} for 0 ≤ t ≤ 1. Take inner products with elements {ψ_{I}, I ∈ I} of a complete orthonormal basis for L_{2}[0, 1] and set y_{I} = ⟨ψ_{I}, dY ⟩, θ_{I} = ⟨ψ_{I}, f⟩ and z_{I} = ⟨ψ_{I}, dW ⟩. This gives an infinite sequence version of (1), to be used in Section 1.1 below. To recover precisely (1), consider an initial segment of cardinality n of the index set I. A discrete orthogonal wavelet transform of model (i) yields an approximation to this initial segment (after calibrating ϵ = σn^{−1/2}, cf. [10]).

(c) Redundant regression. Suppose that there are given vectors (or signals) x_{1}, … x_{p} ∈ ℝ^{n}, and that it is thought useful, for reaons of parsimony, interpretability or otherwise, to represent μ in terms a few of the x_{i}. Collecting x_{i} as columns of a “design matrix” X = [x_{1} ⋯ x_{n}], one obtains the standard, homoscedastic Gaussian linear regression model y = Xβ + ϵz. In traditional parametric regression analysis, it is supposed that p < n and that μ ∈ span{x_{i}}. However, we specifically consider two “non-parametric” cases: a) p = n and x_{i} orthogonal (i.e. equivalent to (1)), and b) p > n and not orthogonal - here the x_{i} might be a class of basic signals from a (possibly highly) redundant dictionary D and we seek a parsimonious representation of μ in terms of as few elements of D as possible.

Assessing error. An estimator μ̂ = μ̂(y) is a function of observed data y: we wish to quantify and compare the quality of estimation as μ̂ varies. Simplest to work with is mean squared error (MSE):

r_{ϵ}(μ̂, μ) = E_{μ}|μ̂ − μ|^{2} = ∫ |μ̂ − μ|^{2}ϕ_{ϵ}(y − μ)dy.  (2)

Here ϕ_{ϵ}(z) denotes the probability density function of ϵz. The notation r_{ϵ}(μ̂, μ), mnemonic for “risk”, hints at the possible and frequently desirable use of more general error norms ∥μ̂ − μ∥ or loss functions L(μ̂, μ).

The error μ̂−μ is usually decomposed into a zero-mean stochastic component μ̂ − E_{μ}μ̂ and a deterministic component, the bias, E_{μ}μ̂ − μ. For quadratic error measures, these components are uncorrelated, so that the MSE is the sum of variance and squared bias terms. In particular, for a linear estimator μ̂_{L}(y) = Ly,

r(μ̂_{L}, μ) = ϵ^{2}tr LL^{t} + |Lμ − μ|^{2}.  (3)

The quality of approximation of μ by the operator L is thus balanced against the complexity of L, as measured by the variance term, which for example becomes ϵ^{2}m in the case of orthogonal projection onto a subspace of dimension m. Already visible here is the important role that approximation theory plays in analysing the deterministic component of error. For non-linear estimators that, implicitly or explicitly, involve a choice among linear estimators, the analysis of the stochastic term is facilitated by the concentration of measure phenomenon (Section 4).

Models and Estimators. A model is a subset M of the full parameter space ℝ^{n}. A family of models {M_{α}, α ∈ A} is one device commonly used to represent imperfect and partial information about the unknown μ. Often there is a natural estimator μ̂_{α} associated with each model and in this paper we simplistically conflate choice of model with choice of the associated estimator.

Examples 2. (a) Spheres and linear shrinkage. For positive α, let M_{α} be the sphere |μ| = α: this might correspond to prior information about the signalto-noise ratio. Natural corresponding estimators are given by linear shrinkage: μ̂_{α} = γy where γ = γ(α) is obtained by minimizing the MSE in (3), namely nγ^{2} + (1 − γ)^{2}|μ|^{2}, on M_{α} to obtain the Wiener filter γ(α) = α^{2}/(n + α^{2}) ∈ (0, 1).

(b) Subspaces and projections. In the regression setting of Example 1(c) above, to each subset J ⊂ {1, …, p} of the full variable list is associated a linear model M_{J} = span {x_{j}, j ∈ J}. The corresponding estimators are orthogonal projections P_{J} on M_{J}: these are the least squares estimators on the assumption that μ ∈ M_{J}.

Ideal Risk Given a family A of models (or corresponding estimators), and for a given unknown μ, the best attainable MSE is given by the ideal risk

R_{ϵ}(μ, A) = inf_{α} R(μ̂_{α}, μ).

Thus, in example (a), the ideal linear shrinkage risk is

R_{ϵ}(μ, LS) = nϵ^{2}|μ|^{2}/(nϵ^{2} + |μ^{2}|).  (4)

Outline of paper. Of course, μ is not known, and without access to an oracle who divulges the best α, the ideal risk is not attainable by an estimator depending on the data y alone. Nevertheless, it acts as a useful benchmark, and we seek estimators that in an appropriate sense optimally mimick the ideal risk. Such estimators turn out to be non-linear, and in particular, not members of the family μ̂_{α}. For three settings and estimators, oracle inequalities are presented in Theorems 3, 5 and 8 – we emphasize that the inequalities are non-asymptotic and uniform in character, holding for all n, ϵ and for all μ ∈ ℝ^{n}.

Corresponding lower bounds (although asymptotic in n) show that without some restriction on, or further information about μ, the inequalities cannot be improved, and thus represent in some sense the necessary “price” for searching over a class of models/estimators of a given size.

Oracle inequalities are neither the beginning nor the end of a theory, but when available, are informative tools. For example, Theorems 3, 5 and 8 may also be used to derive asymptotic (i.e. low noise ϵ) results within a framework of adaptive minimax estimation: this class of applications is considered in a connected sequence of “illustrations” in the continuous Gaussian white noise model, which we now introduce.

1.1 Illustration: Asymptotic Minimax Estimation.

The continuous Gaussian white noise model is that of Example 1(b). Because of Parseval’s inequality ∫_{0}^{1}(f̂ − f)^{2} = ∑_{I}(θ̂_{I} − θ_{I})^{2} = ∥θ̂ − θ∥^{2}, estimation error can equally well be measured in the sequence domain. To evaluate estimators, we use the minimax principle - although inherently conservative and not universally accepted, we find that it leads to clear structures and informative results. Thus, estimators are assessed by their worst case risk over a given Θ. The minimax risk measures the best attainable such maximum risk, within a class E of estimators: R_{E}(Θ, ϵ) = inf_{θ̂∈E} sup_{θ∈Θ} r_{ϵ}(θ̂, θ). The symbols E = N, L, D, … refer to specific estimator classes: all non-linear, all linear, all threshold rules etc. Finally, estimator θ̂ is called asymptotically E− minimax if

sup_{θ∈Θ} r_{ϵ}(θ̂, θ) = R_{E}(Θ, ϵ)(1 + o(1)), ϵ → 0.

In order to describe a flexible and scientifically meaningful class of parameter spaces Θ, we employ a dyadic sequence notation, in which I = (j, k), with j = 0, 1, … and k = 1, …, 2^{j}. The primary motivation comes from orthonormal bases of wavelets {ψ_{jk}}, which, under suitable regularity and decay conditions on the wavelets, and with suitable modifications to handle intervals, form unconditional bases for many function spaces of interest ([22, 15, 3]). Thus their norms may be characterized in terms of conditions on |θ_{I}| . For example, let χ_{I} denote the indicator function of the interval [(k − 1)2^{−j}, k2^{−j}]: the sequence of (quasi-)norms ∥ · ∥_{α,p}, defined for 0 < α < ∞, 0 < p ≤ ∞ by

∥θ∥_{α,p}^{p} = ∫_{0}^{1}[∑_{I}(2^{aj}|θ_{I}|χ_{I})^{2}]^{p/2}, a = α + 1/2,

are equivalent, (for p > 1 and α ∈ ℕ) to the traditional Sobolev norms ∥f∥_{W_{p}^{α}}^{p} =

∫_{0}^{1} |f^{(α)}|^{p} + |f|^{p}. In the Hilbertian case p = 2, these take the simpler form

∥θ∥_{α,2}^{2} = ∑_{j≥0} 2^{jα}|θ_{j}|^{2}, |θ_{j}|^{2} = ∑_{k} |θ_{jk}|^{2}.

As parameter spaces, we thus use norm balls: Θ_{α,p}(C) = {(θ_{I}) : ∥θ∥_{α,p} ≤ C}, which are analogs of size restrictions on derivatives, but measured in L_{p} norms.

In practice, the values of (α, p, C) will not be known, and rather than seeking a minimax estimator for a single such Θ_{α,p}(C), we look for estimates with an adaptive minimaxity property. Thus, suppose that a scale of spaces S = {Θ_{ν}(C) : ν ∈ V, C > 0} is given, where ν is an order parameter, such as (α, p) above, and C a scale parameter. Then θ̂ is adaptively E−minimax if (i) the definition of θ̂ is independent of (ν, C), and (ii) θ̂ is asymptotically E−minimax for all (ν, C).

2 Linear Shrinkage and Orthogonal Invariance

A celebrated result in parametric statistics, due to Stein [24], is the inadmissibility of the maximum likelihood estimator μ̂^{0}(y) = y in model (1) as soon as n ≥ 3. Indeed, [17] showed that adaptive linear shrinkage

μ̂^{JS+}(y) = (1 − γ̂)_{+}y,  γ̂ = (n − 2)ϵ^{2}/|y|^{2},

is everywhere better than μ̂^{0}, in the sense that for all μ ∈ ℝ^{n}, r_{ϵ}(μ̂^{JS+}, μ) < r_{ϵ}(μ̂^{0}, μ) ≡ nϵ^{2}. Here a_{+} = max(a, 0). The result was and remains surprising because it can seem counterintuitive that combining data from statistically completely independent problems, represented by each coordinate in (1), leads to better MSE properties.

A simple proof was later given by Stein [25], using his unbiased estimate of risk to show that μ̂^{JS}(y) = (1 − γ̂)y, necessarily worse than μ̂^{JS+}, satisfies

r(μ̂^{JS}, μ) = E_{μ}{n − (n − 2)^{2}|y|^{−2}} < n.  (5)

(where, for simplicity, ϵ = 1 here.) Using in (5) the fact that the distribution of |y|^{2} can be represented as the mixture of central chi-squared distributions χ_{n+2P}^{2} with P distributed as a Poisson variate with mean |μ|^{2}/2, and applying Jensen’s inequality, one obtains our first oracle inequality.

Theorem 3 ([7]). In model (1), suppose n ≥ 3. For all μ ∈ ℝ^{n},

E|μ̂^{JS} − μ|^{2} ≤ 2ϵ^{2} + \frac{(n − 2)ϵ^{2}|μ|^{2}}{(n − 2)ϵ^{2} + |μ|^{2}} .  (6)

In view of (4), this implies

r_{ϵ}(μ̂^{JS+}, μ) ≤ 2ϵ^{2} + R_{ϵ}(μ, LS).  (7)

Thus, the classical James-Stein estimator comes within an additive penalty of 2ϵ^{2} of mimicking the ideal linear shrinkage estimator. This performance is impressive when calibrated against the minimax risk R_{N}(ℝ^{n}, ϵ), in this problem nϵ^{2}.

However it should be noted that this inequality is orthogonally invariant, and makes no use of the particular basis in which the unknown signal μ is represented.

2.1 Illustration: Levelwise shrinkage in the dyadic sequence model.

In the dyadic sequence model of Section 1.1, group coefficients by level j : y_{j} = (y_{jk})_{k=1}^{2^{j}}. Form a levelwise James Stein estimator θ̂^{LJS} by applying James-Stein shrinkage to y_{j}: θ̂_{j}^{LJS} = θ̂^{JS+}(y_{j}), at least for levels j below a cutoff J = log_{2} ϵ^{−2}, above which θ̂_{j}^{LJS} simply estimates zero. [Recall the calibration n = ϵ^{−2} of Example 1(b).] The MSE of the θ̂^{LJS} may then also be represented levelwise:

E∥θ̂^{LJS} − θ∥^{2} = ∑_{j<J(ϵ)} E|θ̂^{JS+}(y_{j}) − θ_{j}|^{2} + ∑_{j≥J(ϵ)} |θ_{j}|^{2}.

The oracle inequality (7) may be applied to each level j in the first sum, while the geometric weights 2^{aj} used to define Θ_{α,2} imply that the second sum is negligible for small ϵ. For the scale of Hilbert spaces S_{2} = {Θ_{α,2}(C) : α > 0, C > 0} :

Theorem 4 ([7]). θ̂^{LJS} is adaptively minimax over S_{2}.

This recovers and extends a notable result of Efroimovich & Pinsker [14], originally formulated in the Fourier basis. In fact, one verifies relatively easily that θ̂ is adaptively minimax among linear estimates (from the ideal linear shrinkage risk) and then appeals to the celebrated theorem of Pinsker [23]), which shows that for the ellipsoids occurring in S_{2}, linear minimax rules are actually asymptotically minimax among all non-linear estimates.

This levelwise application of an oracle inequality is shows how the dyadic sequence model allows a “lifting” of results from a symmetric and “parametric” setting (an exchangeable multivariate normal law at each level) to a non-parametric, infinite-dimensional model. Other examples of this type may be found in [7, 9].

3 Orthogonal Regression and Thresholding

To this point, we have considered only orthogonally invariant estimators. However, a basic principle is that sparsity of representation of a signal in a given basis leads to better estimation, and to exploit such sparsity, non-linear estimators are needed.

Thus, assume the orthonormal basis leading to coefficients (1) is chosen so that {μ_{i}} contains few large coefficients, although of course it is not known in advance which among the co-ordinates are important.

In this orthogonal regression setting, the least squares subset selection estimators have a simple co-ordinatewise representation: the j−th component of μ̂_{J}(y) equals y_{j} if j ∈ J and 0 otherwise. Thus, the least squares estimators have the form of diagonal projections (DP below). The mean squared error of μ̂_{J} is then the sum of terms which measure either variance or bias:

r(μ̂_{J}, μ) = ∑_{j∈J} ϵ^{2} + ∑_{j∉J} μ_{j}^{2}.

The ideal risk for among all such diagonal projection estimators can therefore be found by minimizing termwise:

R_{ϵ}(μ, DP) = inf_{J} r(μ̂_{J}, μ) = ∑_{j} μ_{j}^{2} ∧ ϵ^{2}.

To quantify sparsity, order the squared magnitudes of the components of μ via μ_{(1)}^{2} ≥ μ_{(2)}^{2} ≥ … ≥ μ_{(n)}^{2} and define compression numbers c_{j}^{2} = ∑_{k>j} μ_{(k)}^{2}. The number of large coefficients is measured by N(ϵ) = #{j : |μ_{j}| > ϵ}, and we have

R_{ϵ}(μ, DP) = ϵ^{2}N(ϵ) + c_{N(ϵ)}^{2},

which shows an intimate connection between ideal risk and the compressibility of the signal in this basis.

Various forms of thresholding estimator can be introduced: here we consider soft thresholding:

μ̂_{j}^{ST} (y) = sgn(y_{j})(|y_{j}| − λ)_{+}.

The key points are that the estimator acts co-ordinatewise and that there is a threshold zone [−λ, λ] in which the data is interpreted as noise and “discarded”.

Theorem 5 ([6]). If λ = \sqrt{2 log n}, then for all μ ∈ ℝ^{n},

r_{ϵ}(μ̂^{ST} , μ) ≤ (2 log n + 1)[ϵ^{2} + R_{ϵ}(μ, DP)].  (8)

Since the logarithmic penalty is of small order relative to n, the result shows that sparsity, as measured by ideal risk, implies good estimation. The bound is valid for all sample sizes and all μ. There has been much work on the choice of threshold λ - the choice given here is attractive for its conservatism: since for independent and identically distributed N(0, 1) variates z_{i}, P (max_{1≤j≤n} |z_{j}| > \sqrt{2 log n}) → 0, it follows that P (μ̂^{ST} = 0|μ = 0) → 1. For more on these issues and numerical examples, see [11]. Smaller choices of λ, even depending on the data y, lead to better mean squared error in exchange for less conservatism [7]. Natural extensions of Theorems 5 and 6 to correlated noise exist [19]

Optimality. Absent extra restrictions on μ, the factor 2 log n is optimal:

Theorem 6 ([6]). As n → ∞,

inf_{μ̂} sup_{μ∈ℝ^{n}} \frac{r(μ̂, μ)}{ϵ^{2} + R_{ϵ}(μ, DP)} ≥ (2 log n)(1 + o(1)).

The lower bound arises from the difficulty of distinguishing rare true signal components from the also infrequent extremes of the white Gaussian noise z_{i}. Indeed, suppose ϵ = 1 and that the values μ_{i} are drawn independently from a two point prior distribution with masses of probability 1 − δ_{n} at 0 and δ_{n} at μ̄_{n}. Choosing δ_{n} = log n/n and μ̄_{n} ∼ (2 log δ_{n}^{−1})^{1/2}, it turns out that the posterior distribution of μ_{i}, having observed even a value of y_{i} > μ̄_{n}, is still concentrated on 0 : P (μ = 0|y = μ̄_{n} + z) ≈ 1, for z large and fixed, as n → ∞. Hence, with probability δ_{n}, the estimator is forced to make an error of order μ̄_{n}^{2} ∼ 2 log n.

3.1 Illustration: thresholding in the dyadic sequence model.

Return to the dyadic sequence model, and apply soft thresholding at λ = ϵ\sqrt{2 log ϵ^{−2}} to the first n = ϵ^{−2} coefficients. In other words, θ̂_{I}^{T} (y) = η_{ST} (y_{I}, λ) for all I with j < J(ϵ). Applying the thresholding oracle inequality (8) to the first n co-ordinates,

r_{ϵ}(θ̂^{T} , θ) ≤ c · log ϵ^{−2} · [ϵ^{2} + R_{ϵ}(θ, DP)] + ∑_{j≥J(ϵ)} |θ_{j}|^{2}  (9)

In contrast with the scale S_{2} of Section 2.1, consider now a broader scale of Sobolevtype parameter spaces: S = {Θ_{α,p}(C) : α > 1/p − 1/2, p > 0, C > 0}. For such spaces there is a bound relating ideal to minimax risk. First, the geometric weights in the definition imply ([12]) that for Θ = Θ_{α,p}(C) and on setting r = 2α/(2α+1),

R_{ϵ}(Θ, DP) := sup_{Θ} R_{ϵ}(θ, DP) = sup_{θ∈Θ} ∑ θ_{I}^{2} ∧ ϵ^{2} ≤ c_{α}C^{2(1−r)}ϵ^{2r}.

Second, the minimax risk over Θ is minorized by that over any inscribed hypercube of dimension m and side length ϵ : R_{N}(Θ, ϵ) ≥ c_{0}mϵ^{2}. Optimizing over the dimension m and combining with the previous display, we obtain the basic ideal to minimax risk inequality:

R_{ϵ}(Θ, DP) ≤ c_{α}C^{2(1−r)}ϵ^{2r} ≤ c′_{α}R_{N}(Θ, ϵ).  (10)

In combination with the oracle inequality and negligibility of the tail sum in (9), this yields an adaptive near-minimaxity property for thresholding:

Theorem 7. For all Θ_{α,p}(C) ∈ S,

sup_{Θ_{α,p}(C)} r_{ϵ}(θ̂^{T} , θ) ≤ c_{α,p} · log ϵ^{−2} · R_{N}(Θ_{α,p}(C), ϵ).

The term near-minimaxity refers to the logarithmic term in the upper bound, which is negligible with respect to the algebraic rate ϵ^{2r}. In fact, this logarithmic term can also be removed by a lower, data-dependent choice of threshold [7, 18].

Important here is that in contrast to the linear adaptivity of Theorem 4, this result applies for all p > 0, and in particular for p < 2. These latter spaces contain spatially inhomogeneous functions with localized discontinuities or other singularities. The ability of an estimator to adapt to such functions is in practice more important than the attractive, but limited adaptation of the levelwise James-Stein estimator, and its cousins, the spatially homogeneous kernel methods, even with bandwidth selected from data. This is discussed further in [11].

4 Redundant Dictionaries & Complexity Penalized Model Selection

In seeking a sparse representation for a signal, one may build build rich dictionaries D = {x_{1}, …, x_{p}} in various ways: for example by combining many orthonormal bases (as in libraries of wavelet and cosine packets, [4]), or by considering redundant discretizations of continuously parametrized families, or by allowing products (interactions) of many simple elements, such as B-splines with knots at individual data locations (e.g. [16]). In all these cases, the dictionary size p greatly exceeds that data size n, and estimation methods will have to allow for the effects of searching over such a vast domain (in principle, 2^{p} models).

Recalling Examples 1(c) and 2(b), the data may be represented in the form y = Xβ + ϵz, where we now assume that span(X) = ℝ^{n}. Thus, the models of interest correspond to subsets J ⊂ {1, …, p}, M_{J} = span {x_{j} : j ∈ J}, and μ̂_{J} = P_{J}y, orthogonal projection on M_{J}. The risk of individual projection estimators is given by (3), so the ideal risk of subset selection from dictionary D becomes

R_{ϵ}(μ, SS(D)) = min_{J} r_{ϵ}(μ̂_{J}, μ) = min_{J} |μ − P_{J}μ|^{2} + ϵ^{2}rank(P_{J}).

To obtain an estimator that mimicks ideal risk, we use the penalized least squares principle. This balances the fit of the estimate, which in the absence of any penalty could be made arbitrarily close to the data, against some measure of complexity of the estimate:

μ̂_{P} = argminμ̃|y − μ̃|^{2} + ϵ^{2}P (μ̃).

In the orthonormal basis setting, μ̂_{P} can be evaluated explicitly when the penalty P has an additive form: for example, P (μ) = c ∑ μ_{i}^{2} implies linear shrinkage, P (μ) = 2λ ∑ |μ_{i}| implies soft thresholding, and P (μ) = λ^{2} ∑ I{μ_{i} ≠ 0} implies μ̂_{P,i}(y) = y_{i}I{|y_{i}| > λ}, or hard thresholding. For the redundant linear model y = Xβ + ϵz, we generalize the third case by setting

P (μ) = λ^{2}N(μ), N(μ) = min{|J| : μ = ∑_{j∈J} β_{j}x_{j}}.

The resulting penalized least squares estimator may expressed in terms of the residual sums of squares RSS_{J} = |y − μ̂_{J}|^{2} of the possible models:

min_{μ̃} |y − μ̃|^{2} + λ^{2}ϵ^{2}N(μ̃) = min_{J} RSS_{J} + λ^{2}ϵ^{2}rank(P_{J}).

Hence we call this the Complexity Penalized Residual Sum of Squares (CPRSS) estimate. Certain choices of the factor λ lead to well known estimators: λ^{2} = 2 (AIC), log p (BIC), 2 log n (RIC) (For details and references see [8]).

Theorem 8 ( [8]). Let ζ > 1, β > 0 and λ = λ_{p} = ζ[1 + \sqrt{2(1 + β) log(p + 1)}]. Then for all n, p ≥ n,and μ ∈ ℝ^{n},

r_{ϵ}(μ̂_{CP RSS}, μ) ≤ L_{p}[(2 + γ_{p})ϵ^{2} + R_{ϵ}(μ, SS(D))],  (11)

where L_{p} = (1 − ζ^{−1})^{−1}λ_{p}^{2}, and γ_{p} = γ(p, β) → 0 as p → ∞.

The penalty factor λ_{p}^{2} is slightly larger than 2 log p, where p is the cardinality of the dictionary. We emphasize that the result holds for all μ, n and p ≥ n, and in particular the inequality depends only on p, not n! Building on the remarkable [1], Birgé & Massart are conducting a thorough study of penalties P (μ) for which such oracle inequalities and improvements hold. While the constant L_{p} in (11) is certainly not optimal, there is a lower bound similar to Theorem 6:

Theorem 9 ( [8]). For each fixed r ∈ ℕ, there exists a sequence of dictionaries D_{n} with p(n) = |D_{n}| ≍ n^{r} such that as n → ∞,

inf_{μ̂} sup_{μ∈ℝ^{n}} \frac{E|μ̂ − μ|^{2}}{ϵ^{2} + R_{ϵ}(μ, SS(D))} ≥ [2 log p(n)](1 + o(1)].

Role of Concentration Inequalities. The stochastic part of the proof of Theorem 8 depends on an early example (due to Cirelson-Ibragimov-Sudakov [2, 21]) of what are now in probability called concentration (or deviation) inequalities. Suppose f : ℝ^{n} → ℝ is Lipschitz with ∥f∥_{Lip} = L. If Z ∼ N_{n}(0, I), then

P {f(Z) ≥ Ef(Z) + t} ≤ exp{−t^{2}/2L^{2}}.

The key points are the Gaussian tail behaviour of f(Z) and the fact that it does not depend on dimension n - hence the dimension-free aspect of Theorem 8. This inequality can then be applied to all projections onto model subsets of cardinality |J| = ℓ, and then summed over ℓ. Thus, since f(z) = ∥P_{J}z∥ has ∥f∥_{Lip} = 1, and since Ef(Z) ≤ \sqrt{ℓ}, we have, on setting t = \sqrt{2ℓ(1 + β) log p},

P \{ sup_{|J|=ℓ} ∥P_{J}z∥ ≥ \sqrt{ℓ} + t\} ≤ \binom{p}{ℓ}p^{−ℓ(1+β)} ≤ \frac{1}{p^{ℓβ}ℓ!}.

Documenta Mathematica · Extra Volume ICM 1998 · III · 267–278

4.1 Illustration: Minimaxity for non-standard function classes

The penalized least squares formalism can be applied in situations where no unconditional basis exists. To give a simple example, consider again the model dY_{t} = fdt + ϵdW_{t}, where now t ∈ [0, 1]^{2}, and the horizon model for edges in images, studied earlier by, for example, Korostelev and Tsybakov [20]. It is supposed that f takes only the values 0 and 1, and further that the boundary is such that f(t_{1}, t_{2}) = I{t_{2} ≤ θ(t_{1})}. The boundary, or horizon, is supposed to be Hölder continuous: more specifically, we say that f ∈ Hölder_{s}(B) if ∥θ∥_{∞} + ∥θ^{(r)}∥_{β} ≤ B, where r ∈ ℕ, β = s − r ∈ (0, 1] and ∥g∥_{β} = sup |g(t) − g(t′)|/|t − t′|^{β}…. ]

Dictionaries and minimax risk. While D is often conceptually infinite, in practice one must work with a family of finite subdictionaries D_{ϵ} with cardinality m(ϵ) being at most a polynomial function of ϵ^{−2} : m(ϵ) ≤ β_{1}ϵ^{−2β_{2}}. [8] defines a notion of universal dictionary for a scale S = {F_{ν}(C)}offunctionclasses, which has as consequence the same type of ideal to minimax risk inequality as used in the orthobasis case (compare (10)): for all F_{ν}(C) ⊂ S and ϵ < ϵ(ν, C), there exists r = r(ν) such that

R_{ϵ}(F_{ν}(C), D_{ϵ}) ≤ K_{ν}C^{2(1−r)}ϵ^{2r} ≤ K′_{ν}R_{N}(F_{ν}(C), ϵ).

This may then be combined with the oracle inequality of Theorem 8 to obtain adaptive near-minimaxity.

Thus, in the horizon example, we start with a continuum trapezoid dictionary, parametrized by γ = (a, b, c, d), representing a function taking value 1 on the trapezoid in [0, 1]^{2} with abscissae a < b and corresponding ordinates c, d. Thus D_{T rap} = {T_{γ} : γ ∈ [0, 1]^{4}, b ≥ a}. To obtain finite subdictionaries, discretize the unit interval into I_{N} = {i/N : 0 ≤ i ≤ N} and set D_{N} = {T_{γ} : γ ∈ I_{N}^{2} × I_{N^{2}}^{2}}. Choose N(ϵ) = ϵ^{−2}, and set D_{ϵ} = D_{N(ϵ)}. It can be verified [8] that D_{T rap} is universal for S = {Hölder_{s}(B) : 0 < s ≤ 2, 0 < B}, with ν = s/2, C = B^{1/2}.

Corollary 10 ([8]). On Hölder_{s}(B), for 0 < s ≤ 2, and setting r = s/(s + 1),

r_{ϵ}(f̂_{CP RSS}, f) ≤ c_{0} · log_{2} ϵ^{−2} · B^{1−r}ϵ^{2r}.

A key remark is that this adaptively (near minimax) rate of convergence is better than the rate attainable using a two dimensional tensor product wavelet basis when s > 1.

Nevertheless, a serious practical defect of Theorem 8 is the combinatorial search implicit in the definition of μ̂_{CP RSS}. The development of fast algorithms suitable for specific cases is an active direction of current research [5, 13].

Acknowledgements The results described here were obtained in collaboration with David Donoho, and some also in joint work with Gérard Kerkyacharian, Dominique Picard, and Bernard Silverman. Work on this manuscript and talk was supported in part by NSF, NIH and the Guggenheim Foundation.

References

[1] A. Barron, L. Birgé, and P. Massart. Risk bounds for model selection via penalization. Probability Theory and Its Applications, 00:in press, 1999.

[2] B.S. Cirel’son, I.A. Ibragimov, and V.N. Sudakov. Norm of gaussian sample function. In Proceedings of the 3rd Japan-U.S.S.R. Symposium on Probability Theory, Lecture Notes in Mathematics, 550, pages 20–41, 1976.

[3] A. Cohen, I. Daubechies, and P. Vial. Wavelets and fast wavelet transform on an interval. Applied Computational and Harmonic Analysis, 1:54–81, 1993.

[4] R.R. Coifman, Y. Meyer, and M.V. Wickerhauser. Wavelet analysis and signal processing. In B. Ruskai et. al., editor, Wavelets and their Applications, pages 153–178. Jones and Bartlett, 1992.

[5] R.R. Coifman and M.V. Wickerhauser. Entropy-based algorithms for bestbasis selection. I.E.E.E. Transactions on Information Theory, 38:713–718, 1992.

[6] D. L. Donoho and I. M. Johnstone. Ideal spatial adaptation via wavelet shrinkage. Biometrika, 81:425–455, 1994.

[7] D. L. Donoho and I. M. Johnstone. Adapting to unknown smoothness via wavelet shrinkage. J. Amer. Statist. Assoc., 90:1200–1224, 1995.

[8] D. L. Donoho and I. M. Johnstone. Empirical atomic decomposition. Technical report, Stanford University, 1995.

[9] D. L. Donoho and I. M. Johnstone. Minimax estimation via wavelet shrinkage. Annals of Statistics, 1998. To appear.

[10] D. L. Donoho and I. M. Johnstone. Asymptotic minimaxity of wavelet estimators with sampled data. Statistica Sinica, 00:in press, 1999.

[11] D. L. Donoho, I. M. Johnstone, G. Kerkyacharian, and D. Picard. Wavelet shrinkage: Asymptopia? Journal of the Royal Statistical Society, Series B, 57:301–369, 1995. With Discussion.

[12] D. L. Donoho, I. M. Johnstone, G. Kerkyacharian, and D. Picard. Universal near minimaxity of wavelet shrinkage. In Pollard D., Torgersen E., and Yang G.L., editors, Festschrift for L. Le Cam, pages 183–218. Springer Verlag, 1997.

[13] D.L. Donoho. Wedgelets: Nearly minimax estimation of edges. Technical report, Dept. of Statistics, Stanford University, 1997.

[14] S.Yu. Efroimovich and M.S. Pinsker. A learning algorithm for nonparametric filtering. Automat. i Telemeh., 11:58–65, 1984. (in Russian).

[15] M. Frazier, B. Jawerth, and G. Weiss. Littlewood-Paley Theory and the study of function spaces. NSF-CBMS Regional Conf. Ser in Mathematics, 79. American Mathematical Society, Providence, RI, 1991.

[16] J.H. Friedman. Multivariate adaptive regression splines. Annals of Statistics, 19:1–67, 1991. (with discussion).

[17] W. James and C. Stein. Estimation with quadratic loss. In Proceedings of Fourth Berkeley Symposium on Mathematical Statistics and Probability Theory, pages 361–380. University of California Press, 1961.

[18] I. M. Johnstone. Wavelet shrinkage for correlated data and inverse problems: adaptivity results. Statistica Sinica, 00:in press, 1999.

[19] I. M. Johnstone and B. W. Silverman. Wavelet threshold estimators for data with correlated noise. Journal of the Royal Statistical Society, Series B., 59:319–351, 1997.

[20] A.P. Korostelev and A.B. Tsybakov. Minimax Theory of Image Reconstruction: Lecture Notes in Mathematics. Springer Verlag: New York, 1993.

[21] M. Ledoux. Isoperimetry and gaussian analysis. In P. Bernard, editor, Lectures on Probability Theory and Statistics, Ecole d’Eté de Probabilities de Saint Flour, 1994. Springer Verlag, 1996.

[22] P.G. Lemarié and Y. Meyer. Ondelettes et bases Hilbertiennes. Revista Matematica Iberoamericana, 2:1–18, 1986.

[23] M.S. Pinsker. Optimal filtering of square integrable signals in gaussian white noise. Problems of Information Transmission, 16:120–133, 1980. originally in Russian in Problemy Peredatsii Informatsii 16 52-68.

[24] C. Stein. Efficient nonparametric estimation and testing. In Proc. Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 1, pages 187–195. University of California Press, Berkeley, CA., 1956.

[25] Charles Stein. Estimation of the mean of a multivariate normal distribution. Annals of Statistics, 9:1135–1151, 1981.

Iain M. Johnstone

Department of Statistics Stanford University Stanford CA 94305 U.S.A.