Determinantal point processes based on orthogonal polynomials for sampling minibatches in SGD
Year: 2022

Abstract:
Stochastic gradient descent (SGD) is a cornerstone of machine learning. When the number N of data items is large, SGD relies on constructing an unbiased estimator of the gradient of the empirical risk using a small subset of the original dataset, called a minibatch. Default minibatch construction involves uniformly sampling a subset of the desired size, but alternatives have been explored for variance reduction. In particular, experimental evidence suggests drawing minibatches from determinantal point processes (DPPs), tractable distributions over minibatches that favour diversity among selected items. However, like in recent work on DPPs for coresets, providing a systematic and principled understanding of how and why DPPs help has been difficult. In this work, we contribute an orthogonal polynomial-based determinantal point process paradigm for performing minibatch sampling in SGD. Our approach leverages the specific data distribution at hand, which endows it with greater sensitivity and power over existing data-agnostic methods. We substantiate our method via a detailed theoretical analysis of its convergence properties, interweaving between the discrete data set and the underlying continuous domain. In particular, we show how specific DPPs and a string of controlled approximations can lead to gradient estimators with a variance that decays faster with the batchsize than under uniform sampling. Coupled with existing finite-time guarantees for SGD on convex objectives, this entails that, for a large enough batchsize and a fixed budget of item-level gradients to evaluate, DPP minibatches lead to a smaller bound on the mean square approximation error than uniform minibatches. Moreover, our estimators are amenable to a recent algorithm that directly samples linear statistics of DPPs (i.e., the gradient estimator) without sampling the underlying DPP (i.e., the minibatch), thereby reducing computational overhead. We provide detailed synthetic as well as real data experiments to substantiate our theoretical claims. ∗Alphabetical order †Corresponding author 35th Conference on Neural Information Processing Systems (NeurIPS 2021).

text: ∗Alphabetical order †Corresponding author
35th Conference on Neural Information Processing Systems (NeurIPS 2021).

heading: 1 Introduction
text: Consider minimizing an empirical loss
min θ∈Θ
1 N · N∑ i=1 L(zi, θ) + λ(θ), (1)
with some penalty λ : Θ 7→ R+. Many learning tasks, such as regression and classification, are usually framed that way [1]. When N 1, computing the gradient of the objective in (1) becomes a bottleneck, even if individual gradients ∇θL(zi, θ) are cheap to evaluate. For a fixed computational budget, it is thus tempting to replace vanilla gradient descent by more iterations but using an approximate gradient, obtained using only a few data points. Stochastic gradient descent (SGD; [2]) follows this template. In its simplest form, SGD builds an unbiased estimator at each iteration of gradient descent, independently from past iterations, using a minibatch of random samples from the data set. Theory [3] and practice suggest that the variance of the gradient estimators in SGD should be kept as small as possible. It is thus natural that variance reduction for SGD has been a rich area of research; see for instance the detailed references in [4, Section 2].
In a related vein, determinantal point processes (DPPs) are probability distributions over subsets of a (typically large or infinite) ground set that are known to yield samples made of collectively diverse items, while being tractable both in terms of sampling and inference. Originally introduced in electronic optics [5], they have been turned into generic statistical models for repulsion in spatial statistics [6] and machine learning [7, 8]. In ML, DPPs have also been shown to be efficient sampling tools; see Section 2. Importantly for us, there is experimental evidence that minibatches in (1) drawn from DPPs and other repulsive point processes can yield gradient estimators with low variance for advanced learning tasks [4, 9], though a conclusive theoretical result has remained elusive. In particular, it is hard to see the right candidate DPP when, unlike linear regression, the objective function does not necessarily have a geometric interpretation.
Our contributions are as follows. We combine continuous DPPs based on orthogonal polynomials [10] and kernel density estimators built on the data to obtain two gradient estimators; see Section 3. We prove that their variance isOP (p−(1+1/d)), where p is the size of the minibatch, d is the dimension of data; see Section 4. This provides theoretical backing to the claim that DPPs yield variance reduction [4] over, say, uniformly sampling without replacement. In passing, the combination of analytic tools –orthogonal polynomials–, and an essentially discrete subsampling task –minibatch sampling– sheds light on new ways to build discrete DPPs for subsampling. Finally, we demonstrate our theoretical results on simulated data in Section 5.
A cornerstone of our approach is to utilise orthogonal polynomials to construct our sampling paradigm, interweaving between the discrete set of data points and the continuum in which the orthogonal polynomials reside. A few words are in order regarding the motivation for our choice of techniques. Roughly speaking, we would like to use a DPP that is tailored to the data distribution at hand. Orthogonal Polynomial Ensembles (OPEs) provide a natural way of associating a DPP to a given measure (in this case, the probability distribution of the data points), along with a substantive body of mathematical literature and tools that can be summoned as per necessity. This makes it a natural choice for our purposes.
Notation. Let data be denoted by D := {z1, . . . ,zN}, and assume that the zi’s are drawn i.i.d. from a distribution γ on Rd. Assume γ is compactly supported, with support D ⊂ [−1, 1]d bounded away from the border of [−1, 1]d. Assume also that γ is continuous with respect to the Lebesgue measure, and that its density is bounded away from zero on its support. While our assumptions exclude learning problems with discrete labels, such as classification, we later give experimental support that our estimators yield variance reduction in that case too. We define the empirical measure γ̂N := N −1 · ∑N i=1 δzi , where δzi is the delta measure at the point zi ∈ Rd. Clearly, γ̂N → γ in P(Rd); under our operating assumption of compact support, this amounts to convergence in P(D). For simplicity, we assume that no penalty is used in (1), but our results will extend straightforwardly. We denote the gradient of the empirical loss by ΞN = ΞN (θ) := N−1 · ∑N i=1∇θL(zi, θ). A minibatch is a (random) subset A ⊂ [N ] of size |A| = p N such that the random variable ΞA = ΞA(θ) :=
∑ i∈A wi∇θL(zi, θ), (2)
for suitable weights (wi), provides a good approximation for ΞN .

