Wasserman Zhou 2010 - Orthogonal Series Density Estimation

Summary

For $r=1$ and a density in a Sobolev ellipsoid of order $\gamma$, the minimax squared-$L_2$ rate is $n^{-2\gamma /(2\gamma +1)}$. An exponential mechanism built on the $L_2$ distance between orthogonal-series estimates achieves only $n^{-\gamma /(2\gamma +1)}$ (with $k\asymp \sqrt{n}$, $\gamma >1$; Thm. 6.2), while perturbing the series coefficients with Laplace noise attains the minimax rate for any $k\ge n$ (Thm. 6.3). Again perturbation beats the exponential mechanism, though the authors do not know whether their exponential-mechanism bound is tight.

Overview

The section develops “an exponential scheme based on density estimation” and compares it with a perturbation approach. For simplicity $r=1$. The mathematical structure: expand $p$ in an orthonormal basis, estimate the leading coefficients, and either (a) sample a sanitized database whose own series estimate is close to the data’s estimate (exponential mechanism) or (b) add Laplace noise to the coefficient estimates and sample from the resulting density.

Main Content

Set-up (pp. 15-16)

Let $\{1,{\psi }_1,{\psi }_2,\dots \}$ be an orthonormal basis for $L_2(0,1)=\{f:{\int }_0^1{f}^2<\infty \}$ and suppose $p\in L_2(0,1)$:

$$p(x)=1+\sum _{j=1}^{\infty }{\beta }_j{\psi }_j(x),{\beta }_j={\int }_0^1{\psi }_j(x)p(x)dx.$$

Assume the basis is uniformly bounded:

$$c_0\equiv \underset{x}{\sup }\underset{j}{\sup }|{\psi }_j(x)|<\infty .\text{\hspace{1em}(13)}$$

(Lemma 9.3 uses the cosine basis ${\psi }_j(x)=\sqrt{2}\cos (\pi jx)$, so $c_0=\sqrt{2}$.)

Sobolev ellipsoid (p. 16)

For $\gamma >1/2$ and $C>0$, $B(\gamma ,C)=\{\beta =({\beta }_1,{\beta }_2,\dots ):{\sum }_{j=1}^{\infty }{\beta }_j^2{j}^{2\gamma }\le C^2\}$ and $\mathcal{P}(\gamma ,C)=\{p(x)=1+{\sum }_{j\ge 1}{\beta }_j{\psi }_j(x):\beta \in B(\gamma ,C)\}$.

The minimax rate in $L_2$ norm for $\mathcal{P}(\gamma ,C)$ is $n^{-2\gamma /(2\gamma +1)}$ (Efromovich 1999): ${\inf }_{\hat{p}}{\sup }_{P\in \mathcal{P}(\gamma ,C)}\mathbb{E}\int (\hat{p}-p{)}^2\ge c_1{n}^{-2\gamma /(2\gamma +1)}$. It is achieved by the projection estimator

$$\hat{p}(x)=1+\sum _{j=1}^{m_n}{\hat{\beta }}_j{\psi }_j(x),m_n=n^{1/(2\gamma +1)},{\hat{\beta }}_j=n^{-1}\sum _{i=1}^n{\psi }_j(X_i).\text{\hspace{1em}(14)}$$

For $u\in L_2(0,1)$ write $\Vert u{\Vert }_{{\ell }_2}=({\int }_0^1|u(x){|}^{2}dx{)}^{1/2}$.

Exponential mechanism on the series loss

Take the loss

$$\xi (X,Z)=(\int ({\hat{p}}_X(x)-{\hat{p}}_Z^{\ast }(x){)}^{2}dx{)}^{1/2}=\Vert {\hat{p}}_X-{\hat{p}}^{\ast }{\Vert }_{{\ell }_2},\text{\hspace{1em}(15)}$$

where ${\hat{p}}_X$ is (14) from the data and, from the released sample $Z$,

$${\hat{p}}^{\ast }(x)=1+\sum _{j=1}^{m_k}{\hat{\beta }}_j^{\ast }{\psi }_j(x),m_k=k^{1/(2\gamma +1)},{\hat{\beta }}_j^{\ast }=k^{-1}\sum _{i=1}^k{\psi }_j(Z_i).\text{\hspace{1em}(16)}$$

Lemma 6.1 (p. 17)

Under this scheme $\Delta \le \frac{2c_0^2{m}_n}{n}$. Hence the exponential-mechanism density satisfies

