Compressive Statistical Learning with Random Feature Moments

TL;DR

This paper introduces a framework called compressive statistical learning that compresses large datasets into low-dimensional sketches for efficient learning, demonstrated on PCA, clustering, and Gaussian mixture models.

Contribution

It presents a novel approach to resource-efficient large-scale learning by using random feature moments to create informative data sketches.

Findings

Sketch size controls generalization error

Framework effectively applied to PCA, clustering, and Gaussian mixture models

Reduces computational resources needed for large-scale learning

Abstract

We describe a general framework -- compressive statistical learning -- for resource-efficient large-scale learning: the training collection is compressed in one pass into a low-dimensional sketch (a vector of random empirical generalized moments) that captures the information relevant to the considered learning task. A near-minimizer of the risk is computed from the sketch through the solution of a nonlinear least squares problem. We investigate sufficient sketch sizes to control the generalization error of this procedure. The framework is illustrated on compressive PCA, compressive clustering, and compressive Gaussian mixture Modeling with fixed known variance. The latter two are further developed in a companion paper.

Tables

Table 1. Table 1: Summary of the application of the framework on our three main examples (detailed in Section  4 and the companion paper (Gribonval et al., 2020 ) ) in 𝒵 = ℝ d 𝒵 superscript ℝ 𝑑 \mathcal{Z}=\mathbb{R}^{d} . 𝕊 k − 1 subscript 𝕊 𝑘 1 \mathbb{S}_{k-1} denotes the ( k − 1 ) 𝑘 1 (k-1) -dimensional simplex (i.e. the sphere with respect to the ℓ 1 superscript ℓ 1 \ell^{1} -norm in the non-negative orthant of ℝ k superscript ℝ 𝑘 \mathbb{R}^{k} ), and ‖ x ‖ 𝚺 = x T ​ 𝚺 − 1 ​ x subscript norm 𝑥 𝚺 superscript 𝑥 𝑇 superscript 𝚺 1 𝑥 \mathopen{}\mathclose{{}\left\|x}\right\|_{{\boldsymbol{\Sigma}}}=x^{T}{\boldsymbol{\Sigma}}^{-1}x the Mahalanobis norm associated to the positive definite covariance matrix 𝚺 𝚺 {\boldsymbol{\Sigma}} . The order of the sketch size is indicated up to universal numerical multiplicative factor and logarithmic dependencies on the parameters δ 𝛿 \delta and ζ 𝜁 \zeta from Theorem  5.7 . The displayed average kernels are up to a multiplicative constant.

Task	PCA	$k$-med./means ($p=1$/$p=2$)	Gaussian Mixture Model.
Hypothesis $h$	subspace	$k$ cluster centers	mixture ${\pi }_h$ of $k$ Gaussians
	$h\subset {ℝ}^d$	$c_1,\dots ,c_k\in {ℝ}^d$	-means $c_l\in {ℝ}^d$
	$dimh=k$		-covar. ${𝚺}_l\in {ℝ}^{d\times d}$
			-mixture parameters ${\alpha }_l$
Loss $\ell (x,h)$	${\Vert x-P_{h}x\Vert }_2^2$	${\min }_{1\le l\le k}{\Vert x-c_l\Vert }_2^p$	$-\log {\pi }_h\left(x\right)$
Feature function	quadratic polyn.	weighted Fourier features	Fourier features
$\Phi \left(x\right)$	$\frac{1}{\sqrt{m}}{\left(x^T{𝐋}_{j}x\right)}_{j=1}^m$	$\frac{1}{\sqrt{m}}{\left(\frac{e^{ȷ{\omega }_j^{T}x}}{w\left({\omega }_j\right)}\right)}_{j=1}^m$	$\frac{1}{\sqrt{m}}{\left(e^{ȷ{\omega }_j^{T}x}\right)}_{j=1}^m$
Sampling law $\Lambda$	$\mathbb{P}\left(𝐋\right)\propto e^{-{\Vert 𝐋\Vert }_F^2}$	$\mathbb{P}\left(\omega \right)\propto w^2\left(\omega \right)e^{-\frac{s^2{\Vert \omega \Vert }^2}{2}}$	$\mathbb{P}\left(\omega \right)\propto e^{-\frac{s^2{\Vert \omega \Vert }_{{𝚺}^{-1}}^2}{2}}$
Average kernel $\kappa (x,x^{\prime })$	${\Vert x{x}^T-x^{\prime }{x}^{\prime T}\Vert }_F^2$	$\exp \left(-\frac{{\Vert x-x^{\prime }\Vert }_2^2}{2s^2}\right)$	$\exp \left(-\frac{{\Vert x-x^{\prime }\Vert }_{𝚺}^2}{2s^2}\right)$
Proxy
$R(h,𝐲)$	Low-rank  recovery	$\underset{\alpha \in {𝕊}_{k-1}}{\min }{\Vert \sum _{l=1}^k{\alpha }_l\Phi \left(c_l\right)-𝐲\Vert }_2$	${\Vert 𝒜\left({\pi }_h\right)-𝐲\Vert }_2$
from  (27)
Restrictions on	N/A	${\min }_{c_l\ne c_{l^{\prime }}}{\Vert c_l-c_{l^{\prime }}\Vert }_2\ge 2\epsilon$	${\min }_{c_l\ne c_{l^{\prime }}}{\Vert c_l-c_{l^{\prime }}\Vert }_{𝚺}\ge 2\epsilon$
hypothesis class		${\max }_l{\Vert c_l\Vert }_2\le R$	${\max }_l{\Vert c_l\Vert }_{𝚺}\le R$
$\mathscr{H}$ when optim-		$\epsilon :=4s\sqrt{\log \left(ek\right)}$	$\epsilon :=4\sqrt{\left(2+s^2\right)\log \left(ek\right)}$
mizing proxy		$w\left(\omega \right):=1+\frac{s^2{\Vert \omega \Vert }_2^2}{d}$	known covariance ${𝚺}_l=𝚺,\forall l$
Sketch size $m$	$kd$	$k^{2}d\log \left(ekdR/\epsilon \right){\log }^2\left(ek\right)$	$k^{2}d\log \left(ekdR\right){\log }^2\left(ek\right)$
			when $s=\sqrt{d}$, cf (Gribonval et al., 2020)(Table 1) for other values of $s$.

$x\in 𝒵$	sample and sample space
$𝐲$	sketch vector (3)
$\Phi$	sketching function (3), (43),
$𝒜$	sketching operator (5)
$\pi ,\tau$	probabilities on sample space
$\mu$, $\nu$	measures on sample space
$⟨\pi ,f⟩$	${𝔼}_{X\sim \pi }f\left(X\right)$
$⟨\mu ,f⟩$	$\int f\left(x\right)𝑑\mu \left(x\right)$ (App. A.2)
$h$	hypothesis
$\mathscr{H}$	class of hypotheses
$\ell (\cdot ,h)$	loss function
$ℛ,\Delta {ℛ}_h$	risk (1), excess risk (Def. 2.4)
$h^{\star }=h_{\pi }^{\star }$	best hypothesis (1)
$R$	proxy for the risk (4), (27)
$\hat{h}$	learned hypothesis (4)
$P_h$	projection function for comp.-type task (Def. 3.1)
$ℒ=ℒ\left(\mathscr{H}\right)$	class of loss functions (9)
$\Delta ℒ=\Delta ℒ\left(\mathscr{H}\right)$	class of loss differences (17)
$ℱ={\left\{{\varphi }_{\omega }\right\}}_{\omega \in \Omega }$	class of features (Def. 5.1)
$\Lambda$	probability distribution of feature parameters $\omega$ (Def. 5.1)
$\kappa (x,x^{\prime })$	psd kernel (Def. 5.2)
$\kappa (\pi ,{\pi }^{\prime })$	kernel mean embedding (44)
$C_{𝒜}$,$C_{\kappa }$, $C_{{\kappa }_{\Phi }}$	kernel constants (11); (47); (48)
${\Vert \mu \Vert }_{𝒢}$	${sup}_{f\in 𝒢}\left|⟨\mu ,f⟩\right|$ ((8), App. A.2)
${\Vert \mu \Vert }_{\kappa }$	MMD norm (45), (58)
$\parallel \cdot {\parallel }_{ℒ}$, $\parallel \cdot {\parallel }_{\Delta ℒ}$	task-driven norms (9); (19)
$D_h\left(\pi \parallel {\pi }^{\prime }\right)$	excess-risk divergence (18)
$d_h(\pi ,𝔖)$,	bias term wrt. model (25)
$d_{ℱ}(\pi ,{\pi }^{\prime })$	feature-based metric (54)
$𝔖$	model set (of probabilities)
${𝔖}_h$	probabilities s.t. $h$ optimal (26)
${𝔖}_h^{\mathrm{𝙲𝚃}}$, ${𝔖}^{\mathrm{𝙲𝚃}}\left(\mathscr{H}\right)$	compression-type model set (31)
${𝔖}_h^{\mathrm{𝙼𝙻}}$, ${𝔖}^{\mathrm{𝙼𝙻}}\left(\mathscr{H}\right)$	max. likelihood model set (28)
$𝒮={𝒮}_{\kappa }\left(𝔖\right)$	normalized secant set (49)
$\Vert ℰ\Vert$	radius of a set of measures (51)
$c_{\kappa }\left(t\right)$	concentration function (52)
$N(\parallel \cdot \parallel ,A,\epsilon )$	covering numbers (Def. 5.6)
$ℬ$	Closed ball (Def. A.2)

Equations

$$\mathcal{R}(\pi ,h):={\mathbb{E}}_{X\sim \pi }\ell (X,h),$$

$$h_{\pi }^{\star }\in \arg \underset{h\in \mathcal{H}}{\min }\mathcal{R}(\pi ,h).$$

$$\mathcal{R}(\pi ,\hat{h})-\mathcal{R}(\pi ,h^{\star })\le {\eta }_n={\eta }_n(\zeta ),$$

$$\mathbf{y}:=\mathtt{Sketch}(\mathbf{X}):=\frac{1}{n}\sum _{i=1}^n\mathrm{\Phi }(x_i)\in {\mathbb{R}}^m\text{or}{\mathbb{C}}^m;$$

$$\hat{h}\in \arg \underset{h\in \mathcal{H}}{\min }R(\mathbf{y},h)$$

$$\mathbf{y}=\mathtt{Sketch}(\mathbf{X})=\frac{1}{n}\sum _{i=1}^n\mathrm{\Phi }(x_i)=\mathcal{A}({\hat{\pi }}_n),$$

$$\mathcal{A}(\pi ):={\mathbb{E}}_{X\sim \pi }\mathrm{\Phi }(X).$$

$$\tilde{\pi }=\Delta [\mathbf{y}];$$

$$\hat{h}\in \arg \underset{h\in \mathcal{H}}{\min }\mathcal{R}(\tilde{\pi },h).$$

$$\underset{h\in \mathcal{H}}{\sup }|\mathcal{R}(\pi ,h)-\mathcal{R}(\tilde{\pi },h)|\le \frac{1}{2}{\eta }_n,$$

$$\mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathcal{G}}:=\sup_{f\in\mathcal{G}}\mathopen{}\mathclose{{}\left\lvert\mathbb{E}_{X\sim\pi}f(X)-\mathbb{E}_{X^{\prime}\sim\pi^{\prime}}f(X^{\prime})}\right\rvert,$$

$$\mathcal{L}(\mathcal{H}):=\mathopen{}\mathclose{{}\left\{\ell(\cdot,h):h\in\mathcal{H}}\right\}.$$

$$\mathopen{}\mathclose{{}\left\|\widetilde{\pi}-\pi}\right\|_{\mathcal{L}}\lesssim d(\pi,\mathfrak{S})+\mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\hat{\pi}_{n})}\right\|_{2}$$

$$\mathopen{}\mathclose{{}\left\|{\tau}^{\prime}-{\tau}}\right\|_{\mathcal{L}}\leq C_{\mathcal{A}}\mathopen{}\mathclose{{}\left\|\mathcal{A}({\tau}^{\prime})-\mathcal{A}({\tau})}\right\|_{2}\quad\forall{\tau},{\tau}^{\prime}\in\mathfrak{S}.$$

$$\Delta[\mathbf{y}]:=\mathrm{argmin}_{{\tau}\in\mathfrak{S}}\mathopen{}\mathclose{{}\left\|\mathcal{A}({\tau})-\mathbf{y}}\right\|_{2},$$

$$d^{\circ}(\pi,\mathfrak{S}):=\inf_{{\tau}\in\mathfrak{S}}\Big{\{}\mathopen{}\mathclose{{}\left\|\pi-{\tau}}\right\|_{\mathcal{L}}+C_{\mathcal{A}}\mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}({\tau})}\right\|_{2}\Big{\}}.$$

$$\mathcal{R}(\pi,\hat{h})-\mathcal{R}(\pi,h^{\star})\leq 4d^{\circ}(\pi,\mathfrak{S})+4C_{\mathcal{A}}\mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\hat{\pi}_{n})}\right\|_{2}$$

$$\mathcal{R}(\pi ,\hat{h})-\mathcal{R}(\pi ,h^{\star })$$

$$\displaystyle\leq\underbrace{\mathbb{E}_{X\sim\pi}[\mathopen{}\mathclose{{}\left\lvert X-\mathrm{Mean}(\pi)}\right\rvert-\mathopen{}\mathclose{{}\left\lvert X-\mathrm{Med}(\pi)}\right\rvert]}_{=:\mathcal{B}(\pi)}+\mathopen{}\mathclose{{}\left\lvert\mathbf{y}-\mathrm{Mean}(\pi)}\right\rvert.$$

$$d^{\circ }(\pi ,\mathfrak{S})$$

$$\displaystyle\geq\inf_{h\in\mathbb{R}}\Big{\{}\mathbb{E}_{X\sim\pi}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\lvert X-h}\right\rvert}\right]+\mathopen{}\mathclose{{}\left\lvert\mathrm{Mean}(\pi)-h}\right\rvert\Big{\}}$$

$$\displaystyle\geq\mathbb{E}_{X\sim\pi}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\lvert X-\mathrm{Mean}(\pi)}\right\rvert}\right]$$

$$=\mathcal{B}(\pi )+\mathcal{R}(\pi ,h^{\star }).$$

$$\Delta\mathcal{L}(\mathcal{H}):=\mathcal{L}(\mathcal{H})-\mathcal{L}(\mathcal{H})=\mathopen{}\mathclose{{}\left\{g:x\mapsto g(x)=\ell(x,h)-\ell(x,h^{\prime}),h,h^{\prime}\in\mathcal{H}}\right\}.$$

$$\Delta {\mathcal{R}}_{h_0}(\pi ,h):=\mathcal{R}(\pi ,h)-\mathcal{R}(\pi ,h_0)={\mathbb{E}}_{X\sim \pi }[\ell (X,h)-\ell (X,h_0)],$$

$$D_{h_{0}}(\pi\|\pi^{\prime}):=\sup_{h\in\mathcal{H}}\mathopen{}\mathclose{{}\left(\Delta\mathcal{R}_{h_{0}}(\pi,h)-\Delta\mathcal{R}_{h_{0}}(\pi^{\prime},h)}\right).$$

$$D_{h_0}(\pi \Vert {\pi }^{\prime })\le D_{h_0}(\pi \Vert {\pi }^{\prime \prime })+D_{h_0}({\pi }^{\prime \prime }\Vert {\pi }^{\prime }).$$

$$\underset{h_0\in \mathcal{H}}{\sup }{D}_{h_0}(\pi \Vert {\pi }^{\prime })$$

$$\displaystyle=\mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\Delta\mathcal{L}}\leq 2\mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathcal{L}}.$$

$$\mathopen{}\mathclose{{}\left\|{\tau}-{\tau}^{\prime}}\right\|_{\Delta\mathcal{L}}\leq C_{\mathcal{A}}\mathopen{}\mathclose{{}\left\|\mathcal{A}({\tau})-\mathcal{A}({\tau}^{\prime})}\right\|_{2}+\eta,\qquad\forall{\tau},{\tau}^{\prime}\in\mathfrak{S},$$

$$\mathbf{y}$$

