Orthogonal-Series and Wavelet Density Estimation

Intuition Beginner

You have a sample of numbers and you want to draw the curve they came from — the density that says where values are common and where they are rare. One way you have already met is the histogram, or its smooth cousin the kernel estimate: slide a little bump over each data point and add them up. This unit takes a different route. Instead of building the curve out of bumps, build it out of waves.

Think of any musical chord as a sum of pure tones: a low hum for the overall shape, then higher tones adding finer wiggles. A density can be written the same way. Pick a fixed library of standard wave shapes — say cosine curves of increasing frequency. The true density is some particular recipe: so much of the slow cosine, so much of the next, and so on. Each amount in the recipe is just an average you can estimate from the data.

Keep the first few terms, throw away the rest, and you have a smooth estimate. How many terms you keep is the one knob: too few and you miss real features, too many and you start tracing the noise.

The catch is that fixed waves spread across the whole line. If your density is smooth in most places but has one sharp spike, the cosine recipe has to use many high-frequency terms everywhere just to build that one spike, and those terms add ripples far away from it. Kernels have the same trouble: one bandwidth for the whole curve cannot be both wide where the curve is flat and narrow where it is sharp.

Wavelets fix this. A wavelet is a little wave that is localized — it lives in one small region and is silent elsewhere — and it comes in a whole family of sizes and positions. You spend coarse wavelets on the broad shape and fine wavelets only near the spike, where detail is actually needed. Better still, you can look at all the estimated wavelet amounts and simply delete the small ones, which are almost surely just noise. Keeping only the few large amounts gives a clean estimate that automatically uses fine detail where the function wiggles and stays smooth where it does not.

Visual Beginner

The table contrasts how the two families spend their "budget" of terms, and what the deletion rule does.

Approach	What the building blocks look like	Where detail goes	What you tune
Cosine series	full-width waves, low to high frequency	spread evenly over the whole curve	how many terms to keep (the cutoff)
Kernel / bump	one bump shape, one width	the same everywhere	the single bandwidth
Wavelets, keep all coarse	small localized waves at many scales and places	a fixed amount at each scale	the finest scale
Wavelets, delete small ones	same localized waves	concentrated only where the curve wiggles	the deletion level (threshold)

The single picture to carry forward: a series estimate decides how many waves to keep; a thresholded wavelet estimate decides which waves to keep, deleting the small noisy ones. Choosing which, not just how many, is what lets wavelets put detail exactly where the curve needs it.

Worked example Beginner

We estimate the first few cosine amounts of a density from a tiny sample, by direct averaging. Use the cosine library on the interval from $0$ to $1$. Its members are the flat function ${\varphi }_0(x)=1$ and the cosines ${\varphi }_j(x)=\sqrt{2}\cos (\pi jx)$ for $j=1,2,\dots$. The amount of ${\varphi }_j$ in the true density $f$ is the average value of ${\varphi }_j$ when $x$ is drawn from $f$, written ${\beta }_j$. So we can estimate ${\beta }_j$ by plugging each data point into ${\varphi }_j$ and averaging.

Take four data points: $0.1,0.2,0.7,0.9$. We estimate ${\beta }_1$, the amount of the first cosine ${\varphi }_1(x)=\sqrt{2}\cos (\pi x)$.

Step 1. Evaluate ${\varphi }_1$ at each point. With $\sqrt{2}\approx 1.414$: at $0.1$, $\cos (0.1\pi )=\cos (18^{\circ })\approx 0.951$, so ${\varphi }_1\approx 1.345$. At $0.2$, $\cos (36^{\circ })\approx 0.809$, so ${\varphi }_1\approx 1.144$. At $0.7$, $\cos (126^{\circ })\approx -0.588$, so ${\varphi }_1\approx -0.831$. At $0.9$, $\cos (162^{\circ })\approx -0.951$, so ${\varphi }_1\approx -1.345$.

Step 2. Average the four values: $(1.345+1.144-0.831-1.345)/4=0.313/4\approx 0.078$.

Step 3. Read the meaning. The estimate ${\hat{\beta }}_1\approx 0.08$ is small and positive, saying the data lean only slightly toward the low end where $\cos (\pi x)$ is positive. The estimate of ${\beta }_0$ is always exactly $1$, because ${\varphi }_0=1$ averages to $1$ for any sample.

Step 4. Build a two-term estimate: $\hat{f}(x)=1+0.078\cdot \sqrt{2}\cos (\pi x)$. This is a gentle tilt away from flat.

What this tells us: every amount in the recipe is just a sample average, so building a series density estimate is the same averaging idea behind the empirical distribution, applied to wave shapes instead of step thresholds. Keep more terms to add finer features; stop early to stay smooth.

Check your understanding Beginner

Formal definition Intermediate+