heading: 2 Background and relevant literature
text: Stochastic gradient descent. Going back to [2], SGD has been a major workhorse for machine learning; see e.g. [1, Chapter 14]. The basic version of the algorithm, applied to the empirical risk minimization (1), is to repeatedly update
θt+1 ← θt − ηtΞA(θt), t ∈ N, (3) where ηt is a (typically decreasing) stepsize, and ΞA is a minibatch-based estimator (2) of the gradient of the empirical risk function (1), possibly depending on past iterations. Most theoretical analyses assume that for any θ, ΞA(θ) in the t-th update (3) is unbiased, conditionally on the history of the Markov chain (θt) so far. For simplicity, we henceforth make the assumption that ΞA does not depend on the past of the chain. In particular, using such unbiased gradients, one can derive a nonasymptotic bound [3] on the mean square distance of θt to the minimizer of (1), for strongly convex and smooth loss functions like linear regression and logistic regression with `2 penalization. More precisely, for ηt ∝ t−α and 0 < α < 1, [3, Theorem 1] yields
E‖θt − θ?‖2 ≤ f(t) ( E‖θ0 − θ?‖2 + σ2
L
) + C σ2
tα , (4)
where C,L > 0 are problem-dependent constants, f(t) = O(e−tα), and σ2 = E[‖ΞA(θ?)‖2|D] is the trace of the covariance matrix of the gradient estimator, evaluated at the optimizer θ? of (1). The initialization bias is thus forgotten subexponentially fast, while the asymptotically leading term is proportional to σ2/tα. Combined with practical insight that variance reduction for gradients is key, theoretical results like (4) have motivated methodological research on efficient gradient estimators [4, Section 2], i.e., constructing minibatches so as to minimize σ2. In particular, repulsive point processes such as determinantal point processes have been empirically demonstrated to yield variance reduction and overall better performance on ML tasks [4, 9]. Our paper is a stab at a theoretical analysis to support these experimental claims.
The Poissonian benchmark. The default approach to sample a minibatch A ⊂ [N ] is to sample p data points from D uniformly, with or without replacement, and take wi = 1/p constant in (2). Both sampling with or without replacement lead to unbiased gradient estimators. A similar third approach is Poissonian random sampling. This simply consists in starting from A = ∅, and independently adding each element of D to the minibatch A with probability p/N . The Poisson estimator ΞA,Poi is then (2), with constant weights again equal to 1/p. When p N , which is the regime of SGD, the cardinality of A is tightly concentrated around p, and ΞA,Poi has the same fluctuations as the two default estimators, while being easier to analyze. In particular, E[ΞA,Poi|D] = ΞN and Var[ΞA,Poi|D] = OP (p−1); see Appendix S1 for details.
DPPs as (sub)sampling algorithms. As distributions over subsets of a large ground set that favour diversity, DPPs are intuitively good candidates at subsampling tasks, and one of their primary applications in ML has been as summary extractors [7]. Since then, DPPs or mixtures thereof have been used, e.g., to generate experimental designs for linear regression, leading to strong theoretical guarantees [11, 12, 13]; see also [14] for a survey of DPPs in randomized numerical algebra, or [15] for feature selection in linear regression with DPPs.
When the objective function of the task has less structure than linear regression, it has been more difficult to prove that finite DPPs significantly improve over i.i.d. sampling. For coreset construction, for instance, [16] manages to prove that a projection DPP necessarily improves over i.i.d. sampling with the same marginals [16, Corollary 3.7], but the authors stress the disappointing fact that current concentration results for strongly Rayleigh measures (such as DPPs) do not allow yet to prove that DPP coresets need significantly fewer points than their i.i.d. counterparts [16, Section 3.2]. Even closer to our motivation, DPPs for minibatch sampling have shown promising experimental performance [4], but reading between the lines of the proof of [4, Theorem 1], a bad choice of DPP can even yield a larger variance than i.i.d. sampling!
For such unstructured problems (compared to, say, linear regression) as coreset extraction or lossagnostic minibatch sampling, we propose to draw inspiration from work on continuous DPPs,
where faster-than-i.i.d. error decays have been proven in similar contexts. For instance, orthogonal polynomial theory motivated [10] to introduce a particular DPP, called a multivariate orthogonal polynomial ensemble, and prove a faster-than-i.i.d. central limit theorem for Monte Carlo integration of smooth functions. While resorting to continuous tools to study a discrete problem such as minibatch sampling may be unexpected, we shall see that the assumption that the size of the dataset is large compared to the ambient dimension crucially allows to transfer variance reduction arguments.

heading: 3 DPPs, and two gradient estimators
text: Since we shall need both discrete and continuous DPPs, we assume that X is either Rd or D, and follow [17] in introducing DPPs in an abstract way that encompasses both cases. After that, we propose two gradient estimators for SGD that build on a particular family of continuous DPPs introduced in [10] for Monte Carlo integration, called multivariate orthogonal polynomial ensembles.

heading: 3.1 DPPs: The kernel machine of point processes
text: A point process on X is a distribution on finite subsets S of X ; see [18] for a general reference. Given a reference measure µ on X , a point process is said to be determinantal (DPP) if there exists a function K : X × X → R, the kernel of the DPP, such that, for every n ≥ 1,
E [∑ 6= f(xi1 , . . . , xin) ] = ∫ (X )n f(x1, . . . , xn) ·Det [ K(xi, x`) ]n i,`=1 dµ⊗n(x1, . . . , xn) (5)
for every bounded Borel function f : Xn → R, where the sum in the LHS of (5) ranges over all pairwise distinct n-uplets of the random finite subset S. A few remarks are in order. First, satisfying (5) for every n is a strong constraint on K, so that not every kernel yields a DPP. A set of sufficient conditions on K is given by the Macchi-Soshnikov theorem [5, 19]. In words, if K(x, y) = K(y, x), and if K further corresponds to an integral operator
f 7→ ∫ K(x, y)f(y)dµ(y), f ∈ L2(X , µ),
that is trace-class with spectrum in [0, 1], then the corresponding DPP exists. Second, note in (5) that the kernel of a DPP encodes how the points in the random configurations interact. A strong point in favour of DPPs is that, unlike most interacting point processes, sampling and inference are tractable [17, 6]. Third, (5) yields simple formulas for the mean and variance of linear statistics of a DPP. Proposition 1 (See e.g. [20, 21]). Let S ∼ DPP(K,µ) and Φ : X → R be a bounded Borel function,
E [∑ x∈S Φ(x) ] = ∫ Φ(x)K(x, x)dµ(x), (6)
Var [∑ x∈S Φ(x) ] = ∫∫ ‖Φ(x)− Φ(y)‖22 |K(x, y)| 2dµ(x)dµ(y)
+ ∫ ‖Φ(x)‖22 ( K(x, x)− ∫ |K(x, y)|2dµ(y) ) dµ(x). (7)
Since the seminal paper [7], the case X = D has been the most common in machine learning. Taking µ to be γ̂N , any kernel is given by its restriction to D, usually given as an N × N matrix K = K|D. Equation (6) with n = p, A a subset of size p of D, and Φ the indicator of A yields P(A ⊂ S) = N−pDetKA. This is the usual way finite DPPs are introduced [7], except maybe for the factor N−p, which comes from using γ̂N as the reference measure instead of Nγ̂N . In this finite setting, a careful implementation of the general DPP sampler of [17] yields m DPP samples of average cardinality Tr(K) in O(Nω +mNTr(K)2) operations [22]. We now go back to a general X and fix p ∈ N. A canonical way to construct DPPs generating configurations of p points almost surely, i.e. S = {x1, . . . , xp}, is the following. Consider p orthonormal functions φ0, . . . , φp−1 in L2(µ), and take for kernel
K(x, y) = p−1∑ k=0 φk(x)φk(y). (8)
In this setting, the (permutation invariant) random variables x1, . . . , xp with joint distribution
1 p! Det
[ K(xi, x`) ]p i,`=1 p∏ i=1 dµ(xi) (9)
generate a DPP {x1, . . . , xp} with kernel K(x, y), called a projection DPP. For further information on determinantal point processes, we refer the reader to [17, 7].

heading: 3.2 Multivariate orthogonal polynomial ensembles
text: This section paraphrases [10] in their definition of a particular projection DPP on X = Rd, called a multivariate orthogonal polynomial ensemble (OPE). Fix some reference measure q(x)dx on [−1, 1]d, and assume that it puts positive mass on some open subset of [−1, 1]d. Now choose an ordering of the monomial functions (x1, . . . , xd) 7→ xα11 · · ·x αd d ; in this work we use the graded lexical order. Then apply the Gram-Schmidt algorithm in L2(q(x)dx) to these ordered monomials. This yields a sequence of orthonormal polynomial functions (φk)k∈N, the multivariate orthonormal polynomials w.r.t. q. Finally, plugging the first p multivariate orthonormal polynomials φ0, . . . , φp−1 into the projection kernel (8), we obtain a projection DPP with kernel denoted as K(p)q , referred to as the multivariate OPE associated with the measure q(x)dx.

heading: 3.3 Our first estimator: reweight, restrict, and saturate an OPE kernel
text: Let h > 0 and
γ̃(z) = 1
Nhd N∑ i=1 k ( z − zi h ) (10)
be a single-bandwidth kernel density estimator of the pdf of the data-generating distribution γ; see [23, Section 4.2], In particular, k is chosen so that ∫ k(x)dx = 1. Note that the approximation kernel k is unrelated to any DPP kernel in this paper. Let now q(x) = q1(x1) . . . qd(xd) be a separable pdf on [−1, 1]d, where each qi is Nevai-class3. Let K(p)q be the multivariate OPE kernel defined in Section 3.2, and form a new kernel
K (p) q,γ̃(x, y) :=
√ q(x)
γ̃(x) K(p)q (x, y)
√ q(y)
γ̃(y) .
The form of K(p)q,γ̃ is reminiscent of importance sampling [24], which is no accident. Indeed, while the positive semidefinite matrix
K (p) q,γ̃ |D :=
( K
(p) q,γ̃(zi, zj) ) 1≤i,j≤N
(11)
is not necessarily the kernel matrix of a DPP on ({1, . . . , N}, γ̂N ), see Section 3.1, we built it to be close to a projection of rank p. More precisely,∫ K
(p) q,γ̃(zk, y)K (p) q,γ̃(y,z`)dγ̂N (y) =
√ q(zk)
γ̃(zk)
[ 1
N N∑ n=1 K(p)q (zk, zn)K (p) q (zn, z`) q(zn) γ̃(zn)
]√ q(z`)
γ̃(z`) .
IfN is large compared to p and d, so that in particular γ̃ ≈ γ, the term within brackets will be close to∫ K (p) q (zk, z)K (p) q (z, z`)q(z)dz = K (p) q (zk, z`), so that K (p) q,γ̃ |D is almost a projection in L2(γ̂N ).
Let us actually consider the orthogonal projection matrix K̃ with the same eigenvectors as K(p)q,γ̃ |D, but with the p largest eigenvalues replaced by 1, and the rest of the spectrum set to 0. By the Macchi-Soshnikov theorem, K̃ is the kernel matrix of a DPP; see Section 3.1. We thus consider a minibatch A ∼ DPP(K̃). Coming from a projection DPP, |A| = p almost surely, and we define the gradient estimator
ΞA,DPP := ∑ i∈A ∇θL(zi, θ) K̃ii . (12)
3See [10, Section 4] for details. It suffices that each qi is positive on [−1, 1]d.
In Section 4, we shall prove that ΞA,DPP is unbiased, and examine under what assumptions its variance decreases faster than 1/p.
On the computational cost of ΞA,DPP. The bottleneck is computing the p largest eigenvalues of matrix (11), along with the corresponding eigenvectors. This can be done once before running SGD, as a preprocessing step. Note that storing the kernel in diagonalized form only requires O(Np) storage. Each iteration of SGD then only requires sampling a rank-p projection DPP with diagonalized kernel, which takes O(Np2) elementary operations [22]. In practice, as the complexity of the model underlying z, θ 7→ ∇L(z, θ) increases, the cost of computing p individual gradients shall outweigh this O(Np2) overhead. For instance, learning the parameters of a structured model like a conditional random field leads to arbitrarily costly individual gradients, as the underlying graph gets more dense [25]. Alternately, (12) can be sampled directly, without sampling the underlying DPP. Indeed the Laplace transform of (12) is a Fredholm determinant, and it is shown in [26] that Nyström-type approximations of that determinant, followed by Laplace inversion, yield an accurate inverse CDF sampler.
Finally, we stress the unusual way in which our finite DPP kernel K̃ is constructed, through a reweighted continuous OPE kernel, restricted to the actual dataset. This construction is interesting per se, as it is key to leveraging analytic techniques from the continuous case in Section 4.

heading: 3.4 Our second estimator: sample the OPE, but smooth the gradient
text: In Section 3.3, we smoothed the empirical distribution of the data and restricted a continuous kernel to the dataset D, to make sure that the drawn minibatch would be a subset of D. But one could actually define another gradient estimator, directly from an OPE sample A = {w1, . . . ,wp} ∼ DPP(K(p)q , q). Note that in that case, the “generalized minibatch" A ⊂ [−1, 1]d is not necessarily a subset of the dataset D. Defining a kernel density estimator of the gradient,
∇̂θL(z, θ) := 1
Nhd N∑ i=1 ∇θL(zi, θ) · k ( z − zi h ) ,
we consider the estimator
ΞA,s = ∑ wj∈A ∇̂θL(wj , θ) q(wj)K (p) q (wj ,wj) .
On the computational cost of ΞA,s. Since each evaluation of this estimator is at least as costly as evaluating the actual gradient ΞN (θ), its use is mostly theoretical: the analysis of the fluctuations of ΞA,s is easier than that of ΞA,DPP, while requiring the same key steps. Moreover, the computation of all pairwise distances in ΞA,s could be efficiently approximated, possibly using random projection arguments [27], so that the limited scope of ΞA,s might be overcome in future work. Note also that, like ΞA,DPP, inverse Laplace sampling [26] applies to ΞA,s.

heading: 4 Analysis of determinantal sampling of SGD gradients
text: We first analyze the bias, and then the fluctuations, of the gradient estimators introduced in Section 3. By Proposition 1 with Φ = ∇θL(·, θ), it comes
E[ΞA,DPP|D] = ∫ D ∇θL(z, θ) K (p) q,γ̃(z, z) −1 · K (p) q,γ̃(z, z)dγ̂N (z) = 1 N N∑ i=1 ∇θL(zi, θ),
so that we immediatly get the following result. Proposition 2. E[ΞA,DPP|D] = ΞN .
Thus, ΞA,DPP is unbiased, like the classical Poissonian benchmark in Section 2. However, the smoothed estimator ΞA,s from Section 3.4 is slightly biased. Note that while the results of [3], like (4), do not apply to biased estimators, SGD can be analyzed in the small-bias setting [28]. Proposition 3. Assume that k in (10) has compact support and q is bounded onD. Then E[ΞA,s|D] = ΞN +OP (ph/N).
Proof. Using the first part of Proposition 1 again, it comes
E[ΞA,s|D] = EA∼DPP(K(p)q ,q)  ∑ wj∈A ∇̂θL(wj , θ) q(wj)K (p) q (wj ,wj)  =
∫ D
∇̂θL(w, θ)
(((( (((q(w)K (p) q (w,w)
·((((( (( K(p)q w,w)q(w)dw = ∫ D ( 1 Nhd N∑ i=1 ∇θL(zi, θ)k ( z − zi h )) dw
= 1
N N∑ i=1 (∫ D 1 hd · k ( w − zi h ) dw ) · ∇θL(zi, θ) = 1 N N∑ i=1 ∇θL(zi, θ) +OP (ph/N), (13)
where, in the last step, we have used the fact that h−d · k ( h−1(w − zi) ) is a probability density. The OP (ph/N) error term arises because we integrate that pdf on D rather than Supp(k). But since k is a compactly supported kernel, for points zi that are within a distance O(h) from the boundary of D, we incur an error of OP (1). By Proposition 1, the expected number of such points zi is∫ Dh K (p) q (w,w)q(w)dw, where Dh is the h-neighbourhood of the boundary of D. By a classical asymptotic result of Totik, see [10, Theorem 4.8], w 7→ K(p)q (w,w) is O(p) on Dh; whereas q is a bounded density, implying that ∫ Dh
q(w)dw = O(Vol(Dh)) = O(h). Putting together all of these, we obtain the OP (ph/N) error term in (13).
The rest of this section is devoted to the more involved task of analyzing the fluctuations of ΞA,DPP and ΞA,s. For the purposes of our analysis, we discuss certain desirable regularity behaviour for our kernels and loss functions in Section 4.1. Then, we tackle the main ideas behind the study of the fluctuations of our estimators in Section 4.2. Details can be found in Appendix S3.

heading: 4.1 Some regularity phenomena
text: Assume that ( 1p ·K (p) q (z, z))−1∇θL(z, θ) is bounded on the domain D uniformly in N (note that p possibly depends on N ). Furthermore, assume that z 7→ (p−1 ·K(p)q (z, z)q(z))−1∇θL(z, θ) · γ̃(z) is 1-Lipschitz, with Lipschitz constant OP (1) and bounded in θ. Such properties are natural in the context of various known asymptotic phenomena, in particular asymptotic results on OPE kernels and the convergence of the KDE γ̃ to the distribution γ. The detailed discussion is deferred to Appendix S2; we record here in passing that the above asymptotic phenomena imply that
∇θL(z, θ) 1 p ·K (p) q (z, z)q(z) · γ̃(z) −→ ∇θL(z, θ) · γ(z) d∏ j=1 √ 1− z[j]2. (14)
Our desired regularity properties may therefore be understood in terms of closeness to this limit, which itself has similar regularity. At the price of these assumptions, we can use analytic tools to derive fluctuations for ΞA,DPP by working on the limit in (14). For the fluctuation analysis of the smoothed estimator ΞA,s, we similarly assume that (q(w)· 1pK (p) q (w,w))−1∇̂θL(w, θ) isOP (1) and 1-Lipschitz with an OP (1) Lipschitz constant that is bounded in θ. These are once again motivated by the convergence of (q(w) · 1pK (p) q (w,w))−1∇̂θL(w, θ) to∇θL(w, θ) ∏d j=1 √ 1−w[j]2.

heading: 4.2 Reduced fluctuations for determinantal samplers
text: For succinctness of presentation, we focus on the theoretical analysis of the smoothed estimator ΞA,s, leaving the analysis of ΞA,DPP to Appendix S3, while still capturing the main ideas of our approach.
Proposition 4. Under the assumptions in Section 4.1, Var[ΞA,s|D] = OP (p−(1+1/d)).
Proof Sketch. Invoking (7) in the case of a projection kernel, Var[ΞA,s|D] = ∫∫ ∥∥∥∥∥ ∇̂θL(z, θ)q(z)K(p)q (z, z) − ∇̂θL(w, θ)q(w)K(p)q (w,w) ∥∥∥∥∥ 2
2
|K(p)q (z,w)|2dq(z)dq(w)
= 1
p2 ∫∫ ∥∥∥∥∥ ∇̂θL(z, θ)q(z) · 1pK(p)q (z, z) − ∇̂θL(w, θ) q(w) · 1pK (p) q (w,w) ∥∥∥∥∥ 2
2
|K(p)q (z,w)|2dq(z)dq(w)
.Mθ · 1
p2
∫ ∫ ‖z −w‖22 |K (p) q (z,w)|2dq(z)dq(w), (15)
where we used the 1-Lipschitzianity of ∇̂θL(z,θ) q(z)K (p) q (z,z) , withMθ = OP (1) the Lipschitz constant.
We control the integral in (15) by invoking the renowned Christoffel-Darboux formula for the OPE kernel K(p)q [29]. For clarity of presentation, we outline here the main ideas for d = 1; the details for general dimensions are available in Appendix S3. Broadly speaking, since K(p)q is an orthogonal projection of rank p in L2(q), we may observe that ∫∫ |K(p)q (z,w)|2dq(z)dq(w) = p; so that without the ‖z − w‖22 term in (15), we would have a OP (p−1) behaviour in total, which would be similar to the Poissonian estimator. However, it turns out that the leading order contribution to∫ ∫ ‖z − w‖22 |K (p) q (z,w)|2dq(z)dq(w) comes from near the diagonal z = w, and this turns out to be neutralised in an extremely precise manner by the ‖z −w‖22 factor that vanishes on the diagonal. This neutralisation is systematically captured by the Christoffel-Darboux formula [29], which implies
K(p)q (x, y) = (x− y)−1Ap · (φq(x)φq−1(y)− φq(y)φq−1(x)),
where Ap = O(1) and φq, φq−1 are two orthonormal functions in L2(q). Substituting this into (15), a simple computation shows that ∫∫ ‖z − w‖22 |K (p) q (z,w)|2dq(z)dq(w) = O(1). This leads to a variance bound Var[ΞA,s|D] = OP (p−2), which is the desired rate for d = 1. For general d, we use the fact that q = ⊗d i=1 qi is a product measure, and apply the Christoffel-Darboux formula for each qi, leading to a variance bound of Var[ΞA,s|D] = OP (p−(1+1/d)), as desired.
The theoretical analysis of Var[ΞA,DPP|D] follows the broad contours of the argument for Var[ΞA,s|D] as above, with additional difficulties introduced by the spectral truncation from K(p)q,γ̃ to K̃; see Section 3. This is addressed by an elaborate spectral approximation analysis in Appendix S3. In combination with the ideas expounded in the proof of Proposition 4, our analysis in Appendix S3 indicates a fluctuation bound of Var[ΞA|D] = OP (p−(1+1/d)).

heading: 5 Experiments
text: In this section, we compare the empirical performance and the variance decay of our gradient estimator ΞA,DPP to the default ΞA,Poi. We do not to include ΞA,s, whose interest is mostly theoretical; see Section 4. Moreover, while we focus on simulated data to illustrate our theoretical analysis, we provide experimental results on a real dataset in Appendix S4. Throughout this section, the pdf q introduced in Section 3.2 is taken to be q(w) ∝ ∏d j=1(1 + w[j])
αj (1 − w[j])βj , with αj , βj ∈ [−1/2, 1/2] tuned to match the first two moments of the jth marginal of γ̂N . All DPPs are sampled using the Python package DPPy [30], which is under MIT licence.
Experimental setup. We consider the ERM setting (1) for linear regression Llin((x, y), θ) = 0.5(〈x, θ〉 − y)2 and logistic regression Llog((x, y), θ) = log[1 + exp(−y〈x, θ〉)], both with an additional `2 penalty λ(θ) = (λ0/2)‖θ‖22. Here the features are x ∈ [−1, 1]d−1, and the labels are, respectively, y ∈ [−1, 1] and y ∈ {−1, 1}. Note that our proofs currently assume that the law γ of z = (x, y) is continuous w.r.t. Lebesgue, which in all rigour excludes logistic regression. However, we demonstrate below that if we draw a minibatch using our DPP but on features only, and then deterministically associate each label to the drawn features, we observe the same gains for logistic regression as for linear regression, where the DPP kernel takes into account both features and labels.
For each experiment, we generate N = 1000 data points xi with either the uniform distribution or a mixture of 2 well-separated Gaussian distributions on [−1, 1]d−1. The variable yi ∈ R is generated as 〈xi, θtrue〉+ εi, where θtrue ∈ Rd−1 is given, and ε ∈ RN is a white Gaussian noise vector. In the linear regression case, we scale yi by ‖y‖∞ to make sure that yi ∈ [−1, 1]. In the logistic regression case, we replace each yi by its sign. The regularization parameter λ0 is manually set to be 0.1 and the stepsize for the t-th iteration is set as 1/t0.9, so that (4) applies. In each experiment, performance metrics are averaged over 1000 independent runs of each SGD variant.
Performance evaluation of sampling strategies in SGD. Figure 1 summarizes the experimental results of ΞA,Poi and ΞA,DPP, with p = 5 and p = 10. The top row shows how the norm of the complete gradient ‖ΞN (θt)‖ decreases with the number of individual gradient evaluations t × p, called here budget. The bottom row shows the decrease of ‖θt − θ?‖. Note that using t× p on all x-axes allows comparing different batchsizes. Error bars on the bottom row are ± one standard deviation of the mean. In all experiments, using a DPP consistently improves the performance of Poisson minibatches of the same size, and the DPP with batchsize 5 sometimes even outperforms Poisson sampling with batchsize 10, showing that smaller but more diverse batches can be a better way of spending a fixed number of gradient evaluations. This is particularly true for mixture data (middle column), where forcing diversity with our DPP brings the biggest improvement. Maybe less intuitively, the gains for the logistic regression in d = 11 (last column) are also significant, while the case of discrete labels is not covered yet by our theoretical analysis, and the moderately large dimension makes the improvement in the decay rate of the variance minor. This indicates that there is variance reduction beyond the change of the rate.
Variance decay. For a given dimension d, we want to infer the rate of decay of the variance σ2 = E[‖ΞA(θ?)‖2|D], to confirm the OP (p−(1+1/d)) rate discussed in Section 4. We take L = Llin as an example, with N = 1000 i.i.d. samples D from the uniform distribution on [−1, 1]d. For d = 1, 2, 3, we show in Figure 2 the sample variance of 1000 realizations of the variance of ‖ΞA,Poi(θ?)‖2 (white dots) and ‖ΞA,DPP(θ?)‖2 (black dots), conditionally on D. Blue and red dots indicate standard 95% marginal confidence intervals, for indication only. The slope found by maximum likelihood in a linear regression in the log-log scale is indicated as legend. The experiment confirms that σ2 is smaller for the DPP, decays as p−1−1/d, and that this decay starts at a batchsize p that increases with dimension.

heading: 6 Discussion
text: In this work, we introduced an orthogonal polynomial-based DPP paradigm for sampling minibatches in SGD that entails variance reduction in the resulting gradient estimator. We substantiated our proposal by detailed theoretical analysis and numerical experiments. Our work raises natural questions and leaves avenues for improvement in several directions. These include the smoothed estimator ΞA,s, which calls for further investigation in order to be deployed as a computationally attractive procedure; improvement in the dimension dependence of the fluctuation exponent when the gradients are smooth enough, like [31, 32] did for [10]; sharpening of the regularity hypotheses for our theoretical investigations to obtain a more streamlined analysis. While our estimators were motivated by a continuous underlying data distribution, our experiments suggest notably good performance in situations like logistic regression, where the data is at least partially discrete. Extensions to account for discrete settings in a principled manner, via discrete OPEs or otherwise, would be a natural topic for future research. Another natural problem is to compare our approach with other, non-i.i.d., approaches for minibatch sampling. A case in point is the method of importance sampling, where the independence across data points suggests that the variance should still be OP (1/p) as in uniform sampling. More generally, incorporating ingredients from other sampling paradigms to further enhance the variance reducing capacity of our approach would be of considerable interest. Finally, while our results already partly apply to more sophisticated gradient estimators like PolyakRupert averaged gradients [3, Section 4.2], it would be interesting to introduce repulsiveness across consecutive SGD iterations to further minimize the variance of averaged estimators. In summary, we believe that the ideas put forward in the present work will motivate a new perspective on improved minibatch sampling for SGD, more generally on estimators based on linear statistics (e.g. in coreset sampling), and beyond.
Acknowledgments and Disclosure of Funding
RB acknowledges support from ERC grant BLACKJACK (ERC-2019-STG-851866) and ANR AI chair BACCARAT (ANR-20-CHIA-0002). S.G. was supported in part by the MOE grants R-146-000250-133 and R-146-000-312-114.

Review summaries:
This paper studies unbiased estimators of SGD minibatch via DPP with orthogonal polynomials. The estimator can be constructed by sampling DPP whose kernel consists of the orthogonal projection of some kernel proposition. Samples from the projection DPP are guaranteed fixed sizes (corresponds to the rank) and the variance of minibatch SGD can be achieved by O(p^{-(1 + 1/d)}) for batch size p and feature dimension d. This is smaller than the uniform batch sampling with O(p^{-1}) which leads us that DPP can play a role of variance reduction. Empirical results on both synthetic and real-world datasets show the superior performance of DPP sampling than uniform one in terms of the norm of gradients, distance to the optimal solution as well as the variance of gradient estimations.

The authors introduce the first complete theoretical analysis of minibatch sampling using determinantal point processes in the context of (batched) stochastic gradient descent.  The core of this analysis is built upon continuous DPP results (specifically, on orthogonal polynomial ensembles, or OPEs) in which it has been shown that DPPs achieve faster-than-iid error decay. Assuming knowledge of a KDE $\tilde \gamma$ for the distribution over data points, the authors introduce two batched estimators of the full gradient: 1) A weighted sum of gradients at points sampled by a DPP with kernel $K$, where $K$ is build by restricting the eigenvalues of an OPE kernel reweighted by $\tilde \gamma$. This kernel (and associated DPP) can be computed and sampled from efficiently. 2) An estimator of more theoretical nature, which can sample batches of points that do not arise in the actual data, and is as least as costly in practice as estimating the full gradient.  The authors then prove that the first estimator is unbiased, and obtain a bound on the bias of the second. They then also show that the variance of the resulting estimator decays at a rate $\le p^{-(1 + \frac1d)}$, where $p$ is an integer parameter controlling the number of polynomials used to build the OPE and $d$ is the dimensionality of the input data.  Finally, the authors show on synthetic data that the resulting samplers improve the empirical convergence speed of SGD as a function of the number of required gradient evaluations, and confirm empirically the theoretical variance decay rate obtained earlier.

The authors introduce a sampling minibatches method in SGD based on determinantal point processes (DPPs) to reduce the estimation variance. Specifically, the authors proposed two gradient estimators using multivariate orthogonal polynomial ensembles, and provide theoretical results on the bias and the variance of the estimators. One of the two proposed estimators is tested on a synthetic dataset.  

This paper presents a framework for variance reduction in stochastic gradient descent (SGD) that involves detrimental point processes (DPPs) based on orthogonal polynomials.  This framework uses a class of continuous DPPs based on orthogonal polynomials, and kernel density estimators built from the data, to perform minibatch sampling in SGD.  Theoretical analysis of this approach proves a bound on the variance for two gradient estimators, leading to variance reduction compared to uniform minibatch sampling without replacement.  This approach is empirically validated using experiments performed on synthetic simulated data.