$$\displaystyle\widetilde{\pi}\in\mathfrak{S}\text{ satisfying }\mathopen{}\mathclose{{}\left\|\mathcal{A}(\widetilde{\pi})-\mathbf{y}}\right\|_{2}$$

$$\hat{h}\text{ satisfying }\mathcal{R}(\tilde{\pi },\hat{h})$$

$$\forall h_{0}\in\mathcal{H}:\;\;\Delta\mathcal{R}_{h_{0}}(\pi,\hat{h})\leq d_{h_{0}}(\pi,\mathfrak{S})+(2+\nu)C_{\mathcal{A}}\mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\hat{\pi}_{n})}\right\|_{2}+\eta+C_{\mathcal{A}}\varepsilon+\varepsilon^{\prime},$$

Full text

Keywords: Kernel mean embedding, random features, random moments, statistical learning, dimension reduction

1 Introduction

Large-scale machine learning faces a number of fundamental computational challenges, triggered both by the high dimensionality of modern data and the increasing availability of very large training collections. Besides the need to cope with high-dimensional features extracted from images, volumetric data, etc., a key challenge is to develop techniques able to fully leverage the information content and learning opportunities opened by large training collections of millions to billions or more items, with controlled computational resources.

Such training volumes can severely challenge traditional statistical learning paradigms based on batch empirical risk minimization. Statistical learning offers a standardized setting where learning problems are expressed as the optimization of an expected loss, or risk, $\mathcal{R}(\pi,h):=\mathbb{E}_{X\sim\pi}{\ell(X,h)}$ over a parameterized family of hypotheses $\mathcal{H}$ (where $\pi$ is the probability distribution of the training collection). This risk is empirically estimated on a training collection, and parameters that empirically minimize it are seeked, possibly with some regularization. Empirical minimization typically requires access to the whole training collection, either in batch mode or iteratively with one or more passes of stochastic gradient. This can become prohibitively costly when the collection is large and each iteration has non-negligible cost. An alternative is to sub-sample the collection, but this may come at the price of neglecting some important items from the collection. Besides online learning (e.g. Mairal et al., 2010), sampling techniques such as coresets (Feldman and Langberg, 2011) or Nyström’s method (e.g. Rudi et al., 2015) have emerged to circumvent computational bottlenecks and preserve the ability to exploit latent information from large collections.

Can we design an alternative learning framework, with the ability to compress the training collection before even starting to learn? We advocate a possible route, compressive statistical learning, which is inspired by the notion of sketching and is endowed with favorable computational features especially in the context of the streaming and distributed data model (Cormode et al., 2011) (see Section 1.3). Rooted both in the generalized method of moments (Hall, 2005) and in compressive sensing (Foucart and Rauhut, 2012), it leverages techniques from kernel methods such as kernel mean embeddings (Gretton et al., 2007; Sriperumbudur et al., 2010) and random Fourier features (Rahimi and Recht, 2007) to obtain innovative statistical guarantees.

As a trivial example, assume $x,h$ belong to $\mathbb{R}^{d}$, and consider the squared loss \ell(x,h)=\mathopen{}\mathclose{{}\left\|x-h}\right\|^{2}, whose risk minimizer is $\mathbb{E}[X]$. In this specific example, keeping only the $d$ empirical averages of the coordinates of $X$ is obviously sufficient. The vision developed in this paper is that, for certain learning problems, all the necessary information can be captured in a sketch: a vector of empirical (generalized) moments of the collection that captures the information relevant to the considered learning task. Computing the sketch is then feasible in one pass, and a near-minimizer of the risk can be computed from the sketch with controlled generalization error.

This paper is dedicated to show how this phenomenon can be generalized: roughly speaking, can the sketch size be taken to be proportional to the number of “intrinsic parameters” of the learning task? Another fundamental requirement for the sketching operation is to be online. When recording the training collection, it should be possible to update the sketch at almost no additional cost. The original training collection can then be discarded and learning can be performed from the sketch only, potentially leading to privacy-preservation. As shown in the companion paper (Gribonval et al., 2020), a sketching procedure based on random generalized moments meets these requirement for clustering and Gaussian mixture estimation.

1.1 Inspiration from compressive sensing

Another classical example of learning task is (centered) Principal Component Analysis (PCA). In this setting, $x\in\mathbb{R}^{d}$, $h$ is an arbitrary linear subspace of dimension $k$, and the loss is \ell(x,h)=\mathopen{}\mathclose{{}\left\|x-P_{h}x}\right\|_{2}^{2} with $P_{h}$ the orthogonal projector onto $h$. The matrix of second moments $\mathbf{\Sigma}_{\pi}:=\mathbb{E}_{X\sim\pi}XX^{T}$ is known to summarize all the information needed to select the best subspace for a training collection. It thus constitutes a natural sketch (of finite dimension $d^{2}$) of the training set.

A much smaller sketch can in fact be computed. Results from compressive sensing and low-rank matrix completion (Foucart and Rauhut, 2012) allow to compress the matrix of second moments to a sketch of dimension of the order of $kd$ (much smaller that $d^{2}$ when $k\ll d$) from which the best rank-$k$ approximation to $\mathbf{\Sigma}_{\pi}$ can be accurately estimated (this rank-$k$ approximation allows to calculate the PCA with appropriate learning guarantees, as we will see in Section 4). This compression operation is made using random linear projections on $\mathbf{\Sigma}_{\pi}$, which can be seen as random second order moments of the training collection.

We propose to generalize such a sketching procedure to arbitrary random generalized moments. Given a learning task and training collection, we study the following questions:

• •

How can we perform learning from a sketch of the training collection?

• •

What statistical learning guarantees can we obtain with such a procedure?

1.2 Contributions

In this paper, we present a general compressive learning framework.

• •

We describe a generic sketching mechanism with random generalized moments and provide a theoretical learning procedure from the sketched data.

• •

We derive general learning guarantees for sketching with random generalized moments.

In the companion paper (Gribonval et al., 2020), we exploit this framework to establish statistical learning guarantees for compressive clustering and compressive Gaussian mixture estimation. We conclude this paper by briefly discussing the potential impact of the proposed framework and its extensions in terms of privacy-aware learning and of the insight it may bring on the information-theoretic properties of certain convolutional neural networks.

1.3 Related work

Sketching and streaming methods.

Sketches are closely linked with the development of streaming methods (Cormode et al., 2011), in which data items are seen once by the user then discarded. A sketch is a small summary of the data seen at a given time, that can be queried for a particular piece of information about the data. As required by the streaming context, when the database is modified, e.g. by inserting or deleting an element, the subsequent update of the sketch must be very fast. In practice, sketches are often applied in various contexts where the data are stored in multiple places. In this heavily distributed framework, a popular class of sketches is that of linear sketches, i.e. structures such that the sketch of the union of two databases is the sum of their sketches – then the sketch of a database distributed over several parts is simply the sum of all their sketches. The sketch presented in this work is indeed a linear sketch (when considered without the normalization constant $1/n$) and as such, updates operations are excessively simple and fast. Sketches have been used for a large variety of operations (Cormode et al., 2011) such as the popular detection of heavy-hitters (Cormode and Muthukrishnan, 2005; Cormode and Hadjieleftheriou, 2009). Closer to our framework, sketches have been used to approximately maintain histograms (Thaper et al., 2002) or quantiles (Gilbert et al., 2002), however these methods are subject to the well-known curse of dimensionality and are unfeasible even in moderate dimension.

Learning in a streaming context.

Various learning algorithms have also been directly adapted to a streaming context. Examples include the Expectation-Maximization algorithm (Andrieu and Doucet, 2003; Cappé and Moulines, 2009), the $k$-means algorithm (Guha and Mishra, 2016; Ailon et al., 2009), or Principal Component Analysis (Ghashami et al., 2016). In each case, the result of the algorithm is updated as new data becomes available. However these algorithms do not fully benefit from the many advantages of sketches. Sketches are simpler to merge in a distributed context, update operations are more immediate, and the learning step can be delocalized and performed on a dedicated machine.

Coresets.

Another popular class of structures that summarize a database for learning is called coresets. Coresets were initially developed for $k$-means (Har-Peled and Mazumdar, 2004) or, more generally, subspace approximation (Feldman et al., 2010; Feldman and Langberg, 2011) and also applied to learning Gaussian Mixture Models (Feldman et al., 2011; Lucic et al., 2017). In a sense, the philosophy behind coresets is situated halfway between sketches and streaming learning algorithms. Like the sketching approaches, coresets methods construct a compressed representation of the database (or “coreset”), but are somehow closer to already approximately performing the learning task. For instance, the coreset described in (Frahling and Sohler, 2005) already incorporates steps of Lloyd’s $k$-means algorithm in its construction. Similar to the $k$-means++ algorithm (Arthur and Vassilvitskii, 2007), many coresets have been developed as (weighted) adaptive subsampling of the data (Feldman et al., 2011; Lucic et al., 2017).

Linear sketches vs Coresets.

It is in general difficult to compare sketching and coresets methods (including the sketching method presented in this paper) in terms of pure performance or theoretical guarantees, since they are very different approaches that can be more or less adapted to certain contexts. We can however outline some differences. Unlike sketches, coresets are not specifically build for the streaming context, and they may require several passes over the data. Nevertheless they can still be adapted to streams of data (as described e.g. in Har-Peled and Mazumdar, 2004; Feldman and Langberg, 2011; Lucic et al., 2017) by using a merge-and-reduce hierarchical strategy: for each batch of data that arrives sequentially, the user builds a coreset, then groups these coresets and builds a coreset of coresets, and so on. This update method is clearly less direct than updating a linear sketch, and more importantly the user must balance between keeping many coresets and letting the size of the overall summary grow with the number of points in the database, or keeping only highest-level coresets at the cost of losing precision in the theoretical guarantees each time the height of the hierarchical structure increases. As a comparison, the sketch presented in the companion paper (Gribonval et al., 2020)for $k$-means does not have these limitations: like with any linear sketch, updates are totally independent of previous events, and for a fixed sketch size the ability to perform the learning task strictly increases with the number of points.

Generalized Method of Moments and Compressive Sensing.

The methodology that we employ to develop the proposed sketching framework is similar to a Generalized Method of Moments (GeMM) (Landau, 1987; Hall, 2005): the parameters $\theta$ of a model are learned by matching a collection of theoretical generalized moments from the distribution $\pi_{\theta}$ with empirical ones from the data. GeMM is often seen as an alternative to Maximum Likelihood estimation, to obtain different identifiability guarantees (Belkin and Sinha, 2015; Hsu and Kakade, 2013; Anderson et al., 2014) or when the likelihood is not available. Traditionally, a finite number of moments is considered, but modern developments give guarantees when an infinite (integral) number of generalized moments are available (Carrasco and Florens, 2000, 2014), in particular generalized moments associated to the (empirical) characteristic function (Carrasco and Florens, 2002; Feuerverger and Mureika, 1977). Our point of view is slightly different: we consider the collection of moments as a compressed representation of the data and as a means to achieve a learning task.

Compared to the guarantees usually obtained in GeMM such as consistency and efficiency of the estimator $\hat{\theta}$, the results that we obtain are more akin to Compressive Sensing and Statistical Learning. For instance, when learning Gaussian Mixture Models, we prove in the companion paper (Gribonval et al., 2020)that learning is robust to modeling error (the true distribution of the data is not exactly a GMM but close to one), which is generally overlooked in GeMM. In the proof technique, this is done by replacing the so-called “global identifiability condition”, (i.e. injectivity of the moment operator), which is a classical condition in GeMM but is already difficult to prove and sometimes simply assumed by practitioners (see Newey and McFadden, 1994, p. 2127) by the strictly stronger Lower Restricted Isometry Property (LRIP) from the Compressive Sensing literature (Donoho, 2006; Candès et al., 2006; Baraniuk, 2007; Foucart and Rauhut, 2012). This is achieved by considering random feature moments (related to random features (Rahimi and Recht, 2007, 2009; Bach, 2017) and kernel mean embeddings (Sriperumbudur et al., 2010)), so in a sense the resulting Compressive Statistical Learning framework could be considered as a Method of Random Feature Moments. While the LRIP is reminiscent of certain kernel approximation guarantees with random features (see e.g. Sriperumbudur and Szabó, 2015; Bach, 2017), it is in fact of a different nature, and none seems to be a direct consequence of the other.

1.4 Outline

Section 2 describes our general framework for compressive statistical learning. We define here statistical learning guarantees, introduce the required notions and state our general Theorem on statistical learning guarantees for compressive learning. An important concept is the notion of Lower Restricted Isometry Property (LRIP) using the notion of a model set (a set of “simple” probability distributions) which is futher discussed in Section 3. To illustrate the proposed framework, we detail in Section 4 a procedure for Compressive PCA, where we do not intend to match the latest developments in the domain of PCA such as stochastic and incremental PCA (Arora et al., 2012; Balsubramani et al., 2013), kernel PCA using Nyström sampling (Sterge et al., 2020) or random features (Ullah et al., 2018; Sriperumbudur and Sterge, 2020); but rather to give a first illustration. Generic techniques to establish the LRIP property for sketches of controlled size are described in Section 5. In the companion paper (Gribonval et al., 2020), we specify a sketching procedure and state the associated learning guarantees for compressive clustering and compressive Gaussian mixture estimation. We discuss in Section 6 possible extensions of the proposed framework as well as the insight it may bring on the information flow across one layer of a convolutive neural network with average pooling. Finally, all proofs are stated in the Appendix.

2 A general compression framework for statistical learning

This section is dedicated to the introduction of our compressive learning framework.

2.1 Statistical learning

Statistical learning offers a standardized setting where many learning problems (supervised or unsupervised) can be expressed as the optimization of an expected risk over a parameterized family of functions. Formally, we consider a training collection $\mathbf{X}=\{x_{i}\}_{i=1}^{n}\in\mathcal{Z}^{n}$ drawn i.i.d. from a probability distribution $\pi$ on the measurable space $(\mathcal{Z},\mathfrak{Z})$. In our examples, $\mathcal{Z}=\mathbb{R}^{d}$ is endowed with the Borel $\sigma$-algebra $\mathfrak{Z}$. One wishes to select a hypothesis $h$ from a hypothesis class $\mathcal{H}$ to perform the task at hand. How well the task can be accomplished with the hypothesis $h$ is typically measured through a loss function $\ell:(x,h)\mapsto\ell(x,h)\in\mathbb{R}$ and the expected risk associated to $h$:

[TABLE]

where (here and in the sequel) we will always assume that we restrict our attention to probability distributions $\pi$ such that $x\mapsto\ell(x,h)$ is measurable and $\pi$-integrable for all $h\in\mathcal{H}$. In the idealized learning problem, one selects a function $h^{\star}_{\pi}$ that minimizes the expected risk (we will assume existence of this minimum for a simpler presentation, although most statements to come can be transformed if needed using a sequence of approximate minimizers)

[TABLE]

We will use the shorthand $h^{\star}$ for $h^{\star}_{{\pi}}$ whenever there is no ambiguity from the context. In practice one has no access to the true risk $\mathcal{R}(\pi,h)$ since the expectation with respect to the underlying probability distribution, $\mathbb{E}_{X\sim\pi}[\cdot]$, is unavailable. Instead, methods such as empirical risk minimization (ERM) produce an estimated hypothesis $\hat{h}$ from the training dataset $\mathbf{X}$ by minimizing the risk $\mathcal{R}(\hat{\pi}_{n},\cdot)$ (or a regularized version) associated to the empirical probability distribution $\hat{\pi}_{n}:=\tfrac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$ of the training samples. One expects to produce, with high probability at least $1-\zeta$ on the draw of the training set, the bound on the excess risk

[TABLE]

where $\eta_{n}$ typically decays as $1/\sqrt{n}$ or better. We will use the following running examples.

Examples:

• •

PCA: as stated in the introduction, the loss function is \ell(x,h)=\mathopen{}\mathclose{{}\left\|x-P_{h}x}\right\|_{2}^{2} where $P_{h}$ is the orthogonal projection onto the subspace hypothesis $h$ of prescribed dimension $k$.

• •

