Sketch and shift: a robust decoder for compressive clustering
Authors: Ayoub Belhadji
Rémi Gribonval
Venue: Accepted by TMLR ()

Abstract: Compressive learning is an emerging approach to drastically reduce the memory footprint of large-scale learning, by first summarizing a large dataset into a low-dimensional sketch vector, and then decoding from this sketch the latent information needed for learning. In light of recent progress on information preservation guarantees for sketches based on random features, a major objective is to design easy-to-tune algorithms (called decoders) to robustly and efficiently extract this information.  To address the underlying non-convex optimization problems, various heuristics have been proposed. In the case of compressive clustering, the standard heuristic is CL-OMPR, a variant of sliding Frank-Wolfe. Yet, CL-OMPR is hard to tune, and the examination of its robustness was overlooked. 
In this work, we undertake a scrutinized examination of CL-OMPR to circumvent its limitations. In particular, we show how this algorithm can fail to recover the clusters even in advantageous scenarios. To gain insight, we show how the deficiencies of this algorithm can be attributed to optimization difficulties related to the structure of a correlation function appearing at core steps of the algorithm. To address these limitations, we propose an alternative decoder offering substantial improvements over CL-OMPR. Its design is notably inspired from the mean shift algorithm, a classic approach to detect the local maxima of kernel density estimators. The proposed algorithm can extract clustering information from a sketch of the MNIST dataset that is 10 times smaller than previously.

Published in Transactions on Machine Learning Research (04/2024)
Sketch and shift: a robust decoder for compressive clustering
Ayoub Bclhadji
ayoob.bbllhadjji@gnaaaa
Univ Lyon, ENS de Lyon, Inria, CNRS, UCBL, LIP UMR 5668, Lyon, France
Rémi Gribonval
remi.grioonval@Bariaa.
Unio Lyon, ENS de Lyon, Inria, CNRS, UCBL, LIP UMR 5668, Lyon, France
Reviewed on OpenReview: https: // openreview. net/ forum? id= G3WWWW
Abstract
Compressive learning is emerging approach to drastically reduce the memory footprint of large-scale learning, by first summariziin large dataset into low-dimensional sketch vector, and then decoding from this sketch the latent information needed for learning. In light of recent progress on information preservation guarantees for sketches based on random features, a major objective is to design easy-to-tune algorithms (called decoders) to robustly and efficiently extract this information. To address the undorlying non-convex optimization problems, various heuristics have been proposed. In the case of compressive clustering, the standard heuristic is CI-OMPR, variant of sliding Frank-Wolfe. Yet, CL-OMPR is hard to tune, and the examination of its robustness was overlooked. In this work, undertake scrutinized examination of CL-OMPR to circumvent its limitations. In particular, we show how this algorithm can fail to recover the clusters even in advantageous scenarios. To gain insight, we show how the deficiencies of this algorithm can be attributed to optimization difficulties related to the structure of a correlation function appearing at core steps of the algorithm. To address these limitations, we propose an alternative decoder offering sub- stantial improvements
CL-OMPR. Its design is notably inspired from the mean shift algorithm, classic approach to detect the local maxima of kernel density estimators. The proposed a algorithm extract clustering information from sketch of the MNIST (resp. of the CIFAR10) dataset that is 10 times smaller than previously and much easier to tune.
1 Introduction
In an era where resources are becoming scarcer, reducing the footprint of learning algorithms is of paramount importance. In this context, sketching is a promising paradigm. This consists conducting a learning task on low dimensional vector called a sketch that captures the essential structure of the initial datasct (Cormode let al. 2011). In other words, the sketch is an informative summary of the initial dataset which may be used to reduce the memory footprint of a learning task.
This article studies decoding in the context of compressive clustering. This is an emerging topic of rescarch in machine learning aiming to scale up the (unsupervised learning) task of clustering by conducting it on a sketch (Gribonval et al. 2021b). In this context, the sk is the mcan of a given feature map over the dataset (Keriven et al. 2017), and decoding consists in the recovery of the cluster centroids from such a sketch. So far, this learning step was conducted using heuristic decoder called CL-OMPR (Keriven et all, 2017). Despite existing proofs of concept for compressive learning with CL-OMPR, in-depth examination of the properties of this decoder has been lacking.
The starting point of our study is scrutinized examination of CL-OMPR. First, using numerical exper- iment, we show how CL-OMPR might fail to recover the clusters even in an advantageous scenario where the datasct generated through a Gaussian mixture model with well-soparated clusters. Second, to provide
The word "sketching" is highly overloaded machine learning. this paper refers vector computed as in equation)

---

Published in Transactions on Machine Learning Research (04/2024)
explanation of this deficiency, we analyze the first step of the algorithm which consists in finding new centroid that maximizes a correlation function. The latter is an approximation of the kernel density esti- mator (KDE) associated to shift-invariant kernel when the feature map corresponds to a Fourier feature map. We show that the main shortcomings of CL-OMPR stem from optimization difficulties related to the correlation function. Third, we propose an alternative decoder which takes inlo account these difficulties and validate on synthetic and real datasets. In particular, in the case of spectral features of MNIST, we demonstrate that our algorithm can extract clustering information from a sketch that is 10 times smaller.
This article is structured as follows. Section 2 recalls definitions and notions related to compressive clustering. In Section 3 we conduct a numerical simulation that highlights the limitations of CL-OMPR. In Section | we propose an alternative decoder that circumvents the shortcomings of CL-OMPR. Section 5 gathers numerical simulations that illustrate the improvement of the alternative decoder upon CL-OMPR.
2 Background
Clustering algorithms group elements into categories, also called clusters, based on their similarity. Numerous clustering algorithms have been proposed in the literature; see e.g. Jain (0110).
Beyond the classical Lloyd-Max approach (Looyd 1982), compressive clustering Keriven et al., 2017 was showed to offer an alternative both well-matched to distributed implementations and streaming scenarios and able to preserve privacy, see e.g. Gribonval et al. (2021b). The approach aims at drastically summarizing a (large) training collection while retaining the information needed to cluster. This objective is similar to that of coresets (see e.g. Feldman & Langberg (2011); Guo et al. (2022)), but the approach is radically different. Corescts summarize a collection of training samples by selecting a limited subset of representative samples (or sometines of transformed samples). In contrast, compressive clustering is somewhat more democratic: instead of selecting few "meritorious" samples, it computes unique vector, called a sketch (of dimension moderately higher than the dimension of cach training sample), that depends equally on all samples. The computation of this sketch is designed to capture the relevant information of the whole training collection to perform given clustering task.
For instance, in Keriven et al. (2018)), the authors compressod 1000 hours of speech data (50 gigabytes) into a sketch of a few kilobytes on a single laptop, using a random Fourier feature map. In this work, the sketch xx of dataset X (.1,...,,NN Rdd was taken to be
1 Σ (xx),
(1)
2x == N
where Φ: : .Rd + Cm is a given feature map. For instance, random Fourier features (RFF) are defined as ((x) = (dw,(()))en (x)))e{m) Cm", where
(x) Hlew
(2)
Vm
and wi,..., 4m are i.i.d. samples from o with o > 0 RRahii & Recht, 2007)
The sketch zx can be seen as the mean of the generalized moment defined by the feature map $ i.e., xx = Ez- where ny := N /N is the empirical distribution of the dataset X. In other words, xx AnN, where is the so-called sketching operator : P(R.) + on P(Rd), the set of probability distributions on Rd, defined by
Ar = / (())dT(x).
Jud
A theorctical analysis of compressive clustering using random Fourier features was conducted in Gribonval lt al.] (2021a); [elladji & Gribonval (2022). In particular, it was shown that the sketch size m must depend on the ambient dimension of the dataset and the number & of centroids in order to recover the centroids with high probability. More preciscly, sufficient sketch size was shown to scale as Okk2d), while numorical investigation in Keriven et al. (2017) showed that a sketch size m = O(kd) is enough in practice recover the cluster conts Alternatively, data-dependent Nyström feature maps were shown to require sketch
2