$$g(z\mid x)=\exp (-\frac{\alpha \Vert {\hat{p}}^{\ast }-\hat{p}{\Vert }_{{\ell }_2}}{\Delta })\le \exp (-\frac{\alpha n\Vert {\hat{p}}^{\ast }-\hat{p}{\Vert }_{{\ell }_2}}{2c_0^2{m}_n})\text{almost surely.}\text{\hspace{1em}(17)}$$

(The printed exponent in (17) is $\alpha \Vert \cdot \Vert /\Delta$ whereas definition (9) has $\alpha \Vert \cdot \Vert /(2\Delta )$; with the factor 2 the privacy parameter is $\alpha$, without it $2\alpha$. This affects only constants, not rates.)

Proof (Sec. 9.11, pp. 32-33)

Let $X=(x,X_2,\dots ,X_n)$ and $Y=(y,X_2,\dots ,X_n)$ so $\delta (X,Y)=1$. By the reverse triangle inequality in ${\ell }_2$: $|\xi (X,Z)-\xi (Y,Z)|=|\Vert u{\Vert }_{{\ell }_2}-\Vert v{\Vert }_{{\ell }_2}|\le \Vert u-v{\Vert }_{{\ell }_2}=\Vert {\hat{p}}_X-{\hat{p}}_Y{\Vert }_{{\ell }_2}$ with $u={\hat{p}}_X-{\hat{p}}^{\ast }$, $v={\hat{p}}_Y-{\hat{p}}^{\ast }$. Pointwise, $|{\hat{p}}_X(x)-{\hat{p}}_Y(x)|=|\frac{1}{n}{\sum }_{j=1}^{m_n}({\psi }_j(X_1)-{\psi }_j(Y_1)){\psi }_j(x)|\le \frac{1}{n}{\sum }_{j=1}^{m_n}(|{\psi }_j(X_1)|+|{\psi }_j(Y_1)|)|{\psi }_j(x)|\le \frac{2c_0^2{m}_n}{n}$. Integrating over $[0,1]$ gives the ${\ell }_2$ bound.

Theorem 6.2 - exponential mechanism rate (p. 17)

Let $Z=(Z_1,\dots ,Z_k)$ be drawn from $g_x(z)$ in (17). Assume $\gamma >1$. If $k\asymp \sqrt{n}$ then

$${\rho }^2(p,{\hat{p}}^{\ast })=O_P\left(n^{-\gamma /(2\gamma +1)}\right).$$

The sanitized estimator converges at a slower rate than the minimax rate $n^{-2\gamma /(2\gamma +1)}$ (the exponent halves).

Proof outline (Sec. 9.11, pp. 33-35)

Follow the proof of Thm. 5.3 with $\rho (F,{\hat{F}}_u)=\Vert p-{\hat{p}}_u{\Vert }_{{\ell }_2}$ and $\rho ({\hat{F}}_X,{\hat{F}}_u)=\Vert {\hat{p}}_X-{\hat{p}}_u{\Vert }_{{\ell }_2}$; the triangle inequalities $\Vert {\hat{p}}_u-{\hat{p}}_X\Vert \ge \Vert {\hat{p}}_u-p\Vert -\Vert {\hat{p}}_X-p\Vert$ and $\le \Vert {\hat{p}}_u-p\Vert +\Vert p-{\hat{p}}_X\Vert$ replace the cdf ones. Condition (10) follows from Markov’s inequality: $P(\Vert {\hat{p}}_X-p\Vert >ϵ)\le \mathbb{E}\Vert {\hat{p}}_X-p{\Vert }^2/{ϵ}^2$ with polynomially decaying MSE. The small ball probability is the bottleneck:

1. By Parseval, $\int (p-{\hat{p}}_k{)}^2={\sum }_{j\le m_k}({\hat{\beta }}_j-{\beta }_j{)}^2+{\sum }_{j>m_k}{\beta }_j^2\le Q+c{k}^{-2\gamma /(2\gamma +1)}$ with $Q={\sum }_{j\le m_k}({\hat{\beta }}_j-{\beta }_j{)}^2$, the second term being the Sobolev tail.
2. Standardise: $U_i=({\psi }_1(X_i)-{\beta }_1,\dots ,{\psi }_{m_k}(X_i)-{\beta }_{m_k}{)}^T$ with covariance ${\Sigma }_k$, $Y_i={\Sigma }_k^{-1/2}{U}_i$, $S=k^{-1/2}{\sum }_i{Y}_i$ so $Q=k^{-1}{S}^T{\Sigma }_{k}S$. If ${\lambda }_k$ is the largest eigenvalue of ${\Sigma }_k$ then $P(Q\le {\delta }^2)\ge P(S^{T}S\le k{\delta }^2/{\lambda }_k)$. Lemma 9.3: $\lambda =\mathrm{limsup}{\lambda }_k<\infty$ (for the cosine basis, ${\Sigma }_{jj}=1+{\beta }_{2j}/\sqrt{2}-{\beta }_j^2$, and off-diagonal entries involve only ${\beta }_{|j-k|},{\beta }_{j+k},{\beta }_{j/2},{\beta }_{2k}$ with $\sum |{\beta }_j|=O(1)$, giving a bounded operator norm).
3. Multivariate Berry-Esseen (Bentkus 2003): ${\sup }_c|P(S^{T}S\le c)-P({\chi }_{m_k}^2\le c)|=O\left(\sqrt{m_k^3/k}\right)=O\left(k^{-(\gamma -1)/(2\gamma +1)}\right)$. This is where $\gamma >1$ is needed.
4. The ${\chi }^2$ tail (Rohde-Duembgen): $P({\chi }_m^2\le m+a)\ge 1-e^{-a^2/(4(m+a))}$. Write the target as ${\rho }^2(p,{\hat{p}}^{\ast })\le {ϵ}_n$ with ${ϵ}_n=c_1{n}^{-\gamma /(2\gamma +1)}$, i.e. norm radius $\sqrt{{ϵ}_n}$, and take $k=\sqrt{n}$. Since $k^{-2\gamma /(2\gamma +1)}=n^{-\gamma /(2\gamma +1)}$ is of the same order as ${ϵ}_n$, the bias term is absorbed and one obtains $P(\Vert p-{\hat{p}}_k{\Vert }_{{\ell }_2}\le \sqrt{{ϵ}_n}/2)\ge \frac{1}{2}-O\left(k^{-(\gamma -1)/(2\gamma +1)}\right)$, bounded away from zero.
5. Plug into the Thm. 5.3-type bound with $\Delta \le 2c_0^2{m}_n/n$. The first term is at most $c_2\exp \{k\log {\sup }_{x}p-\alpha c_{3}n\sqrt{{ϵ}_n}/m_n\}$. Here $k\log \sup p\asymp n^{1/2}$ while $n\sqrt{{ϵ}_n}/m_n\asymp n^{1-\gamma /(2(2\gamma +1))-1/(2\gamma +1)}=n^{3\gamma /(2(2\gamma +1))}$. Because $\frac{3\gamma }{2(2\gamma +1)}>\frac{1}{2}$ iff $\gamma >1$, the negative term dominates and the bound is $c_2\exp \{-\alpha c_4{n}^{3\gamma /(2(2\gamma +1))}\}\to 0$. Hence the theorem (and this is the second place, after Berry-Esseen, where $\gamma >1$ is used).

Perturbation approach

Perturbed series estimator (p. 17)

$\hat{q}(x)=1+{\sum }_{j=1}^{m_n}({\hat{\beta }}_j+{\nu }_j){\psi }_j(x)$ with ${\nu }_1,\dots ,{\nu }_m$ i.i.d. Laplace, density $g(\nu )=\frac{n\alpha }{2c_{0}m}{e}^{-n\alpha |\nu |/(c_{0}m)}$ (as printed). If $\hat{q}(x)<0$ for some $x$, replace $\hat{q}$ by $\hat{q}(x)I(\hat{q}(x)>0)/\int \hat{q}(s)I(\hat{q}(s)>0)ds$ (Hall and Murison 1993). Let $Z=(Z_1,\dots ,Z_k)$ be an i.i.d. sample from $\hat{q}$. In the notation of Lemma 2.6, $R=({\nu }_1,\dots ,{\nu }_m)$, so for any $k$ the release is DP.

(Sensitivity check, own derivation: by the computation in Lemma 6.1 each ${\hat{\beta }}_j$ moves by at most $2c_0/n$ when one record is replaced, so the coefficient vector has ${\ell }_1$-sensitivity at most $2c_{0}m/n$; $\alpha$-DP needs Laplace scale $2c_{0}m/(n\alpha )$, versus the printed scale $c_{0}m/(n\alpha )$. The printed scale is smaller by a factor 2, the same kind of constant-level slip as in (17); rates are unaffected.)