$k$-means clustering: each hypothesis corresponds to a set of $k$ candidate cluster centers, h=\mathopen{}\mathclose{{}\left\{c_{1},\ldots,c_{k}}\right\}, and the loss is defined by the $k$-means cost \ell(x,h)=\min_{1\leq l\leq k}\mathopen{}\mathclose{{}\left\|x-c_{l}}\right\|_{2}^{2}. The hypothesis class $\mathcal{H}$ may be further reduced by defining constraints on the considered centers (e.g., in some domain, or as detailed in the companion paper (Gribonval et al., 2020)with some separation between centers).

• •

Gaussian Mixture Modeling: each hypothesis $h$ corresponds to the collection of weights, means and variances of a mixture of $k$ Gaussians, whose probability density function is denoted $\pi_{h}(x)$. The loss function is based on the maximum likelihood $\ell(x,h)=-\log\pi_{h}(x)$.

2.2 Compressive learning

Our aim, and one of the major achievements of this paper, is to control the excess risk (2) using an estimate $\hat{h}$ obtained from the sole knowledge of a sketch of the training collection. As we will see, the resulting philosophy for large-scale learning is, instead of addressing an ERM optimization problem of size proportional to the number of training samples, to first compute a sketch vector of size driven by the complexity of the task, then to address a nonlinear least-squares optimization problem associated to the Generalized Method of Moments (GeMM) on this sketch.

Taking its roots in compressive sensing (Donoho, 2006; Candès et al., 2006; Foucart and Rauhut, 2012) and the generalized method of moments (Landau, 1987; Hall, 2005), but also on kernel mean embeddings (Smola et al., 2007; Sriperumbudur et al., 2010), random features (Rahimi and Recht, 2007, 2009; Bach, 2017), and streaming algorithms (Gilbert et al., 2002; Cormode and Muthukrishnan, 2005; Cormode et al., 2011), compressive learning relies on the choice of a measurable (nonlinear) feature function $\mathrm{\Phi}:\mathcal{Z}\rightarrow\mathbb{R}^{m}\ \text{or}\ \mathbb{C}^{m}$ and has two main steps:

Compute generalized empirical moments using the feature function on the training collection to summarize it into a single sketch vector

[TABLE] 2. 2.

Produce a hypothesis from the sketch using an appropriate learning procedure: $\hat{h}=\mathtt{Learn}(\mathbf{y})$.

Overall, the goal is to design the feature function $\mathrm{\Phi}(\cdot)$ and the learning procedure $\mathtt{Learn}(\cdot)$ given a learning task (i.e., a loss function) such that the resulting hypothesis $\hat{h}$ has controlled excess risk (2) (if $\mathrm{\Phi}$ is drawn at random according to some specification, we want (2) to hold with high probablity also with respect to the draw of $\mathrm{\Phi}$). To anticipate on what will be developed in Section 3.1 (notably Eq. (27)), let us mention that learning from a sketch will take the form of a minimization problem

[TABLE]

where in a sense $R(\mathbf{y},\cdot)$ will play the role of a proxy for the empirical risk $\mathcal{R}(\hat{\pi}_{n},\cdot)$.

Trivial examples.

• •

Estimation of the mean: Assume $x,h$ belong to $\mathbb{R}^{d}$, and consider the squared loss \ell(x,h)=\mathopen{}\mathclose{{}\left\|x-h}\right\|^{2}, whose risk minimizer is $\mathbb{E}[X]$. In this specific example, it is obviously sufficient to keep only the $d$ empirical averages of the coordinates of $X$, i.e., to use $\mathrm{\Phi}(x):=x$.

• •

PCA: As the principal components are calculated from the eigenvalue decomposition of the matrix of second moments of the samples, we can simply use $\mathrm{\Phi}(x):=xx^{T}$.

A less trivial example is Compressive PCA. Instead of estimating the full matrix $\mathbf{\Sigma}_{\pi}~{}=~{}\mathbb{E}_{X\sim\pi}XX^{T}$, of size $d\times d$, it is known that computing $m$ random gaussian linear measurements of this matrix makes it possible to manipulate a vector $\mathbf{y}$ of dimension $m$ of the order of $kd$ from which one can accurately estimate the best rank-$k$ approximation to $\mathbf{\Sigma}_{\pi}$, that gives the $k$ first principal components. Nuclear norm minimization is typically used to produce this low rank approximation given the vector $\mathbf{y}$. We will describe this procedure in details in Section 4 as a first illustration of our framework.

In the companion paper (Gribonval et al., 2020), for the more challenging examples of Compressive $k$-means and Compressive Gaussian Mixture Modeling, we provide a feature function $\mathrm{\Phi}$ and a method “$\mathtt{Learn}$” (based on a specific proxy (4) corresponding to a non-convex least-squares minimization) that leads to a control of the excess risk.

As described below, these results are achieved by establishing links with the formalism of linear inverse problems and low complexity recovery (i.e., sparse/structured vector recovery, low-rank matrix recovery) and extending theoretical tools to the setting of compressive statistical learning.

2.3 Compressive learning as a linear inverse problem

The most immediate link with linear inverse problems is the following. The sketch vector $\mathbf{y}$ can be seen as a linear function of the empirical probability distribution $\hat{\pi}_{n}:=\tfrac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$:

[TABLE]

where $\mathcal{A}$ is a linear operator from the set of distributions $\pi$ such that $\mathrm{\Phi}$ is integrable with respect to $\pi$, to $\mathbb{R}^{m}$ (or $\mathbb{C}^{m}$), defined by

[TABLE]

This is linear in the sense that111One can extend $\mathcal{A}$ to a linear operator on the space of finite signed measures such that $\mathrm{\Phi}$ is integrable, see Appendix A.2. $\mathcal{A}(\theta\pi+(1-\theta)\pi^{\prime})=\theta\mathcal{A}(\pi)+(1-\theta)\mathcal{A}(\pi^{\prime})$ for any $\pi,\pi^{\prime}$ and $0\leq\theta\leq 1$.

Since for large $n$ we should have $\mathcal{A}(\hat{\pi}_{n})\approx\mathcal{A}(\pi)$, the sketch $\mathbf{y}$ can be viewed as a noisy linear observation of the underlying probability distribution $\pi$. This viewpoint allows to formally leverage the general methodology of linear inverse problems to produce a hypothesis from the sketch $\mathbf{y}$.

Conceptually, we construct the learning-from-sketch procedure $\hat{h}=\mathtt{Learn}(\mathbf{y})$ in two steps:

• •

Define a so-called decoder $\Delta$ that finds a probability distribution $\widetilde{\pi}$ given a sketch $\mathbf{y}$:

[TABLE]

• •

Find a best hypothesis from this estimate:

[TABLE]

As a first coarse analysis of this scheme, notice that if the decoder step is such that a uniform approximation between the risk of $\widetilde{\pi}$ and of $\pi$ holds:

[TABLE]

then we will be able to control the excess risk (2) – our goal. Indeed, using (6) and the triangle inequality, it is easy to show that (7) directly implies (2). (We will see later that (7) can be too coarse and will introduce a more refined analysis based on excess risks in Section 2.5.) In a way, this is very similar to ERM except that instead of using the empirical risk $\mathcal{R}(\hat{\pi}_{n},\cdot)$, we use an estimate of the risk $\mathcal{R}(\widetilde{\pi},\cdot)$ where $\widetilde{\pi}$ is deduced directly from the sketch $\mathbf{y}$.

Remark 2.1**.**

At first sight, the above conceptual view may wrongly suggest that compressive learning replaces statistical learning with the much more difficult problem of non-parametric density estimation. Fortunately, as we will see, this is not the case, thanks to the fact that our objective is never to accurately estimate $\pi$ in the standard sense of density estimation as, e.g., in (Bertin et al., 2011), but only to accurately estimate the risk $\mathcal{R}(\pi,\cdot)$. On practical examples, a natural decoder will be based on best moment matching over a parametric family of probability distributions, which will be expressed more directly as the minimization of a proxy for the risk (4), cf Section 3.1.

2.4 Statistical learning guarantees: a first control of the excess risk

In this section, for simplicity we first focus on how to establish uniform control of the risks of the form (7) using general results from linear inverse problems (we shall introduce in the next section a sharper but also slightly more notation-heavy analysis). To leverage the links between compressive learning and general inverse problems, we further notice that $\sup_{h\in\mathcal{H}}|\mathcal{R}(\pi,h)-\mathcal{R}(\pi^{\prime},h)|$ can be viewed as a metric on probability distributions. Given a class $\mathcal{G}$ of measurable functions $f:\mathcal{Z}\to\mathbb{R}\ \text{or}\ \mathbb{C}$, we use the following notation throughout this work:

[TABLE]

which defines a semi-norm on the space of finite signed measures (see Appendix A.2) on $(\mathcal{Z},\mathfrak{Z})$ such that all $f\in{\mathcal{G}}$ are integrable. In order to be explicit about the integrability assumptions in the results to come, we will call this space the set of $\mathcal{G}$-integrable finite signed measures (resp. probability distributions, when appropriate).

With this notation, we have \sup_{h\in\mathcal{H}}\mathopen{}\mathclose{{}\left\lvert\mathcal{R}(\pi,h)-\mathcal{R}(\pi^{\prime},h)}\right\rvert=\mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathcal{L}(\mathcal{H})} where

[TABLE]

We will usually abbreviate the latter notation by dropping the dependence on $\mathcal{H}$, considered fixed. The desired guarantee (7) then reads \mathopen{}\mathclose{{}\left\|\pi-\Delta[\mathbf{y}]}\right\|_{\mathcal{L}}\leq\eta_{n}/2.

In the usual context of linear inverse problems, producing an accurate estimate from noisy underdetermined linear observations requires some “regularity” assumption. Such an assumption often takes the form of a “low-dimensional” model set that the quantity to estimate is close to.

Example 2.2**.**

In the case of sparse vector recovery (respectively low-rank matrix recovery), one wishes to estimate $\mathbf{x}\in\mathbb{R}^{n}$ (resp. $\mathbf{X}\in\mathbb{R}^{n\times n}$) from $\mathbf{y}\approx\mathbf{A}\mathbf{x}$ (resp. $\mathbf{y}\approx\mathbf{A}\mathrm{vec}(\mathbf{X})$). Guarantees are achieved when $\mathbf{x}$ is close to the set of $k$-sparse vectors (resp. when $\mathbf{X}$ is close to the set of rank-$r$ matrices).

Similarly here, estimating $\pi$ from $\mathbf{y}\approx\mathcal{A}(\pi)$ may require considering some model set $\mathfrak{S}$, whose choice and definition will be discussed in Section 3.

Remark 2.3**.**

While in classical compressive sensing the model set plays the role of prior knowledge on the data distribution that completes the observations, in the examples considered here we will often obtain distribution free excess risk guarantees using models derived from the loss function.

Given a model set222We will always assume that the models $\mathfrak{S}$ under consideration are such that loss and feature functions are integrable with respect to any distribution belonging to $\mathfrak{S}$, i.e. $\mathfrak{S}$ is both $\mathcal{L}$-integrable and \mathopen{}\mathclose{{}\left\{\mathrm{\Phi}}\right\}-integrable, using the terminology introduced after (8). $\mathfrak{S}$ that plays the role of regularizer, and a sketching operator $\mathcal{A}$, an “ideal” decoder $\Delta$ should be robust to two different sources of error: the distribution of the data $\pi$ generally does not belong to $\mathfrak{S}$ but is “close” to it, introducing some modelling error, and the empirical sketch is used instead of the true generalized moments, which adds some noise. Generalizing early formulations for sparsity-regularized inverse problems, a decoder robust to both noise and modelling error is usually reffered to as instance optimal (Cohen et al., 2009; Bourrier et al., 2014). Mathematically, this can be expressed as: for any distribution $\pi$, any draw of the training samples from $\pi$ (embodied by the empirical distribution $\hat{\pi}_{n}$), with $\mathbf{y}=\mathcal{A}(\hat{\pi}_{n})$ and $\widetilde{\pi}=\Delta[\mathcal{A}(\hat{\pi}_{n})]$

[TABLE]

where $\lesssim$ hides multiplicative constants, and $d(\cdot,\mathfrak{S})$ is some measure of distance to the model set $\mathfrak{S}$. In the rest of the paper, we refer to this first term as “bias”. A significant part of later sections will be devoted to the control of the bias and the choice of a good model set. Proving that a decoder satisfies (10) ultimately serves to establish bounds such as (7) to control the excess risk.

It turns out that general results from abstract linear inverse problems (Bourrier et al., 2014) can be adapted to already characterize the existence of a decoder satisfying property (10). By (Bourrier et al., 2014, Section IV-A), if a decoder with the above property exists then a so-called lower Restricted Isometry Property (LRIP) must hold: there is a finite constant $C_{\mathcal{A}}<\infty$ such that

[TABLE]

Conversely, the LRIP (11) implies (Bourrier et al., 2014, Theorem 7) that the following decoder (also known as ideal decoder)

[TABLE]

which corresponds to best moment matching, is instance optimal, i.e., (10) holds for any $\pi$ and $\hat{\pi}_{n}$, with the particular distance

[TABLE]

As a consequence, the LRIP (11) implies a control of the excess risk achieved with the hypothesis $\hat{h}$ selected with (6), where $\widetilde{\pi}=\Delta[\mathbf{y}]$, as

[TABLE]

where we used explicit constants from (Bourrier et al., 2014, Theorem 7). Note that in the above argument, it was never used that the data is distributed i.i.d. from $\pi$. Estimate (14) therefore holds under this form for any fixed data sample (in fact, for any empirical distribution $\hat{\pi}_{n}$, being understood that it determines $\hat{h}$) and any distribution $\pi$. Of course, the data distributional assumption is useful to control the second term in the bound.

2.5 Improved excess risk analysis

The analysis of the previous section has the merits of simplicity, generality, and using existing results from linear inverse problems. However it has some limitations, in particular when the bias term $d^{\circ}(\pi,\mathfrak{S})$ is not close to zero. To emphasize this point, we consider a simple example and compare the excess risk control (for the same sketched learning procedure) obtained through the general bound (14), to a direct computation specific to this example.

Consider the problem of estimating the median of a distribution on $\mathbb{R}$: we assume $\mathcal{Z}=\mathcal{H}=\mathbb{R}$, and consider the absolute value loss \ell(x,h)=\mathopen{}\mathclose{{}\left\lvert x-h}\right\rvert, whose risk minimizer under the distribution $\pi$ is the median $h^{\star}=\mathrm{Med}(\pi)$. As a sketching operator we take simply $\mathrm{\Phi}(x)=x$, resulting in the sketch given by the empirical mean $\mathbf{y}=\mathcal{A}(\hat{\pi}_{n})=\frac{1}{n}\sum_{i=1}^{n}x_{i}$. Finally, as a model we consider the family of 1-point Dirac measures, \mathfrak{S}=\mathopen{}\mathclose{{}\left\{\delta_{x},x\in\mathbb{R}}\right\} (this is the model consisting of all distributions with vanishing optimal risk, see Section 3.2 for a more general discussion). Obviously, it then holds with these choices that the ideal decoder given by (12) is $\widetilde{\pi}=\Delta[\mathbf{y}]=\delta_{\mathbf{y}}$, and further $\hat{h}=\mathbf{y}$.

On the one hand, the excess risk for this sketching/decoding scheme is bounded as follows, by a simple direct calculation, putting $\mathrm{Mean}(\pi):=\mathbb{E}_{X\sim\pi}[X]$:

[TABLE]

On the other hand, it is easy to check that the LRIP (11) holds, with equality, for 1-point Dirac measures with constant $C_{\mathcal{A}}=1$. If we consider the general bound (14), while we recover (up to factor 4) the second term above \mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\hat{\pi}_{n})}\right\|_{2}=\mathopen{}\mathclose{{}\left\lvert\mathbf{y}-\mathrm{Mean}(\pi)}\right\rvert (of order $\mathcal{O}(1/\sqrt{n})$), the first term in (14) is a bias term which is driven by

[TABLE]