---

Published in Transactions on Machine Learning Research (04/2024)
size that depends the effective dimension of the dataset, which may be smaller than O (Chatalic
et al. 2022). In these works, centroid recovery was achieved by addressing the following inverse problem
min aERt Ei= a==1 1. E0
1x
Σ (c))l
(4)
where
θ C Rdd is a compact set. As a heuristic to address the non-convex moment-matching problem, Keriven
et al. 2017 proposed Compressive Learning-OMP (CL-OMPR), which is an adaptation of Orthonormal
Matching P Pursuil (OMP), a widely used algorithm in the field of sparse recovery (Mallat Zhang 1993 Pati
et al. 1993). In a nutshell, CL-OMPR minimizes the residual | -E (c)(c) by adding "atoms" x (<))
in a greedy fashion; see Algorithm 1 for a simplified version of CL-OMPR and Algorithm 3 in Appendix A for the details. The definition of CI-OMPR was indeed inspired by variant of OMP, namely OMP with replacement (OMPR), and adapted to the continuous setting with gradient descent steps. This is similar in spirit to Sliding Franke-Wold (Denoyelle al. 2020), continuous adaptation of Franke-Wole which is itself related to OMP, see e.g. Cherfaoui et al. (2019). Despite promising behaviour in several empirical proofs of concept, CL-OMPR is only a heuristic and our first main contribution in the next section is to highlight its weaknesses, before using the resulting diagnoses to propose a new decoder and demonstrating the resulting improved performance and robustness.
Data: Sketch zx, sketching operator A, parameters k,T 2 k, domain θ
Result: Set of centroids C, vector of weights a
2x; C 5 0;
for i i= . do
Step 1: Find a new centroid and expond the support C
fr(c) == A8.)
11(2,2) == =E 2z
Ve, EE=12
Cm
e+ arg max.ce (c(c)
C' Cuf
end
Step 2: Reduce the support by Hard Thresholding when i
Step 3: Project to find a a + arg minaao20 2| - Ele ac)(l
Step 4: Fine tuning by gradient descent steps: C, + arg mincco, 11ax - E.g a, (())
Step Update the residual: r +- 2x E (C))
end
Algorithm 1: CL-OMPR
3 Illustrating and diagnosing failures of CL-OMPR
Since its introduction by Keriven et al. (2017), CL-OMPR has been popular among methods that cluster using the sketch equation [ (Chatalic al., 2018 2022 In this section, we show how CL-OMPR is in fact not robust in the task of the recovery of centroids. To simplify the analysis, since the goal is show that failures can happen even in favorable scenarios, consider setting where: i) the data is drawn from mixture of well-separated Gaussians, ii) the dataset size N is very large, and iii) the sketch size is very large. We show through numerical experiment, that CL-OMPR fails to accurately reconstruct the centroids even in these optimal conditions.
3.1 A numerical experiment
We consider dataset X {II,..., EN) C R., where N 100000 and the x, i.id. draws from mixture of isotropic Gaussians Σk= a,, E,), where k 3, a 02 03 = 1/3, Cı 02, R and Σ Σx 22 E 22x2 as shown in Figure Ib We consider a sketching operator defined through random Fourier features associated to i.i.d. Gaussian Trequencies Wi,...,4M drawn from -T, with
3

---

