STORM: SKETCH TOWARD ONLINE RISK MINIMIZA-
Year: 2021

Abstract:
Empirical risk minimization is perhaps the most influential idea in statistical learning, with applications to nearly all scientific and technical domains in the form of regression and classification models. The growing concerns about the high energy cost of training and the increased prevalence of massive streaming datasets have led many ML practitioners to look for approximate ERM models that can achieve low cost on memory and latency for training. To this end, we propose STORM, an online sketching-based method for empirical risk minimization. STORM compresses a data stream into a tiny array of integer counters. This sketch is sufficient to estimate a variety of surrogate losses over the original dataset. We provide rigorous theoretical analysis and show that STORM can estimate a carefully chosen surrogate loss for regularized least-squares regression and a margin loss for classification. We perform an exhaustive experimental comparison for regression and classification training on real-world datasets, achieving an approximate solution with size even less than a data sample.

heading: 1 INTRODUCTION
text: Empirical risk minimization (ERM) has been independently studied by the machine learning, statistics, and computer science communities. However, conventional ERM algorithms are becoming more expensive in terms of latency, memory, and energy due to the need for accurate inference over rapidly increasing data volumes (He et al., 2016; Vaswani et al., 2017). The net effect is an increase in storage costs of large data matrices (Guhaniyogi & Dunson, 2015), network latency for personalized inference/recommendations, and carbon emissions. There is a growing need for approximation algorithms that can minimize the cost, memory, and energy associated with training and inference (Lin et al., 2020; Motamedi et al., 2016). Additionally, there are distinct advantages of training on a small memory local/personalized device, as data is generally collected on these devices. This limits the data transfer to big servers that consumes substantial resources (Bill Allcock, 2001; Xu, 2013) and minimizes the privacy vulnerabilities (Mendez Mena et al., 2018). Existing compression methods generally fall into two categories: sketching and sampling.
Linear Sketches: A large family of projection-based sketches uses ideas from compressed sensing, matrix factorization, and dimensionality reduction. A common theme is to approximate the large N × d data matrix with a smaller representation. For example, multiplication with a random ±1 matrix yields a sketch that attains the information-theoretic lower bound for compressed linear regression (within a constant factor) (Clarkson & Woodruff, 2009). Random Gaussian, Hadamard, Haar, and Bernoulli sign matrices have also been used in this framework with varying degrees of theoretical and practical success (Dobriban & Liu, 2019). The “sketch-and-solve” technique has also been extended to regression with the `1 objective and several other linear problems (Sohler & Woodruff, 2011). This framework cannot accommodate regularized objectives because the model is computed directly from a set of precomputed sufficient statistics.
Recent work by Shi & Phillips (2021) has shown that a matrix sketch based on the frequent directions method (Liberty, 2013; Ghashami et al., 2016) can be used for ridge regression. This algorithm uses a sketch to approximate the most costly part of the ridge regression objective, XTX + λI, in a streaming setting. The memory is O(ld), where X ∈ Rn×d and l > 2. That is, the memory is linear in the dimensions, with a minimum size of (2d+ 1)× floating precision that can be very large for high-dimensional data.
Coresets and Sampling: A reasonable strategy is to extract a core subset (or coreset) of points that efficiently represent the complete dataset. The coreset framework supports the approximation of loss function sums within a (1± ) error margin. There are coresets for important problems such as `p-regression (Dasgupta et al., 2009) and logistic regression (Munteanu et al., 2018). Classification objectives also admit small coresets (Har-Peled et al., 2007; Tukan et al., 2020), as do many other regularized ERM problems (Reddi et al., 2015). The problem is that many coresets require superlinear construction time. For example, many of the algorithms in the previous paragraph require k-means clustering, leverage scores, or multiple cycles through the dataset (Tukan et al., 2020; Munteanu et al., 2018; Har-Peled et al., 2007). Coresets are also tricky to merge since it is not straightforward to decide which points to discard. As a result, composable and streaming coresets are an active area of research (Indyk et al., 2014; Karnin & Liberty, 2019). Although leverage scores can be approximated online (Cohen et al., 2016), the process is computationally expensive and not well suited to edge devices. In many practical problems, random sampling is the only feasible sampling baseline that supports regularized ERM.
In this work, we propose a new sketching framework (STORM) that supports both regularization and streaming algorithms. Given an N -point dataset D = {(x1, y1)...(xN , yN )} of d-dimensional examples observed in the streaming setting (Fiat, 1998), the goal is to construct a small structure S that can estimate a regression and classification loss over D. STORM consists of two main components, 1) STORM loss for regression and classification, 2) A tiny array of integer count values called STORM sketch. This array is sufficient to perform empirical risk minimization on STORM loss. In practice, our data structure allows for optimizations that are not possible with other methods. Specifically, our contributions are as follows:
1. We propose an online sketching algorithm (STORM) that can estimate a special class of loss functions using only integer count values. Our sketch is appropriate for distributed data settings because it is small, embarrassingly parallel, and mergeable via addition.
2. We characterize the set of STORM-approximable loss functions. The price of efficient loss estimation is a restriction on the kinds of losses we can use. While we cannot directly estimate popular regression and classification losses, we construct surrogate losses for `2 regression and large margin classification.
3. We show how to optimize models over STORM sketches. The count-based nature of our sketch makes gradient descent difficult, and black-box optimization techniques are of limited use in high-dimensional spaces. We propose an optimization algorithm that is custom-tailored to our sketch.
4. We conduct a rigorous empirical evaluation with ridge regression and classification. We show that our approach has better train/test loss performance and lower variance than sampling baselines in the low-memory, resource-constrained computational regime. For some datasets, we can achieve strong performance with memory smaller than even one sample.

heading: 2 BACKGROUND

heading: 2.1 EMPIRICAL RISK MINIMIZATION
text: In the standard statistical learning framework, we are given a training dataset D = {(xi, yi) ∈ D × Y}ni=1 of examples and asked to select a function h : X → Y that can predict y given x. In this paper, we will present our algorithms for the data space X = Rd and output space Y = R, but most of our results hold generally for other metric spaces too. The learning problem is to select a hypothesis hθ (a function parameterized by θ) that yields good predictions, as measured by a loss function loss(hθ(x), y). In empirical risk minimization (ERM), we select θ to minimize the average loss on the training set.
θ̂ = arg min θ n∑ i=1 loss(hθ(xi), yi)
In the case of Linear regression, hθ(x) = 〈θ,x〉. For most applications, θ is found using the least-squares or `2 loss: loss(hθ(x)) = ||hθ(x) − y||22. The unconstrained `2 loss is smooth and
strongly convex, with favorable convergence criteria. The parameter θ is typically found using gradient descent or via a closed-form solution from the matrix formulation of the model y = Xθ. For our discussion, it will be important to express the loss in terms of the concatenated vector [xi, yi]:
θ̂ = arg min θ n∑ i=1 (〈[xi, yi], [θ,−1]〉)2 = arg min θ ‖y− Xθ‖22

heading: 2.2 LOCALITY SENSITIVE HASHING
text: A Locality sensitive hash (LSH) function (l(x)) provides the same hash value for similar points with high probability. The notion of similarity is usually based on the distance measure of the metric space X . For instance, there are LSH families for the Jaccard (Broder, 1997), Euclidean (Datar et al., 2004; Dasgupta et al., 2011) and angular distances (Charikar, 2002). If we allow asymmetric hash constructions (i.e. l1(x1) = l2(x2), l1 6= l2), then x1 and x2 may collide based on other properties such as their inner product (Shrivastava & Li, 2014). Asymmetric LSH functions use one hash function (l1) for data and a different function (l2) for queries. To accommodate asymmetric functions, we use the following definition of LSH, Our definition is a strict generalization of the original definition from Indyk (1998) Indyk & Motwani (1998), which can be recovered by setting probability thresholds.
Definition 1. We say that a hash family F is locality-sensitive with collision probability k(·, ·) if for any two points x1 and x2 in Rd, l(x1) = l(x2) with probability k(x1,x2) under a uniform random selection of l(·) ∈ F .
The signed random projection (SRP) family for the angular distance (Goemans & Williamson, 1995; Charikar, 2002) is an example of a symmetric LSH family. The SRP family is the set of functions l(x) = sign(w>x), where w ∼ N (0, Id). The SRP collision probability is
k(x1,x2) = 1− 1
π cos−1 ( 〈x1,x2〉 ‖x1‖2‖x2‖2 ) The inner product hash Shrivastava & Li (2014) is an example of an asymmetric LSH. To get this hash function, replace x1 with [ x1, 0, √ 1− ‖x1‖22 ] and x2 with [ x2, √ 1− ‖x2‖22, 0 ]
in the SRP function. This procedure essentially uses different hash functions for x1 and x2, since we augment x1 (the data) in one way and x2 (the query) in a different way. The collision probability k(x1,x2) is now a monotone function of the inner product 〈x1,x2〉. This construction implicitly assumes that x1 and x2 are inside the unit sphere, but other hash functions do not require this assumption (Shrivastava & Li, 2015; Neyshabur & Srebro, 2015; Huang et al., 2018).1
RACE Sketch: The RACE sketch is a 2D integer array S ∈ ZR×B that we index using LSH functions (Coleman & Shrivastava, 2020; Luo & Shrivastava, 2018). Each row of S is called an Array-of-Counts (ACE) data structure (Luo & Shrivastava, 2018), and the entire sketch is a set of R independent ACE repetitions (RACE). To construct a RACE sketch from a dataset D, we begin with an empty R×B integer array and a set of R independent LSH functions {l1(x)...lR(x)}. Each LSH function has the collision probability k(θ,x). When an element x arrives from the data stream, we hash x to get R different hash values. We increment row i at location li(x) and repeat this process for all elements in the dataset.
To query the sketch with a query θ, we return the average value of S[i, li(θ)] over the R rows. RACE sketches are interesting because they estimate the sum of LSH collision probabilities over the dataset with low variance for any query θ.
Theorem 1. Unbiased RACE EstimatorColeman & Shrivastava (2020) Let X = S[i, li(θ)]. Then
E[X] = ∑ x∈D k(x, θ) var(X) ≤ (∑ x∈D √ k(x, θ) )2 1Regardless, we often scale and center the dataset when using inner product hash functions in practice.

heading: 3 STORM
text: While collision probabilities have been well-studied in the context of improving near-neighbor search (Gionis et al., 1999), we study them from a new perspective. Our goal is to compose useful loss functions from LSH collision probabilities and perform ERM on sketches. In this context, we query the sketch with the parameter θ to estimate the empirical risk. This requires some extensions to the RACE sketch.
Intuition To understand our approach, consider the task presented in Figure 1. We partition the data space into B random regions using Gaussian random hyperplanes. The count of points in each region gives an estimate of its density. We repeat this partitioning many times so that we have a good estimate of densely-populated and sparsely-populated regions. In the regression settings, the optimal model hyperplane passes through the dense regions. Our objective is to find a θ vector normal to the optimal model hyperplane. This motivates to look for the sparsest region. STORM sketch (inspired from RACE) provides counts of data points in each partition. We will see that we can optimize over these counts to find a θ that resides in the sparsest region. The main challenge is to design an LSH function l(θ) that has the minimal count when θ is a good model. Next, we will talk about an appropriate hash function, STORM loss for regression and classification and, the STORM sketch for approximating the loss.

heading: 3.1 STORM SURROGATE LOSS FOR LINEAR REGRESSION
text: Since we query the STORM sketch with θ, we must be able to express loss(hθ(x), y) as loss(θ, [x, y]). This condition is not limiting, but it does restrict our attention to hypothesis classes where θ directly interacts with x and y in a way that can be captured by an LSH function. Linear regression is the simplest example of such a model. Here the θ interacts with the data via the inner product 〈[θ,−1], [x, y]〉. The asymmetric hash function from Section 2.2 seems like a natural choice. However, the linear regression loss function is a monotone function of the absolute value of the inner product.
To obtain a surrogate loss with the correct dependence on 〈[θ,−1], [x, y]〉, observe that the collision probability of SRP (Section 2.2) is monotone decreasing if we hash [−x,−y] instead of [x, y]. Thus, we can obtain a function that is monotone in |〈[θ,−1], [x, y]〉| by adding together the collision probability for [x, y] (which is monotone increasing) and for [−x,−y] (which is monotone decreasing). We call this LSH family “paired random projections” (PRP) because the method can be implemented by hashing both [x, y] and [−x,−y] with a given SRP function and updating both locations in the sketch. When we query the vector θ̃ = [θ,−1] with the same SRP function, STORM estimates the following collision probability/ surrogate loss for linear regression
g(θ̃, [x, y]) = 1
2
( 1− 1
π cos−1(〈θ̃, [x, y]〉)
)p + 1
2
( 1− 1
π cos−1(−〈θ̃, [x, y]〉)
)p ,
where we are free to choose any number of random hyperplanes - p ≥ 2. The proposed loss is convex, unbiased and consistent. More info in section 4.

heading: 3.2 STORM SURROGATE LOSS FOR CLASSIFICATION
text: Consider linear hyperplane classifiers of the form hθ(x) = sign(〈θ, [x,−1]〉). The most popular losses used to find θ are classification-calibrated margin losses. A loss function is classificationcalibrated, or Bayes consistent, if the optimal hypothesis under the loss is the same as the Bayes optimal hypothesis. With these factors in mind, we propose Classification Random Projections (CRP) to estimate the following loss function:
Theorem 2. Consider labels y ∈ {−1, 1}. The loss function g(θ, [x, y]) for the linear hyperplane classifier is a classification-calibrated margin loss and is STORM-approximable.
g(θ, [x, y]) = 2p (
1− 1 π
cos−1(−y〈θ,x〉) )p
By optimizing θ in a similar fashion to as with regression, we can train linear classifiers with STORM sketches. The LSH function that implements g(θ̃, [x, y]) is the asymmetric inner product (SRP) hash, but with the argument to the hash function multiplied by y ∈ {−1,+1}. See the appendix for details and proof of theorems.

heading: 3.3 STORM SKETCH
text: We use asymmetric LSH functions to hash data points x and θ with different functions. Most LSH functions have positive definite collision probabilities, but asymmetric LSH functions allow much more flexibility. Our Paired Random Projections (PRP) hashing scheme premultiply the data vector by a {±1} to estimate a surrogate regression loss. Classification Random Projections (CRP) LSH premultiplies x with y. Figure 2 shows our sketch process for PRP. Here we insert elements into the sketch with two SRP hashes, but only use one SRP for the query. In Section 3.1 and Section 3.2, we show that this sketch estimates a surrogate loss for linear regression and classification.
Once we have a sketch and an appropriate hash function, we want to optimize the model parameter θ to minimize the STORM estimate of the empirical risk. Proposed optimization method is discussed in details in the next section.

heading: 3.4 OPTIMIZATION
text: Gradient descent is traditionally used to find the best θ for ERM problems. Unfortunately, our sketch consists of discrete integer counts that have neither an analytic gradient nor a subgradient. However, following the intuition from Figure 1, we know that our surrogate loss tries to find the hyperparameter θ̃ in the sparsest region. This allows us to develop a specialized optimization heuristic that exploits the hyperplanes of the proposed hash to locate the minimum-risk θ more effectively. Our surrogate objective is: θ̂ = arg minθ g([θ,−1], [x, y]) = arg minθ E(S[L([θ,−1])]).
Algorithm 2 Hyperplane-based stochastic gradient descent Input: STORM sketch S ∈ ZR×B Result: Model θ Initialise θ̃0 ← 0d+1 for i = 1..R do Zi = 2∗ binary(argmin S[i])−1
end for Aggregate all labels Z ′is to Z and store it. for iterations t = 1...T do
random sample i Get random hyper-planes Wi from the given seed. Update θ̃jt = θ̃ j t−1 − η∇‖tanh(W>i θ̃)− Zi‖22−2λθ̃
Project last dimension of θ̃ to −1 end for Select best θ̃j using a validation set.
Figure 3: Left: Shows the algorithm for optimization on the extracted labels of least count bins. Right: Figure that shows the further compression on the PRP sketch (works the same for CRP) by storing just the labels (Z) of least count buckets.
where S is the sketch and L = {l1..lp} is the hash function. This implies that θ lies in a region where most of the STORM partitions have lower count values (small values of S[l(θ̃)]). One could do a brute-force minimization over all possible partitions and select a model θ̃ in the winning partition, but this is computationally expensive. Instead, we use the stochastic gradient descent approach and transform our objective as follows: ĝ(θ) = ∑R i=1 S[i, Li([θ,−1])]
The gradient is∇ĝ(θ) = ∑R i=1∇S[i, Li([θ,−1])].
It turns out that the individual minimization problem (for each ith sample), arg minθ S[i, Li([θ,−1])] can also be written as arg minθ‖Li([θ,−1])− Zi‖22. Here Zi ∈ {0, 1}p represents the binary label vector of the minimum count bucket in row i.
Because we use SRP as the hash function, Li(θ̃) = sign(W>i θ̃). We now have the objective function in terms of 12 |sign(W > i θ̃) − Zi‖22. This can be made diffrentiable by using tanh instead of sign function. Hence, the gradient update becomes
θ̃t+1 = θ̃t − η∇‖tanh(W>i θ̃)− Zi‖22
This becomes the SGD over a standard single layer neural net with tanh activation. (Refer to appendix for more details). Also, we can simply add any regularization term with this. In case of Ridge regression we have the following gradient update
θ̃t+1 = θ̃t − η∇‖tanh(W>i θ̃)− Zi‖22−2λθ̃
It is interesting to note that, we are learning over randomly generated hyperplanes (Wi) (these only need a single seed value for generation) and very cheap binary vectors Z ∈ {0, 1}R×p. The sketch size required is Rp or R logB. Although this method (Algorithm 2) finds approximate gradients, we find that the procedure works well in practice. Since we can regenerate the hyperplanes from random seeds, we only need to store the hash codes (Z) of the least populated buckets with this approach.

heading: 4 THEORY: SKETCHABLE SURROGATE LOSS
text: We begin with a formal description of the families of losses that STORM can approximate.
Formally, a loss g(θ, [x, y]) is STORM-approximable if there is a hash function and a STORM sketch such that querying the STORM sketch with θ yields an unbiased estimator of g(θ, [x, y]).
Using the compositional properties of LSH, we can approximate any sum and/or product of LSH collision probabilities using one or more STORM sketches. For example, we estimate sums by
inserting elements with more than one LSH function (or using a second sketch). We can estimate a product of collision probabilities by combining multiple LSH outputs into a single hash code. Theorem 3 formalizes this argument (we defer all proofs to the appendix). Theorem 3. The set of STORM-approximable functions SL contains all LSH collision probabilities and is closed under addition, subtraction, and multiplication.
Given the number of known LSH families, the flexibility of asymmetric LSH and the closure of SL under addition and multiplication, Theorem 3 suggests that SL is an expressive space of functions. It remains for us to show that SL contains useful losses for machine learning problems. Theorem 4. When p ≥ 2, the PRP collision probability g([θ,−1], [x, y])is a convex surrogate loss for the linear regression objective such that
arg min θ g([θ,−1], [x, y]) = arg min θ ‖y − xθ‖
Check Figure 4 (a). A detail proof is given in the Appendix. We analyse, how the sketch size affects the PRP/CRP loss estimation from the STORM sketch. Consider the following Axiom Let Z1, ...ZR be R i.i.d. random variables with mean E[Z] = µ and variance ≤ σ2. Divide the R variables into g groups so that each group contains m = R/g elements, and take the empirical average within each group. The median-of-means estimate µ̂ is the median of the g group means. If g = 8 log(1/δ) and m = R/g, then the following statement holds with probability 1− δ.
|µ̂− µ|≤ 6 σ√ R
√ log 1/δ
This proof is given in Alon et al. (1999) as the proof of Theorem 2.1 (which is a slightly more general version of the statement above). From the variance result from Theorem 1, and a median of means estimate of the PRP/CRP loss from STORM sketch, we state the followingTheorem 5. Let ĝ(θ) is the median of means estimate on the STORM sketch, with R rows, then with probability more than 1− δ
|ĝ(θ)− g(θ)|≤ 6 ∑ x∈D √ k(x, θ)
√ R
√ log(1/δ)
Where ∑ x∈D √ k(x, θ) is the collision probability of the hash function in use.

heading: 4.1 PARAMETERS
text: Number of hash functions (p): The parameter p is the number of PRP/CRP hash hyperplanes that creates 2p partitions. Larger p is good for more resolution of the space, but it also leads to many buckets with few/zero counts. (Figure 4(a) shows that higher p makes the loss surface flatter near the optimal region). Figure 4 (b) shows that setting p ≈ 4 maximizes the gradient of the surrogate loss near the minima, and we find that in practice too (Refer Appendix). It is also efficient because for p = 4, B is only 16. Also, we need fewer rows in the sketch to optimize g(θ̃, [x, y]) when it is strongly convex.
Number of repetitions (R): The square root of collision probability √ k(x, θ) ≤ 1 always. Hence,
in the worst case the value of ∑ x∈D √ k(x, θ) is N . For |ĝ(θ) − g(θ)|≤ C √ log(1/δ), we need number of sketch rowsR < O(N2) (Using the Theorem 5). The worst-caseO(N2) situation happens when every point in the dataset has maximum loss. Practical modeling situations typically have a small fraction of points for which the training loss is small. In experiments, we vary the number of rows and see the change in Mean Square Error (MSE).

heading: 5 EXPERIMENTS
text: Experiment setup: We use PRP and CRP with p = 4 (most optimal, check Figure 4(b)) to create the sketch and we vary the number of repetitions to show the memory tradeoff. We have used a variety of datasets. Refer table 1
All models are trained and evaluated with a fixed 80/20 train/test split. We average over 100 runs for our baselines and sketches, where each run has an independently constructed sketch or random sample. Thus, our average is over the random LSH functions used to construct the sketch and the stochastic gradient descent instances. We use the hyperplane optimization heuristic for regression and derivative-free (Conn et al., 2009) gradient descent for classification. To run the code we used our personal laptop with i7-8750H processor, 2.2 GHz base frequency and 9MB cache size.
Ridge Regression: We evaluate our method on ridge regression, which penalizes θ to avoid the double descent and over-fitting problem (Nakkiran, 2019). To implement this model, we simply add ‖θ‖22 to the STORM objective function. The gradient updates take place as described in Section 3.3, but with an added regularization term. For regularized ERM in our computational setting, random sampling and Frequent direction sketch is the only reasonable baseline. Coresets for ridge regression exist, but cannot be constructed in the streaming setting. Also, Frequent direction based sketch takes space much more than STORM and random samples for higher-dimensional data (near d > 100). Hence a direct comparison is made for the low dimensional- Diabetes and Boston housing dataset. Figure 5 shows the results of our comparison focused mostly on sub-kB region. We find that STORM achieves a lower MSE in the low-memory regime. In particular, for the gas, soil, fruit, and indoor datasets, STORM produces a reasonable model using less memory than even a single sample. Finally, the STORM results have a substantially smaller variance (typically 2x), which is a critical factor in practice. Additional plots and sanity check on a 2D synthetic data is available in the appendix.
STORM for Non-linear Objectives: STORM supports regularized linear regression and classification. With this formulation, a non-linear model can be trained with feature mapping such as random Fourier features (Rahimi & Recht, 2008). To achieve a non-linear boundary the mapping is done to a higher dimension space. Practically this is 3− 6 times the original dimensions. Since STORM does well with high-dimensional datasets, often requiring less memory than a single sample, we expect good performance using high–dimensional feature transformations. Figure 5 last row demonstrate the mean square error (MSE) with the size of sketch variables required for linear regression.
Linear classification: We run the linear classification on datasets from the libsvm repository. We compare against a linear SVM trained on a random sample of the dataset. Figure 5 shows this comparison for the “Australian credit” (Dua & Graff, 2017) and the popular MNIST (Lecun et al., 1998) dataset.

heading: 6 DISCUSSION
text: STORM is a scalable framework to solve regularized ERM problems. It is a memory-efficient framework that supports regularization and streaming, as long as the loss function can be described using LSH collision probabilities. The sketch provides compression ranging from 1/50 to 1/10000 of the original size and can take size effectively even less than a sample. Additionally, our approach is appealing because it performs well in high dimensions, which can be a prerequisite for feature transformations such as random Fourier features (Rahimi & Recht, 2008). We expect that STORM can be extended to work for important tasks such as kernel classifiers and kernel ridge regression. This will greatly expand the modeling options for applications like distributed learning via ERM on edge devices.

heading: A APPENDIX
text: A.1 PROOF OF THEOREMS
Theorem 1. Unbiased RACE Estimator Let X = S[i, li(θ)]. Then
E[X] = ∑ x∈D k(x, θ) var(X) ≤ (∑ x∈D √ k(x, θ) )2
Proof: The proof is detailed in Coleman(2019) Coleman et al. (2019) and Luo(2018) Luo & Shrivastava (2018).
Theorem 3. The set of STORM-approximable functions SL contains all LSH collision probabilities and is closed under addition, subtraction, and multiplication.
It is straightforward to see that STORM can approximate ∑ x∈D k(x, y) as long as there is an LSH function with the collision probability k(x, y). To prove the theorem, it is sufficient to show that given two LSH collision probabilities k1(x, y) and k2(x, y), STORM sketches can approximate the following two functions
f1(y) = ∑ x∈D k1(x, y)± k2(x, y) f2(y) = ∑ x∈D k1(x, y)k2(x, y)
Note that one can always write a product of (weighted) sums ( ∑ n wnkn(x, y)) ( ∑ m wmkm(x, y))
as the sum of (weighted) products ∑ n,m wnwmkn(x, y)km(x, y). Therefore, the previous two situations ensure that the set is closed under addition, subtraction and multiplication.
Addition and Subtraction: Because of the distributive property of addition,∑ x∈D k1(x, y)± k2(x, y) = ∑ x∈D k1(x, y)± ∑ x∈D k2(x, y)
One can then construct a STORM sketch S1 for the k1(x, y) summation and a second STORM sketch S2 for k2(x, y) summation. We can estimate any linear combination of k1(x, y) and k2(x, y) by with a weighted sum of S1 and S2. Multiplication: To approximate sums over the product k1(x, y)k2(x, y), we rely on LSH hash function compositions. Suppose we have an LSH function l1(x) with collision probability k1(x, y) and l2(x) with k2(x, y). Consider the hash function l(x) = π(l1(x), l2(x)) where π(a,b) is an injective (or unique) mapping from Z2 → Z. An example of such a mapping is the function π(a,b) = pa1p b 2 where p1 and p2 are coprime. Since the mapping is injective, this means that l(x) = l(y) only when l1(x) = l1(y) and l2(x) = l2(y). Therefore, Pr[l(x) = l(y)] = Pr[l1(x) = l1(y) ∩ l2(x) = l2(y)] Make the choice of the LSH functions l1 and l2 independently, so that the probability factorizes
= Pr[l1(x) = l1(y)]Pr[l2(x) = l2(y)] = k1(x, y)k2(x, y)
Therefore, one can construct a STORM sketch for the product using the LSH function l(x) = π(l1(x), l2(x)).
Theorem 4. When p ≥ 2, the PRP collision probability g([θ,−1], [x, y])is a convex surrogate loss for the linear regression objective such that
arg min θ g([θ,−1], [x, y]) = arg min θ ‖y − xθ‖
Proof: For the surrogate ERM problem to have the same solution as the linear regression ERM problem, it is sufficient to show two things: that the surrogate loss is convex and that the global minima of the surrogate loss and the linear regresssion loss appear in the same location. The surrogate loss (L(θ,x, y)) is
g([θ,−1], [x, y]) = 1 2
( 1− 1
π cos−1(〈[θ,−1], [x, y]〉)
)p + 1
2
( 1− 1
π cos−1(−〈[θ,−1], [x, y]〉)
)p
For the sake of notation, we will put a = [θ,−1], b = [x, y], and
f(a,b) = ( 1− 1
π cos−1(〈[θ,−1], [x, y]〉)
) L(a,b) = 1
2 f(a,b)p +
1 2 f(a,−b)p
We will use the fact that
∇af(a,b) = ∇af(a,−b) = b π √ 1− |〈a,b〉|2
Location of Minima: The minimum of the surrogate loss is same as the minimum for least squares linear regression. Using the chain rule, the gradient of the surrogate loss is
∇aL(a,b) = p ( f(a,b)p−1 − f(a,−b)p−1 ) 2π √ 1− |〈a,b〉|2 b
When p = 1, the gradient is always zero. When p ≥ 2, the derivative is zero when 〈a,b〉 = 〈[θ,−1], [x, y]〉 = 0 because that is where f(a,b) = f(a,−b). Thus, the surrogate loss has the same minimizer as the least squares loss. Convexity: At index (i, j), the Hessian of the surrogate loss is
[∇2aL(a,b)](i,j) = ∂
∂aj [∇aL(a,b)]i =
∂ ∂aj bi p ( f(a,b)p−1 − f(a,−b)p−1 ) 2π √ 1− |〈a,b〉|2
Thus, the Hessian
∇2aL(a,b) = b>∇a
( pf(a,b)p−1 − pf(a,−b)p−1
2π √ 1− |〈a,b〉|2
)
The gradient
∇a
( pf(a,b)p−1 − pf(a,−b)p−1
2π √ 1− |〈a,b〉|2
)
=
( p(p− 1)f(a,b)p−2
2π √ 1− |〈a,b〉|2 ∇af(a,b) + p(p− 1)f(a,−b)p−2 2π √ 1− |〈a,b〉|2 ∇af(a,−b)
)
+ ( pf(a,b)p−1 − pf(a,−b)p−1 ) ∇a
( 1
2π √ 1− |〈a,b〉|2 ) Simplifying, we obtain
bTb
2π
( p(p− 1)f(a,b)p−2 − p(p− 1)f(a,−b)p−2
π(1− |〈a,b〉|2) + 〈a,b〉pf(a,b)
p−1 + pf(a,−b)p−1
(1− |〈a,b〉|2)3/2
)
Which gives the following expression for the Hessian
∇2θL(a,b) = b Tbρ(a,b)
where
ρ(a,b) = 1
2π
( p(p− 1)f(a,b)p−2 − p(p− 1)f(a,−b)p−2
π(1− |〈a,b〉|2) + 〈a,b〉pf(a,b)
p−1 + pf(a,−b)p−1
(1− |〈a,b〉|2)3/2 ) for p ≥ 2 ρ(a,b) > 0 ∀〈a,b〉 ∈ (−1, 1). Hence, the Hessian is positive semidefinite and the function is convex in a. Also note that the function is convex in θ, since the restriction of the last dimension of a to −1 is the restriction of a convex function to a convex set. Refer Figure 6(a). Theorem 2. Consider labels y ∈ {−1, 1}. The loss function g(θ, [x, y]) for the linear hyperplane classifier is a classification-calibrated margin loss and is STORM-approximable.
g(θ, [x, y]) = 2p (
1− 1 π
cos−1(−y〈θ,x〉) )p
Proof. First, we show that the loss is classification-calibrated. Then, we show that the loss can be estimated using STORM.
Loss is Classification-Calibrated: A necessary and sufficient condition for a convex2 loss function φ(t) to be classification-calibrated is for ddtφ(t) < 0 at t = 0. Here, t = yh(x), where h(x) is the model. For a linear hyperplane classifier, h(x) = 〈θ,x〉. The loss is therefore
φ(t) = 2p (
1− 1 π
cos−1(−t) )p
φ(t) is convex when p ≥ 2 for the same reasons discussed in the proof of Theorem 4. Note that the simple asymmetric LSH for the inner product that we have used throughout the paper3 requires t ∈ [−1, 1]. The derivative is
d dt φ(t) = 2pp
( 1− 1
π cos−1(−t) )p−1 −1 π √
1− t2 2For non-convex φ(t), the sufficient conditions are more complicated. 3There are other asymmetric inner product LSH functions without this requirement, and in practice one
usually scales the data to make the condition true.
At the origin, ddtφ(t) is −1/π. Therefore, the loss is classification calibrated. Figure 7(b) compares our STORM surrogate classification loss against popular margin losses.
Loss is STORM-Approximable: Consider the asymmetric LSH function for the inner product where we premultiply x by −y. The collision probability under this LSH function is
k(θ,x) = ( 1− 1
π cos−1(〈θ,−yx〉) ) as desired.
Theorem 5. Let ĝ(θ) is the median of means estimate on the STORM sketch, with R rows, then with probability more than 1− δ
|ĝ(θ)− g(θ)|≤ 6 ∑ x∈D √ k(x, θ)
√ R
√ log(1/δ)
Proof. Consider the following axiom-
Let Z1, ...ZR be R i.i.d. random variables with mean E[Z] = µ and variance ≤ σ2. Divide the R variables into g groups so that each group contains m = R/g elements, and take the empirical average within each group. The median-of-means estimate µ̂ is the median of the g group means. If g = 8 log(1/δ) and m = R/g, then the following statement holds with probability 1− δ.
|µ̂− µ|≤ 6 σ√ R
√ log 1/δ
This proof is given in Alon et al. (1999) as the proof of Theorem 2.1 (which is a slightly more general version of the statement above).
Now, et Z(q) be the median-of-means estimate constructed using the R unbiased estimators of the sketch with R rows. Then with probability 1− δ,
|ĝ(θ)− g(θ)|≤ 6 ∑ x∈D √ k(x, θ)
√ R
√ log(1/δ)
From Theorem 1, we know that σ ≤ ∑ x∈D √ k(x, θ) Substituting this variance bound into the given Lemma proves the theorem. Required number of repetitions: The square root of collision probability √ k(x, θ) ≤ 1 always.
Hence, in the worst case the value of ∑ x∈D √ k(x, θ) is N . For |ĝ(θ)− g(θ)|≤ C √ log(1/δ), we need number of sketch rows R = O(N2).
SGD convergence and sketch size for hyperplane based method: Assuming that 1) ĝ(θ) is bounded below, 2) The gradient is L-Lipschiz continuous and 3) the gradients bounded as E [ ‖∇ĝi(θ)‖2 ] ≤ σ2, with constant step size η we will get the error (ĝ(θt) − ĝ(θ∗)) bounded
by O(1/t) +O(η). For a simplistic assumption that we are drawing a new sample every time the total number of steps is equivalent to the number of repetitions R in the sketch. Hence the error at the end of the iteration will be bounded in terms of sketch size by O(1/R) +O(η).
A.2 EXPERIMENTS
A.2.1 EXPERIMENT ON 2D SYNTHETIC DATASET
We plot the estimated parameter θ̂ from the STORM on 2D synthetic data (Train: 1600 points, Test: 400 points) generated from Gaussian distribution. Refer to Figure 8.
A.2.2 NON LINEAR REGRESSION ON SYNTHETIC DATA:
A.2.3 COMPARISON WITH ZERO ORDER GRADIENT DECENT
We perform derivative-free optimization Conn et al. (2009) here on the original STORM sketch. Here, the black-box access to the loss function of interest (or its sharp approximation) is sufficient. For the descent algorithm, we initialize θ̃ to zeros in d+ 1 dimensions. The additional dimension
is due to querying the sketch with [θ,−1] rather than [θ]. We compute the approximate gradient by querying the sketch with equidistant points in a ball around θ̃ and update the parameter. After each iteration, we project the last dimension of θ̃ back onto the constraint θ̃d+1 = −1. This takes R×B size, with an additional factor of sketch cell bit precision. Hence it is more memory intensive than the hyper-plane based stochastic gradient descent algorithm. Refer the figure below.
A.2.4 NON-REGULARISED LINEAR REGRESSION
In Figure 11, we observe a double-descent phenomenon for our sampling baselines, explaining the peak near the intrinsic dimensionality of the problem. This sample-wise double descent behavior was recently proved for linear regression by Nakkiran (2019). STORM does not experience the double descent curve in practice because the entire dataset (not just a subsample) is used to minimize the loss. We perform favorably against baselines in memory regimes affected by double descent and STORM performs competitively in other memory regimes. We also observe that the θ found using STORM converges to the optimal θ under least-squares ERM. The datsets in use are taken from UCI repository.

Review summaries:
This paper considers the online risk minimization setting. To reduce the memory cost and accelerate training, this paper proposes a new sketching algorithm (i.e. STORM) based on the locality sensitive hashing (LSH) technique, which mainly includes two steps -- sketching and training. In particular, it provides two STORM surrogate losses for linear regression and classification, respectively. Besides, it also provide theoretical analysis to justify the use of sketching to estimate the surrogate loss with high probability. Experimental results illustrate the effectiveness and efficiency of STORM over other baselines.

The paper considers approximately computing ERM problems for real data, with an emphasis on regression and classification, in the *very low space* regime (most experiments are below 1Kb). This is done my maintaining a standard locality sensitive hashing table, but with specialized hash functions that allow for computing a surrogate for common regression and classification ERM objectives.   Due to the extreme space constraint, optimization is done with a heuristically characterized SGD. A significant amount of experiments is presented.

In this paper, the authors focus on reducing the cost of conventional ERM training and propose a new framework, STORM, which integrates the approaches of locality sensitive hashing and sketches. Besides, the authors use specific optimization algorithms to solve that this framework can’t be directly applied to popular classification and regression tasks. The experimental results prove that the proposed procedure can achieve better results compared to the baselines in the operational environment of low performance and low memory.