The inequality in the third line above is obtained by noticing that \mathopen{}\mathclose{{}\left\lvert x-\mathrm{Mean}(\pi)}\right\rvert\leq\mathopen{}\mathclose{{}\left\lvert x-h}\right\rvert+\mathopen{}\mathclose{{}\left\lvert\mathrm{Mean}(\pi)-h}\right\rvert for each $x,h\in\mathbb{R}$. Hence, using the general bound (14) instead of the specific direct calculation, we get an additional, unwanted term corresponding to the optimal risk $\mathcal{R}(\pi,h^{\star})$ which is nonzero as soon as $\pi\not\in\mathfrak{S}$, and can become arbitrary large (even if $\mathcal{B}(\pi)=0$, e.g. if $\pi$ is symmetric around its mean).

One reason for this lack of sharpness is that the analysis in the previous section concentrated first on (uniform) control of the risk difference \mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathcal{L}}, to deduce only as a second step a control on the excess risk. It is known from the statistical learning literature that it is generally sharper to directly analyze the excess risk; and correspondingly consider the excess loss class (see, for instance, Bartlett et al., 2005, Section 5 and Koltchinskii, 2006, Section 7):

[TABLE]

However, the risk minimizer $h^{\star}$ depends on the distribution $\pi$; for this reason we will consider a family of excess losses and risks with respect to some reference hypothesis $h_{0}$.

Definition 2.4**.**

The excess risk relative to a reference hypothesis $h_{0}$ is defined as:

[TABLE]

and the associated excess risk divergence with respect to $h_{0}$ is:

[TABLE]

Observe that $\Delta\mathcal{R}_{h_{0}}(\pi,h_{0})=0$ for each $\pi$, hence the latter quantity is nonnegative (although no absolute value is involved in its definition), but not symmetric in general. Yet, it satisfies an (oriented) triangle inequality: for any $\pi,\pi^{\prime},\pi^{\prime\prime}$

[TABLE]

It is therefore a hemimetric (see Definition B.1 in the Appendix).

The excess risk divergence will play a role similar to that of \mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathcal{L}}, and satisfies in particular

[TABLE]

With this setting we have the following result, which can be seen as a refinement of (14).

Theorem 2.5**.**

Consider a loss class $\mathcal{L}(\mathcal{H})$, a feature function $\mathrm{\Phi}$, and a model set $\mathfrak{S}$ such that every probability distribution ${\tau}\in\mathfrak{S}$ is both $\mathcal{L}$-integrable and \mathopen{}\mathclose{{}\left\{\mathrm{\Phi}}\right\}-integrable. Assume that the sketching operator $\mathcal{A}$ associated to $\mathrm{\Phi}$ satisfies the following LRIP inequality:

[TABLE]

for some finite constants $C_{\mathcal{A}}>0$ and $\eta\geq 0$.

Consider any training collection $\mathbf{X}=\{x_{i}\}_{i=1}^{n}\in\mathcal{Z}^{n}$, and denote $\hat{\pi}_{n}:=\tfrac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$. Define

[TABLE]

Then, for any probability distribution $\pi$ that is both $\mathcal{L}$-integrable and \mathopen{}\mathclose{{}\left\{\mathrm{\Phi}}\right\}-integrable:

[TABLE]

where

[TABLE]

Similarly to  (14), the above estimate holds regardless of any distributional assumptions on the training collection $\mathbf{X}$. Nevertheless, estimate (24) is primarily of interest when $\mathbf{X}$ is drawn i.i.d. according to $\pi$ and with $h_{0}=h^{\star}_{\pi}$, in which case the left-hand side is the excess risk with respect to the optimum risk, which is what one generally aims at controlling. However, in some situations it may be also helpful to consider excess risk with respect to other reference hypotheses $h_{0}$; this can include situations where $h^{\star}_{\pi}$ itself is not well-defined if the infimum of the risk is not attained.

As compared to (14), we observe that this result is more general, as it allows for a $(\nu,\varepsilon)$-approximate decoder (22), an $\varepsilon^{\prime}$-approximate ERM (23), an $\eta$-approximate LRIP condition (20); more importantly, the main bound (24) involves the sharper excess risk divergence rather than the loss norm \mathopen{}\mathclose{{}\left\|.}\right\|_{\mathcal{L}}. It may also be useful to consider $\pi=\hat{\pi}_{n}$, to predict the quality compared to the empirical risk minimizer.

Moreover, inequality (19) implies that the lower LRIP condition (11) considered in the previous section implies the relaxed LRIP condition (20) (with $\eta=0$, and up to a factor 2 in the constant), so establishing (11) is sufficient in order to obtain the improved inequality (24).

The proof of Theorem 2.5 follows the structure outlined in the previous section, but requires to formally extend the result of (Bourrier et al., 2014, Section IV-A) (leading from the LRIP (11) to instance optimality (10)) to the case of a hemimetric, $\eta$-approximate LRIP and $\varepsilon$-approximate decoder. These technical aspects are relegated to Appendix B.

Discussion:

• •

Computing the sketch (21) is highly parallelizable and distributable. Multiple sketches can be easily aggregated and updated as new data become available.

• •

As discussed in Remark 2.1, while (22) may appear as a general nonparametric density estimation problem, in all the examples considered in this paper and the companion one (Gribonval et al., 2020), it is indeed a nonlinear parametric least-squares fitting problem and the existence of the minimizer follows in practice from compactness arguments.

• –

For Compressive PCA (Section 4) it is a low-rank matrix reconstruction problem. Provably good algorithms to estimate its solution have been widely studied.

• –

For Compressive $k$-means and Compressive Gaussian Mixture Modeling (cf the companion paper (Gribonval et al., 2020)), the resulting optimization problem has been empirically addressed through the CL-OMPR algorithm (Keriven et al., 2015, 2018). Algorithmic success guarantees are an interesting challenge. This is however beyond the scope of this paper. We note that the classic (non-compressed) $k$-means problem by minimization of the empirical risk is known to be NP-hard (Garey et al., 1982; Aloise et al., 2009) and that guarantees for approaches such as K-means++ (Arthur and Vassilvitskii, 2007) are only in expectation and with a logarithmic sub-optimality factor.

• •

In Section 3 we discuss choices of the model set $\mathfrak{S}$ that are driven by the learning task only and make the minimization problem (23) trivial to solve. With these choices the combined solution of (22)-(23) is explicitly turned into the minimization of a proxy for the risk, as in (4).

• •

The second term in the bound (24) of the excess risk, $\eta_{n}$, is the empirical estimation error \mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\hat{\pi}_{n})}\right\|_{2}. It is easy to show that it decays as $1/\sqrt{n}$ when the data is drawn i.i.d. according to $\pi$, this will be done explicitly for the considered examples.

For large collection size $n$ drawn i.i.d. according to $\pi$, the term \mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\hat{\pi}_{n})}\right\|_{2} becomes small and (24) shows that compressive learning will benefit from accurate excess risk guarantees provided the model $\mathfrak{S}$ and the feature function $\mathrm{\Phi}$ (or equivalently the sketching operator $\mathcal{A}$) are chosen so that:

the LRIP (20) holds; ideally for a “small” value of $m$, as we also seek to design compact sketches and, eventually, tractable algorithms to learn from them. 2. 2.

the distance $d_{h^{\star}_{\pi}}(\pi,\mathfrak{S})$ is “small”; this vague notion will be exploited in Section 3 below to guide our choice of model set $\mathfrak{S}$, and will be made more concrete on examples.

We illustrate the improvement obtained for the bias term with respect to the coarser analysis on the toy example of median estimation considered previously. In that setting we have $h^{\star}_{\pi}=\mathrm{Med}(\pi)$, and for ${\tau}=\delta_{x}\in\mathfrak{S}$:

[TABLE]

so that

[TABLE]

The inequality in the third line is obtained by using $x=\mathrm{Mean}(\pi)$. Note that the presence of the last term is unavoidable (since \mathopen{}\mathclose{{}\left\lvert x-\mathrm{Med}(\pi_{0})}\right\rvert+\mathopen{}\mathclose{{}\left\lvert x-\mathrm{Mean}(\pi)}\right\rvert\geq\mathopen{}\mathclose{{}\left\lvert\mathrm{Med}(\pi_{0})-\mathrm{Mean}(\pi)}\right\rvert), and that it is still larger than $\mathcal{B}(\pi)$ which we recall is the only bias term appearing in the direct calculation (15). (A situation where it is much larger is the following: assume $\pi=(\frac{1}{2}+\varepsilon)\delta_{0}+(\frac{1}{2}-\varepsilon)\delta_{1}$, for $0<\varepsilon<1/2$. Then $\mathrm{Med}(\pi)=0$, $\mathrm{Mean}(\pi)=\frac{1}{2}-\varepsilon$, \mathbb{E}_{X\sim\pi}\mathopen{}\mathclose{{}\left\lvert X-\mathrm{Mean}(\pi)}\right\rvert=(\tfrac{1}{2}+\varepsilon)(\tfrac{1}{2}-\varepsilon)+(\tfrac{1}{2}-\varepsilon)(\tfrac{1}{2}+\varepsilon)\ =(1+2\varepsilon)(\tfrac{1}{2}-\varepsilon), \mathbb{E}_{X\sim\pi}\mathopen{}\mathclose{{}\left\lvert X-\mathrm{Med}(\pi)}\right\rvert=\mathbb{E}_{X\sim\pi}X=\tfrac{1}{2}-\varepsilon, hence \mathopen{}\mathclose{{}\left\lvert\mathrm{Med}(\pi)-\mathrm{Mean}(\pi)}\right\rvert=\tfrac{1}{2}\ -\varepsilon while $\mathcal{B}(\pi)=2\varepsilon(1/2-\varepsilon)$.) In this sense, even using the improved excess risk analysis, the general bound (24) can lack some tightness. It is nevertheless much sharper than the bound (16), in particular $d_{h^{\star}_{\pi}}(\pi,\mathfrak{S})=0$ as soon as $\mathrm{Mean}(\pi)=\mathrm{Med}(\pi)$, while $d_{h^{\star}_{\pi}}(\pi,\mathfrak{S})>0$ in general in this case, see (16).

3 Task-driven model sets

An important ingredient of the proposed framework is the model set $\mathfrak{S}$ and we now discuss its choice. If prior knowledge on the data distribution is available, it is of course possible to choose $\mathfrak{S}$ to incorporate such knowledge into compressive statistical learning. However, it is more common in statistical learning to seek “distribution-free” statistical guarantees. Therefore, it may be more desirable to derive a model set entirely or mostly from the learning task itself, to the extent possible.

To begin with, we show that for any model set $\mathfrak{S}$, the abstract two-step learning mechanism (cf steps (22)-(23) in Theorem 2.5) can be written as the minimization of a more explicit proxy (4) for the empirical risk. This rewriting exploits a partition of $\mathfrak{S}$ into certain submodels $\mathfrak{S}_{h}$ driven by the learning task, i.e. by the loss family $\{\ell(\cdot,h)\}_{h\in\mathcal{H}}$. For certain learning tasks, we further show that the submodels $\mathfrak{S}_{h}$ and the corresponding proxy (4) have a simple expression provided we choose a “natural”, task-driven, model set $\mathfrak{S}$. Finally, for so-called “compression type” learning tasks, choosing such a task-driven model set allows to control the bias term in the excess risk (24) by a function of the optimal risk $\mathcal{R}(\pi,h^{\star})$.

3.1 Learning from a sketch without explicit density estimation

Consider a family $\mathfrak{S}$ of $\mathcal{L}$-integrable distributions and assume for each $\pi\in\mathfrak{S}$ the risk admits a minimizer, i.e., there is $h\in\mathcal{H}$ such that $\mathcal{R}(\pi,h)=\inf_{h^{\prime}\in\mathcal{H}}\mathcal{R}(\pi,h^{\prime})$. When this holds the model set $\mathfrak{S}$ can be decomposed as $\mathfrak{S}=\cup_{h\in\mathcal{H}}\mathfrak{S}_{h}$ where for each hypothesis $h\in\mathcal{H}$ we define

[TABLE]

and the hypothesis $\hat{h}$ selected using steps (22)-(23) in Theorem 2.5 (with $\varepsilon^{\prime}=0$) is equivalently obtained as a near-minimizer of the following proxy for the risk

[TABLE]

in the sense that \inf_{{\tau}\in\mathfrak{S}}\mathopen{}\mathclose{{}\left\|\mathcal{A}({\tau})-\mathbf{y}}\right\|_{2}=\inf_{h}\inf_{{\tau}\in\mathfrak{S}_{h}}\mathopen{}\mathclose{{}\left\|\mathcal{A}({\tau})-\mathbf{y}}\right\|_{2}. With this expression in hand, it is possible to directly cast the estimation of $\hat{h}$ as (4). To turn this into a concrete proxy for the risk it is helpful to consider a model set $\mathfrak{S}$ such that $\mathfrak{S}_{h}$ has a simple characterization.

3.2 Choosing a model set: with or without prior knowledge ?

Learning tasks such as maximum likelihood estimation directly involve a natural model set for which $\mathfrak{S}_{h}$ as in (26) is easily characterized. Consider the loss $\ell(x,h)=-\log\pi_{h}(x)$ with \mathopen{}\mathclose{{}\left\{\pi_{h},h\in\mathcal{H}}\right\} a parameterized family of distributions with $\pi_{h^{\prime}}\neq\pi_{h}$ for $h^{\prime}\neq h$. In this setting, it is natural to consider the following model set:

[TABLE]

which is nothing more than a statistical model in the usual sense. Moreover, up to a constant additive term, it then holds $\mathcal{R}(\pi,h)=\textnormal{KL}(\pi\|\pi_{h})$, where $\textnormal{KL}(\cdot\|\cdot)$ is the Kullback-Leibler divergence, so that for any $h\in\mathcal{H}$:

[TABLE]

see (Cover and Thomas, 1991, Chapter 9). As a conclusion, the proxy (27) reads R(\mathbf{y},h)=\mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi_{h})-\mathbf{y}}\right\|_{2}.

For many other learning tasks, the choice of the model set $\mathfrak{S}$ results from a tradeoff between several needs. On the one hand, results from compressed sensing suggest that given a model set $\mathfrak{S}$ that has proper “low-dimensional” properties, it is possible to choose a small sketch size $m$ and design the sketching operator $\mathcal{A}$ such that the LRIP (20) holds, and the ideal decoder $\Delta$ in (12) — or its relaxed version in (22) — is guaranteed to stably recover probability distributions in $\mathfrak{S}$ from their compressed version obtained with $\mathcal{A}$. This calls for the choice of a “small” model set. On the other hand, and perhaps more importantly, the model set should not be “too small” in order to ensure that the obtained control of the excess risk is nontrivial.

Ideally, in the common case of compression-type tasks, as defined below, the bias term in the excess risk (24) should be small when the true optimum risk is small, and even vanish when the true optimum risk vanishes, i.e. when $\inf_{h\in\mathcal{H}}\mathcal{R}(\pi,h)=0$.

Definition 3.1**.**

We call the learning task a compression-type task if the loss can be written as $\ell(x,h)=d^{p}(x,P_{h}x)$, where $d$ is a metric on $\mathcal{Z}$, $p>0$, and $P_{h}:\mathcal{Z}\to\mathcal{Z}$ is a “projection function”, i.e.,

[TABLE]

Typical examples of compression-type tasks are PCA, $k$-means, and $k$-medians. For PCA, $P_{h}$ is the orthogonal projector onto subspace $h$. For $k$-means and $k$-medians, $P_{h}$ maps $x\in\mathcal{Z}=\mathbb{R}^{d}$ to the closest center $c_{i}$ from $h=(c_{1},\ldots,c_{k})$, with ties broken arbitrarily. In other words, given an arbitrary Voronoi partition corresponding to $k$ disjoint sets $W_{j}$ such that $\cup_{j}W_{j}=\mathbb{R}^{d}$ and $d(x,c_{j})=\min_{l}d(x,c_{l})$ for each $x\in W_{j}$, $P_{h}x=c_{j}$ if and only if $x\in W_{j}$. Manifold learning tasks where $P_{h}$ is a projection onto a manifold parameterized by $h$ (with ties broken arbitrarily) would also fit under this framework.

For a compression-type task, a natural model set is the family of $\mathcal{L}$-integrable probability distributions

[TABLE]

We consider a few examples:

• •