Published in Transactions on Machine Learning Research (04/2024)
1.00
CL-MMP (m=:30)
CLIOMP (m=1000)
Proposed algortthm (m-30)
0.75
Proposed a algerithm (m=1000)
Lloyd's algorithmm
0.50
0.25
0.00
-0.2
-0.50
-0.75
0.01
0.01
0.05
0.1
0.2
-1.00
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
(a) MSE as a function of a
(b) The synthetic dataset
Figure 1: Average MSE of Lloyd's algorithm, CT-OMPR for two sketch sizes, and the proposed Algorithm 2 on synthetic dataset with three well-separated clusters in dimension 2. For well chosen values of the MSE is of the order of 0.01, to be compared with inter-cluster squared distances of the order of 0.25, and intra-cluster variances of the order of 0.05.
o > 0. Recall that, in the context k-means clustering, the mean squares error (MSE) of a set of centroids C {cI,.... ckl is defined by
1 MSE(C; x): N KEIM
Σ min I c/
(5)
This the quantity directly optimized by Lloyd's algorithm, while CL-OMPR instead addresses the opti- mization problem equation as proxy, both with constraint G; E e [-1, here. Figure la compares the MSE of Lloyd's algorithm, compressive clustering using CL-OMPR, and compressive clustering using Algorithm 2 for sketch sizes {30, 1000), averaged over 50 realizations of the sketching operator, function of the "bandwidth" (of random Fourier features). As in previous work, observe that the performance of CL-OMPR improves for large values of m. Moreover, we observe that, for performance CL-OMPR can be close to that of Lloyd-max, but the range of for which this is the case is narrow (a zs 0.3), indicating that the parameter a of the sketched clustering pipeline hard to tune. On the contrary, the performance of Algorithm 2 is the same as Lloyd's algorithm in a wider range of the bandwidth (a [0.03, 0.3]), even moderate sketch sizes m 30 for which CL-OMPR is unable to approach this performance for any value of
1000, the
3.2 Diagnoses via links with kernel density estimation
As we now show, compressive clustering bears strong links with kernel density estimation, and some of the difficulties in tuning observed above are reminiscent of classical kernel size tuning issues (Chen 2017). We also highlight more specific difficulties of CL-OMPR as decoder, which will help us design a better decoder (Algorithm in the next section.
To relate compressive clustering to kernel density estimation, we examine the correlation function, which appears in Step of CL-OMPR (Algorithm D.
Definition 1. Given a sketching operator P(TRd) - Cm", and given r e Cm define the correlation function fr Rd IR associated to and A
fx(x) : A6.).
(6)
Figure 2 shows the correlation function fax associated r 2x = Any, in the first step of CL-OMPR, with the sketching operator defined in Section 3.1 for {100, 1000), and (0.01, 0.05, 0.1, 0.2, 1) where
4

---

Published in Transactions on Machine Learning Research (04/2024)
TN -i=1dd /2,/N corresponds to the dataset represented by Figure 1b We observe that, for 1, the correlation function has single global maximum, which does not reflect the dataset's "clustered" geometry. As for a E (0.1,0.2), the correlation function exhibits multiple local maxima, the most significant ones being aligned with the cluster locations: the correlation function retains the datasct's geometry. In the case a 0.05, the correlation function continues depict the dataset's 's geometry while introducing a number of spurious local maxima, especially for small sketch size m 100. Finally, for a = 0.01, the correlation function heavily distorted by 'noise' especially for m 100. In light of these observations, we expect that recovering the clusters from the sketch will be infeasible for high and small values of o, whereas it might be possible but poses challenges from optimization perspective for intermediate d's.
(a) 0 0.01
(b) 0.05
(c) 0.1
(d) 0.2
(e) o=
(f) 0.01
(g) 0.05
(h) 0.1
(i) 0.2
(i) a=
Figure 2: The correlation function for 100 (top) versus 1000 (bottom).
An important fact is that the correlation function approximates the kernel density estimator
fKDE(0) = N Σ (z, 24),
(7)
where No (I Ba-Ni 41,) AL, A8y).
Indeed, 2x ATN, then
N
Eo (x) =
N Eu 21, Re(A,, = fKDE(*).
1
Σ
(8)
As consequence, step 1 in CL-OMPE (Algorithm I be seen search for local maxima of fKDE through its approximation fzx- In other words, compressive clustering is related to density-based clus- tering methods, to which it brings the ability to work in memtry-constrained scenarios, with distributed implementations, in streaming contexts, or with privacy constraints.
Classical difficulties of density-based clustering. Density-based clustering methods are attractive thanks their capacity to identify clusters of arbitrary shapes: they define clusters as high density regions in the feature space (Fukunaga & Hostetler, 1975). Thus, some of the difficulties that the practitioner of compressive clustering may face are inherent to density-based clustering methods such the selection of the kernel or the selection of the parameter of the kernel (Comaniciu & Meer 2002): a very large value of a yields a smeared kernel and thus a smooth KDE for which the whole dataset belongs to a single cluster, while a very low value of yields a trivial KDE where every point in the dataset belongs to its own cluster.
Specific difficultics of CL-OMPR. In addition, CL-OMPR (Algorithm 1 comes with own difficultics. First, Step performed by randomly initializing a candidate centroid and performing gradient ascent on fr. Such gradient. ascent fail for several reasons:
5

---

Published in Transactions on Machine Learning Research (04/2024)
If the data is well clustered, then gradients are vanishingly small beyond distance of roughly away from the true cluster centroids. Gradient ascent will essentially not move the randomly initialized centroid c. This is illustrated on Figure 3 in Section 4.1
For moderate sketch sizes m, spurious local maxima of Sax appear due to Gibbs-like phenomena: the approximation of the smooth FKDE of (6) by fax gives rise to oscillations, see Figure 2
Second, the residual is updated in Step 5 by removing from ax the sketch of an estimated [Cl-mixture of Diracs corresponding to the term Σ/e (, = AEE aqdeii see equation 3. In the most favorable case where each cluster is very localized, this is indced likely to lead to a new correlation function fr with no local maximum around each of the already found centroids. However, the picture is quite different in the more realistic case where clusters have non-negligible (and possibly non-isotropic) intra-cluster variance. In the latter case, the corrclation function f, of the residual updated CL-OMPR indcod often still has a local maximum quite close to an already found centroid, leading the algorithm to repeatedly select the same dominant clusters and missing the other ones. This is illustrated in Figure in Section 4.2 where we propose fix. Finally, observed that due to Step 2 of Algorithm D after first increasing the number of selected centroids from 1 to k, CL-OMPR repeats T - k times the addition of a new one followed by the removal of one, thus never exploring more than h + candidates at a time. As we will, a variant allowing more exploration is beneficial.
4 Towards a robust decoder for compressive clustering
In this soction, we propose a decoder that overcome the limitations of CL-OMPR mentioned in the end of Section 3.2 For this, we introduce three ingredients: i) several approches to handle the non-couvex optimization problem that appears in Step of CL-OMPR, which consists in seeking the centroids among the local maxima of the correlation function fex: ii) reducing the support C through hard-thresholding is executed after selecting all the candidate centroids, iii) allowing to fit a ICl-mixture of Gaussians instead of Diracs to take into consideration clusters with non-negligible (and possibly nonisotropic) intra-cluster (co)variance: considering the selected e from Step 1 as the mean of a new Gaussian, this requires estimating the covariance matrix of this added Gaussian.
Algorithm 2 proceeds in two steps. The first one, which is reminiscent of Orthogonal Matching Pursuit PPai et al. 1993 Mallat & Zhang 1993). consists seeking an initial support Ō O of 2 candidate centroids. This is achieved by ileralively seeking a candidate point o; that locally maximizes the correlation function fr, then seeking the vector a R4 that minimizes residual, which expression depends on the fitted model: mixture of Diracs, or a mixture of Gaussians. The latter requires the estimation of the local covariance matrix Σ; using a procedure Est EstimateSigma, that we will motivate in Section 4.2 and detail in Appendix B This procedure either outputs a valid (i.e., positive definite) covariance matrix, or zero, in which case the corresponding component is set to a Dirac. Observe that updating the residual in every iteration crucial (but not sufficient- this motivates to consider Gaussian mixtures) order to obtain candidate centroids that are different at each iteration. This step depends the fitted model. When fitting mixture of Gaussians, updating the residual boils down to calculating the terms (N Σ,), which are equal to the characteristic functions of the corresponding Gaussians when A is build using random Fourier features (Keriven) 2017 Equation 5.7). The second step consists in pruning the support to C' to obtain the targeted number l of centroids.
We next discuss the choice of the procedures GetLocalMaximi and EstimateSigma.
4.1 Detecting local maxima of the correlation function
We describe in this section, three algorithms to implement GetLocalMaximum in Step 1 of Algorithm 2 The pros and cons of these variants are discussed theoretically and will be assessed numerically in Section
Discretized approach: The first option consists in replacing the optimization problem maX fe(c) by max fe(c), where Ô θ is (finite) discretization of @. This approach is recurrent in the literature of
6

---

Published in Transactions on Machine Learning Research (04/2024)
Data: Sketch xx, sketching operator A, number of centroids k number of atoms T 2 k, fitted model 'model', domain @ C IRd
Result: The set of centroids C' = fei....,ca. Ck), the vector of weights a (a1,..., ,gk
r f- 2x; C++
Step Look for on initial support
for i = 1.....7 T do
(x) t Re(r, A5a)
G - GetLocalMaximmnt fr; Q) //Find new centroid
C Cecur (G) / //Expand the support
if model = Dirac then
a arg min OER Cl ||E E/ 2j= Ab,b //Project to find the weights
10l r 2x - 1 ajAbe, //Update the residual
end
else if model Gaussian then
Σ, EstimateSigma(fxx ctt/sssiamm///mmttttttet the covariance matrix if Σ, = 0 then
T < (c,, /.e//Define new Gaussian component
end
else if Σ, + 0 then
n. t Se, //Revert to Dirac component
arg min
DER_ || E/ An //Find the weights
CI 2x Σ 11 aN E,) //Update the residual
end
Step 2: Reduce the inilial support when k if [C| > k then Select indices i1,-..,44 is of the largest elements of oz Oip Oii C+ (c....., , //Reduce the support
a (aig...., ,...) //Reduce a
end
Algorithm 2: Proposed docoder
40
40
48
(a) A comparison of the dynamics
(b) [ lla (left) vs. [| fax lla (right) in log scale
Figure 3: A comparison between plain gradient ascent and sketched mean shift in the identification of a local maximum of fzx; dynamic ranges of V l and [[ llz.
superresolution (Tang et al. 2013 Traonmilin al. 2020 where θ is usually a subset of low-dimensional vector space. We show later that this approach yields significant improvement over CL-OMPR. However, in high-dimensional settings θ c Rd, the cardinality of é must grow exponentially with the dimension d.
Sketched mean shift approach: Instead of discretizing, we may use a gradient-based method to identify local maxima of the correlation function. In this context, a simple gradient ascent is ineffective in practice

---

Published in Transactions on Machine Learning Research (04/2024)
because the gradient V fr is vanishingly weak beyond the immediate vicinity of some points, resulting in slow gradient dynamic; see Figure for an illustration in the case of r zx. As an alternative, we propose to perform reweighted gradient ascent iterations
c cllcc [f((c]
" Vf((c)),
(9)
where Ill is the projection onto θ with respect to the Euclidean distance. This boils down to the ascent algorithm based on the gradient of log fr. when f. takes positive values. As illustrated in Figure [a in the case of the correlation function fax, this algorithm exhibits improved dynamics compared to plain ascent of the gradient of f,: the trajectory from the same starting point hardly changes when using plain gradient ascent iteralions, but does reach a centroid when performing the reweighted gradient ascent of equation Figure 3b compares the norm of Vfex the norm of / we observe that the former is very small most regions while the latter reaches much larger values in general, except nearby stationary points where it vanishes. Now, equation g is reminiscent of the mean shift algorithm (Fukunaga & Hostetler 1975 Cheng 1995p mentioned in Section 3.2 The latter boils down to gradient ascent applied to the logarithm of the KDE, defined by equation when Ko is the Gaussian kernel, hence the name sketched mean shift. Still, the usual mean shift algorithm requires access to the whole dataset X' every iteration, while equation doesn't, thanks the very principle of the sketching approach which summarizes the dataset. into xx- Observe that making use of the iteration equation requires that the function fr does not vanish on the whole path of optimization. This properly was empirically observed in the numerical simulations presented in Section
The sketched mean shift algorithm might still converge to spurious local maximum, particularly when the sketch size m is low. Thus, the best output of sketched mean shift over L random initializations of c is computed. As we will see in the experiments Section 5 this approach is more efficient than discretization
4.2 The importance of the fitted model
Figure f illustrates the benefit of fitting a k-mixture of Gaussians instead of a k-mixture of Diracs during the iterations of Algorithm 2 with mixtures of Dirac, once the residual r is updated, the corresponding correlation function fr still has a (high) local maximum close to the previously found centroid, leading the algorithm to select a nearly identical cluster. This no longer happens when fitting a mixture of Gaussians, with properly chosen covariance. The implementation of this in Algorithm 2 requires to estimate the local covariance matrix once centroid is selected in Step 1. In Appendix B we propose estimator of this matrix. As we shall see in Section 5 this allows to robustify further the decoder for lower values than if we simply fit k-mixture of Diracs. The theoretical study of this estimator, especially in high-dimensional domains, left for a future work.
(a) Two successive iterations of Algorithm when fitting a mixture of Diracs
(b) Two successive iterations of Algorithm when fitting mixture of Gaussians
Figure 4: The correlation function f. associated to the residual r at the first and second iteration of Algo- rithm using two fitted models: mixture of Diraes and mixture of Gaussians. The red points correspond the selected centroids.

---

Published in Transactions on Machine Learning Research (04/2024)
5 Numerical simulations
In this section we present numerical experiments that illustrate the performance of the newly proposed decoder. In Section 5.1 we conduct numerical experiments on synthetic datasets to study the performance of the decoder in high-dimensional settings. In Section 5.2 we investigate the importance of the fitted model using synthetic datasets. Finally, in Section 5.3 we perform the experiments on real datasets.
5.1 The discretized approach compared to the sketched mean shift approach
In this section, we compare the two variants of Algorithm M proposed in Section 4.1 For this purpose, we conduct numerical simulations on a synthetic dataset X C RR made of three clusters: the dataset X' (IN contains N 100000 vectors in R which are obtained by ii.d. sampling from a mixture
of isotropic Gaussians Σi= (c,((c,E)) where k = 3, 01 = 02 = as = 1/3, c1, C2 E RR, and E1 -.. Ex T2 12 E R222, The centroids 6; and variance o%2 were chosen using a function implemented in the Python package Pycleto to ensure enough separation while essentially fitting the dataset in [1, For all algorithms we use as a domain O the hypercube O [-1,
Performance is defined relative to the performance of Lloyd's algorithm, ie., using normalized version of the MSE called the relative squared error (RSE) and widely used in previous work on compressive clustering. The RSE is defined as
R)) = MSE(C; X)
(10)
MSE(C 2))
where MSE(C; X)) is given by equation 5 and CLooyd is the configuration of centroids given by Lloyd's algorithm (bost of runs with different controid sceds).
Figure 5a compares the RSE of compressive clustering using CL-OMPR, compressive clustering using Al- gorithm based on the discretized approach, and compressive clustering using Algorithm 2 based on the sketched mcan shift approach. By definition, the RSE of Lloyd's algorithm is onc.
For each compressive approach, the RSE averaged on 50 realizations of the sketching operator and presented for two sketch sizes m E {2000, 1000}. For the three compressive algorithms, we take T 2k. Moreover, for the two variants of Algorithm 2 we take L = 10000 random initializations that are i.i.d. samples from the uniform distribution on θ 11. We observe that both CL-OMPR and Algorithm 2 based on the discretized approach fail to match or even approach the performance of Lloyd's algorithm (i.e., to achieve an RSE close to one) for any bandwidth of a, while the performance of Algorithm 2 based on the sketched mean shift approach is very similar to that of Lloyd's algorithm for m = 200, and matches it for m = 1000 in a range of bandwidths (a & [0.110.3)). In other words, CL-OMPR fails to capture the local maxima of fr- The use of weighted gradient descent as in equation D in the second variant allows to better capture the local maxima, hence Algorithm 2 based on the sketched mean shift approach outperforms Algorithm 2 based on the discretized approach. As result from now on when considering Algorithm 2 we systematically focus the sketched mean shift approach.
Now, we examine how the number of initializations L influences the performance of Algorithm 2 Figure compares the RSE of compressive clustering using CL-OMPR, and compressive clustering using Algorithm basod on the sketched mcan shift approach associated to various values of L. We observe that the range of the bandwidth a for which Algorithm 2 matches Lloyd's algorithm (i.e., achieves an RSE close to one) increases with L, for both considered values of m. In particular, we observe that the performance of Algorithm 2 for 1000 and a > 0.4 when m 200) does not depend on L, while it is
large values of (o > 0.3 when highly dependent on the value of L for low values of σ. In light of the 2D illustrations of Section B this is likely due to the existence of spurious local maxima for low values of the bandwidth a. Indeed, spurious local maxima proliferate as a gets smaller; see Figure for an illustration. This numerical experiment shows that the recovery of the centroids from the sketch using Algorithm 2 is possible even for low values of by increasing L. Figure 6 shows the relative squared error (RSE) of compressive clustering using CL-OMPR,
htttpp://ggtaub..o///ceeelll The datasel and the code
ppycle
generate
similar
one
can
be downloaded from
wftte//ofo/mddaaaeeeee WoWWng
9

---

Published in Transactions on Machine Learning Research (04/2024)
Algorithm 2 (Sketched mean shitt)
Algorithm (sketched mean shift
Algerihh 2 (dsccretize)
Algoritim 2 (discretizedi algorithm
Lloye's algorithm
RSE
RSE
(a) Average RSE of CL-OMPR for two sketch sizes (left: 200; right:
1000), and Algorithm
with two
approaches (discretized and sketched mean shift, I. 1000).
CLOMMR
CL-OMPR
2 2 (sketchdd mean shift, L-
Algontrm (sketched mean shift L 1)
Aigor Ihm (shetched mean shit, 10)
Algorithen 2 (slatched mean shift 10)
(sketched shit.
Algorithm (skatched
Aigar tim (Abetched
1000)
Algorithm (sketched mean shift, 1000)
Lioyd's algeritm
Lkkyd's alger tte
RSE
RSE
5
2
a
(b) Average RSE of CL-OMPR for two sketch sizes (left: 200; right: 1000), and Algorithm using sketched mean shift for T. E (10, 100, 10003.
Figure 5: Comparison of CL-OMPR and Algorithm with three synthetic clusters [-1, 1]6
Algorithm 2 based on the sketched mcan shift approach (with thrce values of L 10, L 100 and L 1000), using 50 realizations of the sketching operator and presented for various values of the sketch size and the bandwidth parameter
We observe that the range of the bandwidth for which Algorithm 2 matches the performance of Lloyd's algorithm (this corresponds to reaching an RSE close to 1) increases with m and L. More precisely, as L gets larger, the performance of Algorithm g improves for & 0.3, yet without an improvement for larger values of (a 2 0.5). In comparison, the average RSE of CL-OMPR stays above for every value of the sketch size m and the bandwidth
Remark 1. All experiments in the main body of this paper are conducted in moderate dimension d S 10, which is pretty much the standard dimension of state of the art proofs of concept in sketched clustering. While Algorithm g is designed to overcome the important issues of CL-OMPR scrutinized in this poper, we not expect this algorithm to be particularly good at overcoming issues related to a curse of dimension. In particular, is expected that in high dimension, success will require an exponential growth of L. This is empirically conftrmed with additional experiments in Appendir C.1 Another algorithmic challenge left to
We Figure 6ala threshold (RSF 5) different than threshold (RSE used in Figures
10

---

Published in Transactions on Machine Learning Research (04/2024)
future work is to handle large numbers K of clusters, EL-OMPR. is known to badly scale with K Keriven et ol. 2017).
5000
5000
2000
2000
1000
RSE
FSE
0.05
0.07
0.1
0.2
0.3
0.4
0.5
0.7
0.05
0.07
0.1
0.2
0.3
0.4
0.5
0.7
(a) CL-OMPR
(b) Algorithm
(L 10)
5000
5000
2000
2000
1000
1000
500
500
200
200
0.5
(c) A Agorithm
(T. = 100)
(d) A Algorithm
2
1000)
Figure 6: Comparison of CL-OMPR and Algorithm 2 on three synthetic clusters in [-1, 1]
5.2 On the importance of the fitted model
In this section, we investigate the importance of the fitted model. For this purpose, we compare fitting a k-mixture of Gaussians instcad of a k-mixture Diracs on two synthetic datasets. The first one is the one considered in Section 3 while the second one is formed by clusters [_1, that are not well separated; see Figure 7b We compare Algorithm P using fitting with k-mixture of Diracs to Algorithm 2 using fitting with k-mixture of Gaussians, In the latter case, the local covariance matrix is estimated
detailed in Appendix B
For each dataset, use as baseline the best out of 5 replicas (with different centroid seeds) of Lloyd's algorithm. For each instance of Algorithm 2 (with Diracs, or with Gaussians; both using the sketched mean shift approach with L = 1000), and each considered value of the sketch size m and the bandwidth parameter 0, we compute the average RSE (cf equation 10 over 10 realizations of the sketching operator. Figure shows the obtained RSE for various values of the sketch size m and the bandwidth o.
In the case of the dataset represented by Figure 7a observe that using Gaussians instead of Diracs primarily enlarges the range of the bandwidth for which Algorithm 2 matches Lloyd's performance (RSE close to onc). For the the datasct represented by Figure 7b. we make several observations: i) when fitting mixtures of Diracs, the largest bandwidth where the perfomance of Lloyd's algorithm is matched is smaller than with the dataset of Figure 7a even for large values of the sketch size (m 2 1000), this is due to the increased difficulty in scparating the clusters; ii) fitting mixture of Gaussians improves the performance of Algorithm 2 for low values a: the smallest bandwidth where the perfomance of Lloyd's algorithm is matched is smaller than when fitting mixture of Diracs.
11

---

Published in Transactions on Machine Learning Research (04/2024)
1.00
1.00
0.50
0.50
0.25
0.25
0.00
0.00
0.25
0.50
0.75
0.75
-1.00
-1.00
-1,00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
-1.00
-0.75
-0.50
-0.25
0.00
0.25
0.50
0.75
1.00
(a)
(b)
5000
2000
2000
1000
14
ASE
500
500
11
100
100
0.01
0
(c)
(d)
2000
2000
E
E
RSE
500
500
100
100
0.005 0.01 0.02 0.03
0.05
0.2
0.4
0.5
0.7
0.005 0.01 0.02 0.03
0.05
0.1
0.2
0.3
0.4
0.5
0.7
(f)
Figure Performance of compressive clustering using Algorithm on two datasets (left column: a well- separated dataset; right column: a dataset with closeby clusters). Top: representation of the datasets; middle: average RSE obtained for various sketch sizes and scales using mixture of Diracs model. Bottom: average RSE achieved using a mirture of Gaussians model
5.3 Experiments on MNIST
In this section, we investigate whether the obscrvations of Soction 5.2 hold for real datasots. For this purpose, we perform experiments on spectral features of the MNIST dataset which consist. of = 70000 handwrit- ten digits features with = 10 classes. These spectral features are computed by taking the eigenvectors associated the 10 smallest positive eigenvalues of the normalized Laplacian matrix associated to the
nearest neighbors matrix MMuja & Lowe 2009), associated to SIFT descriptors of each image (Vedaldi &
Further experiments CIFARIO Appendix
C.2 display similar behavior.
12

---

Published in Transactions on Machine Learning Research (04/2024)
5000
5000
2000
14
2000
1000
-18
1000
500
500
200
200
100
100
0.1
(a) CL-OMPR
(b) Algorithm
(L 10)
5000
5000
2000
2000
3.4
12
100
100
0.1
0.2
0.3
0.4
0.5
0.7
0.1
0.2
0.3
0.4
0.5
0.7
(c) Algorithm
100)
(d) Algorithm
2 1000)
Figure
comparison of CL-OMPR and Algorithm on spectral features of MNIST
Fulkerson 2010 The resulting matrix can downloaded from
www.:aa...o../oaaeoooommooooooot /
8 shows the RSE, defined by equation 10 based on Lloyd's algorithm (best. of 5 runs of Lloyd's 's al- gorithm with different centroid soeds), of compressive clustoring using CL-OMPR, and Algorithm with mixtures of D and the sketched mean shift approach (with three values of L: L 10, L 100, and L 1000), based 10 realizations of the sketching operator for various values of the sketch size and the bandwidth a. Similarly to the experiment conducted above on the synthetic dataset, we observe that Algo- rithm 2 ooutperform CL-OMPR for a S 0.5, and get a performance close to Lloyd's algorithm (RSE < 1.5) on a range of a that increases with m and L. In particular, for T. = 1000, Algorithm 2 achieves an RSE smaller than 1.5 for m = 500 (only five times the number degrees of freedom kd 100 needed to describe k 10 centroids in dimension achieved for sketch size 10 times larger: lating the sketch is approximately 0.2 seconds, and the average rumning time for decoding using CL-OMPR is about 3 seconds. Meanwhile, the running time for decoding using Algorithm g depends linearly the number of initializations L; for instance, it averages 12 seconds when L = 10. Improving the running time of Equation 1 is left to future work, as our primary effort in this paper in to demonstrate improved and robustified RSE performance over a wide range of values of and m compared to CL-OMPR.
Figure
10), while the best achievable RSE by CL-OMPR close to and was
5000. For this dataset, the average running time for caleu-
6 Conclusion
By diagnosing CL-OMPR, we were able to propose a robustified decoder which significantly outperforms CL-OMPR on synthetic datasets. In particular, we demonstrated on a simple dataset that CL-OMPR fails, while algorithm matches the performance of Lloyd's algorithm, even though (unlike Lloyd's algorithm) it only has access to a low-dimensional sketch, which dimension is independent of the size of the dataset. Moreover, compared to previous proofs of concept of compressive learning, our algorithm is able to extract clustering information from smaller sketches, therefore further reducing the needed memory footprint. Given the potential of comprossive clustcring mmeoory-constrimee scenarios, with distributed implementations, in streaming contexts, or with privacy constraints), this work thus opens up unprecedented prospects for making
6 Algorithm with Ganssian mixture led worse results, likely due limitations of the covariance estimation procedure of A Appendix dimension d= Improving this procedure for high dimensions challenge left. future work.
13

---

Published in Transactions on Machine Learning Research (04/2024)
the full pipeline of this paradigm both competitive, casy to implement and easy to tune. Considering the choice of the bandwidth scale o, which is currently something of an art. the robustness of our algorithm combined with recently proposed selection criteria (Giffon & Gribonval 2022) open a promising avenue to design a fully turnkey pipeline for compressive clustering. The next upcoming challenge
course validate the versatility of the resulting approach on datasets living in high-dimensional domains, while controlling the number of required intializations of the sketched mean shift control its computational complexity. To achieve this goal, leveraging alternative feature maps would be beneficial (Avron et al,, 2017 Chatalic et al. 2022). On another vein, promising research path involves estimating local covariance matrices from the sketch. Given the strong connection of our approach with the mean shift algorithm, for which the convergence was actively investigated recently (Li et al. 2007 Ghassabeh] 2015 Huang et al. 2018 Yamasaki & Tanaka 2019 2023), we aim to investigate the theoretical properties of sketched mean shift. Finally, the connection of the mean shift algorithm to EM algorithm (Carreira- Pareeiia-eeriaann 2007) suggests the possibility extend our approach in paramcter estimation of more general mixture models such mixture of alpha-stable distributions (Keriven et al., 2018b
Acknowledgments
This project was supported by the AllegroAssai ANR project ANR-19-CHTA-0009. The authors thank the Blaise Pascal Center (CBP) for the computational means. It uses the SIDUS solution developed by Quemener & Corvellec 2013).
14

---

Published in Transactions on Machine Learning Research (04/2024)
References
H. Avron, M. Kapralov, C. Musco, Ch. Musco, A. Velingker, and A. Zandieh. Random fourier features for kernel ridge regression: Approximation bounds and statistical guarantees. In International conference on machine learning, pp. 253 262. PMLR, 2017.
A. Belhadji and R. Gribonval. Revisiting RIP guarantees for sketching operators on mixture models. Novem- ber 2022. URL hette://baassseeeec//hl111177771 preprint.
M. A. Perreina- Gaussian mean-shift is an em algorithm. IEFE Transactions on Pottern Analysis and Machine Intelligence, 29(5):767- 776, 2007.
A. Chatalic, R. Gribonval, and N. Keriven. Large-scale high-dimensional clustering with fast sketching. In 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 4714 4718. IEEE, 2018.
A. Chatalic, L. Carratino, E. De Vito, and L. Rosasco. Mean nyström embeddings for adaptive compressive learning. In International Conference on Artificial Intelligence and Statistics, pp. 9869 9889. PMLR. 2022.
Y.C. Chen. A tutorial on kernel density estimation and recent advances. Biostatistics & Epidemiology, 1 (1):161-187, 2017.
Y. Cheng. Mean shift, mode seeking, and clustering. IEEE transactions pattern analysis and machine intelligence, 17(8):790-799, 1995.
Farah Cherfaoui, Valentin Emiya, Liva Ralaivola, and Sandrine Anthoine. Recovery and convergence rate of the frank wolfe algorithm for the m-exacl-sparse problem. IEEE Transactions on Information Theory, 65(11):7407-744) 2019. doi: 10. /III 201..................
D. Comaniciu and P. Meer. Mean shift: A robust approach toward feature space analysis. IEEF Transactions on pattern analysis and machine intelligence, 24(5):603 619, 2002.
G. Cormode, M. Garofalakis, J. Haas, and Ch. Jermaine. Synopses for massive data: Samples, histograms, wavelets, sketches. Foundations and TrendsS in Databases, 4(1-3):1-294, 2011.
Quentin Denoyelle, Vincent Duval, Gabriel Peyré, and Emmanuel Soubies. The sliding Frank Wolfe algo- rithm and its application to super-resolution microscopy. Inverse problems, 36(1):014001, Jamuary 2020. doi: 10.1088 10.1088/1361-644 ab2229. URL lattps://oosccienceee http:ee......ogooaeeeeeeeeeeeeee 10001008811100//////11
ab2a29 Publisher: IOP Publishing.
Dan Feldman and Michael Langberg. A Unified Framework for Approximating and Clustering Data. STOC, Ls.69: 578, 2011. URL
K. Fukunaga and L. Hostetler. The estimation of the gradient of density function, with applications in pattern recognition. IEEE Transactions information theory, 21(1):32-45 40, 1975.
Y. A. Ghassabeh. A sufficient condition for the convergence of the mean shift algorithm with gaussian kernel. Journal of Multivariate Analysis, 135:1-10, 2015.
L. Giffon and R. Gribonval. Compressive clustering with optical processing unit. GRETSI Colloque Francophone de Traitement du Signal et des Images, 2022.
R. Gribonval, G. Blanchard, N. Keriven, and Y. Traonmilin. Compressive statistical learning with random feature moments. Mathematical Statistics and Learning, 3(2):113-1664 2021a.
R. Gribonval, A. Chatalic, N. Keriven, V. Schellekens, L. Jacques, and P. Schniter. Sketching data sets for large-scale learning: Keeping only what you need. IEEE Signal Processing Magazine, 38(5):12 36, 2021b.
Chengeheng Guo, Bo Zhao, and Yanbing Bai. DeepCore: A Comprehensive Library for Coreset Selection in Deep Learning, June 2022. URL httt:://axxxv http//aa////....0a5///3// arXiv:2204.08499 [cs].
15

---

Published in Transactions on Machine Learning Research (04/2024)
K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. In Proceedings of the IEFE conference computer vision and pattern recognition, pp- 770-778, 2016.
K. Huang, X. Fu, and N. Sidiropoulos. On convergence of epanechnikov mean shift. In Proceedings of the AAAT Conference on Artificial Intelligence, volume 32, 2018.
A. K. Jain. Data clustering: 50 years beyond k-means. Pultern recognition letters, 31(8):651-660 2010.
N. Keriven. Apprentissage de modèles de mélange à large échelle par Sketching. PhD thesis, Rennes 2017.
N. Keriven, N. Tremblay, Y. Traonmilin, and R. Gribonval. Compressive k-means. In 2017 IEEE Interna- tional Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 6369 6373. TEEE, 2017.
N. Keriven, A. Bourrier, R. Gribonval, and P. Pérez. Sketching for large-scale learning of mixture models. Information and Inference: A Journal of the IMA, 7(3):447-3) -008 2018a.
N. Keriven, A. Deleforge, and A. Liutkus. Blind Source Separation Using Mixtures of Alpha- Stable Distri- butions. In 2018 IEEE International Conference Acoustics, Speech and Signal Processing (ICASSP), pp. 771-775, April 2018b. doi: 10.1109 11009IIAASPP2008.84600000 URL https://ieeexpports ieee. org/ document/8462095, ISSN: 2379-190X.
X. Li, Z. Hu, and F. Wu. A note the convergence of the mean shift. Pottern recognition, 4066):1756 2007.
S. Lloyd. Least squares quantization in pem. IEEE transactions on information theory, (8:2):129 137, 1982.
S. Mallat and Z. Zhang. Matching pursuits with time-frequency dictionaries. IEEE Transactions on signal processing, 41 41(12):3397 3337-333373334 1993.
K. Muandet, K. Fukumizu, B. Sriperumbudur, B. Schölkopf, et al. Kernel mean embedding of distributions:
A review and beyond. Foundations and TrendsS in Machine Learning, 10(1-2):1 141, 2017.
M. Muja and D. G. Lowe. Fast approximate nearest. neighbors with automatic algorithm configuration. VISAPP (1), 2(331-340):2, 2009.
Y.C. Pati, R. Rezaiifar, and P.S. Krishnaprasad. Orthogonal matching pursuit: Recursive function approxi- mation with applications to wavelet decomposition. In Proceedings of 27th Asilomar conference on signals, systems and compulers, pp. 40-44. TEEE, 1993.
E. Quemener and M. Corvellec. Sidus the solution for extreme deduplication of an operating system. Linax Journal, 2013(335):3, 2013.
A. Rahimi and B. Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007.
G. Tang, B. N. Bhaskar, and B. Recht. Sparse recovery over continuous dictionaries-jus discretize. In 2013 Asilomar Conference on Signals, Systems and Computers, pp. 1043-1047. IEEE, 2013.
Y. Traonmilin, J.-F. Aujol, and A. Leclaire. Projected gradient descent for non-convex sparse spike estima- tion. IEEE Signal Processing Letters, 27:1110 1114, 2020.
A. Vedaldi and B. Fulkerson. Vlfeat: An open and portable library of computer vision algorithms. In Proceedings of the 18th ACM international conference on Mullimedia, pp. 1469 1472, 2010.
Yamasaki and Tanaka. Properties of mean shift. IEEE bransactions on paltern analysis and machine intelligence, (2233-227322 2286, 2019.
R. Yamasaki and T. Tanaka. Convergence analysis of mean shift. arXiv preprinl arXin:2305.08 084663 2023.
16

---

Published in Transactions on Machine Learning Research (04/2024)
A The CL-OMPR algorithm
Data: Sketch 2x, sketching operator A, parameter k, bounds e,
Result: Centroids C, weights a
r f- 2x; C 0;
for 1 1. ... 2k do
Step 1: Find new centroid
c maxe Re A8
end
Step 2: Espand support
C' C CUCufc
end
Step 3: Enforce sporsily by Hard Thresholding when L k
if [C] h then
8 arg m2n xx E/El Bi
Abe, TT
Select k largest entries Be,..... Biz Reduce the support C .ci....., Cia
end end
Step 4: Project to find a
arg min xx E AnAdes
end
Step 5: Global gradient descent
C, + mincaa (1ex- E1Cl a;A .
end
Update residual: r- 2x E Ad
end
A Algorithm 3: CL-OMPR
B An estimator of the local covariance matrix
Algorithm 2 allows to fit k-mixture of Gaussians, which requires to estimate the local covariance matrix once a centroid is selected in Step 1. In this section, we propose an estimator of this matrix. This provides an implementation of EstimateSigma for Algorithm 2
To begin with an intuition, consider the following setting: the x; are i.i.d. draws from a mixture of isotropic Gaussians D a,N(c, E,), where oi,...,ak [0, 1] such that Σ 04 1 and CI.....Cy Ck and E,,..., Ex st+ with & N We consider sketching operator defined through random Fourier features associated to i.i.d. Gaussian frequencies w1,...,WM WM drawn from (V, with > 0. As the number of samples N goes to too the KDE conyerges to the kernel mean embedding JKME of the Gaussian mixture Σ 1 (c; Ex) (see Section 3.1.2 in (Muandet et al. 2017)). ) which is given by (see Lemma 6.4.1 in [eriien 2017
Rd
k
Σ Us exp ( 1 (x - cF)T(EA + -ZId) (x-)))x))))))))) Rd.
fKME(I) :
(11)
17

---

Published in Transactions on Machine Learning Research (04/2024)
Now, when the clusters c ck are well-separated, we have for i E {1....,k}
log fKME Rota - 2 ((a-aa)(Ea + (o - -c)= + log Us
(12)
In other words, log /KMR behaves locally as quadratic function associaled to the p.s.d. matrix (E; + olla)= Thus, if ft is an estimate the Hessian of - -lg fKME then A - o² is an estimate of E, Since an estimated covariance matrix must be p.s.d., the matrix H o turns out be not p.s.d., then we set as our estimate Σ, == to revert to Dirac component, see Algorithm Thus, to build an estimator of E; we need to estimate the Hessian of log SKME-
The following estimator is simpy the Hessian of - log fax evaluated at the point and may be seen surrogate of the Hessian of - log fRME evaluated at the point
Definition 2. Let e Rd such that Jax (c) > O. Desine the matria A by
A: - (Hess(Jax: (c) -VV Vel
(13)
where Hess(Saxic) is the Hessian matria of the function fax evaluated at
Now, with & the random Fourier feature map with components ow, defined by equation function fex from Definition
the correlation
is
T
(cx
Re mm
2ythy (c),
(14)
j=1
where 21,..., 2m E C are the entries of the sketch zx. Thus
Velex
Re
2j4n sow (c),
(15)
Vm
and the Hessian matrix Hess(faxi c) given by
1 mm
Hess(Jzxic) =
Re
رسرو Bwy (c),
(16)
In other words, the matrix A defined by equation values of the sketch 2x and the frequencies w Wyne
13 be numerically evaluated at given e using the
C Additional numerical experiments
C.1 Growth of with the dimension d
In this scction, bricfly illustrate that the number of initializations L needed to achieve good performance sulfers from the curse of dimensionality. For this, we conduct the same experimenl across different dimensions d E {2, 3, 5, 10, 20): the dataset X' is formed by N = 10, 000 elements in Rd obtained through a generative model. The elements I; are i.i.d. samples from a mixture (.,,(c1, TZ + 0.5. o where éj and 02 are centers of the two normal distributions, dp > 0, and I is the d-dimennsional identity matrix. Figure g shows the RSE equation 10 compressive clustering using the sketched mean shift approach with m = 10d based on 10 realizations of the sketching operator for L E (1, 10, 100, 1000} and d E 2, 10, 20), for two different values of bandwidth This empirically confirms that L should grow with the dimension d.
C.2 The case of CIFAR-10 10
In this section, we investigate whether the observations of Section 5.3 hold for the CIFAR-10 dataset. For this purpose, perform experiments the training set of the CIFAR-10 dataset, which consists of N 60000
18

---

Published in Transactions on Machine Learning Research (04/2024)
1000
1000
100
100
1-
(a) 10d, o 0.1
(b) 10d, = 0.2
Figure 9: The dependency between L and d
images features with = 10 classes. More preciscly, we extracted features from the last fully connected layer (of dimension D - 512) of a trained ResNet18 (He et al. 2016). These extracted features are then reduced dimension of d 10 using linear PCA. The network is trained on the training set of CIFAR-10 for 50 epochs with SGD with momentum 0.9, learning rate 0.1, learning rate decay 0.1, batch-size 512 and weight-decay 5 10
Figure 10 shows the RSE, defined by equation 10 based on Lloyd's algorithm of compressive clustering using CL-OMPR, and Algorithm 2 with mixtures of Diracs and the sketched mean shift approach with L = 250, based on 10 realizations of the sketching operator for various values of the sketch size m and the bandwidth o. Similarly to the experiments conducted in Section 5.3 we observe that A Algorithm performs CL-OMPR for S 0.5, and get performance close Lloyd's algorithm on range of that increases with m. We
have observed running times that of the samee order of magnitude as in Section 5.3
2-5
5000
5000
2000
2000
RSE
200
(a) CL-OMPR
(b) Algorithm 250)
Figure 10: comparison of CL-OMPR and Algorithm on the features of CIFAR-10
19


Reviews:
[{"review_id": "5ypBdCVekr", "reviewer": "Reviewer_XiFY", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n\n1. The paper is very instructive, demonstrating the issues of the existing approach in several aspects (e.g., vanishing gradients, multiple local maxima, lack of exploration). Each of these aspects is then addressed by certain modification: normalization of gradient ascent step by the function value, replacing the mixture of Diracs with mixture of Gaussians. \n\n2. The performance gains presented in Section 5 (especially Figure 5) are substantial. \n\n3. The experiments are mostly performed in a controllable environment with small datasets, which allows to understand the issue and potential gains of the new algorithm. The error bars are included where appropriate. \n\n **Weaknesses:**\n\n1. The illustration/description of experiments in section 5.2 is unclear to me. Figure 7 has several duplicate captions for different plots, which makes it difficult to understand the results. I assume the left plots (a), (c) and (e) correspond to a separable case and the right plots are for non-separable? \n\n2. In section 5, the paper suddenly introduces RSE metric instead of previously used MSE metric. Is this a standard metric in the field? It's unclear to me why Figure 5 uses MSE and the rest of experiments use RSE. Some explanation/motivation would be helpful. \n\n3. It would be useful to include the running time and memory cost of the proposed algorithms. \n\n**Questions:**\n \n1. What is exactly meant by “phantom local maxima”? Cannot find the definition for this term in the paper. \n\n2. What is the connection between CL-OMPR and Frank-Wolfe method (stated in the abstract)? Does the connection also hold for Alg. 2?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a new heuristic algorithm for compressive clustering of data. Through carefully designed synthetic experiments, the authors are able to identify the deficiencies of the baseline method (CL-OMPR) and propose and algorithmic fixes in several dimensions. The modified algorithm is then further tested on synthetic data and MNIST dataset.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "For requested changes, see \"Weaknesses\" section above. \n\n\n**Typos:**  \n\n1. Page 2 \"In Section 3, we show conduct\", \"ths non-convex\", Page 6 \"litterature\", Page 8 \"we we perform\", \n\n2. “In this section, we investigate whether the observations of Section 5.3 hold for real datasets” Should be Section 5.2?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "XBQWmhAsmA", "reviewer": "Reviewer_hypU", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n- an improvement over CL-OMPR is clear for considered datasets\n- natural while interesting ideas on improvement at different steps\n\nWeaknesses\n- no analysis for large datasets (CIFAR10 or CIFAR100 would work here)\n- no theoretical properties proved; only a partial discussion is conducted\n- no connection to the core set literature [1] and the literature on transforming the training sample to a set of pseudoimages [2, 3]\n- the findings from the papers above can be used to enhance experiments with additional checks that the new objects from clustering are beneficial\n- the problem with the underrepresentation of other papers is highlighted by citing mostly papers from a single research group, more diversity would help\n\n\n1. Guo, Chengcheng, Bo Zhao, and Yanbing Bai. \"Deepcore: A comprehensive library for coreset selection in deep learning.\" International Conference on Database and Expert Systems Applications. Cham: Springer International Publishing, 2022.\n2. Inoue, Hiroshi. \"Data augmentation by pairing samples for images classification.\" arXiv preprint arXiv:1801.02929 (2018).\n3. Naveed, Humza, et al. \"Survey: Image mixing and deleting for data augmentation.\" Engineering Applications of Artificial Intelligence 131 (2024): 107791.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper addresses the problem of compressive learning. The goal is to reduce the whole dataset to a set of points (in the embedding space). A previous approach considered is CL-OMPR. It has a number of problems, including numerical instability and a requirement for a high-precision tuning of hyperparameters. Another alternative, LLoyd’s algorithm, is also computationally demanding. The authors propose improving the CL-OMPR approach and call the result Algorithm 2. We see that the proposed approach can recover the core set for artificial datasets and the classic MNIST.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Inclusion of Analysis for Large Datasets:\nConduct a comprehensive analysis using larger datasets such as CIFAR10 or CIFAR100 to evaluate the scalability and performance of the proposed methods in more complex scenarios. This will provide a better understanding of how the model behaves with an increased number of classes and higher-dimensional data.\n\nProving Theoretical Properties:\nDevelop and present formal proofs for the theoretical properties discussed in the paper. A rigorous mathematical foundation will strengthen the validity of the partial discussions and support the empirical findings. \n\nConnection to Core Set and Training Sample Transformation Literature:\nIntegrate references to the core set literature, explicitly drawing parallels and distinctions with the work presented in [1]. Furthermore, engage with the methodologies on transforming the training sample to a set of pseudoimages, as seen in [2, 3]. This should not just be a literature review but should also include a discussion of how these approaches relate to and can potentially enhance the current work.\n\nEnhancement of Experiments:\nUtilize the findings from [1], [2], and [3] to design additional experimental checks. Specifically, verify the efficacy of the new objects obtained from clustering, ensuring they contribute positively to the model's performance for a reduced sample. This validation will provide empirical evidence of the practical benefits derived from integrating insights from the referenced literature.\n\nDiversification of Citations:\nBroaden the range of references to include diverse research groups and papers. This not only addresses the issue of underrepresentation but also enriches the paper with a variety of perspectives and findings. It will demonstrate a more comprehensive engagement with the wider academic community and enhance the paper's credibility.", "broader_impact_concerns": "-", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6ZD2TJvnTd", "reviewer": "Reviewer_8iyN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- This paper is overall clearly written and easy to follow. \n- Compressive clustering is an important research topic, and this paper explains its background and motivation well.\n- This paper carefully analyzes the current approach (CL-OMPR) and proposes its improvement, which is a beneficial contribution to the community.\n\n**Weaknesses**\n- The empirical evaluation of this paper relies on only two datasets (one synthetic and one real-world). Given that the main contribution is tied to empirical results, it is necessary to include experiments conducted on a more diverse set of datasets.\n- Furthermore, each dataset used is low-dimensional, with dimensions being 6 for the synthetic dataset and 10 for the real-world dataset. Since higher-dimensional datasets are more and more common in practical situations, conducting experiments on high-dimensional data is particularly desirable.\n- While this paper explores the dependency of $L$, the limited diversity of datasets raises concerns about understanding how $L$ correlates with dataset properties. For example, if $L$ depends on the dimensionality and larger $L$ is needed for higher dimensional datasets, the effectiveness of the proposal is less convincing.\n- The paper lacks the measurement of the actual runtime for each tested method. Including the runtime information would be valuable in assessing the efficiency of these methods.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies compressive clustering, a task of clustering on sketches and decoding centroids from them. The focus is primarily on the standard heuristic CL-OMPR, revealing its limitations, particularly its failure to detect clearly separated clusters. To address this issue, the paper not only identifies the constraints of CL-OMPR but also introduces a novel algorithm with several improvements. Notably, two approaches are proposed to detect local maxima of the correlation function: a discretized approach and a sketched mean shift approach. The empirical evaluation, conducted on a synthetic dataset and the MNIST dataset, demonstrates the effectiveness of the sketched mean shift approach when compared to the original CL-OMPR.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please address the above listed weaknesses, that is, including a more diverse set of datasets, a more careful treatment of $L$, and reporting runtime.\n\n**Minor issues:**\n- P. 2, just after Eq.(4): \"ths\" -> \"this\"\n- Sometimes 'Lloyd' is spelled as 'LLoyd,' whereas 'Lloyd' should be consistently used. Please correct this inconsistency.", "broader_impact_concerns": "I do not have any concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The reviewers mostly agreed that this was a clear improvement over the baseline method, CL-OMPR, and one reviewer found the paper very instructive in diagnosing the problems with the baseline method and deriving the new algorithm.  There were some requested changes by reviewers, and in the AE's opinion, the authors suitably addressed these requests, including running more experiments and adding more text about alternative approaches (such as coresets).  No reviewer raised issues about novelty nor correctness. The biggest criticisms were that (1) this is a *narrow* paper, designed only for a very specific setup of compressive clustering, and (2) that the method is limited to moderate dimensions.  In response, I think that (1) there is still a sufficient TMLR audience for this paper, and given the novelty, it meets the TMLR criteria, and (2) the dimension limitation is a drawback of this family of methods but there is still a clear use-case for these kinds of methods (as well as room for future improvement).
Decision: Accept (poster)