Let $X_1,\dots ,X_n$ be i.i.d. with an unknown density $f\in L^2([0,1])$, and fix an orthonormal basis $\{{\varphi }_0,{\varphi }_1,{\varphi }_2,\dots \}$ of $L^2([0,1])$ with respect to the Lebesgue inner product $⟨g,h⟩={\int }_0^{1}gh$. The development follows Wasserman ^{[Wasserman, L. — All of Nonparametric Statistics]} Ch. 8-9 and Tsybakov ^{[Tsybakov, A. B. — Introduction to Nonparametric Estimation]} Ch. 1-3. The plug-in averaging that defines the coefficient estimates is the same device that produces the empirical distribution and its functionals 45.05.02.

Definition (Fourier coefficients in a basis). Since $\{{\varphi }_j\}$ is a complete orthonormal system, $f={\sum }_{j\ge 0}{\beta }_j{\varphi }_j$ in $L^2$, with coefficients

$${\beta }_j=⟨f,{\varphi }_j⟩={\int }_0^1{\varphi }_j(x)f(x)dx=\mathbb{E}{\varphi }_j(X_1).$$

The representation of $\Vert f{\Vert }_2^2={\sum }_j{\beta }_j^2$ is Parseval's identity for the basis $\{{\varphi }_j\}$; the $L^2$-theory of such expansions, with Bessel's inequality and Parseval's identity, is developed for the trigonometric system in 02.10.01, and the abstract Hilbert-space version is 02.11.08. The standard choices are the cosine basis ${\varphi }_0=1$, ${\varphi }_j(x)=\sqrt{2}\cos (\pi jx)$, and the Legendre or Fourier bases.

Definition (orthogonal-series / projection estimator). Because each ${\beta }_j$ is an expectation, estimate it by the sample average ${\hat{\beta }}_j=\frac{1}{n}{\sum }_{k=1}^n{\varphi }_j(X_k)$, an unbiased estimator with $\mathbb{E}{\hat{\beta }}_j={\beta }_j$. Truncating at a resolution (cutoff) $J$ gives the orthogonal-series (or projection) density estimator

$${\hat{f}}_J(x)=\sum _{j=0}^J{\hat{\beta }}_j{\varphi }_j(x).$$

The cutoff $J$ is the smoothing parameter, the series analogue of a kernel bandwidth: small $J$ oversmooths, large $J$ undersmooths.

Definition (Sobolev ellipsoid). A density has Sobolev smoothness of order $p$ with radius $L$ if its coefficients lie in the ellipsoid

$$\Theta (p,L)=\{\theta =({\theta }_j{)}_{j\ge 0}:\sum _{j\ge 1}{a}_j^2{\theta }_j^2\le L^2\},a_j\asymp j^p,$$

which for the trigonometric basis encodes $p$ square-integrable weak derivatives. The decay rate $a_j\asymp j^p$ is the smoothness; faster decay of ${\beta }_j$ means a smoother $f$.

Definition (multiresolution analysis). A multiresolution analysis (MRA) of $L^2(\mathbb{R})$ is a nested sequence of closed subspaces $\cdots \subset V_{-1}\subset V_0\subset V_1\subset \cdots$ with ${\bigcap }_j{V}_j=\{0\}$, $\overline{{\bigcup }_j{V}_j}=L^2(\mathbb{R})$, the scaling relation $g\in V_j⟺g(2\cdot )\in V_{j+1}$, and a father wavelet (scaling function) $\phi$ whose integer shifts $\{\phi (\cdot -k){\}}_{k\in \mathbb{Z}}$ form an orthonormal basis of $V_0$. The orthogonal complement $W_j=V_{j+1}\ominus V_j$ is the detail space, spanned by the shifts of a mother wavelet $\psi$.

Definition (wavelet basis and DWT). From $(\phi ,\psi )$ one forms the dilated-translated families

$${\phi }_{j,k}(x)=2^{j/2}\phi (2^{j}x-k),{\psi }_{j,k}(x)=2^{j/2}\psi (2^{j}x-k),$$

and $\{{\phi }_{j_0,k}{\}}_k\cup \{{\psi }_{j,k}:j\ge j_0,k\}$ is an orthonormal basis of $L^2(\mathbb{R})$. The expansion $f={\sum }_k{\alpha }_{j_0,k}{\phi }_{j_0,k}+{\sum }_{j\ge j_0}{\sum }_k{\beta }_{j,k}{\psi }_{j,k}$ has approximation coefficients ${\alpha }_{j_0,k}=⟨f,{\phi }_{j_0,k}⟩$ at the coarse scale and detail coefficients ${\beta }_{j,k}=⟨f,{\psi }_{j,k}⟩$ at finer scales. The map from samples of $f$ to the finite array $(\alpha ,\beta )$ via the cascade filter recursion is the discrete wavelet transform (DWT); the prototype is the Haar MRA, $\phi =1_{[0,1)}$ and $\psi =1_{[0,1/2)}-1_{[1/2,1)}$.

Definition (thresholding). Given empirical detail coefficients ${\hat{\beta }}_{j,k}$ and a threshold $\lambda >0$, the hard-threshold and soft-threshold estimators are