Compressive PCA: the model set $\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$ consists of all distributions which admit a matrix of second moments of rank at most $k$. Given any $\widetilde{\pi}\in\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$, a minimum risk hypothesis according to (6) is any subspace $\hat{h}$ spanned by eigenvectors associated to the $k$ largest eigenvalues of $\mathbf{\Sigma}_{\widetilde{\pi}}$. More details will be given shortly in Section 4.

• •

Compressive $k$-means or $k$-medians: the model set $\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$ consists of mixtures of $k$ Diracs. Given h=\mathopen{}\mathclose{{}\left\{c_{1},\ldots,c_{k}}\right\} and any $\widetilde{\pi}=\sum_{\ell=1}^{k}\alpha_{\ell}\delta_{c_{\ell}}\in\mathfrak{S}^{\mathtt{CT}}_{h}$, a minimum risk hypothesis according to (6) is $\hat{h}=h$. Since $\mathcal{A}(\delta_{c})=\mathrm{\Phi}(c)$, the proxy (27) reads

[TABLE]

with \mathbb{S}_{k-1}:=\mathopen{}\mathclose{{}\left\{\alpha\in\mathbb{R}^{k}:\alpha_{\ell}\geq 0;\sum_{\ell=1}^{k}\alpha_{\ell}=1}\right\} the simplex.

For compressive PCA we exhibit in Section 4 a feature function $\mathrm{\Phi}$ so that $\mathcal{A}$ satisfies the LRIP (20) with respect to the model set $\mathfrak{S}=\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$. The same is done in the companion paper (Gribonval et al., 2020) for compressive $k$-means, compressive $k$-medians and compressive Gaussian mixture modeling.

3.3 Controlling the bias term for compression-type tasks

The bias term — defined in (25) — is a measure of distance to the model set $\mathfrak{S}$. For compression-type tasks, the particular model set $\mathfrak{S}=\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$ in (31) was designed so that this bias term vanishes when $\pi\in\mathfrak{S}$, and we can further bound the bias term $d_{h^{\star}_{\pi}}(\pi,\mathfrak{S}^{\mathtt{CT}}(\mathcal{H}))$ with an increasing function of the true minimum risk, $\mathcal{R}(\pi,h^{*})$. This leads to recovery guarantees providing distribution-free excess risk guarantees. Whether this holds for other learning tasks, or even generically, is a challenging question left to further work.

The following lemmas allow to obtain an upper bound of the bias term in function of the minimum risk in a number of relevant cases. For a probability distribution $\pi$ on $\mathcal{Z}$ and $P_{h}$ as in Definition 3.1, we denote $P_{h}\pi$ the push-forward of $\pi$ through $P_{h}$, i.e., the probability distribution of a random variable $Y=P_{h}X$ where $X\sim\pi$. Given a loss class $\mathcal{L}(\mathcal{H})$ we recall that $h^{\star}_{\pi}=\arg\min_{h\in\mathcal{H}}\mathcal{R}(\pi,h)$.

Lemma 3.2**.**

Consider a compression-type task on the input space $\mathcal{Z}$. Then

• •

$\mathfrak{S}^{\mathtt{CT}}_{h}$* is the set of probability distributions on $X\in\mathcal{Z}$ such that $X\in\mathcal{E}_{h}:=P_{h}\mathcal{Z}$ almost surely. With the model set $\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$ and the loss class $\mathcal{L}(\mathcal{H})$, the bias term (25) satisfies (for any $h_{0}\in\mathcal{H}$)*

[TABLE]

• •

If $d^{p}$ is a metric (in particular if $p\leq 1$) then $D_{h}(\pi\|P_{h}\pi)=0$ for any $h\in\mathcal{H}$ and $\mathcal{L}$-integrable distribution $\pi$.

Remark 3.3**.**

When $d^{p}$ is not a metric there are $u,v,w\in\mathcal{Z}$ such that $d^{p}(u,v)>d^{p}(u,w)+d^{p}(w,v)$. The loss $\ell(x,h):=d^{p}(x,h)$, with \mathcal{H}:=\mathopen{}\mathclose{{}\left\{v,w}\right\}, defines a compression-type task with $P_{h}x=h$ for all $x\in\mathcal{Z}$. Set $\pi=\delta_{u}$. Since $d(u,w)<d(u,v)$ we have $h^{\star}_{\pi}=w$. We also have for $h=v\in\mathcal{H}$

[TABLE]

Hence, one cannot generically obtain $D_{h}(\pi\|P_{h}\pi)=0$, not even with the restriction $h=h^{\star}_{\pi}$.

For $p=2$, $d$ the Euclidean distance on $\mathbb{R}^{d}$, and $h_{0}=h^{\star}_{\pi}$ (which we recall is generally the primary interest case since our main bound (24) then gives a control of the excess risk with respect to the optimum), we still have $D_{h^{\star}_{\pi}}(\pi\|P_{h^{\star}_{\pi}}\pi)=0$ for certain tasks. In light of the above remark, this is a nontrivial property which is established for PCA in Lemma E.1, and for $k$-means in the companion paper (Gribonval et al., 2020). Beyond these specific situations, it is possible (under somewhat generic additional assumptions) to bound the two terms appearing in (33) by (a power of) the risk itself, as established next.

Lemma 3.4**.**

Consider a compression-type task where $(\mathcal{Z},d)$ is a separable metric space. Then

• •

Assume that $\mathcal{Z}$ has $d$-diameter bounded by $B$. Then for any $p>1$, $h\in\mathcal{H}$ and $\mathcal{L}$-integrable distribution $\pi$:

[TABLE]

• •

Assume that \mathrm{\Phi}:(\mathcal{Z},d)\to(\mathbb{R}^{m},\mathopen{}\mathclose{{}\left\|\cdot}\right\|_{2}) (or (\mathbb{C}^{m},\mathopen{}\mathclose{{}\left\|\cdot}\right\|_{2})) is $L$-Lipschitz. Then, for $h\in\mathcal{H}^{\prime}$ and $p\geq 1$:

[TABLE]

For $h\in\mathcal{H}$ and $p\leq 1$, if the space $\mathcal{Z}$ has $d$-diameter bounded by $B$:

[TABLE]

The proofs of the above lemmas are in Appendix D. For Lemma 3.4, optimal transport is exploited through connections between the considered norms and the norm \mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathrm{Lip}(L,d)}=L\cdot\mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathrm{Lip}(1,d)}, where $\mathrm{Lip}(L,d)$ denotes the class of functions $f:(\mathcal{Z},d)\to\mathbb{R}$ that are $L$-Lipschitz. The two lemmas can be combined to express an “explicit” bound on $d_{h^{\star}_{\pi}}$, this is postponed to concrete examples.

4 Illustration with Compressive PCA

As a first simple illustration, this general compressive statistical framework can be applied to the example of PCA, where most of the tools already exist. Our aim is essentially illustrative, and focuses on controlling the excess risk, rather than to compare the results with state-of-the art PCA techniques.

Definition of the learning task.

The risk associated to the PCA learning problem is defined333for simplicity we assume centered distributions $\mathbb{E}_{X\sim\pi}X=0$ and don’t empirically recenter the data. as \mathcal{R}_{k-\mathtt{PCA}}(\pi,h)=\mathbb{E}_{X\sim\pi}\mathopen{}\mathclose{{}\left\|X-P_{h}X}\right\|_{2}^{2} with $P_{h}$ the orthogonal projector onto subspace $h$. It is minimized by any subspace $h^{\star}_{\pi}$ associated with $k$ largest eigenvalues of the matrix $\mathbf{\Sigma}_{\pi}=\mathbb{E}_{X\sim\pi}XX^{T}$.

It is well established (Foucart and Rauhut, 2012) that matrices that are approximately low-rank can be estimated from partial linear observations under a certain Restricted Isometry Property (RIP). This leads to the following natural way to perform Compressive PCA.

Choice of a model set.

The ”natural” model set from (31) is $\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})=\{\pi:\textrm{rank}(\mathbf{\Sigma}_{\pi})\leq k\}$. More generally we can consider as a model set $\mathfrak{S}_{r}:=\{\pi:\textrm{rank}(\mathbf{\Sigma}_{\pi})\leq r\}$, with $r\geq k$, so that $\mathfrak{S}_{r}\supset\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$.

Choice of feature function.

Choose (at random) a linear operator $\mathcal{M}:\mathbb{R}^{d\times d}\to\mathbb{R}^{m}$ satisfying (with high probability) the RIP on low-rank matrices: for any $\mathbf{M}\in\mathbb{R}^{d\times d}$ of rank at most $2r$,

[TABLE]

with \mathopen{}\mathclose{{}\left\|\cdot}\right\|_{F} the Frobenius norm and $\delta<1$. This is feasible with $m$ of the order of $rd$, by taking the Frobenius inner product of $\mathbf{M}$ with $m$ independent random Gaussian matrices (see e.g. Foucart and Rauhut, 2012).

Given these facts one can define the feature function as $\mathrm{\Phi}:\mathcal{Z}=\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}$ by $\mathrm{\Phi}(x):=\mathcal{M}(xx^{T}).$

Sketch computation.

Given sample points $x_{1},\ldots,x_{n}$ in $\mathbb{R}^{d}$, compute the sketch $\mathbf{y}$ as in (3), i.e., compute empirical estimates of random second moments of the distribution $\pi$ of $X$.

Learning from a sketch.

Given a sketch vector $\mathbf{y}$, estimate a solution of the optimization problem over positive semi-definite (p.s.d.) symmetric matrices ($\mathbf{\Sigma}\succcurlyeq 0$)

[TABLE]

This step estimates the rank-${r}$ p.s.d. matrix whose sketch best matches that of the empirical matrix of second moments, in the least squares sense. Compute the eigen-decomposition $\hat{\mathbf{\Sigma}}=\mathbf{U}\mathbf{D}\mathbf{U}^{T}$ and output

[TABLE]

In Appendix E we control the excess risk of PCA by relating the excess risk divergence $D_{h_{0}}(\pi\|\pi^{\prime})$ — with $h_{0}\in\mathcal{H}$ an arbitrary hypothesis — to the Frobenius norm \mathopen{}\mathclose{{}\left\|\mathbf{\Sigma}_{\pi}-\mathbf{\Sigma}_{\pi^{\prime}}}\right\|_{F}, and upper bounding the “bias” term (25) appearing in the generic bound of Theorem 2.5, to obtain the following result:

Theorem 4.1**.**

Consider any probability distribution $\pi$ with finite second moments and any draw of $x_{i}$, $1\leq i\leq n$ (represented by the empirical distribution $\hat{\pi}_{n}$). Applying the above approach yields, for any $s$, $1\leq s\leq r$:

[TABLE]

where $\lambda_{i}(\mathbf{\Sigma}_{\pi})$ are the eigenvalues of $\mathbf{\Sigma}_{\pi}$ ranked in decreasing order (with multiplicity), $c_{\delta}:=2\sqrt{2}\frac{\sqrt{1+\delta}}{\sqrt{1-\delta}}$, $c^{\prime}_{\delta}:=2\sqrt{2}/\sqrt{1-\delta}$. In particular:

[TABLE]

Discussion:

• •

Bias term. The first term in the right hand side of (41) is a bias term that vanishes when the true risk is low. Since it is proportional to the true risk, it leads to the (non-sharp) oracle inequality \mathcal{R}_{k-\mathtt{PCA}}(\pi,\hat{h})\leq C_{\delta}(k,r)\mathcal{R}_{k-\mathtt{PCA}}(\pi,h^{\star}_{\pi})+c^{\prime}_{\delta}\sqrt{k}\mathopen{}\mathclose{{}\left\|\mathcal{M}(\mathbf{\Sigma}_{\pi}-\mathbf{\Sigma}_{\hat{\pi}_{n}})}\right\|_{2}. We show in (Gribonval et al., 2020), (using Lemma 3.4 and (36)) that this type of property also holds for Compressive $k$-medians; for Compressive $k$-means we prove similar properties where the bias term is bounded by the square root of the true risk (using (34),(35)).

It is notable that the bias multiplier $C_{\delta}(k,r)$ is of order $\sqrt{k}$ if we use the natural model set ($r=k$), but drops to a constant independent of $k$ as soon as we choose e.g. the larger model set $\mathfrak{S}_{r}$ with $r=2k$. Thus, there appears to be a clear advantage, in the sense of the obtained bound, in choosing a reconstruction model that is larger than the natural model set $\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$, while not significantly changing the magnitude of the number of required data sketches. At this point it is unclear to us if the inflation of the bias factor for the natural model is unavoidable or is just a technical artefact.

• •

Sample complexity. Regarding the second term, if we further assume that the support of $\pi$ is contained in a Euclidean ball of radius $R$, then by the RIP (37) we have a.s. \mathopen{}\mathclose{{}\left\|\mathcal{M}(xx^{T})}\right\|_{2}\leq\sqrt{1+\delta}\cdot R^{2} hence, by the vectorial Hoeffding’s inequality (see e.g. Pinelis, 1992), we obtain with high probability w.r.t. data sampling that \sqrt{k}\mathopen{}\mathclose{{}\left\|\mathcal{M}(\mathbf{\Sigma}_{\pi})-\mathcal{M}(\mathbf{\Sigma}_{\hat{\pi}_{n}})}\right\|_{2} is of the order of $R^{2}\sqrt{k/n}$.

• •

Root-$n$ consistency in a high-dimensional scenario. As noticed in the previous point, the statistical estimation error term in the sketched learning bound is of order $\sqrt{k/n}$. Consider a high-dimensional situation where $d$ is large and growing with $n$, and assume a polynomial spectral decay $\lambda_{j}(\mathbf{\Sigma}_{\pi})\leq j^{-\alpha}$ with $\alpha>1$. Then by choosing $s=r/2$, the bias term in (40) is of order $\sqrt{kr^{-(2\alpha-1)}}$. As a consequence, it is sufficient to take $r$ of order $\min(d,n^{\frac{1}{2\alpha-1}})$ so as to ensure that the bias term is at most of the same order as the statistical error term. This gives an advantage compared to the standard approach of storing all $d^{2}$ second order moments, as soon as $d>n^{\frac{1}{2\alpha-1}}$. For comparison, standard statistical analysis of PCA based on the uniform bound on the deviations of the empirical risk from its expectation (see e.g. Shawe-Taylor et al., 2005) leads to a control of order444More refined techniques (Blanchard et al., 2007; Reiß and Wahl, 2020) can lead to a convergence rate of the PCA excess risk of order $n^{-1}$ asymptotically, for $k$ fixed, but depending on eigenvalue gaps. For the present discussion we compare ourselves to the simplest analysis available in the standard learning context. $\sqrt{k/n}$. Hence, using the sketched approach we can reduce storage/memory imprint significantly while keeping statistical guarantees of the same order as in the standard setting.

• •

Relation to efficient kernel-PCA methods. Methods in recent literature have been proposed to make kernel PCA more efficient, relying on Nyström subsampling Sterge et al. (2020) or on approximation of the kernel using random features (Ullah et al., 2018; Sriperumbudur and Sterge, 2020). Even when restricting attention to a linear kernel, a direct comparison to our approach proves delicate. The Nyström subsampling method is very much taylored to the dual point of view, which is canonical in kernel methods (i.e. approximation of the kernel Gram matrix): this implicitly posits to store all training points in order to represent the output of the method as a kernel expansion (the computational gain concerns the storage and manipulation of the $(n,n)$ Gram matrix). On the other hand, the random feature approach could be construed as closer in spirit to ours, however rather than storing generalized moments as we do, it is in essence a (random) dimension reduction of the individual input points from kernel space to a finite-dimensional feature space where regular PCA is performed. Also the theoretical works of Ullah et al. (2018); Sriperumbudur and Sterge (2020) concern reconstruction in $L^{2}(\pi)$ space and not in the original kernel space norm (lifting back the PCA projection found in approximate feature space into original kernel space proves a delicate question). In contrast, we consider sketching the data into generalized empirical moments but propose a reconstruction method directly in the relevant space with theoretical guarantees. It is to be noted however, that the results of Ullah et al. (2018) posit a black-box PCA method applied in the finite (but high-)dimensional approximate feature space, and one could apply a sketching approach at this stage. This suggests the interesting proposal that efficient kernel PCA methods could be combined, rather than be in competition with, the sketching approach, though we did not push the idea further.

Practical algorithms for learning and comparison to prior PCA-specific results.