Theorem 6.3 - perturbation attains the minimax rate (p. 17)

Let $Z=(Z_1,\dots ,Z_k)$ be drawn from $\hat{q}$, assume $\gamma >1$ and choose $k\ge n$. Then

$${\rho }^2(p,{\hat{p}}_Z)=O_P\left(n^{-2\gamma /(2\gamma +1)}\right),$$

where ${\hat{p}}_Z$ is the orthogonal series density estimator based on $Z$.

“Hence, again, the perturbation technique achieves the minimax rate of convergence and so appears to be superior to the exponential mechanism. We do not know if this is because the exponential mechanism is inherently less accurate, or if our bounds for the exponential mechanism are not tight enough.”

(Heuristic for why perturbation is cheap (own calculation, not in the paper). With Laplace scale $2c_{0}m/(n\alpha )$ per coefficient, the added integrated variance is $m\cdot 2(2c_{0}m/(n\alpha ){)}^2=O(m^3/(n^2{\alpha }^2))$, of lower order than the sampling variance $m/n$ whenever $m^2\ll n{\alpha }^2$. At $m_n=n^{1/(2\gamma +1)}$ this holds for any $\gamma >1/2$ and fixed $\alpha$, so the private estimator inherits the usual bias-variance balance $m_n^{-2\gamma }+m_n/n\asymp n^{-2\gamma /(2\gamma +1)}$. The paper’s stronger assumption $\gamma >1$ is carried over from Thm. 6.2.)

Examples

Rate comparison for $\gamma =2$ (smooth density, own arithmetic). Exponential mechanism: $n^{-\gamma /(2\gamma +1)}=n^{-2/5}$; perturbation and minimax: $n^{-2\gamma /(2\gamma +1)}=n^{-4/5}$. At $n=10^4$: $0.025$ versus $0.00025$ in rate scale (a factor $n^{2/5}=40$ in the exponential mechanism’s disfavour, ignoring constants).

import numpy as np
rng = np.random.default_rng(1)
 
def private_cosine_density(x, m, alpha, rng=rng, c0=np.sqrt(2)):
    """Perturbed series estimator: beta_hat_j + Laplace noise (scale 2*c0*m/(n*alpha), from the Lemma 6.1 sensitivity)."""
    n = len(x)
    j = np.arange(1, m + 1)
    beta_hat = np.sqrt(2) * np.cos(np.pi * np.outer(j, x)).mean(axis=1)
    noisy = beta_hat + rng.laplace(0, 2 * c0 * m / (n * alpha), size=m)
    grid = np.linspace(0, 1, 2001)
    q = 1 + (noisy[:, None] * np.sqrt(2) * np.cos(np.pi * np.outer(j, grid))).sum(0)
    q = np.clip(q, 0, None); q /= np.trapezoid(q, grid)
    return grid, q
 
x = rng.beta(10, 10, 5000)
grid, q = private_cosine_density(x, m=8, alpha=1.0)
from scipy.stats import beta
print(np.trapezoid((q - beta.pdf(grid, 10, 10))**2, grid))      # integrated squared error

With seed 1 this gives an integrated squared error of about $0.0019$ for $n=5000$, $\alpha =1$ and $m=8$ cosine terms (so $n\alpha =5000$, larger than the $n\alpha \le 100$ settings of the histogram simulation in Wasserman Zhou 2010 - Rate Tables, Simulation and Open Questions, so the two are not directly comparable). It illustrates that for a smooth density a handful of series coefficients carry the signal and need little noise.

Connections

• Marketing measurement. Smooth continuous variables (days-to-conversion, spend, order value) can be released as a handful of Laplace-noised series coefficients instead of a histogram; the number of coefficients $m$ plays the role of the bin count and the same bias-variance-privacy trade-off applies ($m\asymp n^{1/(2\gamma +1)}$ non-privately). Smoother distributions allow fewer coefficients and therefore less noise: a practical argument for parametric or low-dimensional summaries of response curves before privatisation. Related: nonparametric curve estimation in Nonparametric Models Overview.

See Also

• Wasserman Zhou 2010 - Rate Tables, Simulation and Open Questions
• Wasserman Zhou 2010 - Histogram Release Mechanisms
• Wasserman Zhou 2010 - Exponential Mechanism and Small Ball Probabilities