$${\hat{\beta }}_{j,k}^H={\hat{\beta }}_{j,k}1\{|{\hat{\beta }}_{j,k}|>\lambda \},{\hat{\beta }}_{j,k}^S=\mathrm{sign}({\hat{\beta }}_{j,k})(|{\hat{\beta }}_{j,k}|-\lambda {)}_{+}.$$

Hard thresholding keeps a coefficient unchanged or deletes it; soft thresholding additionally shrinks the survivors toward zero by $\lambda$. The soft-threshold operator is exactly the proximal map of $\lambda |\cdot |$, the same shrinkage that produces sparsity in the lasso 45.06.06.

Counterexamples to common slips

• The estimated density need not be a density. The truncated series ${\hat{f}}_J$ can be negative and need not integrate to one (though $\int {\hat{f}}_J={\hat{\beta }}_0$ is forced to a constant by orthonormality with ${\varphi }_0=1$). Post-processing — truncating negatives and renormalizing — restores a bona fide density without changing the rate.
• The cutoff $J$ is not a count of "good" terms. Linear truncation keeps all coefficients up to $J$, including small ones full of noise, and discards all beyond $J$, including occasional large ones. This is precisely the rigidity that nonlinear thresholding removes by deciding coefficient-by-coefficient.
• A high-resolution wavelet basis is not a higher-order kernel. A linear wavelet estimator at fixed finest scale is a kernel-type (linear) estimator and shares the kernel's inability to adapt to inhomogeneous smoothness 26.08.01; the adaptivity comes only from the nonlinear thresholding step, not from the basis alone.
• Soft and hard thresholding differ in bias. Hard thresholding is unbiased on survivors but discontinuous in the data; soft thresholding is continuous but biases every survivor by $\lambda$. The continuity is why soft thresholding has the clean oracle inequality below.

Key theorem with proof Intermediate+

The signature result for the linear (series) estimator is the bias-variance decomposition of its integrated risk and the optimal cutoff and rate it forces over a Sobolev class.

Theorem (risk of the projection estimator and the optimal rate). Let ${\hat{f}}_J={\sum }_{j=0}^J{\hat{\beta }}_j{\varphi }_j$ with ${\hat{\beta }}_j=n^{-1}{\sum }_k{\varphi }_j(X_k)$, and assume ${\sup }_x{\sum }_{j=0}^J{\varphi }_j(x{)}^2\le C_{\varphi }(J+1)$ (a uniform-boundedness condition met by the cosine and trigonometric bases). Then the integrated mean-squared error decomposes as

$$\mathbb{E}\Vert {\hat{f}}_J-f{\Vert }_2^2=\underset{{\text{bias}}^2}{\underbrace{\sum _{j>J}{\beta }_j^2}}+\underset{\text{variance}}{\underbrace{\sum _{j=0}^J\mathrm{Var}({\hat{\beta }}_j)}},\sum _{j=0}^J\mathrm{Var}({\hat{\beta }}_j)\le \frac{C_{\varphi }(J+1)}{n}.$$

If $f$ lies in the Sobolev ellipsoid $\Theta (p,L)$ with $a_j\asymp j^p$, then choosing $J_n\asymp n^{1/(2p+1)}$ gives

$$\mathbb{E}\Vert {\hat{f}}_{J_n}-f{\Vert }_2^2\lesssim n^{-2p/(2p+1)},$$

which is the minimax-optimal rate over $\Theta (p,L)$.

Proof. Orthonormality gives the Pythagorean split. Write ${\hat{f}}_J-f={\sum }_{j\le J}({\hat{\beta }}_j-{\beta }_j){\varphi }_j-{\sum }_{j>J}{\beta }_j{\varphi }_j$; these two pieces are orthogonal in $L^2$ (disjoint index sets), so $\Vert {\hat{f}}_J-f{\Vert }_2^2={\sum }_{j\le J}({\hat{\beta }}_j-{\beta }_j{)}^2+{\sum }_{j>J}{\beta }_j^2$. Taking expectations and using $\mathbb{E}({\hat{\beta }}_j-{\beta }_j{)}^2=\mathrm{Var}({\hat{\beta }}_j)$ (unbiasedness) yields the displayed decomposition. For the variance, ${\hat{\beta }}_j$ is an average of i.i.d. terms, so $\mathrm{Var}({\hat{\beta }}_j)=n^{-1}\mathrm{Var}({\varphi }_j(X_1))\le n^{-1}\mathbb{E}{\varphi }_j(X_1{)}^2=n^{-1}{\int }_0^1{\varphi }_j^{2}f$. Summing over $j\le J$ and using the uniform bound ${\sum }_{j\le J}{\varphi }_j(x{)}^2\le C_{\varphi }(J+1)$,

$$\sum _{j=0}^J\mathrm{Var}({\hat{\beta }}_j)\le \frac{1}{n}{\int }_0^1\left(\sum _{j=0}^J{\varphi }_j^2\right)f\le \frac{C_{\varphi }(J+1)}{n}{\int }_0^{1}f=\frac{C_{\varphi }(J+1)}{n}.$$