One can consider several relaxations of the nonconvex optimization problem (38) in order to perform compressive PCA. Beside convex relaxations using the minimization of the nuclear norm (Foucart and Rauhut, 2012, Section 4.6), Kabanava et al. (2016) showed (in a complex-valued setting) that the rank constraint in (38) can be relaxed when $\mathcal{M}$ is made of random rank-one projections, i.e. when \mathrm{\Phi}(x)=\frac{1}{\sqrt{m}}(|\mathopen{}\mathclose{{}\left\langle a_{j},x}\right\rangle|^{2})_{j=1,m} where $a_{j}\in\mathbb{C}^{d}$ are independent standard complex Gaussian vectors. In this setting, let

[TABLE]

and the corresponding hypothesis $\hat{h}$ obtained through (39). Combining (Kabanava et al., 2016, Theorem 4 with $p=2$) with Equation (72) in Section E and Equation (63) in Section B, we have the following result: if $m\geq Ckd$ where $C$ is a universal constant, then with high probability on the draw of the $a_{j}$, for any $x_{1},\ldots,x_{n}$, we have the control

[TABLE]

where $D_{1},D_{2}$ are positive universal constants that do not depend on $k$.

Hence, provided we use a model set of dimension $2r$ as discussed above, the error control (41) obtained via our general approach matches what can be obtained using directly the PCA-specific study of Kabanava et al. (2016). Two practical advantages of the latter are (a) that (42) is a convex program, and (b) that the sketches are made using rank-one matrices, which are cheaper to store. Still, the guarantees obtained by our general approach is able to match prior results for setting-specific methods. It will further permit the study of the less trivial setting of compressive clustering and compressive Gaussian mixture estimation as shown in the companion paper (Gribonval et al., 2020).

5 Establishing the LRIP for random sketching operators

In this section, we investigate how to establish the LRIP (20) (with $\eta=0$) when the sketching operator $\mathcal{A}$ is associated to random features. The approach uses connections with the notion of kernel mean embedding of probability distributions.

5.1 Random features and kernel mean embeddings

Definition 5.1** (Random feature map).**

Consider $\mathcal{F}:=\{{\phi_{\omega}}\}_{\omega\in\Omega}$ a parameterized family of (real- or complex-valued) measurable functions, $\Lambda$ a probability distribution $\Lambda$ over the parameter set $\Omega$ (often $\Omega=\mathbb{R}^{d}$), and a sketch size $m$. A random feature map is defined by drawing $m$ i.i.d parameters $(\omega_{j})_{j=1}^{m}$ according to $\Lambda$ and setting

[TABLE]

Any draw of the feature function $\mathrm{\Phi}$ defines a positive semi-definite kernel between samples $\kappa_{\mathrm{\Phi}}(x,x^{\prime}):=\langle\mathrm{\Phi}(x),\mathrm{\Phi}(x^{\prime})\rangle_{\mathbb{R}^{m}}$ (or $\langle\mathrm{\Phi}(x),\mathrm{\Phi}(x^{\prime})\rangle_{\mathbb{C}^{m}}$). Compressive learning is deeply connected to kernel mean embeddings of probability distributions, as the related sketching operator $\mathcal{A}$ defines a so-called kernel mean embedding between probability distributions which are $\mathcal{F}$-integrable.

Definition 5.2** (Kernel mean embedding, Mean Map Discrepancy (Gretton et al., 2007; Sriperumbudur et al., 2010)).**

Any positive semi-definite kernel $\kappa(\cdot,\cdot)$ in the sample space is associated to a Mean Map Embedding (a kernel between distributions). By abuse of notation, we keep the notation $\kappa$ for both the expression of the kernel in the sample space and of the corresponding kernel for probability distributions with appropriate integrability

[TABLE]

The associated Maximum Mean Discrepancy (MMD) metric is

[TABLE]

The average kernel $\kappa$ associated to $(\mathcal{F},\Lambda)$, will play a key role in establishing the LRIP. Given $x,x^{\prime}\in\mathcal{Z}$, the expectation of $\kappa_{\mathrm{\Phi}}(x,x^{\prime})=\frac{1}{m}\sum_{j=1}^{m}{\phi_{\omega}}_{j}(x)\overline{{\phi_{\omega}}_{j}}(x^{\prime})$ over the draws of $\omega_{j}$ is

[TABLE]

Similarly, given $\pi,\pi^{\prime}$, the squared MMD \mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\kappa}^{2} with this kernel is the expectation of

[TABLE]

A characterization of the MMD that we will leverage throughout this section is that for any $\pi,\pi^{\prime}$,

[TABLE]

We observe that $\mathcal{A}$ satisfies the LRIP (20) (with $\eta=0$) for a given model set $\mathfrak{S}$ if, and only if, the metric \mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\Delta\mathcal{L}} is dominated by \mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\kappa_{\mathrm{\Phi}}} for $\pi,\pi^{\prime}\in\mathfrak{S}$. Our overall strategy to check that a random feature function $\mathrm{\Phi}$ defined by $\mathcal{F}$ and $\Lambda$ satisfies the LRIP (20) (with $\eta=0$) with controlled sketch dimension $m$ will be to:

prove that the average kernel $\kappa$ defined by (46) satisfies the Kernel LRIP

[TABLE] 2. 2.

in the spirit of compressive sensing theory, use concentration of measure and covering arguments to show that for any $0<\delta<1$, for large enough $m$, with high probability on the draw of $\omega_{j}$,

[TABLE]

so that the kernel LRIP (47) actually holds with $\kappa_{\mathrm{\Phi}}$ instead of $\kappa$ and constant $C_{\kappa_{\mathrm{\Phi}}}:=C_{\kappa}/\sqrt{1-\delta}$.

Remark 5.3**.**

The expression (48) expresses the control of the relative error of approximation of the MMD, restricted to certain distributions. This contrasts with state of the art results on random features that either control uniformly the approximation of the kernel \mathopen{}\mathclose{{}\left\lvert\kappa_{\mathrm{\Phi}}(\cdot,\cdot)-\kappa(\cdot,\cdot)}\right\rvert (see e.g. Sriperumbudur and Szabó, 2015) or the approximation of functions in the RKHS induced by $\kappa$ by those in the RKHS induced by $\kappa_{\mathrm{\Phi}}$ (Bach, 2017). These types of controls are indeed of a different nature compared to ours, and none seems to be a direct consequence of the others.

5.2 Ingredients to verify the Lower Restricted Isometry Property

In sight of the inequalities (47),(48) we need to prove, the analysis will focus on the so-called normalized secant set of the model $\mathfrak{S}$ with respect to the average kernel $\kappa$, defined as follows (see, e.g. Dirksen, 2016; Puy et al., 2017):

Definition 5.4** (Normalized secant set).**

The normalized secant set of a model set $\mathfrak{S}$ with respect to a kernel $\kappa$ is the following subset of the set of finite, signed measures (see Appendix A.2) with appropriate integrability

[TABLE]

Using the secant set, the LRIP (48) is equivalent to

[TABLE]

The radius of $\mathcal{S}_{\kappa}$ with respect to certain function norms will play an important role. Since this notion will come up repeatedly in the analysis, we introduce the following notation, which will be heavily used in the sequel. Given a norm \mathopen{}\mathclose{{}\left\|\cdot}\right\| on measures, the radius of a subset $\mathcal{E}$ of finite signed measures is denoted

[TABLE]

With these definitions we can observe that

[TABLE]

where we recall that the metric \mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\mathcal{G}} is defined in (8). In particular, the constant from (47) can be equivalently rewritten as C_{\kappa}:=\mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\Delta\mathcal{L}}.

Concerning 48, the strategy to establish it will rely on the following two quantities: first, a concentration function $c_{\kappa}(t)$ characterizing the pointwise (i.e. for fixed ${\tau},{\tau}^{\prime}\in\mathfrak{S}$) concentration of \mathopen{}\mathclose{{}\left\|\mathcal{A}({\tau})-\mathcal{A}({\tau}^{\prime})}\right\|_{2}^{2} around its expectation; secondly, certain covering numbers of $\mathcal{S}_{\kappa}$ needed to step from pointwise to uniform concentration.

Classical arguments from compressive sensing (Baraniuk et al., 2008; Eftekhari and Wakin, 2015; Puy et al., 2017; Dirksen, 2016; Foucart and Rauhut, 2012) prove that certain random linear operators satisfy the RIP by relying on pointwise concentration inequalities. Similarly, a first step to establish that the inequalities (48) hold with high probability consists in assuming first a pointwise version of the same, i.e., for any choice of $m$ in (43):

[TABLE]

for some concentration function $t\mapsto c_{\kappa}(t)$ that should ideally be as small as possible. The following result shows that the radius \mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\mathcal{F}} can be used to control such a concentration function.

Lemma 5.5**.**

Consider a family of functions $\mathcal{F}:=\{{\phi_{\omega}}\}_{\omega\in\Omega}$, $m$ parameters $(\omega_{j})_{j=1}^{m}$ drawn i.i.d. according to some distribution $\Lambda$ on $\Omega$, and $\mathcal{A}$ the (random) operator induced (see (5)) by the feature function \mathrm{\Phi}(x):=\tfrac{1}{\sqrt{m}}\mathopen{}\mathclose{{}\left({\phi_{\omega_{j}}}(x)}\right)_{j=1}^{m}. Denoting $\kappa$ the associated average kernel (cf (46)) we have \mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\mathcal{F}}\geq 1. Moreover, if \mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\mathcal{F}}<\infty then (52) holds for all $\mu\in\mathcal{S}_{\kappa}$, with

[TABLE]

The proof is in Appendix C. Observe that the above estimate only depends on the choice of the feature family $\mathcal{F}$, and holds for any feature sampling distribution $\Lambda$. More refined estimates for mixture models, exploiting moments of $\Lambda$ rather than a uniform bound, are provided in the companion paper (Gribonval et al., 2020) and used to obtain concrete estimates for Compressive Clustering and Compressive GMM.

Finally, we will extrapolate pointwise concentration (52) to all pairs ${\tau},{\tau}^{\prime}\in\mathfrak{S}$ using covering numbers of the normalized secant set with respect to an appropriate metric.

Definition 5.6** (Covering number).**

The covering number \mathrm{N}\mathopen{}\mathclose{{}\left(d(\cdot,\cdot),S,\delta}\right) of a set $S$ with respect to a (pseudo)metric555Further reminders on metrics, pseudometrics, and covering numbers are given in Appendix A. $d(\cdot,\cdot)$ is the minimum number of closed balls of radius $\delta$ with respect to $d(\cdot,\cdot)$ with centers in $S$ needed to cover $S$.

As the normalized secant set is a subset of the infinite-dimensional space of finite-signed measures, it is not obvious when its covering numbers are finite. Controlling them can be nontrivial, yet this is feasible on a case by case basis as will be illustrated in the companion paper (Gribonval et al., 2020).

Covering numbers and pointwise concentration can be then combined to give rise to the following result (whose proof is in Appendix C) where the logarithm of the covering numbers somehow captures an intrinsic dimension of the considered learning task:

Theorem 5.7**.**

Consider \mathcal{F}:=\mathopen{}\mathclose{{}\left\{{\phi_{\omega}}}\right\}_{\omega\in\Omega} a family of functions, $\Lambda$ a probability distribution on $\Omega$, $\mathrm{\Phi}$ the associated random feature function and $\kappa$ the corresponding average kernel. Consider the pseudometric on $\mathcal{F}$-integrable probability distributions666In fact, we consider the extension of $d_{\mathrm{\Phi}}$ to finite, $\mathcal{F}$-integrable signed measures, see Appendix A.2.

[TABLE]

Consider a model set $\mathfrak{S}$ and $\mathcal{S}_{\kappa}=\mathcal{S}_{\kappa}(\mathfrak{S})$ its normalized secant set. Assume that $\mathcal{S}_{\kappa}$ has finite covering numbers with respect to the pseudometric $d_{\mathcal{F}}$. For $0<\delta,\zeta<1$, if

[TABLE]

then, with probability at least $1-\zeta$ on the draw of $(\omega_{j})_{j=1}^{m}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\Lambda$, the operator $\mathcal{A}$ induced by $\mathrm{\Phi}$ (cf (43) and (5)) satisfies

[TABLE]

When (56) holds, the LRIP (20) with $\eta=0$ holds with constant C_{\mathcal{A}}:=\frac{\mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\Delta\mathcal{L}}}{\sqrt{1-\delta}} and $\eta=0$.

5.3 Summary and applications

To briefly summarize the results in this section, in order to establish the LRIP property with respect to a given model $\mathfrak{S}$ in the context of a sketching operator $\mathrm{\Phi}$ associated to a family of random features $\mathcal{F}$ and feature sampling distribution $\Lambda$ we proceed as follows. After identifying the associated average kernel (44), the key quantities to estimate relative to the normalized secant $\mathcal{S}_{\kappa}(\mathfrak{S})$ are its radius \mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\Delta\mathcal{L}} (which serves as a measure of compatibility between the kernel, the learning task, and the model set ${\mathfrak{S}}$), the pointwise concentration function $c_{\kappa}(.)$ from (52), and the covering numbers of $\mathcal{S}_{\kappa}$ with respect to the distance $d_{\mathcal{F}}$ from (54).

Even though the above ingredients and results may look quite abstract at this stage, we can turn them into concrete estimates on several examples. The resulting guarantees are summarized in Table 1 for the examples developed in detail in the companion paper (Gribonval et al., 2020)(compressive $k$-Means, $k$-medians and GMM).

The random sketching results developed in the present section can also be used to revisit the illustrative PCA example from Section 4. Namely, while we have directly lifted from existing literature the RIP property (37) for random Gaussian sketching matrices applied to low-rank covariance matrices, the arguments used there to establish this property follow in essence the canvas of this section (pointwise concentration of the random operator to its average, then unifom concentration via appropriate covering number arguments). In this context, the squared MMD with respect to the averaged kernel is precisely the Frobenius norm between covariance matrices. The additional ingredient needed to complete the analysis is to relate the PCA excess risk to the Frobenius norm of differences of low rank matrices (see (72) in the technical Appendix E, which can be reinterpreted as a bound on \mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\Delta\mathcal{L}} in the PCA setting).

6 Conclusion and perspectives

The principle of compressive statistical learning is to learn from large-scale collections by first summarizing the collection into a sketch vector made of empirical (random) moments, before solving a nonlinear least squares problem. The main contribution of this paper is to set up a general mathematical framework for compressive statistical learning and to demonstrate on an example (compressive PCA) that the excess risk of this procedure can be controlled, as well as the sketch size. The companion paper (Gribonval et al., 2020) completes the illustration of the framework by considering two more examples: compressive clustering and compressive Gaussian mixture estimation — with fixed known covariance.

Sharpened estimates?

Our demonstration of the validity of the compressive statistical learning framework for certain tasks is, in a sense, qualitative, and we expect that many bounds and constants are sub-optimal. A number of non-sharp oracle inequalities have been established in the course of our endeavor. A particular question is to obtain more explicit and/or tighter control of the bias term $d_{h^{\star}_{\pi}}(\pi,\mathfrak{S}^{\mathtt{CT}}(\mathcal{H}))$, and to understand whether Lemma 3.2, which relates this bias term to the optimal risk, can be tightened and/or extended to other loss functions. In the same vein, as fast convergence rates for the excess risk can be established for certain classical statistical learning tasks under appropriate conditions (see e.g. (Levrard, 2013) for the case of $k$-means), it is natural to wonder whether the same holds for compressive statistical learning.

Links with neural networks.

From an algorithmic perspective, the sketching techniques we have explicitly characterized in this paper have a particular structure which is reminiscent of a one-layer (random) neural network with subsequent averaging over multiple data points. Indeed, when the sketching function $\mathrm{\Phi}$ corresponds to random Fourier features, its computation for a given vector $x$ involves first multiplication by the matrix $\mathbf{W}\in\mathbb{R}^{m\times d}$ whose rows are the selected frequencies $\omega_{j}\in\mathbb{R}^{d}$, then pointwise application of the $e^{\jmath\cdot}$ nonlinearity. Here we consider random Fourier moments, hence a subsequent averaging operation is performed. As we have seen, this draws a link with the MMD, as is done e.g. in the so-called MMD-GANS (Li et al., 2015; Binkowski et al., 2018), where the so-called discriminator is a neural net trained to compute MMDs over batches of samples.