For the bias, membership in $\Theta (p,L)$ means ${\sum }_j{a}_j^2{\beta }_j^2\le L^2$ with $a_j\asymp j^p$, so the tail obeys ${\sum }_{j>J}{\beta }_j^2\le a_J^{-2}{\sum }_{j>J}{a}_j^2{\beta }_j^2\le c{L}^2{J}^{-2p}$. The risk is therefore bounded by $c{L}^2{J}^{-2p}+C_{\varphi }(J+1)/n$. Minimizing $J^{-2p}+J/n$ over $J$ balances the two terms at $J^{-2p}\asymp J/n$, i.e. $J^{2p+1}\asymp n$, giving $J_n\asymp n^{1/(2p+1)}$ and common value of order $n^{-2p/(2p+1)}$. $\square$

Bridge. The decomposition's engine is orthonormality, which turns the global $L^2$ error into a coordinatewise sum and so reduces nonparametric estimation to an infinite collection of one-dimensional mean estimation problems in the Fourier coefficients. This builds toward the Master-tier Gaussian sequence model, where ${\hat{\beta }}_j\approx {\beta }_j+(\sigma /\sqrt{n})z_j$ becomes the exact model $y_j={\theta }_j+\epsilon z_j$, and the projection estimator is then literally "keep the first $J$ coordinates," whose risk is the same bias-tail-plus-variance sum; the optimal-rate calculation appears again there as the Pinsker problem of choosing a coordinatewise filter. The foundational reason the curse of dimensionality is escaped is that smoothness forces the coefficients to decay, so only $\asymp n^{1/(2p+1)}$ of them carry signal above the noise floor $1/n$, and this is exactly the effective-dimension bookkeeping that the cutoff $J$ performs. The linear estimator here generalises to the nonlinear thresholded estimator, which keeps coordinates by size rather than by index: where the projection estimator keeps a fixed prefix and is optimal over a single smoothness class, thresholding keeps a data-chosen set and is near-optimal over a whole scale of classes at once. Putting these together, the cutoff and the threshold are the two faces of one question — which coordinates carry signal — and the central insight is that for inhomogeneously smooth functions the answer is not a prefix, which is precisely where kernels 26.08.01 and linear series estimators are beaten by wavelet thresholding.

Exercises Intermediate+

Lean formalization Intermediate+

Mathlib hosts the Hilbert-space substrate (HilbertBasis, the L2 inner-product API, Parseval) but none of the statistical estimation objects, so no module is wired in (lean_status: none). The intended statement of the projection-estimator risk decomposition and of the soft-threshold map, once a sampling-and-empirical-coefficient API exists, would read roughly as follows.

-- Intended shape; not part of the current Babel Bible Lean build.
-- Requires: an orthonormal `HilbertBasis ℕ ℝ (L2 [0,1])` and an empirical-coefficient map.
variable {n : ℕ} (φ : ℕ → (L2 (Set.Icc (0:ℝ) 1)))
  (hφ : Orthonormal ℝ φ) (f : L2 (Set.Icc (0:ℝ) 1))

-- Soft-thresholding: the proximal map of λ‖·‖₁, shared with the lasso unit.
noncomputable def softThreshold (lam y : ℝ) : ℝ :=
  Real.sign y * max (|y| - lam) 0

-- Bias-variance decomposition of the projection estimator's integrated risk.
-- E ‖f̂_J - f‖² = Σ_{j>J} βⱼ² + Σ_{j≤J} Var(β̂ⱼ).
theorem projection_risk_decomposition (J : ℕ) (β : ℕ → ℝ)
    (hβ : ∀ j, β j = inner (f) (φ j)) :
    True := by
  trivial  -- placeholder: the empirical-coefficient/Var API is absent in Mathlib

the Mathlib gap analysis records what is missing: the orthogonal-series and wavelet density estimators, the bias-variance risk decomposition and the $n^{-2p/(2p+1)}$ Sobolev rate, the multiresolution-analysis predicate and the Haar (and smooth) wavelet HilbertBasis, the soft/hard-threshold maps, and the Donoho-Johnstone oracle inequality with the universal threshold $\sigma \sqrt{2\log n}$.

Advanced results Master

The series-estimator risk and the thresholding theory unify in the Gaussian sequence model, the canonical reduction of nonparametric estimation. Under the orthonormal basis the empirical coefficients satisfy ${\hat{\beta }}_j\approx {\beta }_j+(\sigma /\sqrt{n})z_j$ with approximately independent $z_j\sim N(0,1)$, which is the model $y_j={\theta }_j+\epsilon z_j$, $\epsilon =\sigma /\sqrt{n}$. In this model the projection estimator is diagonal truncation and the threshold estimators are diagonal nonlinearities, and the entire rate theory is the analysis of coordinatewise risk.