This suggests that our analysis could help analyze the tradeoffs between reduction of the information flow (dimension reduction) across multiple layers of such networks and the preservation of statistical information (Shwartz-Ziv and Tishby, 2017). For example, this could explain why the pooled output of a layer is rich enough to cluster the input patches. Given the focus on drastic dimension reduction, this seems very complementary to the work on the invertibility of deep networks and pooling representations with random Gaussian weights (Bruna Estrach et al., 2014; Giryes et al., 2016; Gilbert et al., 2017). Finally, we mention the recent popularity of networks with random weights in statistical physics (Gabrié et al., 2018) and in analyzing the initialization point of optimization algorithms with a kernel characterization (Jacot et al., 2018; Bietti and Mairal, 2019), for which information-preservation (non-degeneracy during training) is also an essential feature.

Privacy-aware learning via sketching?

The reader may have noticed that, while we have defined sketching in (3) as the empirical average of (random) features $\mathrm{\Phi}(x_{i})$ over the training collection (or in fact the training stream), the essential feature of the sketching procedure is to provide a good empirical estimator of the sketch vector $\mathcal{A}(\pi)=\mathbb{E}_{X\sim\pi}\mathrm{\Phi}(X)$ of the underlying probability distribution. A consequence is that one can envision other sketching mechanisms, in particular ones more compatible with privacy-preservation constraints (Duchi et al., 2014). For example, one could average $\mathrm{\Phi}(x_{i}+\xi_{i})$, or $\mathrm{\Phi}(x_{i})+\xi_{i}$, or $\mathbf{D}_{i}\mathrm{\Phi}(x_{i})$, etc., where $\xi_{i}$ is a heavy-tailed random vector drawn independently from $x_{i}$, and $\mathbf{D}_{i}$ is a diagonal “masking” matrix with random Bernoulli \mathopen{}\mathclose{{}\left\{0,1}\right\} entries. An interesting perspective is to characterize such schemes in terms of tradeoffs between differential privacy and ability to learn from the resulting sketch. Preliminary results in this direction have been recently achieved Schellekens et al. (2019); Chatalic et al. (2021).

Recipes to design sketches for other learning tasks through kernel design?

Given the apparent genericity of the proposed compressive statistical learning framework, a particular challenge is to extend it beyond the learning tasks considered in this paper and its companion (Gribonval et al., 2020). Kernel versions of these tasks (kernel PCA, kernel $k$-means, or spectral clustering) appear as the most likely immediate extensions. They are expected to lead to sketching architectures reminiscent of two-layer convolutional neural networks with additive pooling. Compressive supervised classification and compressive regression are also natural candidate tasks in a learning setting, but seem more challenging.

Given a learning task, the main bottleneck is to find an adequate sketching function $\mathrm{\Phi}(\cdot)$. As illustrated on Figure 1, this primarily relies on the quest for a task-compatible kernel, i.e., one satisfying the Kernel LRIP (47). Subsequent technical steps would rely on the identification of an integral representation of this kernel using random features with the right concentration properties, and establishing that the associated secant set has finite covering dimension with respect to the feature-based metric (54). On a case by case basis, one may have to identify the analog of the separation conditions apparently needed for compressive $k$-means, see (Gribonval et al., 2020).

Vice-versa, one could wonder which family of learning tasks is compatible with a given kernel. In other words, how “universal” is a kernel, and how much can be learned from a single sketched representation of a database ? We expect that tasks such as compressive ranking, which involve pairs, triples, etc. of training samples, may require further extensions of the compressive statistical learning framework, to design sketches based on $U$-statistics rather than plain moments. These would lead to sketches linear in the product probability $\pi\otimes\pi$ instead of $\pi$. The investigation of such extended scenarios is expected to benefit from analogies with the lifting techniques used in phaseless reconstruction, see e.g. (Candès et al., 2013).

Acknowledgements

This work was supported in part by the European Research Council, PLEASE project (ERC-StG-2011-277906), the german DFG (FOR-1735 “Structural inference in statistics”, SFB-1294 “Data Assimilation”), the Franco-German University through the binational Doktorandenkolleg CDFA 01-18, and the ANR (ANR-19-CHIA-0021-01, project BISCOTTE; ANR-19-CHIA-0009, project AllegroAssai). Rémi Gribonval is very grateful to Michael E. Davies for many enlightening discussions around the idea of compressive statistical learning since this project started several years ago. The authors also wish to warmly thank Bernard Delyon and Adrien Saumard, as well as Gabriel Peyré and Lorenzo Rosasco for their constructive feedback on early versions of this manuscript.

Appendix

Appendix A Notations, definitions

In this section we group all notations and some useful classical results.

A.1 Metrics and covering numbers

Definition A.1**.**

A pseudometric $d$ over a set $X$ satisfies all the axioms of a metric, except that $d(x,y)=0$ does not necessarily imply $x=y$. Similarly, a semi-norm \mathopen{}\mathclose{{}\left\|\cdot}\right\| over a vector space $X$ satisfies the axioms of a norm except that \mathopen{}\mathclose{{}\left\|x}\right\|=0 does not necessarily imply $x=0$.

Definition A.2** (Ball, $\delta$-covering, Covering number).**

Let $(X,d)$ be a pseudometric space. For any $\delta>0$ and $x\in X$, we denote $\mathcal{B}_{X,d}(x,\delta)$ the closed ball of radius $\delta$ centered at the point $x$:

[TABLE]

Let $Y\subseteq X$ be a subset of $X$. A subset $Z\subseteq Y$ is a $\delta$-covering of $Y$ if $Y\subseteq\bigcup_{z\in Z}\mathcal{B}_{X,d}(z,\delta)$. The covering number \mathrm{N}\mathopen{}\mathclose{{}\left(d,Y,\delta}\right)\in\mathbb{N}\cup\{+\infty\} is the smallest $k$ such that there exists a $\delta$-covering of $Y$ made of $k$ elements $z_{i}\in Y$.

A.2 Finite signed measures

The space $\mathfrak{M}$ of finite signed measures on the measurable sample space $(\mathcal{Z},\mathfrak{Z})$ is a linear space that contains the set of probability distributions on $(\mathcal{Z},\mathfrak{Z})$. By the Hahn-Jordan theorem, any finite signed measure $\mu\in\mathfrak{M}$ can be decomposed into a positive and a negative part, $\mu=\mu_{+}-\mu_{-}$, where both $\mu_{+}$ and $\mu_{-}$ are non-negative finite measures on $(\mathcal{Z},\mathfrak{Z})$, hence $\mu_{+}=\alpha\pi_{+}$ and $\mu_{-}=\beta\pi_{-}$ for some probability distributions $\pi_{+},\pi_{-}$, and non-negative scalars $\alpha,\beta\geq 0$. A real-valued, measurable function $f$ on $(\mathcal{Z},\mathfrak{Z})$ is said integrable with respect to $\mu$ when it is integrable both with respect to $\mu_{+}$ and $\mu_{-}$. Noticing that the expectation of an integrable function $f$ is linear in the considered probability distribution, we adopt the inner product notation for expectations:

[TABLE]

being understood that we implicitly assume that $f$ is integrable with respect to $\pi$ when using this notation. This extends to finite signed measures: given a decomposition of $\mu\in\mathfrak{M}$ as $\mu=\alpha\pi-\beta\pi^{\prime}$ with $\pi,\pi^{\prime}$ two probability distributions and $\alpha,\beta\geq 0$, we denote

[TABLE]

which can be checked to be independent of the particular choice of decomposition of $\mu$. With these notations, given a class $\mathcal{G}$ of measurable functions $g:\mathcal{Z}\to\mathbb{R}\ \text{or}\ \mathbb{C}$ we can define

[TABLE]

and check that this is a semi-norm on the linear subspace \mathopen{}\mathclose{{}\left\{\mu\in\mathfrak{M}:\forall f\in\mathcal{G},f\text{ integrable w.r.t. }\mu}\right\} as claimed when we introduced (8). Similarly, pseudometrics similar to (54) can be extended to finite signed measures as

[TABLE]

When the functions in $\mathcal{G}$ are smooth these quantities can be extended to tempered distributions.

The total variation norm is defined on $\mathfrak{M}$ as \mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\textnormal{TV}}=\mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\mathcal{B}} with \mathcal{B}=\mathopen{}\mathclose{{}\left\{f:f\ \text{is continuous and}\ \mathopen{}\mathclose{{}\left\|f}\right\|_{\infty}\leq 1}\right\} (see e.g. Sriperumbudur et al., 2010) and yields a Banach structure on $\mathfrak{M}$ (see e.g. Halmos, 2013).

The mean kernel $\kappa$ (cf (44)) can naturally be extended from probability distributions to finite signed measures. Let $\mu_{1},\mu_{2}\in\mathfrak{M}$ and $\pi_{1},\pi_{1}^{\prime},\pi_{2},\pi_{2}^{\prime},\alpha_{1},\alpha_{2},\beta_{1},\beta_{2}$ such that $\mu_{1}=\alpha_{1}\pi_{1}-\beta_{1}\pi_{1}^{\prime}$ and $\mu_{2}=\alpha_{2}\pi_{2}-\beta_{2}\pi_{2}^{\prime}$ (decompositions as differences of probability measures). Provided that $\kappa(\cdot,\cdot)$ is well-defined on the corresponding probability distributions, we can define

[TABLE]

which can be checked to be independent of the particular choices of decomposition.

By linearity of the integral and the definition of the kernel for probability distributions, we obtain a seminorm \mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\kappa} associated to the mean kernel:

[TABLE]

that coincides with the metric of the mean kernel (45) for probability distributions.

Appendix B Proof of Theorem 2.5

In this section, we start with a suitable generalization of (Bourrier et al., 2014, Section IV-A), working with some relaxed assumptions on the considered metrics.

Definition B.1** (hemimetric).**

A function $d(\cdot\|\cdot):\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ is a hemimetric if

[TABLE]

A hemimetric is a pseudometric if it is symmetric: $d(x\|y)=d(y\|x)$ for any $x,y\in\mathcal{X}$.

Hemimetrics on basic sets such as $\mathcal{X}=\mathbb{R}^{n}$ will be denoted $d(x\|y)$. Hemimetrics between probability distributions will be denoted $D(\pi\|\pi^{\prime})$. The notation $d(x,y)$ will preferentially be used for pseudometrics on basic sets, while $D(\pi,\pi^{\prime})$ will denote pseudometrics between probability distributions.

Definition B.2** (relaxed lower restricted isometry property (rLRIP)).**

A function $\Psi:\mathcal{X}\to\mathcal{Y}$ satisfies the lower restricted isometry property on the subset $\Sigma\subset\mathcal{X}$ with respect to the hemimetric $d_{\mathcal{X}}(\cdot\|\cdot)$ on $\mathcal{X}$ and the pseudometric $d_{\mathcal{Y}}(\cdot,\cdot)$ on $\mathcal{Y}$ with constant $\eta\geq 0$ iff

[TABLE]

Lemma B.3**.**

Assume that $\Psi$ satisfies the rLRIP on $\Sigma$ with respect to $d_{\mathcal{X}}(\cdot\|\cdot)$ and $d_{\mathcal{Y}}(\cdot,\cdot)$ with constant $\eta$. Consider $\varepsilon,\nu\geq 0$ and a decoder $\Delta:\mathcal{Y}\to\Sigma\subset\mathcal{X}$ such that

[TABLE]

Then $\Delta$ satisfies the instance optimality property: $\forall x^{*}\in\mathcal{X},y\in\mathcal{Y}$,

[TABLE]

where

[TABLE]

Proof.

The proof follows very closely (Bourrier et al., 2014) and is adaptated to the fact that $d_{\mathcal{X}}(\cdot\|\cdot)$ is a hemimetric. Consider $x^{*}\in\mathcal{X}$, $y\in\mathcal{Y}$ and $\hat{x}=\Delta(y)$. Consider any $z\in\Sigma$ and write

[TABLE]

As this holds for any $z\in\Sigma$, taking the infimum yields the result. ∎

Remark B.4**.**

Conversely, when $d_{\mathcal{X}}(\cdot\|\cdot)$ is a pseudometric, if some decoder satisfies (61) for each $x^{*}\in\mathcal{X},y\in\mathcal{Y}$ with $\mathcal{D}^{*}(\cdot\|\cdot)$ some function such that $\mathcal{D}^{*}(x\|\Sigma)=0$ for each $x\in\Sigma$, then the rLRIP (59) holds with constant $2(\eta+\varepsilon)$ with respect to $\widetilde{d_{\mathcal{Y}}}(\cdot,\cdot)=(2+\nu)d_{\mathcal{Y}}(\cdot,\cdot)$. Indeed, for $x,x^{\prime}\in\Sigma$, as $\mathcal{D}^{*}(x\|\Sigma)=\mathcal{D}^{*}(x^{\prime}\|\Sigma)=0$, by (61) and the symmetry of $d_{\mathcal{X}}(\cdot\|\cdot)$, we have with $y:=\Psi(x)$, $\hat{x}:=\Delta(y)$:

[TABLE]

The triangle inequality yields $d_{\mathcal{X}}(x\|x^{\prime})\leq d_{\mathcal{X}}(x\|\hat{x})+d_{\mathcal{X}}(\hat{x}\|x^{\prime})\leq\widetilde{d_{\mathcal{Y}}}(\Psi(x),\Psi(x^{\prime}))+2(\eta+\varepsilon)$.

To prove Theorem 2.5, given a fixed $h_{0}\in\mathcal{H}$, we apply Lemma B.3 with

[TABLE]

By (19)-(20) the rLRIP (59) is satisfied by $\Psi$ on $\Sigma:=\mathfrak{S}$. By (22) $y:=\mathbf{y}=\mathcal{A}(\hat{\pi}_{n})$ and $\Delta(\mathbf{y})=\widetilde{\pi}$ satisfy (60) with $\varepsilon$ replaced by $C_{\mathcal{A}}\varepsilon$. Since the definition (62) with yields $\mathcal{D}^{*}(\pi\|\Sigma)=d_{h_{0}}(\pi\|\mathfrak{S})$ as in (25), we get by (61) with $x^{*}:=\pi$:

[TABLE]

On the other hand, using (23), we obtain

[TABLE]

and thus, by definition of $D_{h_{0}}(\pi\|\widetilde{\pi})$ (cf (18)):

[TABLE]

Appendix C Proof of Lemma 5.5 and

Theorem 5.7

To establish Lemma 5.5 we use Bernstein’s inequality for bounded random variables, which is for example a consequence of Massart 2007, Corollary 2.10.

Lemma C.1** (Bernstein’s inequality).**

Let $X_{i}\in\mathbb{R}$, $i=1,\ldots,N$ be $i.i.d.$ bounded random variables such that $\mathbb{E}X_{i}=0$, $|X_{i}|\leq M$ and $Var(X_{i})\leq\sigma^{2}$ for all $i$’s. Then for all $t>0$ we have

[TABLE]

Proof of Lemma 5.5.

First, observe that for any $\mathcal{F}$-integrable probability distributions $\pi,\pi^{\prime}$

[TABLE]

and that

[TABLE]

Specializing to ${\tau},{\tau}^{\prime}\in\mathfrak{S}$ we get 1\leq C:=\mathopen{}\mathclose{{}\left\|{\tau}-{\tau}^{\prime}}\right\|_{\mathcal{\mathcal{F}}}/\mathopen{}\mathclose{{}\left\|{\tau}-{\tau}^{\prime}}\right\|_{\kappa}\leq\mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\mathcal{F}} and $-1\leq Z(\omega)\leq C^{2}-1$, hence \mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\mathcal{F}}\geq 1 and |Z(\omega)|\leq\max\mathopen{}\mathclose{{}\left(1,C^{2}-1}\right)\leq C^{2}\leq\mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\mathcal{F}}^{2}. As $\mathbb{E}_{\omega\sim\Lambda}Z(\omega)=0$ we obtain

[TABLE]

Applying Lemma C.1 with the independent random variables $Z(\omega)$ we obtain for each $t>0$:

[TABLE]

Lemma C.2**.**

Consider a family of functions $\mathcal{F}:=\{{\phi_{\omega}}\}_{\omega\in\Omega}$, $m$ parameters $(\omega_{j})_{j=1}^{m}$ drawn i.i.d. according to some distribution $\Lambda$ on $\Omega$, and $\mathcal{A}$ the (random) operator induced (see (5)) by the feature function \mathrm{\Phi}(x):=\tfrac{1}{\sqrt{m}}\mathopen{}\mathclose{{}\left({\phi_{\omega_{j}}}(x)}\right)_{j=1}^{m}. and $\kappa$ the associated kernel. Assume that (52) holds with concentration function $c_{\kappa}(t)$ and consider $\mathcal{S}\subset\mathcal{S}_{\kappa}$ and $d_{\mathcal{F}}$ the metric defined in (54). For any $\delta>0$ such that

[TABLE]

we have, with probability at least $1-2N\exp(-m/c_{\kappa}(\delta/2))$:

[TABLE]

Proof of Lemma C.2.

Consider \mu=({\tau}-{\tau}^{\prime})/\mathopen{}\mathclose{{}\left\|{\tau}-{\tau}^{\prime}}\right\|_{\kappa} with ${\tau},{\tau}^{\prime}\in\mathfrak{S}$. By definition of the concentration function, for any $t>0$ and $m\geq 1$

[TABLE]

This establishes a pointwise concentration result when $\mu$ is on the normalized secant set $\mathcal{S}_{\kappa}$. We now use a standard argument to extend this to a uniform result on $\mathcal{S}$. Let $\mu_{i}$, $1\leq i\leq N$ be the centers of a $\delta/2$-covering (with respect to the metric $d_{\mathcal{F}}$) of $\mathcal{S}$. Using (67) with $t=\delta/2$, the probability that there is an index $i$ such that \Big{|}\mathopen{}\mathclose{{}\left\|\mathcal{A}(\mu_{i})}\right\|_{2}^{2}-1\Big{|}\geq\delta/2 is at most $\zeta=2N\exp(-m/c_{\kappa}(\delta/2))$. Hence, with probability at least $1-\zeta$, we have: for any $\mu\in\mathcal{S}$, with $i$ an index chosen so that $d_{\mathcal{F}}(\mu,\mu_{i})\leq\delta/2$:

[TABLE]

Proof of Theorem 5.7.

Denote $\zeta=2N\exp(-m/c_{\kappa}(\delta/2))$ with N:=\mathrm{N}\mathopen{}\mathclose{{}\left(d_{\mathcal{F}},\mathcal{S}_{\kappa},\delta/2}\right). By Lemma C.2, the assumptions imply that with probability at least $1-\zeta$ on the draw of $\omega_{j}$, $1\leq j\leq m$, we have

[TABLE]

This implies (48). Since \mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\Delta\mathcal{L}}<\infty, the LRIP (20) holds wrt \mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\Delta\mathcal{L}} with C_{\mathcal{A}}=\frac{\mathopen{}\mathclose{{}\left\|\mathcal{S}_{\kappa}}\right\|_{\Delta\mathcal{L}}}{\sqrt{1-\delta}}. ∎

Appendix D Proof of Lemma 3.2 and Lemma 3.4

If $\pi\in\mathfrak{S}^{\mathtt{CT}}_{h}$ then $0=\mathcal{R}(\pi,h)=\mathbb{E}_{X\sim\pi}\ell(X,h)=\mathbb{E}_{X\sim\pi}d^{p}(X,P_{h}X)$ hence $d(X,P_{h}X)=0$ almost surely, i.e., $X=P_{h}X\in P_{h}\mathcal{Z}=\mathcal{E}_{h}$ almost surely. The converse is trivial. The bound (33) follows directly since for any $h\in\mathcal{H}$, $P_{h}\pi\in\mathfrak{S}^{\mathtt{CT}}(h)\subset\mathfrak{S}^{\mathtt{CT}}(\mathcal{H})$. This establishes the first claim of Lemma 3.2.

Let $h_{0}\in\mathcal{H}$ be fixed. By (29), with $Y\sim P_{h_{0}}\pi$, we have $P_{h_{0}}Y=Y$ hence $\ell(Y,h_{0})=d^{p}(Y,P_{h_{0}}Y)=0$ and for any $h\in\mathcal{H}$:

[TABLE]

For the second claim of Lemma 3.2, by (30), since $d^{p}$ is a metric we have for any $x\in\mathcal{Z}$

[TABLE]

It follows using (68) that

[TABLE]

As this holds for any $h$, and as equality is reached for $h=h_{0}$, we get $D_{h_{0}}(\pi\|P_{h_{0}}\pi)=0$.

In particular when $p\in(0,1]$ we have $(a+b)^{p}\leq a^{p}+b^{p}$ for any $a,b\geq 0$ hence for $u,v,w\in\mathcal{Z}$, by the triangle inequality, $d^{p}(u,v)\leq[d(u,w)+d(w,v)]^{p}\leq d^{p}(u,w)+d^{p}(w,v)$, showing that $d^{p}$ is a metric.

For the claims of Lemma 3.4, we will exploit optimal transport through connections between the considered norms and the norm \mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\mathrm{Lip}(L,d)}=L\cdot\mathopen{}\mathclose{{}\left\|\cdot}\right\|_{\mathrm{Lip}(1,d)}, where $\mathrm{Lip}(L,d)$ denotes the class of functions $f:(\mathcal{Z},d)\to\mathbb{R}$ that are $L$-Lipschitz.

For $p\geq 1$, and $\mathcal{Z}$ with $d$-diameter bounded by $B$, since \mathopen{}\mathclose{{}\left\lvert a^{p}-b^{p}}\right\rvert\leq\max(pa^{p-1},pb^{p-1})\mathopen{}\mathclose{{}\left\lvert a-b}\right\rvert for any $a,b\geq 0$, we have

[TABLE]

by the triangle inequality, hence $\mathcal{L}(\mathcal{H})\subset\mathrm{Lip}(pB^{p-1},d)$. Using (19) this implies that for any $\pi,\pi^{\prime}$ we have in general in the above considered setting:

[TABLE]

It is well-known that the 1-Wasserstein distance between two distributions can be equivalently defined in terms of optimal transport (so-called “earth mover’s distance”) but also as

[TABLE]

as soon as $(\mathcal{Z},d)$ is a separable metric space, see, e.g., (Dudley, 2002, Theorem 11.8.2). By the transport characterization of the Wasserstein distance, considering the transport plan that sends $x$ to $P_{h}x$, where $h\in\mathcal{H}^{\prime}$, we conclude

[TABLE]

by Jensen’s inequality (since $p\geq 1$ here), yielding the claim (34).

For the final claim, we have

[TABLE]

where f_{\mathbf{u}}(x):=\mathopen{}\mathclose{{}\left\langle\mathrm{\Phi}(x),\mathbf{u}}\right\rangle. Moreover, for \mathopen{}\mathclose{{}\left\|\mathbf{u}}\right\|_{2}\leq 1 and any $x,x^{\prime}$, since $\mathrm{\Phi}$ is assumed $L$-Lipschitz:

[TABLE]

i.e., $f_{\mathbf{u}}(\cdot)$ is $L$-Lipschitz with respect to $d(\cdot,\cdot)$. It follows that for any $\pi,\pi^{\prime}$, \mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\pi^{\prime})}\right\|_{2}\leq L\mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\mathrm{Lip}(1,d)}=L\mathopen{}\mathclose{{}\left\|\pi-\pi^{\prime}}\right\|_{\textnormal{Wasserstein}_{1}(d)}.

The claim (35) when $p\geq 1$ follows by (69). When $p\leq 1$ and the space $\mathcal{Z}$ has $d$-diameter bounded by $B$, as $d(X,P_{h}X)=d^{1-p}(X,P_{h}X)d^{p}(X,P_{h}X)\leq B^{1-p}d^{p}(X,P_{h}X)$, we obtain (36) as follows

[TABLE]

Appendix E Proof of Theorem 4.1 on Compressive PCA

For Compressive PCA, we recall that $k$ is the number of PCA components we want to estimate. The hypothesis class $\mathcal{H}$ is the set of linear subspaces of dimension $k$ of the input space $\mathbb{R}^{d}$, which is in one-to-one correspondance with the space $\mathfrak{P}_{k}$ of orthoprojectors $\mathbf{P}$ of rank $k$. In the remainder of this section we therefore use directly $\mathfrak{P}_{k}$ as the hypothesis class, for notational convenience. We recall that for $r\geq k$, we consider the model $\mathfrak{S}_{r}$ consisting of probability distributions having their second moment matrix of rank at most $r$.

Observe that for any $\mathbf{P}\in\mathfrak{P}_{k}$:

[TABLE]

where the inner product \mathopen{}\mathclose{{}\left\langle,}\right\rangle_{F} is the Frobenius product, and the minimum risk is

[TABLE]

where $\mathbf{P}^{*[\ell]}_{\pi}$ ($1\leq\ell\leq d$) denotes an orthoprojector onto the $\ell$ first eigenvectors of $\mathbf{\Sigma}_{\pi}$, and $\lambda_{i}(\mathbf{M})$ denote the eigenvalues, with multiplicity and ordered in nonincreasing sequence, of a matrix $\mathbf{M}$.

We follow the improved risk analysis of Section 2.5. In the above setting, the excess risk divergence (18) with respect to an arbitrary reference hypothesis $\mathbf{P}_{0}\in\mathfrak{P}_{k}$ is given by

[TABLE]

Lower RIP.

We start with the following bound on the excess risk divergence, holding without restriction for any distributions $\pi,\pi^{\prime}$ with existing second moments and any $\mathbf{P}_{0}\in\mathfrak{P}_{k}$:

[TABLE]

Proof of (72).

By the so-called Ky Fan Theorem Fan (1949), for a symmetric matrix $\mathbf{M}\in\mathbb{R}^{d\times d}$, and a positive integer $\ell\leq d$, one has

[TABLE]

As a result, we obtain

[TABLE]

∎

Since the right-hand side of (72) does not depend of $\mathbf{P}_{0}$, we get from (19) in particular that for any $\pi,\pi^{\prime}$ with finite second moments

[TABLE]

Hence, since $\mathcal{M}$ is a linear operator having a RIP (37) on matrices of rank lower than $2r$, $\mathcal{M}$ induces (in the way described in Section 4) a sketching operator $\mathcal{A}:\pi\mapsto\mathcal{A}(\pi):=\mathcal{M}(\mathbf{\Sigma}_{\pi})$ that has the lower RIP described in (20) with constants $C_{\mathcal{A}}=\frac{\sqrt{2k}}{\sqrt{1-\delta}},\eta=0$ on model $\mathfrak{S}_{r}$.

Ideal decoder and generic excess risk control.

The ideal decoder (12) writes

[TABLE]

and is equivalent to (38), i.e.

[TABLE]

Formally, $\Delta[\mathbf{y}]$ can then be taken as any distribution having second moment matrix $\hat{\mathbf{\Sigma}}$. This last step can naturally be shunted since whatever the choice of such a representative distribution, the associated estimated hypothesis $\hat{h}$ is directly given by (39).

By Theorem 2.5, this decoder is instance optimal yielding (24) with $\nu=\eta=\varepsilon=\varepsilon^{\prime}=0$, i.e., the excess risk of $\hat{h}$ of the procedure of Section 4 when the true data distribution is $\pi$ is controlled by

[TABLE]

with the bias term defined in (25).

Control of the bias term.

The next lemma improves over Lemma 3.2 in the special case of PCA.

Lemma E.1**.**

Consider a probability distribution $\pi$ with finite second moments, $\mathbf{P}^{*}:=\mathbf{P}^{*[k]}_{\pi}$, and ${\tau}_{r}\in\mathfrak{S}_{r}$ (any) probability distribution with covariance $\mathbf{\Sigma}_{\pi}^{[r]}:=\mathbf{P}^{*[r]}\mathbf{\Sigma}_{\pi}\mathbf{P}^{*[r]}$. We have $D_{\mathbf{P}^{*}}(\pi\|{\tau}_{r})=0$.

Proof.

Since $\mathbf{\Sigma}_{\pi}-\mathbf{\Sigma}_{\pi}^{[r]}$ is a nonnegative matrix and $r\geq k$, (71) yields

[TABLE]

As a result the “bias” term in the bound (73) is upper bounded as

[TABLE]

From now on we assume $d>k$ (otherwise the bias is trivially 0). We now use the following lemma to bound \mathopen{}\mathclose{{}\left\|\mathcal{M}(\mathbf{\Sigma}_{\pi}-\mathbf{\Sigma}_{\pi}^{[r]})}\right\|_{2}.

Lemma E.2**.**

Let $\mathbf{\Sigma}\in\mathbb{R}^{d\times d}$ be symmetric p.s.d. and $\mathcal{M}:\mathbb{R}^{d\times d}\to\mathbb{R}^{m}$ satisfy the upper RIP in (37), i.e. {\mathopen{}\mathclose{{}\left\|\mathcal{M}(\mathbf{M})}\right\|_{2}^{2}}/{\mathopen{}\mathclose{{}\left\|\mathbf{M}}\right\|_{F}^{2}}\leq 1+\delta for all matrices of rank at most $2r$. For $1\leq\ell\leq d$, consider $\mathbf{\Sigma}^{[\ell]}$ a best rank $\ell$ approximation to $\mathbf{\Sigma}$. Then for any $1\leq s\leq\min(\ell,2r)$ we have

[TABLE]

where $\sigma_{j}:=\sigma_{j}(\mathbf{\Sigma}):=\sum_{i=j+1}^{d}\lambda_{i}(\mathbf{\Sigma})$.

Proof.

This proof follows mainly the ideas of Candès (2008). By definition of $\mathbf{\Sigma}^{[\ell]}$ there is an eigendecomposition $\mathbf{\Sigma}=U\Lambda U^{T}$, where $\Lambda$ is a diagonal matrix containing the eigenvalues of $\mathbf{\Sigma}$ with multiplicity in decreasing order, such that $\mathbf{\Sigma}^{[\ell]}=U\Lambda^{[\ell]}U^{T}$ where $\Lambda^{[\ell]}$ contains the first $\ell$ eigenvalues. Decompose $\Lambda$ into blocks, $\Lambda=\sum_{j\geq 0}\Lambda_{j}$, where $\Lambda_{0}$ contains the first $\ell-s$ eigenvalues, and $\Lambda_{j}$, $j\geq 1$ are the next blocks of $s$ eigenvalues in decreasing order (the last block is of size $\leq s$); that is, for $j\geq 1$ the block $\Lambda_{j}$ contains eigenvalues of indices $m$ such that $\ell+(j-2)s<m\leq\ell+(j-1)s$. Let us also denote $\Lambda_{j}^{+}$, for $j\geq 1$, the blocks of eigenvalues of size $s$ starting one index later, that is, $\Lambda_{j}$ contains eigenvalues of indices $m$ such that $\ell+(j-2)s+1<m\leq\ell+(j-1)s+1$. Let $S_{j}=U\Lambda_{j}U^{T}$, so that $\mathbf{\Sigma}-\mathbf{\Sigma}^{[\ell]}=\sum_{j\geq 2}S_{j}$. As $s\leq 2r$, $\mathrm{rank}(S_{j})\leq 2r$, $j\geq 2$ hence by the upper RIP property,

[TABLE]

∎

The above lemma (with $\ell:=r$) allows us to control \mathopen{}\mathclose{{}\left\|\mathcal{A}(\pi)-\mathcal{A}(\pi_{r})}\right\|=\mathopen{}\mathclose{{}\left\|\mathcal{M}(\mathbf{\Sigma}_{\pi}-\mathbf{\Sigma}_{\pi}^{[r]})}\right\|_{2} using $\sigma_{r-s+1}$ for $1\leq s\leq r$. This, in combination with (74) and (73), establishes (40). Moreover, choosing $s:=r-k+1$, we get $\sigma_{r-s+1}=\sigma_{k}=\mathcal{R}_{k-\mathtt{PCA}}(\pi,\mathbf{P}^{*[k]}_{\pi})$, leading to (41).This proves Theorem 4.1.

Table of notations