Theorem 1 (Pinsker: sharp asymptotic minimaxity of linear filters). Over the Sobolev ellipsoid $\Theta (p,L)$ in the sequence model with noise $\epsilon$, the minimax $L^2$ risk is attained, to leading order, by a linear filter ${\hat{\theta }}_j=c_j{y}_j$ with Pinsker weights $c_j=(1-a_j/\mu {)}_{+}$, and equals $P(p,L){\epsilon }^{4p/(2p+1)}(1+o(1))$ for an explicit Pinsker constant $P(p,L)$ ^{[Tsybakov, A. B. — Introduction to Nonparametric Estimation]}. With $\epsilon =\sigma /\sqrt{n}$ this is the rate $n^{-2p/(2p+1)}$ with the exact optimal constant, and it shows that linear estimators are not merely rate-optimal but constant-optimal over a single Sobolev class — the projection estimator's only deficiency is its fixed grid, which Pinsker's tapering refines.

Theorem 2 (Donoho-Johnstone oracle inequality and the universal threshold). In the sequence model $y_i={\theta }_i+\sigma z_i$, $i=1,\dots ,N$, the soft-threshold estimator with the universal threshold ${\lambda }_N=\sigma \sqrt{2\log N}$ satisfies, for all $\theta$,

$$\mathbb{E}\Vert {\hat{\theta }}^S-\theta {\Vert }_2^2\le (2\log N+1)({\sigma }^2+\sum _{i=1}^N\min ({\theta }_i^2,{\sigma }^2)),$$

so it pays at most a $2\log N+1$ factor over the oracle diagonal-projection risk ^{[Donoho, D. L. & Johnstone, I. M. — Ideal spatial adaptation by wavelet shrinkage]}. The threshold is universal because ${\max }_i|\sigma z_i|\le \sigma \sqrt{2\log N}$ with probability tending to one: at this level pure noise is, with high probability, killed entirely, so the estimator is noise-free under the null $\theta =0$. The same paper proves hard thresholding obeys a comparable bound; soft thresholding is preferred for its continuity and the resulting clean risk algebra.

Theorem 3 (Besov-ball adaptivity and near-minimaxity). Translating the sequence-model oracle inequality back through the wavelet isometry, the thresholded wavelet estimator with $\lambda =\sigma \sqrt{2\log n}$ is simultaneously near-minimax over a full scale of Besov balls $B_{p,q}^s(L)$ — including spaces of inhomogeneous and spatially varying smoothness — attaining the minimax rate up to a $\log n$ factor without knowing the smoothness $s$, $p$, $q$ ^{[Donoho, D. L. & Johnstone, I. M. — Ideal spatial adaptation by wavelet shrinkage]}. The Besov scale is the natural home for wavelet thresholding because the Besov norm is equivalent to a weighted ${\ell }_p$ norm of the wavelet coefficients; functions with a jump or a cusp lie in a Besov ball for $p<2$ but in no Sobolev ball of comparable index, and over such balls no linear estimator — kernel, spline, or fixed-cutoff series — can attain the minimax rate. The nonlinearity of thresholding is essential, not incidental.

Theorem 4 (wavelet density estimation). For i.i.d. $X_1,\dots ,X_n\sim f$, the linear wavelet density estimator $\hat{f}={\sum }_k{\hat{\alpha }}_{j_0,k}{\phi }_{j_0,k}+{\sum }_{j_0\le j<J}{\sum }_k{\hat{\beta }}_{j,k}{\psi }_{j,k}$ with ${\hat{\alpha }}_{j,k}=n^{-1}{\sum }_i{\phi }_{j,k}(X_i)$, ${\hat{\beta }}_{j,k}=n^{-1}{\sum }_i{\psi }_{j,k}(X_i)$ attains the $L^2$ rate $n^{-2p/(2p+1)}$ over $H^p$ when $2^J\asymp n^{1/(2p+1)}$; its thresholded version ${\hat{\beta }}_{j,k}^S$ with ${\lambda }_j\asymp \sqrt{j/n}$ attains the Besov-class rates adaptively ^{[Wasserman, L. — All of Nonparametric Statistics]}. The empirical coefficients are again pure plug-in averages of basis functions — the density-estimation specialization of the empirical-measure functionals of 45.05.02 — so the variance bookkeeping is identical to the projection theorem, scale by scale.

Synthesis. The foundational reason orthogonal-series, kernel, and wavelet estimation are one subject is that each reduces, under an orthonormal basis, to estimating the coordinates of $f$ in the Gaussian sequence model, and the only question is which coordinates carry signal above the noise floor $\sigma /\sqrt{n}$. This is exactly the bias-variance trade-off the cutoff $J$ resolves for linear estimators and the threshold $\lambda$ resolves for nonlinear ones, and the central insight is that the optimal answer is a prefix of coordinates when smoothness is homogeneous (where Pinsker's linear filter is even constant-optimal) but a data-selected sparse set when smoothness is inhomogeneous (where only thresholding attains the rate).

The thresholding step is dual to the lasso 45.06.06: soft thresholding is the proximal map of $\lambda \Vert \cdot {\Vert }_1$ in an orthonormal design, so the oracle inequality $R(\hat{\theta })\le (2\log N+1)({\sigma }^2+{\sum }_i\min ({\theta }_i^2,{\sigma }^2))$ is the sequence-model twin of the lasso oracle inequality, and the universal threshold $\sigma \sqrt{2\log N}$ is the sequence-model twin of the lasso penalty $\sigma \sqrt{2\log p/n}$ — both arise from the Gaussian maximal inequality controlling the largest noise coordinate. Putting these together, wavelets generalise the global waves of Fourier analysis 02.10.01 by localizing them at every scale, which is what lets the coefficient array become sparse for functions with isolated features, and the bridge from kernels 26.08.01 to wavelets is precisely the move from a single bandwidth — one resolution everywhere — to a multiresolution analysis whose coefficients can be selected one at a time, restoring the spatial adaptivity that a fixed-resolution method structurally cannot have.

Full proof set Master

The projection-estimator risk decomposition and optimal rate, the unbiasedness and variance of ${\hat{\beta }}_j$, the soft- and hard-threshold derivations, and the Haar orthonormality are proved in full in the Key theorem and Exercises. The remaining Master claims are recorded here.

Proposition 1 (oracle inequality for soft thresholding, single coordinate). In the model $y=\theta +\sigma z$, $z\sim N(0,1)$, the soft-threshold estimator $\hat{\theta }={\eta }_{\lambda }(y)=\mathrm{sign}(y)(|y|-\lambda {)}_{+}$ with $\lambda =\sigma \sqrt{2\log N}$ obeys $\mathbb{E}(\hat{\theta }-\theta {)}^2\le (2\log N+1)({\sigma }^2{N}^{-1}+\min ({\theta }^2,{\sigma }^2))$. Summing over $N$ independent coordinates gives Theorem 2.

Proof (sketch following Donoho-Johnstone). Two regimes. For the noise-dominated regime $\theta =0$: $\mathbb{E}{\eta }_{\lambda }(\sigma z{)}^2={\sigma }^2\mathbb{E}{\eta }_{\lambda /\sigma }(z{)}^2$, and a Gaussian tail computation gives $\mathbb{E}{\eta }_t(z{)}^2\le (t^2+1)\cdot 2\bar{\Phi }(t)$ where $\bar{\Phi }$ is the standard normal survival function; at $t=\sqrt{2\log N}$, $2\bar{\Phi }(t)\le 1/(t\sqrt{\pi }N)\le N^{-1}$, so the per-coordinate null risk is $\le {\sigma }^2(2\log N+1)N^{-1}$, the source of the additive ${\sigma }^2$ term after summation. For the signal regime, soft thresholding is a non-expansive map with $|{\eta }_{\lambda }(y)-\theta |\le |y-\theta |+\lambda$ pointwise and ${\eta }_{\lambda }$ keeps $\hat{\theta }$ within $\lambda$ of $y$, so a case split on $|\theta |\le \sigma$ versus $|\theta |>\sigma$ bounds $\mathbb{E}(\hat{\theta }-\theta {)}^2$ by $(2\log N+1)\min ({\theta }^2,{\sigma }^2)$ up to the additive null term; the $2\log N$ factor is exactly ${\lambda }^2/{\sigma }^2$. Combining the two regimes coordinatewise gives the stated bound. The constant $2\log N$ is unimprovable in the worst case: it is the price of controlling all $N$ noise coordinates simultaneously through the maximal inequality $\mathbb{E}{\max }_i|\sigma z_i|\sim \sigma \sqrt{2\log N}$. $\square$

Proposition 2 (the universal threshold from the Gaussian maximal inequality). For $z_1,\dots ,z_N$ i.i.d. $N(0,1)$, $\mathrm{Pr}({\max }_i|z_i|>\sqrt{2\log N}(1+\delta ))\to 0$ for every $\delta >0$, so at $\lambda =\sigma \sqrt{2\log N}$ the soft-threshold estimator sets every coordinate of pure noise to zero with probability tending to one.

Proof. By the Gaussian tail $\bar{\Phi }(t)\le \frac{1}{t\sqrt{2\pi }}{e}^{-t^2/2}$ and a union bound, $\mathrm{Pr}({\max }_i|z_i|>t)\le 2N\bar{\Phi }(t)\le \frac{2N}{t\sqrt{2\pi }}{e}^{-t^2/2}$. At $t=(1+\delta )\sqrt{2\log N}$, $e^{-t^2/2}=N^{-(1+\delta {)}^2}$, so the bound is $\le \frac{2}{(1+\delta )\sqrt{4\pi \log N}}{N}^{1-(1+\delta {)}^2}\to 0$ since $(1+\delta {)}^2>1$. Thus ${\max }_i|\sigma z_i|\le \lambda$ eventually, and ${\eta }_{\lambda }$ annihilates each such coordinate. The matching lower bound ${\max }_i|z_i|/\sqrt{2\log N}\to 1$ a.s. shows $\sqrt{2\log N}$ is the exact threshold at which noise is just killed; any smaller multiple leaves spurious survivors. $\square$

Proposition 3 (Besov norm = weighted ${\ell }_p$ of wavelet coefficients). For a sufficiently regular wavelet basis, the Besov seminorm $\Vert f{\Vert }_{B_{p,q}^s}$ is equivalent to the sequence norm $\left({\sum }_j{2}^{jq(s+1/2-1/p)}\right({\sum }_k|{\beta }_{j,k}{|}^p{)}^{q/p}{)}^{1/q}$ on the wavelet coefficients.

Proof (statement; full proof in the cited monograph). The equivalence is the wavelet characterization of Besov spaces: the analysis/synthesis operators between $f$ and its coefficient array $({\beta }_{j,k})$ are bounded both ways on $B_{p,q}^s$ provided the wavelet has more than $s$ vanishing moments and matching smoothness, which the cascade construction supplies for Daubechies wavelets of high enough order. The consequence used in Theorem 3 is structural: a thresholding rule, which acts coordinatewise on $({\beta }_{j,k})$, is natural for Besov balls because the ball is itself a coordinatewise (weighted-${\ell }_p$) constraint, whereas a Sobolev ball is a weighted-${\ell }_2$ ellipsoid for which linear filtering is natural. The mismatch — Besov for $p<2$ versus the ${\ell }_2$ geometry of linear estimators — is the precise reason linear estimators fail to achieve the minimax rate over Besov balls and thresholding succeeds. This is stated without independent proof here; see ^{[Donoho, D. L. & Johnstone, I. M. — Ideal spatial adaptation by wavelet shrinkage]} and the Härdle-Kerkyacharian-Picard-Tsybakov monograph. $\square$

Proposition 4 (variance of the linear wavelet density estimator). With ${\hat{\beta }}_{j,k}=n^{-1}{\sum }_i{\psi }_{j,k}(X_i)$ and the wavelet normalization $\Vert {\psi }_{j,k}{\Vert }_2=1$ on compact support of width $2^{-j}$, $\mathrm{Var}({\hat{\beta }}_{j,k})\le n^{-1}\int {\psi }_{j,k}^{2}f\le n^{-1}\Vert f{\Vert }_{\infty }$, and summing the $O(2^j)$ nonvanishing coefficients at scale $j$ over scales up to $J$ gives total variance $\le C{2}^J/n$.

Proof. The single-coefficient variance bound is the wavelet instance of Exercise 3: $\mathrm{Var}({\hat{\beta }}_{j,k})=n^{-1}\mathrm{Var}({\psi }_{j,k}(X_1))\le n^{-1}\mathbb{E}{\psi }_{j,k}(X_1{)}^2=n^{-1}\int {\psi }_{j,k}^{2}f$, and since $\int {\psi }_{j,k}^2=1$ with $f\le \Vert f{\Vert }_{\infty }$ on the support, this is $\le n^{-1}\Vert f{\Vert }_{\infty }$. At scale $j$ the support partition of $[0,1]$ has $O(2^j)$ cells, so $O(2^j)$ coefficients are nonzero, contributing $O(2^j/n)$ to the variance; summing the geometric series over $j\le J$ gives $O(2^J/n)$. Balancing this against the Sobolev bias tail $O(2^{-2Jp})$ reproduces $2^J\asymp n^{1/(2p+1)}$ and the rate $n^{-2p/(2p+1)}$, matching the projection theorem with the resolution $2^J$ playing the role of the cutoff $J+1$. $\square$

Connections Master

The empirical distribution and Glivenko-Cantelli theory 45.05.02 is the parent device: every coefficient estimate ${\hat{\beta }}_j=n^{-1}{\sum }_k{\varphi }_j(X_k)$ is the empirical measure $P_n$ applied to the basis function ${\varphi }_j$, so the orthogonal-series and wavelet density estimators are plug-in functionals of $P_n$, and their variance bookkeeping is the law-of-large-numbers control of those averages. Where Glivenko-Cantelli gives uniform consistency of $P_n$ over a single function class, this unit organizes the averages into a basis and trades uniform control for the coordinatewise bias-variance trade-off that yields rates.

Kernel density estimation and the bandwidth trade-off 26.08.01 is the linear-estimator counterpart: a kernel estimate is, in the Fourier domain, a multiplier (a linear filter) exactly like the projection estimator's hard cutoff or Pinsker's taper, and it shares the single-resolution limitation that thresholded wavelets overcome. The bandwidth, the cutoff $J$, and the finest scale $2^J$ are three names for one smoothing parameter; the move from kernels to thresholded wavelets is the move from one resolution to a selectable multiresolution.

The lasso, sparsity, and oracle inequalities 45.06.06 shares the exact analytic engine: soft thresholding is the proximal map of $\lambda \Vert \cdot {\Vert }_1$, so wavelet shrinkage is the lasso solved in the orthonormal wavelet design, where the KKT conditions decouple coordinatewise into soft thresholds. The Donoho-Johnstone oracle inequality and the lasso oracle inequality are the same theorem in the Gaussian sequence model, with the universal threshold $\sigma \sqrt{2\log n}$ and the lasso penalty $\sigma \sqrt{2\log p/n}$ both produced by the Gaussian maximal inequality bounding the largest noise coordinate.

Fourier series and the $L^2$-isometry of coefficients 02.10.01 supplies the orthonormal-expansion machinery — Bessel's inequality, Parseval's identity $\Vert f{\Vert }_2^2={\sum }_j{\beta }_j^2$, and completeness — on which the whole risk decomposition rests, since the Pythagorean split of the integrated error into a coefficient-variance sum and a squared-bias tail is Parseval applied to ${\hat{f}}_J-f$. Wavelets extend that global trigonometric system into a localized multiresolution basis, replacing one frequency index by a scale-and-position index.

Historical & philosophical context Master

Orthogonal-series density estimation dates to Čencov (1962) and Schwartz (1967), who proposed estimating the Fourier-coefficient sequence of a density by sample averages and truncating; the method sat alongside the kernel estimators of Rosenblatt (1956) and Parzen (1962) as the two classical routes to nonparametric density estimation. The sharp minimax theory of linear estimators over Sobolev ellipsoids was completed by Pinsker (1980), whose exact asymptotic constant for the $L^2$ risk remains the benchmark against which adaptive procedures are measured ^{[Tsybakov, A. B. — Introduction to Nonparametric Estimation]}.

Wavelets entered statistics after the multiresolution-analysis framework of Mallat and Meyer (1986-1989) and Daubechies's (1988) compactly supported orthonormal wavelets gave a basis that is simultaneously localized in space and frequency and computable by a fast cascade transform. The decisive statistical step was taken by David Donoho and Iain Johnstone in a sequence of papers beginning with the 1994 Biometrika article on ideal spatial adaptation ^{[Donoho, D. L. & Johnstone, I. M. — Ideal spatial adaptation by wavelet shrinkage]}, which introduced the universal threshold $\sigma \sqrt{2\log n}$, the oracle inequality against the ideal diagonal-projection risk, and the proof that simple coordinatewise thresholding is near-minimax over a whole scale of Besov classes at once — a single estimator, ignorant of the unknown smoothness, that matches the best estimator tailored to each class up to a logarithmic factor. With Kerkyacharian and Picard they established the Besov-ball adaptivity that explains why thresholding succeeds precisely where linear methods fail: the Besov norm is a coordinatewise constraint on wavelet coefficients, and functions with localized features inhabit Besov but not Sobolev balls. The convergence of this line with the ${\ell }_1$-penalization of Tibshirani's lasso, both governed by soft thresholding and the Gaussian maximal inequality, is recorded in the high-dimensional-statistics literature of the following two decades.

Bibliography Master

@article{DonohoJohnstone1994,
  author  = {Donoho, David L. and Johnstone, Iain M.},
  title   = {Ideal spatial adaptation by wavelet shrinkage},
  journal = {Biometrika},
  volume  = {81},
  number  = {3},
  year    = {1994},
  pages   = {425--455}
}

@article{DJKP1995,
  author  = {Donoho, David L. and Johnstone, Iain M. and Kerkyacharian, G\'erard and Picard, Dominique},
  title   = {Wavelet shrinkage: asymptopia?},
  journal = {Journal of the Royal Statistical Society, Series B},
  volume  = {57},
  number  = {2},
  year    = {1995},
  pages   = {301--369}
}

@book{Wasserman2006,
  author    = {Wasserman, Larry},
  title     = {All of Nonparametric Statistics},
  publisher = {Springer},
  series    = {Springer Texts in Statistics},
  address   = {New York},
  year      = {2006}
}

@book{Tsybakov2009,
  author    = {Tsybakov, Alexandre B.},
  title     = {Introduction to Nonparametric Estimation},
  publisher = {Springer},
  series    = {Springer Series in Statistics},
  address   = {New York},
  year      = {2009}
}

@book{HKPT1998,
  author    = {H\"ardle, Wolfgang and Kerkyacharian, G\'erard and Picard, Dominique and Tsybakov, Alexander},
  title     = {Wavelets, Approximation, and Statistical Applications},
  publisher = {Springer},
  series    = {Lecture Notes in Statistics 129},
  address   = {New York},
  year      = {1998}
}

@article{Pinsker1980,
  author  = {Pinsker, Mark S.},
  title   = {Optimal filtering of square-integrable signals in {Gaussian} noise},
  journal = {Problems of Information Transmission},
  volume  = {16},
  number  = {2},
  year    = {1980},
  pages   = {120--133}
}

@article{Daubechies1988,
  author  = {Daubechies, Ingrid},
  title   = {Orthonormal bases of compactly supported wavelets},
  journal = {Communications on Pure and Applied Mathematics},
  volume  = {41},
  number  = {7},
  year    = {1988},
  pages   = {909--996}
}

@article{Cencov1962,
  author  = {\v{C}encov, Nikolai N.},
  title   = {Evaluation of an unknown distribution density from observations},
  journal = {Soviet Mathematics Doklady},
  volume  = {3},
  year    = {1962},
  pages   = {1559--1562}
}