Gradient-Based Empirical Risk Minimization Using Local Polynomial Regression

1. Introduction

Empirical risk minimization (ERM) is one of the mainstays of contemporary machine learning. Indeed, training tasks, such as classification, regression, or representation learning using deep neural networks, can all be formulated as specific instances of ERM. In this paper, we consider gradient-based iterative optimization methods that are used to perform ERM, such as batch gradient descent (GD) (Nesterov 2004), stochastic gradient descent (SGD) (Nemirovski et al. 2009), and their refinements (e.g., Nesterov 1983, Dekel et al. 2012, Johnson and Zhang 2013). Because the primary computational bottleneck in such settings is gradient computation, a significant portion of this literature is concerned with analyzing the convergence rates of different iterative algorithms and determining which ones have smaller oracle complexities—that is, smaller numbers of gradient queries required to obtain approximate minimizers—under various levels of convexity, smoothness, sampling, and other assumptions (see, e.g., Bubeck 2015, Jain and Kar 2017, Netrapalli 2019, and the references therein). We focus on the widely studied setting of strongly convex loss functions. In this case, neglecting any implicit dependence on dimension, the (first-order) oracle complexity of GD is known to be $O(n \log ({ϵ}^{-1}))$ (Nesterov 2004, theorem 2.1.15), where n is the number of training samples in the ERM objective function and $ϵ>0$ is the desired approximation accuracy of the solution. In comparison, it is well-known that SGD has an oracle complexity of $O({ϵ}^{-1})$ (Nemirovski et al. 2009, Bottou et al. 2018), again neglecting any implicit dependence on dimension. It turns out that by carefully accounting for dependence on the dimension p of the parameter of interest in such an ERM, the oracle complexities of GD and SGD are $O(n \log (p{ϵ}^{-1}))$ and $O(p{ϵ}^{-1})$, respectively (see Propositions A.1 and A.2 in Appendix A). Typically, when performing ERM, the objective is to learn the true risk minimizer as accurately as possible. However, the empirical risk is statistically inaccurate, generally, up to a scaling of order $n^{-1/2}$. Hence, the typical regime of interest is $ϵ=\mathrm{\Omega }(n^{-1/2})$ when n is large. For small p (i.e., $p\ll \sqrt{n}$), this makes SGD, at least theoretically, preferable to GD. This is corroborated by the ubiquitous use of SGD and its variants for ERM in practice. Our goal is to improve upon the oracle complexities of GD, SGD, and their variants by exploiting additional information about the loss functions that has heretofore been neglected in the literature.

Recall that the general unconstrained ERM problem is usually abstracted as (cf. Bubeck 2015, section 1, and Hastie et al. 2009):

where each $f_i:{\mathbb{R}}^p\to \mathbb{R}$ represents the loss function corresponding to a different sample of training data. (Note that this formulation, modulo the extra scaling $n^{-1}$, is also utilized in distributed optimization, where each f_i represents the local objective function corresponding to an individual agent in a network of n agents (Nedić 2015).) As suggested by the form of (1), typical analysis of iterative gradient-based algorithms, such as GD or SGD, makes no assumptions about f_i as a function of i. However, the loss function often varies quite smoothly with the data in ERM. For example, given training data $\{(x^{(i)},y_i)\in {\mathbb{R}}^p\times \mathbb{R}:i\in [n]\}$ (with $[n]\triangleq \{1,\dots ,n\}$), the classical least-squares formulation of linear regression with regularization is (see, e.g., Hastie et al. 2009, chapter 3):where the objective is to find the optimal coefficients $\theta \in {\mathbb{R}}^p$, and $\mathcal{R}(\theta )$ is some regularization term—for example, ridge or lasso. Letting $f(x,y;\theta )={(y-{\theta }^{\text{T}}x)}^2+\mathcal{R}(\theta )$ for $x,\theta \in {\mathbb{R}}^p$ and $y\in \mathbb{R}$, we see that each $f_i(\theta )=f(x^{(i)},y_i;\theta )$ in formulation (1) for this problem. Clearly, f is very smooth as a function of the data (x, y), regardless of the regularization $\mathcal{R}(\theta )$. This observation that loss functions are often smooth in the data holds in more complex learning scenarios, too. For instance, when training a deep neural network with a smooth loss—for example, cross-entropy loss—and smooth activation functions—for example, logistic (sigmoid) functions—the induced loss function for the ERM objective function is infinitely differentiable in both the data and the network’s weight parameters (Bishop 2006, chapter 5). Because such smoothness in the data has not been exploited by modern analysis of iterative gradient-based optimization algorithms for (1) (see Section 1.3), the main thrust of this work is to speed up algorithms for ERM using this additional smoothness in data when the data dimension is small relative to n.

To this end, in the remainder of this work, we consider the following problem formulation of ERM. For $d,n\in \mathbb{N}$, assume that we are given n samples of d-dimensional training data $\mathcal{D}=\{x^{(i)}\in {[h^{\prime },1-h^{\prime }]}^d:i\in [n]\}$, where $h^{\prime }>0$ is an arbitrarily small parameter (to be determined later; see (25)), and each training sample belongs to the hypercube ${[h^{\prime },1-h^{\prime }]}^d$.¹ Note that when $d\ge 2$, the first d − 1 elements of a sample can be perceived as a $(d-1)$-dimensional feature vector, and the last element of the sample can be construed as the corresponding continuous-valued label (as in the linear regression example above). For a given $p\in \mathbb{N}$, fix the (parametrized) loss function $f:{[0,1]}^d\times {\mathbb{R}}^p\to \mathbb{R},\hspace{0.17em}f(x;\theta )$, where $\theta =({\theta }_1,\dots ,{\theta }_p)\in {\mathbb{R}}^p$ is a p-dimensional parameter vector. Now, define the ERM objective function $F:{\mathbb{R}}^p\to \mathbb{R}$ as follows:

which is the empirical risk—that is, the empirical expectation of the loss function with respect to the training data $\mathcal{D}$. With this objective function, we consider the unconstrained minimization problem:where $F_{*}\in \mathbb{R}$ denotes the minimum value. In contrast to the canonical formulation (1), our formulation in (3) and (4) explicitly states the dependence of the loss function on the data. The next subsection delineates the main idea of our new inexact gradient descent algorithm to solve this problem, which exploits the smoothness of loss functions in the data. This is followed by discussions of our main contributions, related literature, and an outline of the rest of the paper.

1.1. Main Idea: Local Polynomial Interpolation-Based Gradient Descent

To present the main idea of our algorithm, it is useful to introduce the inexact gradient descent framework, where inexact gradients are used with bounded ${\mathcal{\ell }}^2$-norm error guarantees (cf. Friedlander and Schmidt 2012, section 1.1).² Given the objective function $F(\theta )$, a sequence of step sizes $\{{\alpha }_t>0:t\in \mathbb{N}\}$, and an initial parameter vector ${\theta }^{(0)}\in {\mathbb{R}}^p$, an inexact gradient descent algorithm proceeds via the iterations

where ${\theta }^{(t)}$ is the parameter vector generated by the tth iteration, and for each $t\in \mathbb{N},\hspace{0.17em}{\widehat{\nabla F}}^{(t)}$ denotes an approximation of the true gradient ${\nabla }_{\theta }F({\theta }^{(t-1)})$. Moreover, the inexact gradient approximations possess ${\mathcal{\ell }}^2$-norm error guarantees of the formfor each $t\in \mathbb{N}$, where $\{{\mathrm{\Delta }}_t\ge 0:t\in \mathbb{N}\}$ is the sequence of error bounds. Note that for standard GD, ${\mathrm{\Delta }}_t=0$ for $t\in \mathbb{N}$. On the other hand, for stochastic optimization algorithms, Δ_t bounds the expected values of squared ${\mathcal{\ell }}^2$-norm errors. These error guarantees are useful for the convergence analysis of inexact gradient descent methods (e.g., Friedlander and Schmidt 2012, lemma 2.1).

In this work, we propose an inexact gradient descent method that uses local polynomial interpolation to produce the approximations of true gradients over iterations. Specifically, we first construct a (finite) uniform grid $\mathcal{G}\subset {[0,1]}^d$ over the data space ${[0,1]}^d$. Then, for every training data sample $x^{(j)}\in \mathcal{D}$, we precompute the (scalar-valued) interpolation weights $\{w_y^{*}(x^{(j)}):y\in \mathcal{G}\}$ such that

for all i and any $\theta \in {\mathbb{R}}^p$—that is, linear combinations of partial derivatives evaluated at grid points with the aforementioned weights are “good” estimators of the partial derivatives at training samples (this will be rigorized in Section 3). These weights are calculated using basic linear algebraic operations so that local polynomial regression analysis ensures that the above approximation holds; we defer the details of local polynomial regression to Sections 3.1 and 3.2.1. Notably, these linear weights need only be computed once for a specific choice of bounded kernel function, and calculating them does not require any gradient computation of the loss function.

After these precomputation stages, we run inexact gradient descent using a pertinent choice of fixed step size ${\alpha }_t=\alpha >0$ for an appropriately chosen number of iterations $T\in \mathbb{N}$. In particular, at iteration $t\in \mathbb{N}$, due to the finite sum structure of the ERM objective function in (3), we construct the following approximation of the true gradient ${\nabla }_{\theta }F({\theta }^{(t-1)})$ using the previous iteration’s parameter vector ${\theta }^{(t-1)}$:

Constructing this local polynomial interpolation-based approximation ${\widehat{\nabla F}}^{(t)}$ requires us to compute $|\mathcal{G}|$ gradients of the loss function $\{{\nabla }_{\theta }f(y;{\theta }^{(t-1)}):y\in \mathcal{G}\}$. Then, we compute the standard inexact gradient descent update

as shown in (5). We refer to this new approach as the Local Polynomial Interpolation-based Gradient Descent (LPI-GD) algorithm. This algorithm and its analysis are presented in detail in Section 3.

We note that the inexact gradients ${\widehat{\nabla F}}^{(t)}$ are “good” approximations of ${\nabla }_{\theta }F({\theta }^{(t-1)})$, and we can derive ${\mathcal{\ell }}^2$-norm error guarantees Δ_t as long as (7) is a “good” approximation. Intuitively, when the inexact gradients in LPI-GD are good approximations of the true gradients and the grid size $|\mathcal{G}|$ is much smaller than the number of training samples $|\mathcal{D}|=n$, we expect LPI-GD to outperform GD and related methods in oracle complexity. As we will illustrate in the sequel, sufficient smoothness of loss functions f in the data yields accurate approximations of true gradients in LPI-GD and small grid sizes $|\mathcal{G}|$ when the data dimension d is also small with respect to n. This small data dimension requirement, at least theoretically, restricts the scope of problems where LPI-GD can be applied. However, LPI-GD could be useful in settings where the data dimension is inherently small or can be reduced using dimensionality-reduction techniques (as discussed in Section 1.3).

1.2. Main Contributions

Having conveyed the main idea of our approach, we now delineate the main contributions of this work. Note that our theoretical results hold under the smoothness and strong convexity assumptions stated in Section 2.2.

1. We propose a new inexact gradient descent algorithm, LPI-GD, to compute approximate solutions of (4) in Section 3.2.1 (see Algorithm 1). As noted in Section 1.1, the main innovation of this algorithm is to use the smoothness of the loss function f in the data to learn the gradient oracle of f at every iteration by performing local polynomial regression based on oracle queries over a set of virtual data points. Interestingly, our estimators for gradients are linear, and the corresponding weights must only be computed once in the algorithm.

2. We derive the iteration and oracle complexities of this algorithm in Proposition 2 and Theorem 2, respectively, in Section 3.2.2. In particular, we show that the oracle complexity scales like $\tilde{O}({(p{ϵ}^{-1})}^{d/(2\eta )})$ (neglecting subdominant factors), where $\eta >0$ denotes the Hölder exponent determining the smoothness of the gradient of the loss function with respect to the data (see Section 2.2). Although we focus on strongly convex loss functions to leverage simple convergence results on inexact gradient descent from the literature, analogous analysis could be carried out for larger classes of loss functions—for example, when $f(x;·):{\mathbb{R}}^p\to \mathbb{R}$ is nonconvex.

3. Furthermore, we show in Propositions 3 and 4 in Section 3.2.2 that if the data dimension is sufficiently small—for example, $d=O(\log \log (n))$ (as in healthcare data analytics or control applications; see Section 1.3)—and the loss function f exhibits modest smoothness in the data—that is, $\eta =\mathrm{\Theta }(d)$—then the oracle complexity for our method beats the oracle complexities (or, more precisely, the oracle complexity bounds) of GD, SGD, and their variance reduced and mini-batch variants for any $p=O(\text{poly}(n))$ and any $ϵ=\mathrm{\Theta }(\text{poly}(n^{-1}))$. In particular, when we let $ϵ=\mathrm{\Theta }(n^{-1/2})$ have the smallest scaling of interest (as explained earlier) and $p=O(n^{\beta })$ for any $\beta >0$, Table 1 depicts that the oracle complexities of GD, SGD, and LPI-GD are $O(n \log (n)),\hspace{0.17em}O(n^{\beta +(1/2)})$, and $O(n^{\delta })$ for any small $\delta >0$ (when the loss function is smooth enough), respectively. This demonstrates that smoothness of loss functions in the data can be successfully exploited in gradient-based algorithms for ERM to significantly lower their oracle complexity when the data dimension is small.

4. Finally, in order to analyze the convergence rate of our algorithm, we also generalize the technique of local polynomial regression, or, more precisely, local polynomial interpolation (because there is no noise in our problem), in nonparametric statistics to the multivariate setting in Theorem 1 in Section 3.1. Although it is implied in many standard expositions that univariate local polynomial regression can be generalized to multiple variables (see, e.g., Fan and Gijbels 1996, Tsybakov 2009, or Wasserman 2019), the multivariate case requires significantly more care than the univariate case. Consequently, a detailed and rigorous nonasymptotic analysis of local polynomial interpolation is difficult to find in the literature. Therefore, we believe that Theorem 1 may be of independent interest in statistics.

Table 1. Oracle Complexities of GD, SGD, and Our LPI-GD Algorithm for Strongly Convex Loss Functions

1.3. Related Literature

We next discuss relevant existing work to contextualize the contributions of this work. There is an enormous literature concerning the analysis of convergence rates of exact and inexact gradient-based algorithms that perform the optimization in (1). It is not feasible to survey this entire literature here, and, instead, we only mention some closely related works. For example, in our strongly convex setting of interest, the convergence analysis of GD can be found in the standard text Nesterov (2004), and corresponding convergence analysis of SGD can be found in Nemirovski et al. (2009) (also see Jain and Kar 2017, Bottou et al. 2018, Netrapalli 2019, and the references therein). Moreover, various refinements and improvements of these basic iterative algorithms have been analyzed in the literature. Notable examples include the accelerated or momentum-based GD (Polyak 1964, Nesterov 1983), mini-batch SGD methods (Dekel et al. 2012), and variance reduced SGD methods (Le Roux et al. 2012, Johnson and Zhang 2013, Shalev-Shwartz and Zhang 2013, Schmidt et al. 2017, Schmidt et al. 2017). We also refer readers to Bubeck (2015) for a unified treatment of many of the aforementioned approaches. More generally, inexact gradient descent methods for ERM have been analyzed extensively as well—see, for example, Friedlander and Schmidt (2012), Devolder et al. (2013), So and Zhou (2017), and the references therein for work on the strongly convex case. (Note that SGD can be perceived as an example of an inexact GD method, which uses unbiased estimates of the true gradient of the ERM objective function at every iteration.) This paper and most of the aforementioned references focus on obtaining upper bounds on oracle complexity. Conversely, several authors have also established fundamental lower bounds; cf. Agarwal et al. (2009), Carmon et al. (2021), and the references therein. As mentioned earlier, to our knowledge, all of this literature concerns the formulation in (1). Hence, these methods do not exploit smoothness of loss functions in data by learning gradients in every iteration as in our proposed method, LPI-GD.

On the other hand, motivated by the utility of gradient information for various tasks like variable selection and determining coordinate covariations (cf. Mukherjee and Zhou 2006 and Mukherjee and Wu 2006), the problem of learning gradients in the context of supervised learning problems, such as classification and regression, has received a lot of attention from the machine learning community. Specifically, most of this literature aims to simultaneously learn a classification or regression function along with its gradient from training data. For example, the text Fan and Gijbels (1996) presents the classical theory of gradient estimation using local polynomial fitting, and De Brabanter et al. (2013) present a gradient estimation approach that uses local polynomial regression, but does not estimate the regression function. In a related vein, Delecroix and Rosa (1996) study kernel methods for gradient estimation. Various other methods for learning gradients have also been analyzed, such as reproducing kernel Hilbert space methods (Mukherjee and Wu 2006, Mukherjee and Zhou 2006), regression splines (Zhou and Wolfe 2000), and nearest-neighbor methods (Ausset et al. 2021). However, the vast majority of these techniques have not been considered in the context of optimization for machine learning.

Two recent exceptions where gradient estimation has been applied to optimization are Wang et al. (2018) and Ausset et al. (2021). In particular, Wang et al. (2018) illustrate that learning gradients leads to desirable convergence rates for derivative-free optimization algorithms, which only use zeroth order oracles (i.e., these oracles only output function values, not gradient values). Similarly, Ausset et al. (2021) portray the effectiveness of their nearest-neighbor-based gradient estimator for derivative-free optimization. It is worth noting that both of these works analyze optimization of a single objective function (rather than a sum, as in (1)). Hence, the results derived in these papers are not particularly insightful in our ERM setting in (3) and (4), where we seek to learn gradients based on a first-order oracle by exploiting the smoothness of $f(x;\theta )$ in the data x.

Lastly, as we mentioned in Section 1.2, we establish in Propositions 3 and 4 that the oracle complexity of our LPI-GD algorithm beats GD, SGD, and their variants for a range of values of p and ϵ, if the data dimension d is sufficiently small—for example, $d=O(\log \log (n))$. Because we utilize multivariate local polynomial regression in our algorithm, a “small d assumption” is theoretically unavoidable due to the curse of dimensionality (cf. Fan and Gijbels 1996, section 7.1). This restriction on d can be a limitation for a variety of ERM applications where d is very large, and LPI-GD is a suitable choice primarily when data can be restricted to live in small dimensions. We close this discussion of related literature by delineating some application domains where LPI-GD may be suitable. Indeed, among the multitude of contexts where ERM is used, certain contexts have problems where the number of data attributes is inherently small. For example, in the healthcare data analytics space, electronic health records often have very few features per patient—for example, temperature, blood pressure, initial diagnosis, and medications—but machine learning techniques are still widely used to build models for diagnosing diseases or predicting risk (cf. Choi et al. 2018 and its references). As another example, in the problem of system identification for robotics and control applications, neural networks have been used to build nonlinear state-space models (cf. Nechyba and Xu 1994, Nelles 2001, and their references). The number of inputs of these networks is often very small when the model order is known or assumed to be low. On the other hand, in certain applications with high-dimensional “raw data,” data scientists believe that the useful information in the data lives in a low-dimensional manifold. This belief is the basis of a myriad of dimensionality reduction techniques, such as principal component analysis (Pearson 1901, Hotelling 1933), canonical correlation analysis (Hotelling 1936), Laplacian eigenmaps (Belkin and Niyogi 2001) or diffusion maps (Coifman and Lafon 2006), and modal decompositions (or correspondence analysis), cf. Huang et al. (2024), Makur (2019, chapter 4), and the references therein. As a result, certain high-dimensional data analysis pipelines proceed by first preprocessing the data using pertinent dimensionality reduction techniques and then using the extracted low-dimensional features for other learning tasks. Our contributions in this work could be of utility in such high-dimensional settings where a dimensionality reduction procedure is used before ERM. In applications where the data dimension is large and dimensionality reduction cannot be applied, however, LPI-GD is likely not a suitable algorithm.

1.4. Outline

Finally, the outline of the remainder of this paper is as follows. We introduce required notation and definitions, formally state key assumptions, and briefly elaborate on the standard framework of oracle complexity in optimization theory in Section 2. In Section 3, we provide a detailed development of the LPI-GD algorithm and present our main results. We prove our results on local polynomial regression in Section 4 and carry out the convergence analysis of LPI-GD in Section 5. Finally, we conclude our discussion and propose future research directions in Section 6. Furthermore, we present known oracle complexity results for GD, SGD, and their variants in Appendix A and discuss the key assumptions used in this work in Appendix B (cf. Jadbabaie et al. 2023, section II-C, appendix I).

2. Formal Setup

In this section, we present relevant notation in Section 2.1, important assumptions and definitions in Section 2.1, and the notion of oracle complexity in Section 2.3. These notions will be used throughout this work.

2.1. Notation

We briefly introduce some notation that is utilized in this work. Let $\mathbb{N}\triangleq \{1,2,3,\dots \}$ and ${\mathbb{Z}}_{+}\triangleq \{0,1,2,3,\dots \}$ denote the sets of natural numbers and nonnegative integers, respectively. As mentioned earlier, for any $n\in \mathbb{N}$, let $[n]\triangleq \{1,\dots ,n\}$ be the set of positive integers less than n + 1. In the sequel, we analyze a d-variate problem setting with $d\in \mathbb{N}$. So, we briefly introduce the standard multi-index notation. For any d-tuple of nonnegative integers $s=(s_1,\dots ,s_d)\in {\mathbb{Z}}_{+}^d$, we let $|s|=s_1+\cdots +s_d,\hspace{0.17em}s!=s_1!\cdots s_d!$, and $x^s=x_1^{s_1}\cdots x_d^{s_d}$ for $x\in {\mathbb{R}}^d$. Moreover, given a continuously differentiable function $f:{\mathbb{R}}^d\to \mathbb{R},\hspace{0.17em}f(x_1,\dots ,x_d)$, we write ${\nabla }^{s}f={\partial }^{|s|}f/(\partial x_1^{s_1}\cdots \partial x_d^{s_d})$ to denote its sth partial derivative, and ${\nabla }_{x}f={\left[\frac{\partial f}{\partial x_1}\hspace{0.17em}\cdots \hspace{0.17em}\frac{\partial f}{\partial x_d}\right]}^{\text{T}}$ to denote its gradient. We will sometimes find it convenient to index certain vectors or matrices using the set $\{s\in {\mathbb{Z}}_{+}^d:|s|\le l\}$. In these cases, to help with calculations, we will assume that these vectors or matrices are written out with the indices in lexicographical order (i.e., for any $r,s\in {\mathbb{Z}}_{+}^d$ with $r\ne s$, we say that r < s if r_j < s_j for the first position $j\in [d]$, where $r_j\ne s_j$).

Furthermore, for any vector $x\in {\mathbb{R}}^n$ with $n\in \mathbb{N}$, we let x_i be the ith entry of x, where $i\in [n]$ (or more generally, i belongs to an appropriate index set), and we let $\Vert x{\Vert }_p$ be the usual ${\mathcal{\ell }}^p$-norm of x with $p\in [1,+\infty ]$. Similarly, for any matrix $A\in {\mathbb{R}}^{m\times n}$ with $m,n\in \mathbb{N}$, we let $A_{i,j}$ be the (i, j)th element of A, where $i\in [m]$ and $j\in [n]$ (or, more generally, i and j belong to appropriate index sets), and we let $\Vert A{\Vert }_{\text{F}}$ be the Frobenius norm of A. In addition, $A^{\text{T}}$ represents the transpose of A, and when m = n, $A^{-1}$ represents the inverse of A if it is invertible, $\text{tr}(A)$ represents the trace of A, $\det (A)$ represents the determinant of A, and ${\lambda }_{\mathsf{\text{min}}}(A)$ represents the minimum eigenvalue of A. We also let I_n be the n × n identity matrix. For any pair of symmetric matrices $A,B\in {\mathbb{R}}^{n\times n}$, we use $⪰$ to denote the Löwner partial order and write $A⪰B$ to mean that A − B is positive semidefinite.

Finally, we introduce some useful miscellaneous notation. We let $\lceil ·\rceil$ denote the ceiling function, $\mathrm{𝟙}\{·\}$ denote the indicator function that equals one if its input proposition is true and equals zero otherwise, $\exp (·)$ and $\log (·)$ denote the natural exponential and logarithm functions with base e, respectively, and $\mathbb{E}[·]$ denote the expectation operator, where the underlying probability measure can be easily inferred from context. Throughout this work, we will also utilize standard Bachmann-Landau asymptotic notation—that is, $O(·),\hspace{0.17em}\mathrm{\Theta }(·),\hspace{0.17em}\omega (·)$, etc.—where the underlying limiting variable (that usually tends to infinity) will be clear from context. Moreover, we will let $\text{poly}(n)$ and $\text{polylog}(n)$ denote polynomial and polylogarithmic scaling in the variable n, respectively.

2.2. Smoothness Assumptions and Approximate Solutions

In order to perform tractable convergence analysis of gradient-based iterative algorithms that solve (4), we must formally define the problem class by specifying any smoothness conditions on the loss function, the appropriate notion of “approximate solutions,” and the oracle model of computational complexity used (Nesterov 2004, section 1.1). This section and Section 2.3 precisely define our problem class.

Recall our formulation where we are given n samples of d-dimensional training data $\mathcal{D}$ and a loss function $f:{[0,1]}^d\times {\mathbb{R}}^p\to \mathbb{R},\hspace{0.17em}f(x,\theta )$, and we consider the ERM problem in (3) and (4). Let $ϵ>0$ be any small constant representing the desired level of approximation accuracy. We begin by imposing the following assumptions on our loss function (also see Jadbabaie et al. 2023, section II-B):

1. Smoothness in parameter: There exists $L_1>0$ such that for all fixed $x\in {[0,1]}^d$, the gradient of f with respect to θ, denoted ${\nabla }_{\theta }f(x;·):{\mathbb{R}}^p\to {\mathbb{R}}^p$, exists and is L₁-Lipschitz continuous as a function of the parameter vector:

   $$\forall {\theta }_1,{\theta }_2\in {\mathbb{R}}^p,\hspace{1em}{\Vert {\nabla }_{\theta }f(x;{\theta }_1)-{\nabla }_{\theta }f(x;{\theta }_2)\Vert }_2\le L_1{\Vert {\theta }_1-{\theta }_2\Vert }_2.$$(10)

2. Smoothness in data: There exist $\eta >0$ and $L_2>0$ such that for all fixed $\vartheta \in {\mathbb{R}}^p$ and for all $i\in [p]$, the ith partial derivative of f (with respect to θ_i) at ϑ, denoted $g_i(·;\vartheta )\triangleq \frac{\partial f}{\partial {\theta }_i}(·;\vartheta ):{[0,1]}^d\to \mathbb{R}$, belongs to the $(\eta ,L_2)$-Hölder class as a function of the data vector x (cf. Tsybakov 2009, definition 1.2). By definition, this means that each $g_i(·;\vartheta ):{[0,1]}^d\to \mathbb{R}$ is $l=\lceil \eta \rceil -1$ times differentiable, and for every $s=(s_1,\dots ,s_d)\in {\mathbb{Z}}_{+}^d$ such that $|s|=l$, we have

   $$\forall y_1,y_2\in {[0,1]}^d,\hspace{1em}|{\nabla }^s{g}_i(y_1;\vartheta )-{\nabla }^s{g}_i(y_2;\vartheta )|\le L_2{\Vert y_1-y_2\Vert }_1^{\eta -l}.$$(11)

3. Strong convexity: There exists $\mu >0$ such that for all fixed $x\in {[0,1]}^d$, the map $f(x;·):{\mathbb{R}}^p\to \mathbb{R}$ is μ-strongly convex as a function of the parameter vector:

   $$\forall {\theta }_1,{\theta }_2\in {\mathbb{R}}^p,\hspace{1em}f(x;{\theta }_1)\ge f(x;{\theta }_2)+{\nabla }_{\theta }f{(x;{\theta }_2)}^{\text{T}}({\theta }_1-{\theta }_2)+\frac{\mu }{2}{\Vert {\theta }_1-{\theta }_2\Vert }_2^2.$$(12)

   Moreover, as additional notation, we let $\sigma \triangleq L_1/\mu$ represent the condition number of the problem.

Under these assumptions, we define an ϵ-approximate solution to the optimization problem (4) as any parameter vector ${\theta }^{*}\in {\mathbb{R}}^p$ that satisfies

We make two pertinent remarks at this point. Firstly, it can be argued based on the above strong convexity assumption (as shown in Section 5.1) that $F_{*}$ is finite and achievable by a unique global minimizer of F; in particular, the infimum in (4) may be replaced by a minimum. Secondly, both the first and third assumptions are standard in the optimization literature (cf. Nesterov 2004 and Bubeck 2015). As noted earlier, our main idea is the observation that the second assumption, that gradients of loss functions are also smooth in the data, holds in many machine learning applications. For example, in the regularized linear regression example in (2), ${\nabla }_{\theta }f$ is clearly smooth with respect to the data (assuming $\mathcal{R}$ is differentiable). We refer readers to Appendix B, which is taken from (Jadbabaie et al. 2023, section II-C, appendix I), for further discussion of these standard assumptions within the optimization and statistics contexts, where examples of (${\mathcal{\ell }}^2$-regularized) linear and logistic regression satisfying these assumptions are also provided.

2.3. Oracle Complexity

Finally, we briefly introduce the formalism of “first-order oracle complexity,” upon which our analysis will be based. The notion of oracle complexity was first introduced in Nemirovskiĭ and Yudin (1983) to circumvent the difficulties faced by vanilla computational complexity theory style analysis of convex programming (cf. Agarwal et al. 2009, section 1). As elucidated in Nesterov (2004, section 1.1), although “arithmetical complexity,” which accounts for the running time of all basic arithmetic operations executed by an algorithm, is a more realistic measure of computational complexity, the idea of oracle (or “analytical”) complexity—that is, the total number of oracle calls made by an algorithm—is a more convenient abstraction for the analysis of iterative optimization methods. Indeed, the running time of oracle queries is often the main bottleneck in arithmetical complexity.

Observe that any iterative algorithm that approximately solves optimization problems belonging to the class of ERM problems defined in (3) and (4), under the assumptions in Section 2.2, requires two inputs: (a) the training data $\mathcal{D}$, and (b) a “description” of the specific loss function f under consideration. To fulfill the latter requirement, we will assume that such algorithms have access to a local first-order oracle $\mathcal{O}:{[0,1]}^d\times {\mathbb{R}}^p\to {\mathbb{R}}^p$, defined as:

which returns the value of the gradient ${\nabla }_{\theta }f$ with respect to θ for any input query pair $(x,\vartheta )$. (Note that sometimes the oracle is defined to produce values of both the loss function f and its gradient ${\nabla }_{\theta }f$ at $(x,\vartheta )$ (cf. Nesterov 2004, section 1.1.2).) This oracle $\mathcal{O}$ can be construed as a black box “description” of the gradient of the loss function f with respect to θ.

Given any particular iterative algorithm, $\mathcal{A}=\mathcal{A}(\mathcal{D},\mathcal{O})$, that solves our class of ERM problems, its first-order oracle complexity (or just oracle complexity) $\mathrm{\Gamma }(\mathcal{A})$ is formally defined as the minimum number of oracle calls made by the algorithm $\mathcal{A}$ to obtain an ϵ-approximate solution in the sense of (13) (Nesterov 2004, section 1.1). Thus, $\mathrm{\Gamma }(\mathcal{A})$ is typically a function of the notable problem variables $n,d,p,\eta ,ϵ$, and the fixed constants $L_1,L_2,\mu$. In this paper, following the paradigm of the optimization literature, we will focus on determining the first-order oracle complexity of our proposed optimization method.

3. Main Results and Discussion

In this section, we first present a generalization of local polynomial regression to the multivariate case in Section 3.1 and then describe our proposed algorithm LPI-GD and its convergence analysis in Section 3.2.

3.1. Local Polynomial Regression

We commence by introducing the essentials of local polynomial regression specialized to our noiseless setting. Local polynomial regression was originally introduced in Stone (1977) and has experienced much subsequent development. Most references on the subject only present the univariate (i.e., d = 1) framework—for example, Tsybakov (2009, section 1.6), and Wasserman (2019, chapter “Nonparametric Regression”). However, local polynomial regression has also been developed and analyzed in multivariate settings (i.e., for general d). We refer readers to Fan and Gijbels (1996) for a standard reference on local polynomial regression from classical statistics, Cleveland and Loader (1996, section 2) for more classical references and a brief history of smoothing using local regression, and Zhang and Chan (2011) for a succinct, but reasonably complete, review of classical work on multivariate local polynomial regression. These sources and the references therein convey that whereas the univariate local polynomial regression literature contains nonasymptotic analyses (in the spirit of modern high-dimensional statistics, where both dimension and number of samples are finite, e.g., the presentation in Tsybakov 2009, section 1.6), the multivariate local polynomial regression literature primarily presents asymptotic results (with a more classical feel). For example, Masry (1996, 1999) analyzes multivariate local polynomial regression and proves asymptotic normality results (i.e., central limit theorems) for estimators of the regression function under various assumptions. Because we will require a nonasymptotic conditioning analysis of multivariate local polynomial regression in our subsequent analysis, we rigorously carry out the much more involved analysis for the multivariate case by carefully generalizing the univariate development in Tsybakov (2009, section 1.6).

Consider a function $g:{[0,1]}^d\to \mathbb{R}$ that belongs to the $(\eta ,L_2)$-Hölder class, namely, g is $l=\lceil \eta \rceil -1$ times differentiable and satisfies (11) for every $s\in {\mathbb{Z}}_{+}^d$ with $|s|=l$. In order to present the local polynomial regression-based approximation of this function, we first state some required definitions and notation. Let $K:[-1,1]\to [b,c]$ be a bounded kernel function, where $c\ge b>0$ are some fixed constants, and we will use the extension of this kernel (to the domain $\mathbb{R}$), $K:\mathbb{R}\to [0,c]$, in the sequel with the assumption that the support of K is $[-1,1]$. For any fixed $m\in \mathbb{N}$, define the uniform grid

which has cardinality $|{\mathcal{G}}_m|=m^d$. Finally, define the vector(-valued function) $U:{\mathbb{R}}^d\to {\mathbb{R}}^D$ as:which has dimension $D\triangleq \dim (U)={\sum }_{k=0}^l\left(\begin{array}{c}k+d-1\\ d-1\end{array}\right)=\left(\begin{array}{c}l+d\\ d\end{array}\right)$ (via the hockey-stick identity). Then, for any arbitrarily small value of the bandwidth parameter $h\in \left(0,\frac{1}{2}\right)$, we can write the following weighted regression problem:where $\widehat{\mathrm{\Phi }}(x)={[{\widehat{\mathrm{\Phi }}}_1(x)\hspace{0.17em}\cdots \hspace{0.17em}{\widehat{\mathrm{\Phi }}}_D(x)]}^{\text{T}}$. We refer to the first entry $\widehat{\varphi }(x)\triangleq {\widehat{\mathrm{\Phi }}}_1(x)=U{(0)}^{\text{T}}\widehat{\mathrm{\Phi }}(x)$ of $\widehat{\mathrm{\Phi }}(x)$ as the local polynomial interpolator of order l for $g:{[h,1-h]}^d\to \mathbb{R}$ (which is $g:{[0,1]}^d\to \mathbb{R}$ restricted to the domain ${[h,1-h]}^d$).

To intuitively understand the problem in (17), let us temporarily specialize to the setting where d = 1, η = 2, and the kernel $K(u)=\mathrm{𝟙}\{-1\le u\le 1\}$. Then, $g:[0,1]\to \mathbb{R}$ belonging to the $(2,L_2)$-Hölder class means that g is continuously differentiable and $g^{\prime }:[0,1]\to \mathbb{R}$ is L₂-Lipschitz continuous. Moreover, the uniform grid has the form ${\mathcal{G}}_m=\left\{\frac{1}{m},\frac{2}{m},\dots ,\frac{m-1}{m},1\right\}$ and $U:\mathbb{R}\to {\mathbb{R}}^2$ is given by $U(u)={[1\hspace{1em}u]}^{\text{T}}$. Hence, the problem in (17) becomes:

For any fixed x, the problem in (18) performs linear regression to approximate g in a local neighborhood of x using a subset of the samples $\{(y,g(y)):y\in {\mathcal{G}}_m\}$, which live in the neighborhood. Because the mean value theorem and L₂-Lipschitz continuity of $g^{\prime }$ imply the first-order Taylor approximation $g(y)=g(x)+g^{\prime }(x)(y-x)+O({(y-x)}^2)$ for all y in the neighborhood of x (cf. Rudin 1976, theorem 5.10), we can intuitively interpret the (scaled) optimal regression coefficients ${\widehat{\mathrm{\Phi }}}_1(x)$ and ${\widehat{\mathrm{\Phi }}}_2(x)/h$ as estimators of the first-order Taylor approximation coefficients g(x) and $g^{\prime }(x)$, respectively. This is why $\widehat{\varphi }(x)={\widehat{\mathrm{\Phi }}}_1(x)$ is a reasonable approximation of g(x). We refer readers to standard texts, such as Tsybakov (2009), for further intuition.

Returning to our general setting in (17), we proceed to showing that the interpolator $\widehat{\varphi }$ provides a “good” approximation of the function g in the supremum (or Chebyshev) norm. Toward this end, for any $x\in {[h,1-h]}^d$, construct the D × D symmetric matrix

which is also positive semidefinite due to the nonnegativity of the kernel. We next present a key technical proposition, which generalizes (Tsybakov 2009, lemma 1.5) to the d-variate setting and demonstrates that B(x) is actually positive definite. To present this proposition, we also define the (d, l)-dependent constantwhich will be used in subsequent results as well for crisper exposition.

Proposition 1

(Löwner Lower Bound). Suppose that $l\ge d\ge 19$. Then, for all $m\in \mathbb{N}$ and $h\in \left(0,\frac{1}{2}\right)$ such that $mh\ge 4l{(3e)}^d/\mathrm{\Lambda }(d,l)$ and all $x\in {[h,1-h]}^d$, we have

$$B(x)⪰\frac{b^d\mathrm{\Lambda }(d,l)}{2}{I}_D,$$

where the constant $\mathrm{\Lambda }(d,l)>0$ is defined in (20).

Proposition 1 is established in Section 4.1 by exploiting several basic ideas from the theory of orthogonal polynomials and matrix perturbation theory. The proof is much more sophisticated than Tsybakov (2009, lemma 1.5) because explicit quantitative estimates are derived to deal with general d. Because Proposition 1 shows that B(x) is invertible, a straightforward exercise in calculus conveys that the unique solution $\widehat{\mathrm{\Phi }}(x)$ of the weighted regression problem in (17) can be written as:

where the optimal vector-valued weights $w_y(x)\in {\mathbb{R}}^D$ are given by:

Then, defining the first coordinates of these vectors as the real-valued interpolation weights:

the local polynomial interpolator of order l for $g:{[h,1-h]}^d\to \mathbb{R}$ can be written as the following weighted sum of the interpolation points $\{g(y):y\in {\mathcal{G}}_m\}$:

In Section 4.2, we present some important properties of the interpolation weights given in (23) based on the key estimate in Proposition 1.

The ensuing theorem conveys that local polynomial regression can be used to uniformly approximate a function g (in the supremum norm sense) for sufficiently large uniform grids ${\mathcal{G}}_m$.

Theorem 1

(Supremum Norm Interpolation). Fix any constant $\delta \in (0,1)$, and consider the function $g:{[0,1]}^d\to \mathbb{R}$, which belongs to the $(\eta ,L_2)$-Hölder class. Suppose that $l=\lceil \eta \rceil -1\ge d\ge 19$ and that the grid size $|{\mathcal{G}}_m|=m^d$ is sufficiently large so that

$$m\ge \frac{110(2L_2+1)c}{b}\left(\frac{d{(3e)}^d}{\mathrm{\Lambda }{(d,l)}^2}\right){\left(\frac{1}{\delta }\right)}^{1/\eta },$$

where $c\ge b>0$ are the bounds on our kernel, and $\mathrm{\Lambda }(d,l)$ is defined in (20). Then, the supremum norm approximation error of g is bounded by δ, viz.

$$\underset{x\in {[h,1-h]}^d}{\sup }|\widehat{\varphi }(x)-g(x)|\le \delta ,$$

where $\widehat{\varphi }$ is the local polynomial interpolator of order l for g given in (24), and the bandwidth $h\in \left(0,\frac{1}{2}\right)$ is chosen to be

$$h=\frac{4l{(3e)}^d}{m\mathrm{\Lambda }(d,l)}\le \frac{2b}{55(2L_2+1)c}\left(\frac{\mathrm{\Lambda }(d,l)l}{d}\right).$$

Theorem 1 is proved in Section 4.3 using the pertinent properties of the interpolation weights in (23) derived in Section 4.2. We note that the lower bound on m is expected to scale at least exponentially with the dimension d due to the curse of dimensionality (cf. Fan and Gijbels 1996, section 7.1). However, it may be possible to improve the dependence we have above using a different analysis. In the sequel, we will exploit Theorem 1 to perform inexact gradient descent.

3.2. Inexact Gradient Descent Using Local Polynomial Interpolation

We now turn to describing and analyzing a new inexact gradient descent algorithm that performs local polynomial interpolation at every iteration. In particular, we delineate the algorithm in Section 3.2.1 and derive convergence rates in Section 3.2.2.

3.2.1. Description of Algorithm.

Under the setup and assumptions described in Section 1, suppose we seek to solve the optimization problem (4) in the ϵ-approximate sense of (13) based on some training data $\mathcal{D}$ with the parameter $h^{\prime }>0$ given by

which is the upper bound on the bandwidth h in Theorem 1.³ We first present our proposed algorithm, dubbed Local Polynomial Interpolation-based Gradient Descent (LPI-GD), in general, and then specialize the values of different problem variables with judicious choices later on.

Let ${\theta }^{(t)}\in {\mathbb{R}}^p$ denote the “approximate solution” produced by LPI-GD at iteration $t\in {\mathbb{Z}}_{+}$, and let ${\delta }_t\in (0,1)$ denote the maximum allowable supremum norm approximation error at iteration $t\in \mathbb{N}$. Fix any bounded kernel $K:[-1,1]\to [b,c]$ with $c\ge b>0$ as in Section 3.1, and corresponding to each δ_t at iteration $t\in \mathbb{N}$, consider the uniform grid ${\mathcal{G}}_{m_t}$, where $m_t\in \mathbb{N}$ is given by Theorem 1:

and the bandwidth $h_t\in \left(0,\frac{1}{2}\right)$ is also given by Theorem 1:

We initialize ${\theta }^{(0)}\in {\mathbb{R}}^p$ to be any arbitrary vector in ${\mathbb{R}}^p$ at iteration t = 0. Then, at iteration $t\in \mathbb{N}$, we make $|{\mathcal{G}}_{m_t}|=m_t^d$ first-order oracle calls at all virtual data points in the uniform grid ${\mathcal{G}}_{m_t}$ with the approximate solution at iteration t − 1:

For every $i\in [p]$, these oracle calls give us information about the partial derivatives $\{g_i(y;{\theta }^{(t-1)}):y\in {\mathcal{G}}_{m_t}\}$ with respect to θ_i. Hence, akin to (17) and (24) in Section 3.1, we can construct the local polynomial interpolator ${\widehat{\varphi }}_i^{(t)}:{[h^{\prime },1-h^{\prime }]}^d\to \mathbb{R}$ of order $l=\lceil \eta \rceil -1$ for the $(\eta ,L_2)$-Hölder class function $g_i(·;{\theta }^{(t-1)})=\frac{\partial f}{\partial {\theta }_i}(·;{\theta }^{(t-1)}):{[h^{\prime },1-h^{\prime }]}^d\to \mathbb{R}$ as follows:

where the interpolation weights $w_{y,t}^{*}(x)\in \mathbb{R}$ are defined analogously to (23); in particular, they are the first elements of the weight vectors $w_{y,t}(x)\in {\mathbb{R}}^D$, which are defined as in (22), but with ${\mathcal{G}}_m$ replaced by the grid ${\mathcal{G}}_{m_t}$ and h replaced by the bandwidth h_t. Using the learned interpolators in (29) at iteration t, we next approximate the true gradient of the ERM objective function ${\nabla }_{\theta }F({\theta }^{(t-1)})$ at ${\theta }^{(t-1)}$ as follows. For every $i\in [p]$ and every $j\in [n]$, we first estimate the partial derivative $g_i(x^{(j)};{\theta }^{(t-1)})$ at the data sample $x^{(j)}$ by computing ${\widehat{\varphi }}_i^{(t)}(x^{(j)})$. Then, because we have ${\nabla }_{\theta }F({\theta }^{(t-1)})=\frac{1}{n}{\sum }_{j=1}^n{\nabla }_{\theta }f(x^{(j)};{\theta }^{(t-1)})$ due to (3), we can estimate ${\nabla }_{\theta }F({\theta }^{(t-1)})$ by computingwhere ${[{\widehat{\varphi }}_i^{(t)}(x^{(j)}):i\in [p]]}^{\text{T}}\in {\mathbb{R}}^p$ is our approximation of ${\nabla }_{\theta }f(x^{(j)};{\theta }^{(t-1)})$ at the data sample $x^{(j)}$. Finally, we perform the inexact gradient descent updatewhich produces the approximate solution for iteration t. We run the above procedure for $T\in \mathbb{N}$ iterations, and the LPI-GD algorithm returns ${\theta }^{(T)}\in {\mathbb{R}}^p$ as the ϵ-approximate solution for the optimization problem (4). The analysis in Section 3.2.2 elucidates the specific values taken by the constants δ_t and the total number of iterations T in terms of the other problem parameters. In particular, we will see that the approximation errors δ_t, grid sizes m_t, bandwidths h_t, and uniform grids ${\mathcal{G}}_{m_t}$ all will not change over iterations. A pseudocode summary of the LPI-GD algorithm with pertinent δ_t and T (selected as in Section 3.2.2) is presented in Algorithm 1; the algorithm uses the training data $\mathcal{D}$ and a first-order oracle $\mathcal{O}$ to generate an ϵ-approximate solution that solves (4). Notably, the interpolation weights only need to be computed once in the algorithm, and our local polynomial interpolators for gradients are linear estimators, which use these weights.

We remark that Algorithm 1 portrays a “theoretical” algorithm because it assumes that problem parameters, such as the Lipschitz constant L₁, the strong convexity parameter μ, the minimum value $F_{*}$, and the Hölder exponent η and constant L₂, are known to the practitioner. When some or all of these parameters are unknown, we can use appropriate bounds and estimates based on the loss function $f:{[0,1]}^d\times {\mathbb{R}}^p\to \mathbb{R}$. For example, we may use any upper bound on L₁, any lower bound on μ, and any universal lower bound on f (which yields a lower bound on $F_{*}$). Moreover, for smooth functions f, we can choose any desired level of differentiability η and any upper bound on L₂. When such bounds or estimates are also unavailable, the best we can do is to select T, δ_t, m_t, and h_t to have the correct scaling with respect to ϵ, d, p, and η. It is worth mentioning that our objective in this paper is to demonstrate that the oracle complexities of GD, SGD, and their variants can be theoretically beaten by the proposed algorithm. So, we do not delve further into considerations regarding its practical implementation here.

Algorithm 1

(Local Polynomial Interpolation-Based Gradient Descent (LPI-GD) Algorithm to Approximately Solve (4))

Input: n samples of training data $\mathcal{D}=\{x^{(i)}\in {[h^{\prime },1-h^{\prime }]}^d:i\in [n]\}$

Input: first order oracle $\mathcal{O}:{[0,1]}^d\times {\mathbb{R}}^p\to {\mathbb{R}}^p$ defined in (14)

Input: approximation accuracy level $ϵ>0$

Output: ϵ-approximate solution ${\theta }^{*}\in {\mathbb{R}}^p$ satisfying (13)

  Step 1: Initialization

• 1: Set arbitrary initial parameter vector ${\theta }^{(0)}\in {\mathbb{R}}^p$

• 2: Set number of iterations $T\in \mathbb{N}$ according to (32)

• 3: Set supremum norm approximation error $\delta ={(1-{\sigma }^{-1})}^{T/2}$     ⊳ This is the same for all iterations; see (33)

• 4: Set grid size $m\in \mathbb{N}$ according to (26)      ⊳ Substitute m_t = m and ${\delta }_t=\delta$ in (26)

• 5: Set bandwidth $h=4l{(3e)}^d/(m\mathrm{\Lambda }(d,l))$      ⊳ See (27)

    Step 2: Calculation of Interpolation Weights

• 6: Construct uniform grid ${\mathcal{G}}_m$ according to (15)

• 7: for j = 1 to n do     ⊳ Iterate over training data

• 8:    Compute interpolation weights $\{w_y^{*}(x^{(j)}):y\in {\mathcal{G}}_m\}$ via (22) and (23) using kernel $K:[-1,1]\to [b,c]$

• 9: end for

    Step 3: Inexact Gradient Descent

• 10: for t = 1 to T do

• 11:   Make oracle calls $\{\mathcal{O}(y,{\theta }^{(t-1)}):y\in {\mathcal{G}}_m\}$ to get partial derivatives $\left\{\frac{\partial f}{\partial {\theta }_i}(y;{\theta }^{(t-1)}):y\in {\mathcal{G}}_m,\hspace{0.17em}i\in [p]\right\}$

       Step 3a: Local Polynomial Interpolation

• 12:   for j = 1 to n do      ⊳ Iterate over training data

• 13:    Compute local polynomial interpolators $\{{\widehat{\varphi }}_i^{(t)}(x^{(j)}):i\in [p]\}$ according to (29) using the oracle queries $\left\{\frac{\partial f}{\partial {\theta }_i}(y;{\theta }^{(t-1)}):y\in {\mathcal{G}}_m,\hspace{0.17em}i\in [p]\right\}$ and precomputed interpolation weights $\{w_y^{*}(x^{(j)}):y\in {\mathcal{G}}_m\}$

• 14:   end for

• 15:   Construct estimate of gradient ${\widehat{\nabla F}}^{(t)}$ using interpolators $\{{\widehat{\varphi }}_i^{(t)}(x^{(j)}):i\in [p],\hspace{0.17em}j\in [n]\}$ as shown in (30)

       Step 3b: Standard GD Update

• 16:   Update: ${\theta }^{(t)}\leftarrow {\theta }^{(t-1)}-\frac{1}{L_1}{\widehat{\nabla F}}^{(t)}$      ⊳ See (31)

• 17: end for

• 18: return ${\theta }^{*}={\theta }^{(T)}$

Before proceeding to our analysis, it is worth juxtaposing Algorithm 1 with standard batch GD and SGD. At each iteration, GD exactly computes gradients for all data samples using n oracle calls. Because this can be computationally expensive, inexact gradient descent methods attempt to reduce the number of oracle calls per iteration at the cost of increasing the number of iterations. For example, SGD makes only one oracle call per iteration for a randomly chosen data sample and uses this queried gradient value as a proxy for the true gradient. In contrast, at each iteration, our approach evaluates the gradients at certain representative virtual points in the data space and uses these oracle calls to learn a model for the gradient oracle. This learnt model is then used to approximate the true gradient. Intuitively, when the loss function is a sufficiently smooth function of the data and the data dimension d is small, we can learn an accurate model with a reasonably small number of virtual points. It is also worth mentioning that because we perform regressions to learn our models, our proposed method is likely to be robust to noise in the oracle queries. However, we do not analyze this in the present work.

3.2.2. Analysis of Algorithm.

To compute the oracle complexity of the LPI-GD algorithm, we must first compute the minimum number of iterations it takes to obtain an ϵ-approximate solution to (4). The ensuing proposition chooses appropriate values of δ_t, which turn out to be constant with respect to t, and develops an upper bound on the minimum number of iterations T.

Proposition 2

(Iteration Complexity of LPI-GD). Suppose that $l=\lceil \eta \rceil -1\ge d\ge 19$ and the assumptions of Section 2.2 hold. For any (small) accuracy $ϵ>0$, define the number of iterations $T\in \mathbb{N}$ as

$$T=\lceil {\left(\log \left(\frac{\sigma }{\sigma -1}\right)\right)}^{-1}\log \left(\frac{F({\theta }^{(0)})-F_{*}+\frac{p}{2\mu }}{ϵ}\right)\rceil ,$$(32)

and the maximum allowable supremum norm approximation errors $\{{\delta }_t\in (0,1):t\in [T]\}$ via:

$$\forall t\in [T],\hspace{1em}{\delta }_t^2={\left(1-\frac{1}{\sigma }\right)}^T=\mathrm{\Theta }\left(\frac{ϵ}{p}\right),$$(33)

where $F:{\mathbb{R}}^p\to \mathbb{R}$ is the ERM objective function in (3), and $\sigma =L_1/\mu >1$ is the condition number defined in Section 2.2. Then, after T iterations of the LPI-GD algorithm described in Algorithm 1 in Section 3.2.1 with arbitrary initialization ${\theta }^{(0)}\in {\mathbb{R}}^p$ and updates (31), the parameter vector ${\theta }^{(T)}\in {\mathbb{R}}^p$ is an ϵ-approximate solution to the optimization problem in (4) in the sense of (13):

$$F({\theta }^{(T)})-F_{*}\le ϵ.$$

Proposition 2 is proved in Section 5.1 using Theorem 1 and some known ideas from the analysis of inexact gradient descent methods in Friedlander and Schmidt (2012). Using Proposition 2, the ensuing theorem bounds the first-order oracle complexity of the LPI-GD algorithm in Algorithm 1, denoted $\mathrm{\Gamma }(\text{LPI-GD})$, which is the minimum number of oracle calls made by Algorithm 1 to obtain an ϵ-approximate solution (as defined in Section 2.3). To succinctly present the ensuing theorem, we also define the constant

where $\mathrm{\Lambda }(d,l)$ is defined in (20), and $\mu ,L_1,L_2,b,c$ are constants defined in Sections 2.2 and 3.1.

Theorem 2

(Oracle Complexity of LPI-GD). Suppose that $l=\lceil \eta \rceil -1\ge d\ge 19$ and the assumptions of Section 2.2 hold. For any (small) accuracy $0<ϵ\le \frac{(L_1-\mu )p}{2L_1\mu }$, consider the LPI-GD algorithm described in Algorithm 1 in Section 3.2.1 with arbitrary initialization ${\theta }^{(0)}\in {\mathbb{R}}^p$, number of iterations $T\in \mathbb{N}$ as in (32), maximum allowable approximation errors $\{{\delta }_t\in (0,1):t\in [T]\}$ as in (33), and updates (31). Then, the first-order oracle complexity of this algorithm (to obtain an ϵ-approximate solution in the sense of (13)) is upper bounded by

$$\mathrm{\Gamma }(\text{LPI-GD})\le C_{\mu ,L_1,L_2,b,c}(d,l){\left(\frac{p+2\mu (F({\theta }^{(0)})-F_{*})}{ϵ}\right)}^{d/(2\eta )}\log \left(\frac{p+2\mu (F({\theta }^{(0)})-F_{*})}{2\mu ϵ}\right),$$

where the constant $C_{\mu ,L_1,L_2,b,c}(d,l)$ is defined in (34).

Theorem 2 is established in Section 5.2. Unlike the “dimension-free” oracle complexity bounds for batch GD and SGD, our bound in Theorem 2 explicitly depends on d and p. However, even though it may not be immediately obvious, the oracle complexity of LPI-GD scales much better than the oracle complexities of GD or SGD in important regimes. We next elucidate a scaling law for the bound in Theorem 2 in the setting where $n,d,p,\eta ,{ϵ}^{-1}\to \infty$ and all other problem parameters are constant.

As in Appendix A, let ${\theta }_{*}\in {\mathbb{R}}^p$ be the unique global minimizer of the empirical risk $F:{\mathbb{R}}^p\to \mathbb{R}$ in (3) so that $F_{*}=F({\theta }_{*})$. Assume that the parameter θ belongs to an ${\mathcal{\ell }}^{\infty }$-ball with constant radius (as in many applications). Then, in the aforementioned setting, we have using Bubeck (2015, lemma 3.4) that

where the inequality uses the facts that ${\nabla }_{\theta }F:{\mathbb{R}}^p\to {\mathbb{R}}^p$ is L₁-Lipschitz continuous (see the proof of Proposition 2 in Section 5.1) and ${\nabla }_{\theta }F({\theta }_{*})=0$ (which is the zero vector), and the equality holds because $\Vert {\theta }^{(0)}-{\theta }_{*}{\Vert }_2^2=O(p)$ (see Appendix A). Hence, the first-order oracle complexity of the LPI-GD algorithm in Theorem 2 scales as

For comparison, recall from Proposition A.1 in Appendix A (cf. Nesterov 2004, theorem 2.1.15) that when the empirical risk is strongly convex and has Lipschitz continuous gradient, the first-order oracle complexity of GD is bounded by

Similarly, recall from Proposition A.2 in Appendix A (cf. Nemirovski et al. 2009 and Bottou et al. 2018, theorem 4.7) that when the empirical risk is strongly convex and has Lipschitz continuous gradient, and certain universal second moment bounds are also satisfied on stochastic gradients, the first-order oracle complexity of SGD is bounded by

Note that the expressions in (37) and (38) do not depend on $d,\eta$ and neglect all other constant parameters.

To compare the first-order oracle complexities in (36), (37), and (38), we will assume for simplicity that the quantities $d,p,\eta ,{ϵ}^{-1}$ are all functions of the number of training samples n and study the asymptotics as $n\to \infty$. The ensuing proposition presents a broad scaling regime where the oracle complexity of LPI-GD beats the oracle complexities (or, more precisely, the oracle complexity bounds) of both GD and SGD shown above.

Proposition 3

(Comparison with GD and SGD). Fix any constants $\alpha >0,\hspace{0.17em}\beta >0,\hspace{0.17em}\tau >\text{max}\{1,{(\alpha +\beta )}^{-1}\}$, and $\gamma >\text{max}\{1,\tau (\alpha +\beta )/2\}$. Suppose that $ϵ=\mathrm{\Theta }(n^{-\alpha }),\hspace{0.17em}p=\mathrm{\Theta }(n^{\beta }),\hspace{0.17em}d\le \frac{\log \log (n)}{4 \log (e(\gamma +1))}$, and $\eta =\mathrm{\Theta }(d)$ such that $\eta \ge \tau (\alpha +\beta )d/2$ and $d\le l=\lceil \eta \rceil -1\le \gamma d$. Furthermore, suppose that $\mu ,L_1,L_2,b,c$ are $\mathrm{\Theta }(1)$ with respect to n. Then, we have

$$\begin{array}{l}\mathrm{\Gamma }(\text{LPI-GD})=O\left(\log (n) \exp \left(2\sqrt{\log (n)}\right)n^{1/\tau }\right),\\ {\mathrm{\Gamma }}_{*}(\text{GD})=\mathrm{\Theta }(n \log (n)),\\ {\mathrm{\Gamma }}_{*}(\text{SGD})=\mathrm{\Theta }(n^{\alpha +\beta }).\end{array}$$

This implies the following results:

$$\underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{GD})}=\underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{SGD})}=0,$$

as well as the stronger results:

$$\begin{array}{l}\underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{GD}))}\le \frac{1}{\tau },\\ \underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{SGD}))}\le \frac{1}{\tau (\alpha +\beta )}.\end{array}$$

Proposition 3 is proved in Section 5.3. Under its assumptions, the LPI-GD algorithm not only beats the oracle complexities of both GD and SGD, but the stronger bounds on the limits of log-oracle-complexity-ratios show that it beats them in the exponent (because ${\tau }^{-1},{\tau }^{-1}{(\alpha +\beta )}^{-1}<1$). In fact, under the above assumptions, the LPI-GD algorithm beats any $\mathrm{\Theta }(\text{poly}(n))$ oracle complexity when τ is chosen to be large enough (i.e., when there is sufficient smoothness). Furthermore, as explained in Appendix A, comparisons of oracle complexities of GD and SGD in the literature typically neglect the factors of p in (37) and (38) (see our discussion at the outset of Section 1). In contrast, Proposition 3 performs comparisons with these additional factors because they stem from a reasonable assumption that the parameter lives in an ${\mathcal{\ell }}^{\infty }$-ball. However, because the LPI-GD algorithm can beat any $\mathrm{\Theta }(\text{poly}(n))$ oracle complexity, the oracle complexity of LPI-GD would scale better than the oracle complexities of GD or SGD in the regime outlined in Proposition 3, even if we neglected the factors of p in GD and SGD.

We can also calculate interesting regimes where LPI-GD scales better than more refined gradient-based algorithms. In particular, there are a spectrum of algorithms interpolating between GD and SGD that utilize more oracle calls per iteration than SGD, but less than GD. Two notable classes of such algorithms include variance reduced methods (Le Roux et al. 2012, Johnson and Zhang 2013, Shalev-Shwartz and Zhang 2013, Schmidt et al. 2017), which periodically compute full gradients while running SGD to reduce the variance of gradient estimates over iterations, and mini-batch SGD (Dekel et al. 2012, Bottou et al. 2018), where each iteration uses a randomly chosen mini-batch to approximate the full gradient. We next show that the LPI-GD algorithm can also beat mini-batch SGD and stochastic variance reduced gradient (SVRG)—a representative variance reduced SGD method (Johnson and Zhang 2013). Note that the known first-order oracle complexities of SVRG and mini-batch SGD are presented in Propositions A.3 and A.4 in Appendix A, respectively. As before, we assume that the quantities $d,p,\eta ,{ϵ}^{-1}$ are all functions of the number of training samples n and study the asymptotics as $n\to \infty$.

Proposition 4

(Comparison with SVRG and Mini-Batch SGD). As in Proposition 3, fix any constants $\alpha >0,\hspace{0.17em}\beta >0,\hspace{0.17em}\tau >\text{max}\{1,{(\alpha +\beta )}^{-1}\}$, and $\gamma >\text{max}\{1,\tau (\alpha +\beta )/2\}$. Suppose that $ϵ=\mathrm{\Theta }(n^{-\alpha }),\hspace{0.17em}p=\mathrm{\Theta }(n^{\beta }),\hspace{0.17em}d\le \frac{\log \log (n)}{4 \log (e(\gamma +1))}$, and $\eta =\mathrm{\Theta }(d)$ such that $\eta \ge \tau (\alpha +\beta )d/2$ and $d\le l=\lceil \eta \rceil -1\le \gamma d$. Furthermore, suppose that $\mu ,L_1,L_2,b,c$ are $\mathrm{\Theta }(1)$ with respect to n. Then, the existing upper bounds in Propositions A.3 and A.4 on the first-order oracle complexities of SVRG and mini-batch SGD satisfy

$$\begin{array}{l}\mathrm{\Gamma }(\text{SVRG})\le {\mathrm{\Gamma }}_{*}(\text{SVRG})=\mathrm{\Theta }(n \log (n)),\\ \mathrm{\Gamma }(\text{Mini-batch}\hspace{0.17em}\text{SGD})\le {\mathrm{\Gamma }}_{*}(\text{Mini-batch}\hspace{0.17em}\text{SGD})=\mathrm{\Theta }(n^{\alpha +\beta }).\end{array}$$

This implies the following results:

$$\underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{SVRG})}=\underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{Mini-batch}\hspace{0.17em}\text{SGD})}=0,$$

as well as the stronger results:

$$\begin{array}{l}\underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{SVRG}))}\le \frac{1}{\tau },\\ \underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{Mini-batch}\hspace{0.17em}\text{SGD}))}\le \frac{1}{\tau (\alpha +\beta )}.\end{array}$$

Proposition 4 is proved in Section 5.4. It shows that LPI-GD can beat various algorithms interpolating between GD and SGD under the assumptions of Proposition 3.

We close this section by making several remarks. Firstly, we reiterate that although Propositions 3 and 4 illustrate regimes where LPI-GD beats GD, SGD, and their variants in oracle complexity, these results require the assumption that $d=O(\log \log (n))$, which stems from the curse of dimensionality. This small data dimension assumption can be restrictive for applications where d is large. However, because ERM is used in a diverse array of contexts, LPI-GD may be useful in settings where data dimension is inherently small or can be restricted to live in low-dimensional spaces via dimensionality reduction methods (see Section 1.3). We discuss potential avenues to improve this data dimension requirement in the last remark.

Secondly, we note that many refined gradient-based methods exist that lead to improvements in the dependence of oracle complexity bounds on the condition number—for example, SVRG (Johnson and Zhang 2013) and accelerated GD (Nesterov 1983) (which have similar oracle complexity expressions as GD, but with better dependence on condition number). We do not consider or optimize the trade-off with condition number in our analysis; this is, in many ways, an orthogonal question. As indicated at the outset of this paper, our focus has been on beating GD, SGD, and their variants’ oracle complexities with regard to their dependence on n, ϵ, d, and p.

Thirdly, although methods like SGD, SVRG, and mini-batch SGD match or improve upon GD’s oracle complexity (in terms of its trade-off between n and ϵ, neglecting implicit dependencies on p), they are stochastic optimization algorithms that only generate ϵ-approximate solutions on expectation. In contrast, the LPI-GD algorithm can beat GD and is completely deterministic.

Fourthly, because the LPI-GD algorithm has to learn the oracle at every iteration, the computational (or “arithmetical”) complexity of its iterations may seem to be higher than usual iterative methods. However, the motivation for focusing on the concept of oracle complexity in optimization theory is partly that computing gradients can be extremely hard for certain functions. So, the computational complexity of the oracle usually dominates any other reasonable computation done during iterations. We, therefore, adhere to this canonical paradigm in this work (Nesterov 2004). If the precise running time of the first-order oracle $\mathcal{O}$ is known, then the complete computational complexity of the LPI-GD algorithm can be derived from Proposition 2 and Theorem 2 by accounting for the running time of the D × D matrix inversions used to compute the interpolation weights. Note that, as mentioned earlier and shown in Algorithm 1, these weights only need to be precomputed once before running iterations of inexact GD.

Fifthly, we emphasize that Proposition 3 holds in the important regime where $\alpha =\frac{1}{2}$ and $\beta >1$. As noted in Section 1, $\alpha =\frac{1}{2}$ is typically the smallest accuracy of interest because it produces the scaling of the statistical error of the empirical risk. (Table 1 presents the oracle complexities computed in Proposition 3 in this regime.) The latter condition, $\beta >1$, ensures that the number of parameters p can be much larger than the number of training samples n, which corresponds to the popular overparametrized regime considered in the context of deep neural networks (cf. Poggio et al. 2020).

Finally, we note that the dependence of the bounds in Theorem 2 and Propositions 3 and 4 on both p and d could potentially be improved. The dependence on p arises because Theorem 1 holds for functions with codomain $\mathbb{R}$. However, the proof of Proposition 2 in Section 5.1 needs to bound the ${\mathcal{\ell }}^2$-norm error between (30) and the true gradient ${\nabla }_{\theta }F({\theta }^{(t-1)})$ for $t\in \mathbb{N}$. Translating the coordinate-wise guarantee to an ${\mathcal{\ell }}^2$-norm guarantee in this proof introduces the dependence on p (see (60)). If some version of Theorem 1 could be proved for functions with codomain ${\mathbb{R}}^p$, then we could possibly improve the (explicit) dependence of Theorem 2 on p.

On the other hand, the dependence on d could be improved by tightening the minimum eigenvalue lower bound in Lemma 1, as explained in Appendix C. Nevertheless, this eigenvalue lower bound must decay at least exponentially (and, correspondingly, the grid size parameter m must grow at least exponentially) in d, as noted after Theorem 1 in Section 3.1. Thus, it is intuitively straightforward to see that the best scaling of d we can hope for, so that the conclusions of Propositions 3 and 4 remain true, is $d=O(\log (n))$. This logarithmic scaling of data dimension d with number of training samples n is significant. Indeed, for several data science applications of interest, the useful information in a dataset of n high-dimensional feature vectors is contained within the pairwise (Euclidean) distances between different feature vectors. In these scenarios, standard low-dimensional approximate isometric embedding theorems—for example, the Johnson-Lindenstrauss lemma (cf. Dasgupta and Gupta 2003 and references therein)—ensure that the feature vectors can be mapped to a $d=O(\log (n))$ dimensional space such that the pairwise distances are preserved. We leave the development of our results along this direction of improving the conditioning analysis of local polynomial regression as future work. Moreover, we note that to beat a $d=O(\log (n))$ scaling, one would need to apply dimensionality reduction techniques, such as principal component analysis (Pearson 1901, Hotelling 1933) or other methods listed in Section 1.3. Theoretically analyzing dimensionality reduction simultaneously with nonparametric regression can be challenging. So, another possible future research direction is to incorporate and analyze (kernelized) principal component regression (Jolliffe 1982, Agarwal et al. 2021) in the context of gradient descent.

4. Proofs for Local Polynomial Interpolation

In this section, we provide proofs for the results presented in Section 3.1. Specifically, we first establish Proposition 1, then portray some properties of the interpolation weights in (23), and, finally, prove Theorem 1.

4.1. Proof of Proposition 1

To prove Proposition 1, we require the following intermediate lemma, which generalizes Tsybakov (2009, lemma 1.4) to our multivariate setting.

Lemma 1

(Auxiliary Löwner Lower Bound). When $l\ge d\ge 19$, the D × D symmetric matrix

$$\mathcal{B}\triangleq {\displaystyle {\int }_{{[-1,1]}^d}U(u)U{(u)}^{\text{T}}\text{d}u},$$

satisfies the Löwner lower bound

$$\mathcal{B}⪰\mathrm{\Lambda }(d,l)I_D,$$

where $U:{\mathbb{R}}^d\to {\mathbb{R}}^D$ is defined in (16), and the constant $\mathrm{\Lambda }(d,l)>0$ is defined in (20).

Proof.

We break this proof into several claims for clarity. For convenience in the sequel, we begin by establishing some simple facts about $\mathcal{B}$ in Claim 1 and Claim 2. Specifically, Claim 1 shows that $\mathcal{B}$ is positive definite.

Claim 1

(Positive Definiteness of $\mathcal{B}$). $\mathcal{B}$ is positive definite.

Proof of Claim 1.

Following the proof of Tsybakov (2009, lemma 1.4) mutatis mutandis, note that $\mathcal{B}$ is positive semidefinite:

$$\forall v\in {\mathbb{R}}^D,\hspace{1em}{v}^{\text{T}}\mathcal{B}v={\displaystyle {\int }_{{[-1,1]}^d}{(v^{\text{T}}U(u))}^2\text{d}u}\ge 0.$$

Suppose for the sake of contradiction that there exists a nonzero vector $v\in {\mathbb{R}}^D$ such that

$$v^{\text{T}}\mathcal{B}v={\displaystyle {\int }_{{[-1,1]}^d}{(v^{\text{T}}U(u))}^2\text{d}u}=0.$$

Then, we must have $v^{\text{T}}U(u)=0$ almost everywhere (with respect to the Lebesgue measure) on ${[-1,1]}^d$. However, because the nonzero polynomial $u\mapsto v^{\text{T}}U(u)$ is analytic, its set of roots has Lebesgue measure zero (see, e.g., Krantz and Parks 2002, section 4.1). This is a contradiction. So, we actually have $v^{\text{T}}\mathcal{B}v>0$ for all nonzero vectors $v\in {\mathbb{R}}^D$. Thus, $\mathcal{B}$ is positive definite. □

From hereon, our proof of Lemma 1 deviates greatly from that of Tsybakov (2009, lemma 1.4). We next compute the entries of $\mathcal{B}$ in Claim 2.

Claim 2

(Entries of $\mathcal{B}$). For every $r,s\in {\mathbb{Z}}_{+}^d$ with $|r|,|s|\le l$, we have

$${\mathcal{B}}_{r,s}=\mathrm{𝟙}\left\{r+s\hspace{0.17em}is\hspace{0.17em}\mathit{\text{even}}\hspace{0.17em}\mathit{\text{entry-wise}}\right\}\frac{2^d}{r!s!}{\displaystyle \prod _{j=1}^d\frac{1}{r_j+s_j+1}}.$$

Proof of Claim 2.

Because the rows and columns of $\mathcal{B}$ are indexed by $\{s\in {\mathbb{Z}}_{+}^d:|s|\le l\}$, for any $r,s\in {\mathbb{Z}}_{+}^d$ with $|r|,|s|\le l$, we have

$$\begin{array}{l}{\mathcal{B}}_{r,s}={\displaystyle {\int }_{{[-1,1]}^d}\frac{u^{r+s}}{r!s!}\text{d}u}\\ =\frac{2^d}{r!s!}{\displaystyle {\int }_{{[-1,1]}^d}\frac{1}{2^d}{\displaystyle \prod _{j=1}^d{u}_j^{r_j+s_j}}\text{d}u}\\ =\frac{2^d}{r!s!}{\displaystyle \prod _{j=1}^d\mathbb{E}\left[X_j^{r_j+s_j}\right]}\\ =\mathrm{𝟙}\left\{r+s\hspace{0.17em}\text{is}\hspace{0.17em}\text{even}\hspace{0.17em}\text{entry-wise}\right\}\frac{2^d}{r!s!}{\displaystyle \prod _{j=1}^d\frac{1}{r_j+s_j+1}},\end{array}$$

where $X_1,\dots ,X_d\in [-1,1]$ are independent and identically distributed (i.i.d.) uniform random variables, and the kth moment of X₁ is $\mathbb{E}[X_1^k]=0$ for all odd $k\in {\mathbb{Z}}_{+}$ and $\mathbb{E}[X_1^k]={(k+1)}^{-1}$ for all even $k\in {\mathbb{Z}}_{+}$ (cf. Forbes et al. 2011, chapter 40). □

Thus, Claim 1 and Claim 2 demonstrate that $\mathcal{B}$ is a doubly nonnegative matrix with all rational entries. To establish Lemma 1, we need to lower bound the strictly positive minimum eigenvalue of $\mathcal{B}$. To this end, using Claim 2, we first upper bound the trace of $\mathcal{B}$ in Claim 3, compute the determinant of $\mathcal{B}$ in Claim 4, and then lower bound the determinant of $\mathcal{B}$ in Claim 5.

Claim 3

(Upper Bound on Trace of $\mathcal{B}$).

$$\text{tr}(\mathcal{B})\le {(2e)}^d.$$

Proof of Claim 3.

To upper bound $\text{tr}(\mathcal{B})$, observe that

$$\begin{array}{l}\text{tr}(\mathcal{B})={\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}{\mathcal{B}}_{s,s}}\\ =2^d{\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{1}{{(s!)}^2}{\displaystyle \prod _{j=1}^d\frac{1}{2s_j+1}}}\\ \le 2^d{\displaystyle \sum _{k=0}^l\frac{1}{k!}{\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|=k}\frac{k!}{s!}}}\\ =2^d{\displaystyle \sum _{k=0}^l\frac{d^k}{k!}}\\ \le {(2e)}^d,\end{array}$$

where the second equality uses Claim 2, the fourth equality follows from the multinomial theorem, and the fifth inequality uses the Maclaurin series of the exponential function. □

Note that the exponential scaling of $\text{tr}(\mathcal{B})$ in Claim 3 is tight, because $\text{tr}(\mathcal{B})\ge {\mathcal{B}}_{(0,\dots ,0),(0,\dots ,0)}=2^d$. To lower bound $\det (\mathcal{B})$, we will first compute it exactly in the next claim.

Claim 4

(Determinant of $\mathcal{B}$).

$$\det (\mathcal{B})={\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{2^d{4}^{|s|}}{{(s!)}^2}{\displaystyle \prod _{i=1}^d\frac{1}{2s_i+1}}{\left(\begin{array}{c}2s_i\\ s_i\end{array}\right)}^{-2}}.$$

Proof of Claim 4.

To prove this, we introduce the Legendre polynomials, $\{L_k:[-1,1]\to \mathbb{R}|k\in {\mathbb{Z}}_{+}\}$, which are defined by Rodrigues’ formula (Chihara 2011, chapter V, equation (2.1)):

$$\forall t\in [-1,1],\hspace{1em}{L}_k(t)\triangleq \frac{\sqrt{k+\frac{1}{2}}}{2^{k}k!}\frac{{\text{d}}^k}{\text{d}{t}^k}{(t^2-1)}^k,$$(39)

where L_k is a polynomial with degree k.⁴ It is well-known that the family of Legendre polynomials satisfies the orthonormality property (Chihara 2011, chapter I, exercise 1.6):

$$\forall j,k\in {\mathbb{Z}}_{+},\hspace{1em}{\displaystyle {\int }_{-1}^1{L}_j(t)L_k(t)\text{d}t}=\mathrm{𝟙}\{j=k\}.$$(40)

Moreover, for each $k\in {\mathbb{Z}}_{+}$, the leading coefficient (of t^k) of $L_k(t)$ is known to be (Chihara 2011, chapter V, equation (2.7)):

$$\frac{\sqrt{k+\frac{1}{2}}}{2^k}\left(\begin{array}{c}2k\\ k\end{array}\right).$$(41)

Next, for any $r\in {\mathbb{Z}}_{+}^d$, define the Legendre product polynomial:

$$\forall u\in {[-1,1]}^d,\hspace{1em}{L}_r(u)\triangleq {\displaystyle \prod _{j=1}^d{L}_{r_j}(u_j)},$$(42)

which has degree r (with abuse of notation). It is straightforward to verify that these product polynomials also satisfy an orthonormality condition:

$$\forall r,s\in {\mathbb{Z}}_{+}^d,\hspace{1em}{\displaystyle {\int }_{{[-1,1]}^d}{L}_r(u)L_s(u)\text{d}u}={\displaystyle \prod _{j=1}^d{\displaystyle {\int }_{-1}^1{L}_{r_j}(u_j)L_{s_j}(u_j)\text{d}{u}_j}}=\mathrm{𝟙}\{r=s\},$$(43)

where we utilize (40). Moreover, for each $r\in {\mathbb{Z}}_{+}^d$, the leading coefficient (of u^r) of $L_r(u)$ is known to be

$${\displaystyle \prod _{j=1}^d\frac{\sqrt{r_j+\frac{1}{2}}}{2^{r_j}}\left(\begin{array}{c}2r_j\\ r_j\end{array}\right)}=\frac{1}{2^{|r|}}{\displaystyle \prod _{j=1}^d\sqrt{r_j+\frac{1}{2}}\left(\begin{array}{c}2r_j\\ r_j\end{array}\right)},$$(44)

where we utilize (41). Because of (43), we know that for any $r\in {\mathbb{Z}}_{+}^d$, the set $\{L_s:s\in {\mathbb{Z}}_{+}^d,\hspace{0.17em}s\le r\}$ is linearly independent and spans the set of all d-variate polynomials with degree at most r. (Note that we use the lexicographical order over indices in ${\mathbb{Z}}_{+}^d$ here. Furthermore, it is straightforward to verify that all d-variate monomials with nonzero coefficients that appear in L_s have degree at most s.) As a result, defining a vector(-valued) function $\mathcal{L}:{\mathbb{R}}^d\to {\mathbb{R}}^D$ as:

$$\forall u\in {[-1,1]}^d,\hspace{1em}\mathcal{L}(u)\triangleq {\left[L_s(u):s\in {\mathbb{Z}}_{+}^d,|s|\le l\right]}^{\text{T}},$$(45)

which is indexed in lexicographical order, we have the linear relation:

$$\forall u\in {[-1,1]}^d,\hspace{1em}U(u)=R\mathcal{L}(u),$$(46)

for some fixed matrix $R\in {\mathbb{R}}^{D\times D}$ (that does not depend on u). It is evident that this matrix R is lower triangular due to the lexicographical order that we have imposed. Furthermore, it is easy to derive the diagonal entries of R. Indeed, for any index $s\in {\mathbb{Z}}_{+}^d$ with $|s|\le l$, equating the leading coefficients on both sides of (46) produces

$$\frac{1}{s!}=R_{s,s}\frac{1}{2^{|s|}}{\displaystyle \prod _{j=1}^d\sqrt{\frac{2s_j+1}{2}}\left(\begin{array}{c}2s_j\\ s_j\end{array}\right)},$$

using (44). Rearranging this, we get:

$$\forall s\in {\mathbb{Z}}_{+}^d\hspace{0.17em}\text{with}\hspace{0.17em}|s|\le l,\hspace{1em}{R}_{s,s}=\frac{2^{|s|}{2}^{d/2}}{s!}{\displaystyle \prod _{j=1}^d\frac{1}{\sqrt{2s_j+1}}{\left(\begin{array}{c}2s_j\\ s_j\end{array}\right)}^{-1}}.$$(47)

Now, notice that

$$\begin{array}{l}\mathcal{B}={\displaystyle {\int }_{{[-1,1]}^d}U(u)U{(u)}^{\text{T}}\text{d}u}\\ =R\left({\displaystyle {\int }_{{[-1,1]}^d}\mathcal{L}(u)\mathcal{L}{(u)}^{\text{T}}\text{d}u}\right)R^{\text{T}}\\ =R{I}_D{R}^{\text{T}}\\ =R{R}^{\text{T}},\end{array}$$(48)

where the second equality follows from (46), and the third equality follows from (43). We recognize (48) as the Cholesky decomposition of $\mathcal{B}$ (Horn and Johnson 2013, corollary 7.2.9). Thus, using (47) and (48), we can write

$$\begin{array}{l}\det (\mathcal{B})={\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}{R}_{s,s}^2}\\ ={\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\underset{\le 1}{\underbrace{4^{|s|}{\displaystyle \prod _{i=1}^d{\left(\begin{array}{c}2s_i\\ s_i\end{array}\right)}^{-2}}}}·\underset{={\mathcal{B}}_{s,s}}{\underbrace{\frac{2^d}{{(s!)}^2}{\displaystyle \prod _{j=1}^d\frac{1}{2s_j+1}}}}}\\ ={\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{2^d{4}^{|s|}}{{(s!)}^2}{\displaystyle \prod _{i=1}^d\frac{1}{2s_i+1}}{\left(\begin{array}{c}2s_i\\ s_i\end{array}\right)}^{-2}},\end{array}$$

which is the desired result. □

We next establish a lower bound on the determinant of $\mathcal{B}$ using Claim 4.

Claim 5

(Lower Bound on Determinant of $\mathcal{B}$).

$$\det (\mathcal{B})\ge {\left(\frac{\pi d}{4e}\right)}^{dD}\frac{1}{{(l+d)}^{3lD}}.$$

Proof of Claim 5.

Observe from Claim 4 that we can write

$$\begin{array}{l}\det (\mathcal{B})={\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{2^d{4}^{|s|}}{{(s!)}^2}{\displaystyle \prod _{i=1}^d\frac{1}{2s_i+1}}{\left(\begin{array}{c}2s_i\\ s_i\end{array}\right)}^{-2}}\\ \ge {\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{{(\pi /2)}^d}{4^{|s|}{(s!)}^2}\underset{\ge 1}{\underbrace{{\displaystyle \prod _{j=1}^d\frac{4s_j+1}{2s_j+1}}}}}\\ \ge {\left(\frac{\pi }{2}\right)}^{dD}{\left({\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{1}{2^{|s|}s!}}\right)}^2\\ ={\left(\frac{\pi }{2}\right)}^{dD}\frac{1}{4^{{\sum }_{i=0}^{l}i\left(\begin{array}{c}i+d-1\\ d-1\end{array}\right)}}{\left({\displaystyle \prod _{j=0}^l{\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|=j}\frac{1}{s_1!\cdots s_d!}}}\right)}^2\\ \ge {\left(\frac{\pi }{2e^2}\right)}^{dD}{\left(\frac{e^2}{4}\right)}^{{\displaystyle \sum _{i=0}^{l}i\left(\begin{array}{c}i+d-1\\ d-1\end{array}\right)}}{\displaystyle \prod _{j=0}^l{\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|=j}\frac{1}{{(s_1+1)}^{2s_1+1}\cdots {(s_d+1)}^{2s_d+1}}}}\\ \ge {\left(\frac{\pi }{4e}\right)}^{dD}{\displaystyle \prod _{j=0}^l{\displaystyle \prod _{s\in {\mathbb{Z}}_{+}^d:|s|=j}\frac{1}{{(s_1+1)}^{2s_1+1}\cdots {(s_d+1)}^{2s_d+1}}}}\\ \ge {\left(\frac{\pi }{4e}\right)}^{dD}{\displaystyle \prod _{j=0}^l{\left(\frac{d^d}{{(j+d)}^d{(j+1)}^{2j}}\right)}^{\left(\begin{array}{c}j+d-1\\ d-1\end{array}\right)}}\\ \ge {\left(\frac{\pi d}{4e}\right)}^{dD}\frac{1}{{(l+d)}^{{\sum }_{j=0}^l(2j+d)\left(\begin{array}{c}j+d-1\\ d-1\end{array}\right)}}\\ \ge {\left(\frac{\pi d}{4e}\right)}^{dD}\frac{1}{{(l+d)}^{3lD}},\end{array}$$

where the second inequality follows from the following bound (on Catalan numbers; see Dutton and Brigham 1986, theorem):

$$\forall j\in {\mathbb{Z}}_{+},{\left(\begin{array}{c}2j\\ j\end{array}\right)}^2<\frac{4^{2j+1}}{\pi (4j+1)},$$

the fourth equality follows from adding the exponents $|s|$ of two over all $s\in {\mathbb{Z}}_{+}^d$ with $|s|\le l$ and then squaring, the fifth inequality follows from the Stirling’s formula bound (see, e.g., Feller 1968, chapter II, section 9, equation (9.15)):⁵

$$\forall j\in {\mathbb{Z}}_{+},\hspace{1em}j!\le {(j+1)}^{j+\frac{1}{2}}{e}^{-j+1},$$

the sixth inequality follows from the straightforward bound:

$${\displaystyle \sum _{i=0}^{l}i\left(\begin{array}{c}i+d-1\\ d-1\end{array}\right)}\ge l\left(\begin{array}{c}l+d-1\\ d-1\end{array}\right)=\frac{dl}{l+d}\left(\begin{array}{c}l+d\\ d\end{array}\right)\ge \frac{d}{2}\left(\begin{array}{c}l+d\\ d\end{array}\right)=\frac{dD}{2},$$

which uses the fact that $l\ge d$, the seventh inequality uses the following basic facts and calculations:

1. The hypercube maximizes the volume over all d-orthotopes with an isoperimetric constraint:⁶

   $$\underset{\stackrel{s_1,\dots ,s_d\ge 0:}{s_1+\cdots +s_d=j}}{\text{max}}(s_1+1)\cdots (s_d+1)={\left(\frac{j}{d}+1\right)}^d,$$
   where $j\in {\mathbb{Z}}_{+}$ and the maximum is achieved by $s_1=\hspace{0.17em}\cdots \hspace{0.17em}=s_d=\frac{j}{d}$;

2. For any $j\in {\mathbb{Z}}_{+}$ and any $s\in {\mathbb{Z}}_{+}^d$ with $|s|=j$, the vector $(j,0,\dots ,0)\in {\mathbb{Z}}_{+}^d$ majorizes the vector s, and because the map $0\le t\mapsto 2t \log (t+1)$ is convex, Karamata’s majorization inequality yields (cf. Marshall et al. 2011, chapter 1, section A)

   $${\displaystyle \sum _{i=1}^{d}2s_i \log (s_i+1)}\le 2j \log (j+1),$$
   which implies that
   $${\displaystyle \prod _{i=1}^d{(s_i+1)}^{2s_i}}\le {(j+1)}^{2j},$$
   using the monotonicity of the exponential function;

   and the ninth inequality follows from the simple bound:

   $${\displaystyle \sum _{j=0}^l(2j+d)\left(\begin{array}{c}j+d-1\\ d-1\end{array}\right)}\le (2l+d){\displaystyle \sum _{j=0}^l\left(\begin{array}{c}j+d-1\\ d-1\end{array}\right)}=(2l+d)D\le 3lD,$$
   which uses the fact that $l\ge d$. □

Finally, we can complete the proof of Lemma 1 using Claim 1, Claim 3, and Claim 5. Recall the following well-known result (cf. Merikoski and Virtanen 1997, lemma 1):

which utilizes the positive definiteness of $\mathcal{B}$ established earlier in Claim 1.⁷ Combining Claim 3, Claim 5, and (49), we obtainwhere the third inequality holds because $d\ge 19$, the map $1\le t\mapsto t^t$ is strictly increasing, and we utilize the ensuing standard bounds on binomial coefficients:and the final equality follows from (20). Because $\mathcal{B}⪰{\lambda }_{\mathsf{\text{min}}}(\mathcal{B})I_D$, applying the lower bound in (50) completes the proof of Lemma 1. □

We make three pertinent remarks at this point. Firstly, we note that when $d\le l\le \gamma d$ for some constant $\gamma >1$, one can verify that

Hence, $\mathrm{\Lambda }{(d,l)}^{-1}$ scales double exponentially in d:

This observation will be used in the proof of Proposition 3 in Section 5.3. Secondly, we also conjecture that the scaling in (52) can be significantly improved, as shown in Appendix C. Thirdly, we remark that the condition $d\ge 19$, which appears in many of the results in Section 3, originates from the proof of Lemma 1 (see (50)). This condition is simply stated in our results for precision; it becomes immaterial when d scales with n in, for example, Proposition 3.

We are now in a position to establish Proposition 1 using Lemma 1.

Proof of Proposition 1.

We take inspiration from the proof technique of Tsybakov (2009, lemma 1.5), but the details of our analysis are more involved and quite different. Observe that for all $m\in \mathbb{N}$, for all $x\in {[h,1-h]}^d$, and for all $v\in {\mathbb{R}}^D$ with $\Vert v{\Vert }_2=1$, we have

$$\begin{array}{l}{v}^{\text{T}}B(x)v=v^{\text{T}}\left(\frac{1}{{(mh)}^d}{\displaystyle \sum _{y\in {\mathcal{G}}_m}U\left(\frac{y-x}{h}\right)U{\left(\frac{y-x}{h}\right)}^{\text{T}}{\displaystyle \prod _{j=1}^{d}K\left(\frac{y_j-x_j}{h}\right)}}\right)v\\ =v^{\text{T}}\left(\frac{1}{{(mh)}^d}{\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}U(z)U{(z)}^{\text{T}}{\displaystyle \prod _{j=1}^{d}K(z_j)}}\right)v\\ \ge b^d{v}^{\text{T}}\tilde{B}(x)v,\end{array}$$

where the first equality follows from (19), the second equality follows from the substitution $y=x+hz$ and the fact that the kernel K has support $[-1,1]$, the third inequality holds because the kernel K is lower bounded by b > 0 on $[-1,1]$, and we define $\tilde{B}(x)$ as the D × D symmetric positive semidefinite matrix

$$\tilde{B}(x)\triangleq {\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}U(z)U{(z)}^{\text{T}}\frac{1}{{(mh)}^d}}.$$

Then, applying the Courant-Fischer-Weyl min-max theorem (Horn and Johnson 2013, theorem 4.2.6), we obtain

$${\lambda }_{\mathsf{\text{min}}}(B(x))\ge b^d{\lambda }_{\mathsf{\text{min}}}(\tilde{B}(x)).$$(53)

We proceed to lower bounding ${\lambda }_{\mathsf{\text{min}}}(\tilde{B}(x))$. Fix any $x\in {[h,1-h]}^d$, and any $r,s\in {\mathbb{Z}}_{+}^d$ with $|r|,|s|\le l$, and consider the Riemann sum

$$\tilde{B}{(x)}_{r,s}=\frac{1}{r!s!}{\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}\frac{z^{r+s}}{{(mh)}^d}}.$$

Because $x\in {[h,1-h]}^d$, it is easy to verify that the points $\{z\in {[-1,1]}^d:x+hz\in {\mathcal{G}}_m\}$ partition the hypercube ${[-1,1]}^d$ into a grid of d-cells $\{I_z\subseteq {[-1,1]}^d:z\in {[-1,1]}^d,x+hz\in {\mathcal{G}}_m\}$ with side lengths ${(mh)}^{-1}$ and volume ${(mh)}^{-d}$ each, where I_z is the d-cell with tag z. Moreover, because the polynomial function $\gamma :{[-1,1]}^d\to \mathbb{R},\hspace{0.17em}\gamma (u)=u^{r+s}$ is Riemann integrable, if $h=\omega (m^{-1})$ with all other problem parameters fixed, then the above Riemann sum converges to an entry of the matrix $\mathcal{B}$ defined in Lemma 1 (cf. Apostol 1974, section 14.3):

$$\underset{m\to \infty }{\lim }\frac{1}{r!s!}{\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}\frac{z^{r+s}}{{(mh)}^d}}={\displaystyle {\int }_{{[-1,1]}^d}\frac{u^{r+s}}{r!s!}\text{d}u}={\mathcal{B}}_{r,s}.$$

Therefore, we have shown that ${\lim }_{m\to \infty }\tilde{B}{(x)}_{r,s}={\mathcal{B}}_{r,s}$ when $h=\omega (m^{-1})$. We will require an explicit upper bound on the rate of this convergence in the sequel. To this end, let us first compute the Lipschitz constant of γ with respect to the ${\mathcal{\ell }}^{\infty }$-norm. For any $i\in [d]$ and any fixed $(u_1,\dots ,u_{i-1},u_{i+1},\dots ,u_d)\in {[-1,1]}^{d-1}$, notice that for every $u_i,{\tilde{u}}_i\in [-1,1]$,

$$\begin{array}{l}|\gamma (u)-\gamma (u_1,\dots ,u_{i-1},{\tilde{u}}_i,u_{i+1},\dots ,u_d)|\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=|u_1^{r_1+s_1}\cdots u_{i-1}^{r_{i-1}+s_{i-1}}{u}_{i+1}^{r_{i+1}+s_{i+1}}\cdots u_d^{r_d+s_d}||u_i^{r_i+s_i}-{\tilde{u}}_i^{r_i+s_i}|\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le (r_i+s_i)|u_i-{\tilde{u}}_i|.\end{array}$$

Hence, via repeated application of the triangle inequality, we have

$$\begin{array}{l}\forall u,\tilde{u}\in {[-1,1]}^d,\hspace{1em}|\gamma (u)-\gamma (\tilde{u})|=|u^{r+s}-{\tilde{u}}^{r+s}|\le {\displaystyle \sum _{i=1}^d(r_i+s_i)|u_i-{\tilde{u}}_i|}\\ \le (|r|+|s|){\Vert u-\tilde{u}\Vert }_{\infty }.\end{array}$$(54)

Now, observe that for all $r,s\in {\mathbb{Z}}_{+}^d$ with $|r|,|s|\le l$,

$$\begin{array}{l}|\tilde{B}{(x)}_{r,s}-{\mathcal{B}}_{r,s}|=\frac{1}{r!s!}\left|{\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}\frac{z^{r+s}}{{(mh)}^d}-{\displaystyle {\int }_{I_z}{u}^{r+s}\text{d}u}}\right|\\ \le \frac{1}{r!s!}{\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}\left|\frac{z^{r+s}}{{(mh)}^d}-{\displaystyle {\int }_{I_z}{u}^{r+s}\text{d}u}\right|}\\ \le \frac{1}{r!s!}{\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}\frac{1}{{(mh)}^d}\left(\underset{u\in I_z}{\sup } u^{r+s}-\underset{\tilde{u}\in I_z}{\inf } {\tilde{u}}^{r+s}\right)}\\ \le \frac{|r|+|s|}{\mathit{\text{mhr}}!s!}{\displaystyle \sum _{\stackrel{z\in {[-1,1]}^d:}{x+hz\in {\mathcal{G}}_m}}\frac{1}{{(mh)}^d}}\\ \le \frac{|r|+|s|}{\mathit{\text{mhr}}!s!}{\left(2+\frac{1}{mh}\right)}^d\\ \le \frac{3^d(|r|+|s|)}{\mathit{\text{mhr}}!s!},\end{array}$$

where the second inequality follows from the triangle inequality, the fourth inequality follows from (54), the fifth inequality holds because $|\{z\in {[-1,1]}^d:x+hz\in {\mathcal{G}}_m\}|\le {(2mh+1)}^d$ (because the interval $[-1,1]$ can have at most $2mh+1$ points at distance ${(mh)}^{-1}$ from each other), and the sixth inequality holds because we assume that $m\in \mathbb{N}$ and $h\in \left(0,\frac{1}{2}\right)$ are such that $mh\ge 4l{(3e)}^d/\mathrm{\Lambda }(d,l)\ge 1$. This implies that for all such m and h,

$$\begin{array}{l}{\Vert \tilde{B}(x)-\mathcal{B}\Vert }_{\text{F}}\le \sqrt{{\displaystyle \sum _{r,s\in {\mathbb{Z}}_{+}^d:|r|,|s|\le l}\frac{9^d{(|r|+|s|)}^2}{m^2{h}^2{(r!s!)}^2}}}\\ \le \frac{2\hspace{0.17em}l{3}^d}{mh}{\displaystyle \sum _{r\in {\mathbb{Z}}_{+}^d:|r|\le l}\frac{1}{{(r!)}^2}}\\ \le \frac{2\hspace{0.17em}l{3}^d}{mh}{\displaystyle \sum _{j=0}^l\frac{1}{j!}{\displaystyle \sum _{r\in {\mathbb{Z}}_{+}^d:|r|=j}\frac{j!}{r!}}}\\ =\frac{2\hspace{0.17em}l{3}^d}{mh}{\displaystyle \sum _{j=0}^l\frac{d^j}{j!}}\\ \le \frac{2\hspace{0.17em}l{(3e)}^d}{mh},\end{array}$$(55)

where the fourth equality follows from the multinomial theorem, and the fifth inequality uses the Maclaurin series of the exponential function. Next, recall the following consequence of the Wielandt-Hoffman-Mirsky inequality (cf. Horn and Johnson 2013, corollary 6.3.8):

$$|{\lambda }_{\mathsf{\text{min}}}(\tilde{B}(x))-{\lambda }_{\mathsf{\text{min}}}(\mathcal{B})|\le {\Vert \tilde{B}(x)-\mathcal{B}\Vert }_{\text{F}}.$$(56)

Then, combining (55) and (56) produces the lower bound

$${\lambda }_{\mathsf{\text{min}}}(\tilde{B}(x))\ge {\lambda }_{\mathsf{\text{min}}}(\mathcal{B})-\frac{2\hspace{0.17em}l{(3e)}^d}{mh},$$

and employing Lemma 1 to this bound yields

$${\lambda }_{\mathsf{\text{min}}}(\tilde{B}(x))\ge \mathrm{\Lambda }(d,l)-\frac{2\hspace{0.17em}l{(3e)}^d}{mh},$$(57)

where $\mathrm{\Lambda }(d,l)$ is defined in (20).

Finally, because m and h satisfy $mh\ge 4l{(3e)}^d/\mathrm{\Lambda }(d,l)$, or, equivalently,

$$\frac{2\hspace{0.17em}l{(3e)}^d}{mh}\le \frac{\mathrm{\Lambda }(d,l)}{2},$$

using (53) and (57), we get

$${\lambda }_{\mathsf{\text{min}}}(B(x))\ge b^d{\lambda }_{\mathsf{\text{min}}}(\tilde{B}(x))\ge \frac{b^d\mathrm{\Lambda }(d,l)}{2},$$

for all $x\in {[h,1-h]}^d$. Because $B(x)⪰{\lambda }_{\mathsf{\text{min}}}(B(x))I_D$, the above lower bound on ${\lambda }_{\mathsf{\text{min}}}(B(x))$ completes the proof. □

4.2. Properties of Local Polynomial Interpolation Weights

To prove Theorem 1, we will require two useful lemmata regarding the interpolation weights in (23). The first lemma generalizes Tsybakov (2009, proposition 1.12) to the d-variate setting and demonstrates the intuitively pleasing property that polynomials are reproduced exactly by the local polynomial interpolator given in (24).

Lemma 2

(Polynomial Reproduction (Tsybakov 2009)). Consider any d-variate polynomial $g:{[0,1]}^d\to \mathbb{R}$ with total degree at most l (i.e., for every monomial of g with nonzero coefficient, the sum of the exponents of the d variables is at most l). Then, for every $m\in \mathbb{N}$ and every $x\in {[h,1-h]}^d$ such that B(x) is positive definite, we have

$$\widehat{\varphi }(x)={\displaystyle \sum _{y\in {\mathcal{G}}_m}g(y)w_y^{*}(x)}=g(x),$$

where the interpolation weights are given in (23).

Proof.

We follow the proof of Tsybakov (2009, proposition 1.12) mutatis mutandis. Fix any $m\in \mathbb{N}$ and $x\in {[h,1-h]}^d$ such that B(x) is positive definite, and define the vector $q(x)\in {\mathbb{R}}^D$ as

$$q(x)\triangleq {\left[h^{|s|}{\nabla }^{s}g(x):s\in {\mathbb{Z}}_{+}^d,|s|\le l\right]}^{\text{T}}.$$

Then, for every $y\in {\mathcal{G}}_m$, we have

$$\begin{array}{l}g(y)={\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{{\nabla }^{s}g(x)}{s!}{(y-x)}^s}\\ =q{(x)}^{\text{T}}U\left(\frac{y-x}{h}\right),\end{array}$$

where the first equality uses Taylor’s theorem (cf. Apostol 1974, theorem 12.14) and the fact that g is a polynomial with total degree at most l, and U is defined in (16). This implies that the solution to the weighted regression problem in (17) satisfies

$$\begin{array}{l}\widehat{\mathrm{\Phi }}(x)=\underset{\varphi \in {\mathbb{R}}^D}{\text{arg min}}{\displaystyle \sum _{y\in {\mathcal{G}}_m}{(q(x)-\varphi )}^{\text{T}}U\left(\frac{y-x}{h}\right)U{\left(\frac{y-x}{h}\right)}^{\text{T}}(q(x)-\varphi ){\displaystyle \prod _{j=1}^{d}K\left(\frac{y_j-x_j}{h}\right)}}\\ =\underset{\varphi \in {\mathbb{R}}^D}{\text{arg min}}{(q(x)-\varphi )}^{\text{T}}B(x)(q(x)-\varphi )\\ =q(x),\end{array}$$

where the second equality follows from (19), and the third equality holds because B(x) is positive definite. Therefore, by considering the first coordinate, we get

$$\widehat{\varphi }(x)={\displaystyle \sum _{y\in {\mathcal{G}}_m}g(y)w_y^{*}(x)}=g(x),$$

where the first equality uses (24). This completes the proof. □

The second lemma generalizes Tsybakov (2009, lemma 1.3) to the d-variate setting and portrays an ${\mathcal{\ell }}^1$-norm bound on the sequence of interpolation weights in (23).

Lemma 3

(${\mathcal{\ell }}^1$-Norm Bound). Suppose that $l\ge d\ge 19$. Then, for all $m\in \mathbb{N}$ and $h\in \left(0,\frac{1}{2}\right)$ such that $mh\ge 4l{(3e)}^d/\mathrm{\Lambda }(d,l)$ and all $x\in {[h,1-h]}^d$, we have

$${\displaystyle \sum _{y\in {\mathcal{G}}_m}|w_y^{*}(x)|}\le \frac{2}{\mathrm{\Lambda }(d,l)}{\left(\frac{3\sqrt{e}\hspace{0.17em}c}{b}\right)}^d,$$

where $c\ge b>0$ are the bounds on our kernel, and $\mathrm{\Lambda }(d,l)$ is defined in (20).

Proof.

Once again, we follow the proof of Tsybakov (2009, lemma 1.3) mutatis mutandis. First, notice that B(x) satisfies the Löwner lower bound in Proposition 1 because the conditions of the proposition are satisfied. Then, observe that for every $x\in {[h,1-h]}^d$, we have

$$\begin{array}{l}{\displaystyle \sum _{y\in {\mathcal{G}}_m}|w_y^{*}(x)|}\le {\displaystyle \sum _{y\in {\mathcal{G}}_m}{\Vert w_y(x)\Vert }_2}\\ =\frac{1}{{(mh)}^d}{\displaystyle \sum _{y\in {\mathcal{G}}_m}{\Vert B{(x)}^{-1}U\left(\frac{y-x}{h}\right)\Vert }_2{\displaystyle \prod _{j=1}^{d}K\left(\frac{y_j-x_j}{h}\right)}}\\ \le \frac{c^d}{{(mh)}^d}{\displaystyle \sum _{y\in {\mathcal{G}}_m}{\Vert B{(x)}^{-1}U\left(\frac{y-x}{h}\right)\Vert }_2\mathrm{𝟙}\left\{\forall j\in [d],\hspace{0.17em}|y_j-x_j|\le h\right\}}\\ \le \frac{2c^d}{{(\mathit{\text{mhb}})}^d\mathrm{\Lambda }(d,l)}{\displaystyle \sum _{y\in {\mathcal{G}}_m}{\Vert U\left(\frac{y-x}{h}\right)\Vert }_2\mathrm{𝟙}\left\{\forall j\in [d],\hspace{0.17em}|y_j-x_j|\le h\right\}}\\ \le \frac{2}{\mathrm{\Lambda }(d,l)}{\left(\frac{\sqrt{e}\hspace{0.17em}c}{\mathit{\text{mhb}}}\right)}^d{\displaystyle \sum _{y\in {\mathcal{G}}_m}\mathrm{𝟙}\left\{\forall j\in [d],\hspace{0.17em}|y_j-x_j|\le h\right\}}\\ \le \frac{2}{\mathrm{\Lambda }(d,l)}{\left(\frac{\sqrt{e}\hspace{0.17em}c}{b}\right)}^d{\left(\frac{2mh+1}{mh}\right)}^d\\ \le \frac{2}{\mathrm{\Lambda }(d,l)}{\left(\frac{3\sqrt{e}\hspace{0.17em}c}{b}\right)}^d,\end{array}$$

where the first inequality follows from (23) and the Cauchy-Schwarz inequality, the second equality follows from (22), the third inequality holds because the kernel K is upper bounded by c > 0 on the support $[-1,1]$, the fourth inequality follows from Proposition 1, the fifth inequality follows from the bound

$$\begin{array}{l}{\Vert U\left(\frac{y-x}{h}\right)\Vert }_2^2={\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}{\left(\frac{{(y-x)}^s}{h^{|s|}s!}\right)}^2}\\ \le {\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|\le l}\frac{1}{{(s!)}^2}}\\ \le {\displaystyle \sum _{i=0}^l\frac{1}{i!}{\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|=i}\frac{i!}{s!}}}\\ ={\displaystyle \sum _{i=0}^l\frac{d^i}{i!}}\\ \le e^d,\end{array}$$

which uses (16), the fact that $|y_j-x_j|\le h$ for all $j\in [d]$, the multinomial theorem, and the Maclaurin series of the exponential function, the sixth inequality holds because the interval $[-h,h]$ can have at most $2mh+1$ points at distance $m^{-1}$ from each other, and the seventh inequality follows from the assumption that $mh\ge 4l{(3e)}^d/\mathrm{\Lambda }(d,l)\ge 1$. This completes the proof. □

4.3. Proof of Theorem 1

Finally, we are in a position to prove the supremum norm interpolation guarantee in Theorem 1.

Proof of Theorem 1.

Suppose that $m\in \mathbb{N}$ and $h\in \left(0,\frac{1}{2}\right)$ are such that $mh\ge 4l{(3e)}^d/\mathrm{\Lambda }(d,l)$. Then, observe that for every $x\in {[h,1-h]}^d$, we have

$$\begin{array}{l}\widehat{\varphi }(x)-g(x)={\displaystyle \sum _{y\in {\mathcal{G}}_m}g(y)w_y^{*}(x)}-g(x)\\ ={\displaystyle \sum _{y\in {\mathcal{G}}_m}{w}_y^{*}(x)(g(y)-g(x))}\\ ={\displaystyle \sum _{y\in {\mathcal{G}}_m}{w}_y^{*}}(x){\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:1\le |s|\le l}\frac{{\nabla }^{s}g(x)}{s!}{(y-x)}^s}\\ \hspace{1em}\hspace{1em}\hspace{1em}+w_y^{*}(x){\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|=l}\frac{{\nabla }^{s}g(x+\tau (y-x))-{\nabla }^{s}g(x)}{s!}{(y-x)}^s}\\ ={\displaystyle \sum _{y\in {\mathcal{G}}_m}{w}_y^{*}(x){\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|=l}\frac{{\nabla }^{s}g(x+\tau (y-x))-{\nabla }^{s}g(x)}{s!}{(y-x)}^s}},\end{array}$$

where the first equality follows from (24), the second equality holds because

$${\displaystyle \sum _{y\in {\mathcal{G}}_m}{w}_y^{*}(x)}=1,$$

which follows from Lemma 2 applied to the constant unit polynomial, the third equality follows from Taylor’s theorem with Lagrange remainder term (Apostol 1974, theorem 12.14):

$$g(y)-g(x)={\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:1\le |s|\le l}\frac{{\nabla }^{s}g(x)}{s!}{(y-x)}^s}+{\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|=l}\frac{{\nabla }^{s}g(x+\tau (y-x))-{\nabla }^{s}g(x)}{s!}{(y-x)}^s},$$

which holds for some $\tau \in (0,1)$, and the fourth equality holds because

$${\displaystyle \sum _{y\in {\mathcal{G}}_m}{(y-x)}^s{w}_y^{*}(x)}=0,$$

for all $s\in {\mathbb{Z}}_{+}^d$ with $1\le |s|\le l$, which follows from Lemma 2 applied to the polynomials $u\mapsto {(u-x)}^s$. This implies that for every $x\in {[h,1-h]}^d$,

$$\begin{array}{l}|\widehat{\varphi }(x)-g(x)|\le {\displaystyle \sum _{y\in {\mathcal{G}}_m}|w_y^{*}(x)|{\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|=l}\frac{|{\nabla }^{s}g(x+\tau (y-x))-{\nabla }^{s}g(x)|}{s!}|{(y-x)}^s|}}\\ \le L_2{\displaystyle \sum _{y\in {\mathcal{G}}_m}|w_y^{*}(x)|{\Vert y-x\Vert }_1^{\eta -l}{\displaystyle \sum _{s\in {\mathbb{Z}}_{+}^d:|s|=l}\frac{1}{s!}|{(y-x)}^s|}}\\ =\frac{L_2}{l!}{\displaystyle \sum _{y\in {\mathcal{G}}_m}|w_y^{*}(x)|{\Vert y-x\Vert }_1^{\eta }}\\ =\frac{L_2}{l!}{\displaystyle \sum _{y\in {\mathcal{G}}_m}|w_y^{*}(x)|{\Vert y-x\Vert }_1^{\eta }\mathrm{𝟙}\left\{\forall j\in [d],\hspace{0.17em}|y_j-x_j|\le h\right\}}\\ \le \frac{2L_2{(dh)}^{\eta }}{l!\mathrm{\Lambda }(d,l)}{\left(\frac{3\sqrt{e}\hspace{0.17em}c}{b}\right)}^d\\ =\frac{2L_2}{l!\mathrm{\Lambda }(d,l)}{\left(\frac{3\sqrt{e}\hspace{0.17em}c}{b}\right)}^d{\left(\frac{4dl{(3e)}^d}{m\mathrm{\Lambda }(d,l)}\right)}^{\eta },\end{array}$$

where the first inequality follows from the triangle inequality, the second inequality follows from the Hölder class assumption (see (11)) and the fact that $\tau \in (0,1)$, the third equality follows from the multinomial theorem, the fourth equality follows from (22) and (23) because the kernel K has support $[-1,1]$, the fifth inequality follows from Lemma 3, and the sixth equality holds because we set

$$h=\frac{4l{(3e)}^d}{m\mathrm{\Lambda }(d,l)},$$(58)

which minimizes the upper bound in the fifth inequality and achieves equality in the condition imposed on m and h at the outset of the proof. Lastly, to ensure that ${\sup }_{x\in {[h,1-h]}^d}|\widehat{\varphi }(x)-g(x)|\le \delta$, it suffices to impose the condition

$$\frac{2L_2}{l!\mathrm{\Lambda }(d,l)}{\left(\frac{3\sqrt{e}\hspace{0.17em}c}{b}\right)}^d{\left(\frac{4dl{(3e)}^d}{m\mathrm{\Lambda }(d,l)}\right)}^{\eta }\le \delta ,$$

which is equivalent to

$$m\ge \frac{4dl{(3e)}^d}{\mathrm{\Lambda }(d,l)}{\left(\frac{2L_2}{l!\mathrm{\Lambda }(d,l)}\right)}^{1/\eta }{\left(\frac{3\sqrt{e}\hspace{0.17em}c}{b}\right)}^{d/\eta }{\left(\frac{1}{\delta }\right)}^{1/\eta },$$

which, in turn, is implied by the condition

$$m\ge \frac{12\sqrt{e}(2L_2+1)c}{b}\left(\frac{dl{(3e)}^d}{{(l!)}^{1/\eta }\mathrm{\Lambda }{(d,l)}^2}\right){\left(\frac{1}{\delta }\right)}^{1/\eta },$$(59)

where we use the facts that $\eta >l\ge d\ge 19$ (by assumption) and $\mathrm{\Lambda }{(d,l)}^{(\eta +1)/\eta }\ge \mathrm{\Lambda }{(d,l)}^2$ (because $\mathrm{\Lambda }(d,l)\ge 1$), and we slacken the $\mathrm{\Lambda }(d,l)$ term because this does not change its double exponential nature. Moreover, we may further simplify (59) by noting that

$${(l!)}^{1/\eta }\ge \frac{{(2\pi )}^{1/(2\eta )}{l}^{(l/\eta )+(1/(2\eta ))}}{e^{l/\eta }}\ge \frac{l}{e{l}^{\frac{1}{l+1}}}\ge \frac{l}{2e},$$

where we use the Stirling’s formula bound (cf. Feller 1968, chapter II, section 9, equation (9.15)) and the fact that $l^{1/(l+1)}\le 2$ for $l\ge 19$. Applying this bound, we obtain the sufficient condition for (59) presented in the theorem statement:

$$m\ge \frac{110(2L_2+1)c}{b}\left(\frac{d{(3e)}^d}{\mathrm{\Lambda }{(d,l)}^2}\right){\left(\frac{1}{\delta }\right)}^{1/\eta },$$

where we use the fact that $24e\sqrt{e}\le 110$. Finally, plugging this bound into (58) yields the upper bound on h in the theorem statement. This completes the proof. □

We note that this argument is inspired by the proof of Tsybakov (2009, proposition 1.13), but the details are quite different.

5. Proofs of Convergence Analysis

In this section, we prove Propositions 2–4 and Theorem 2 from Section 3.2.2.

5.1. Proof of Proposition 2

To establish Proposition 2, we require the following lemma from the literature, Friedlander and Schmidt (2012, lemma 2.1) (also see So and Zhou 2017, theorem 2, equation (19)), that upper bounds the difference between the ERM objective function (3) evaluated at ${\theta }^{(T)}$, which is the output of our proposed algorithm, and the true infimum in (4). Although Friedlander and Schmidt (2012, lemma 2.1) holds for general inexact gradient descent methods, we unwind the recursion in Friedlander and Schmidt (2012, lemma 2.1) and adapt it to our setting below for convenience.

Lemma 4

(Inexact Gradient Descent Bound (Friedlander and Schmidt 2012)). Suppose the ERM objective function $F:{\mathbb{R}}^p\to \mathbb{R}$ is continuously differentiable and μ-strongly convex, and its gradient ${\nabla }_{\theta }F:{\mathbb{R}}^p\to {\mathbb{R}}^p$ is L₁-Lipschitz continuous. Then, for any $T\in \mathbb{N}$, we have

$$0\le F({\theta }^{(T)})-F_{*}\le {\left(1-\frac{1}{\sigma }\right)}^T(F({\theta }^{(0)})-F_{*})+\frac{1}{2L_1}{\displaystyle \sum _{t=1}^T{\left(1-\frac{1}{\sigma }\right)}^{T-t}{\Vert {\widehat{\nabla F}}^{(t)}-{\nabla }_{\theta }F({\theta }^{(t-1)})\Vert }_2^2},$$

where ${\theta }^{(t)}$ for $t\in {\mathbb{Z}}_{+}$ are the inexact gradient descent updates given by (31) with an arbitrary initialization ${\theta }^{(0)},\hspace{0.17em}{\widehat{\nabla F}}^{(t)}$ for $t\in \mathbb{N}$ are our local polynomial interpolation-based approximations of the true gradients ${\nabla }_{\theta }F({\theta }^{(t-1)})$ shown in (30), and $\sigma =L_1/\mu \ge 1$ is the condition number defined in Section 2.2.

We remark that Friedlander and Schmidt (2012, lemma 2.1) additionally assumes that the infimum $F_{*}$ in (4) can be achieved. However, this condition need not be imposed in Lemma 4, because if F is strongly convex on ${\mathbb{R}}^m$, its infimum $F_{*}$ is finite and achieved by some unique global minimizer ${\theta }_{*}\in {\mathbb{R}}^p$. We next prove Proposition 2 using Lemma 4 and Theorem 1.

Proof of Proposition 2.

We first verify that the conditions of Lemma 4 hold. Recall the ERM objective function $F:{\mathbb{R}}^p\to \mathbb{R}$ given in (3):

$$\forall \theta \in {\mathbb{R}}^p,\hspace{1em}F(\theta )=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}f(x^{(i)};\theta )},$$

where $\{x^{(i)}\in {[h^{\prime },1-h^{\prime }]}^d:i\in [n]\}$ is some fixed training data, and $f:{[0,1]}^d\times {\mathbb{R}}^p\to \mathbb{R}$ is a given loss function satisfying the assumptions in Section 2.2. Specifically, because $f(x^{(i)};·):{\mathbb{R}}^p\to \mathbb{R}$ is μ-strongly convex for all $i\in [n]$, F is also μ-strongly convex due to Nesterov (2004, lemma 2.1.4). Likewise, because ${\nabla }_{\theta }f(x^{(i)};·):{\mathbb{R}}^p\to {\mathbb{R}}^p$ is L₁-Lipschitz continuous for all $i\in [n]$, F is continuously differentiable and ${\nabla }_{\theta }F:R^p\to {\mathbb{R}}^p$ is also L₁-Lipschitz continuous. Indeed, observe that for all ${\theta }_1,{\theta }_2\in {\mathbb{R}}^p$,

$$\begin{array}{l}{\Vert {\nabla }_{\theta }F({\theta }_1)-{\nabla }_{\theta }F({\theta }_2)\Vert }_2\le \frac{1}{n}{\displaystyle \sum _{i=1}^n{\Vert {\nabla }_{\theta }f(x^{(i)};{\theta }_1)-{\nabla }_{\theta }f(x^{(i)};{\theta }_2)\Vert }_2}\\ \le \frac{L_1}{n}{\displaystyle \sum _{i=1}^n{\Vert {\theta }_1-{\theta }_2\Vert }_2}\\ =L_1{\Vert {\theta }_1-{\theta }_2\Vert }_2,\end{array}$$

where we use the triangle inequality.

Next, we apply Lemma 4 and get

$$F({\theta }^{(T)})-F_{*}\le {\left(1-\frac{1}{\sigma }\right)}^T(F({\theta }^{(0)})-F_{*})+\frac{1}{2L_1}{\displaystyle \sum _{t=1}^T{\left(1-\frac{1}{\sigma }\right)}^{T-t}{\Vert {\widehat{\nabla F}}^{(t)}-{\nabla }_{\theta }F({\theta }^{(t-1)})\Vert }_2^2},$$

where ${\theta }^{(t)}$ for $t\in {\mathbb{Z}}_{+}$ are our updates in (31), and ${\widehat{\nabla F}}^{(t)}$ for $t\in \mathbb{N}$ are our approximations of the true gradients in (30). To upper bound the squared ${\mathcal{\ell }}^2$-norms on the right-hand side of this inequality, notice that for any $t\in \mathbb{N}$, we have

$$\begin{array}{l}{\Vert {\widehat{\nabla F}}^{(t)}-{\nabla }_{\theta }F({\theta }^{(t-1)})\Vert }_2\le \frac{1}{n}{\displaystyle \sum _{j=1}^n\sqrt{{\displaystyle \sum _{i=1}^p{\left({\widehat{\varphi }}_i^{(t)}(x^{(j)})-\frac{\partial f}{\partial {\theta }_i}(x^{(j)};{\theta }^{(t-1)})\right)}^2}}}\\ \le \sqrt{{\displaystyle \sum _{i=1}^p\underset{x\in {[h^{\prime },1-h^{\prime }]}^d}{\sup }{|{\widehat{\varphi }}_i^{(t)}(x)-\frac{\partial f}{\partial {\theta }_i}(x;{\theta }^{(t-1)})|}^2}}\\ \le \sqrt{p}\hspace{0.17em}{\delta }_t,\end{array}$$(60)

where the first inequality follows from (30) and the triangle inequality, ${\widehat{\varphi }}_i^{(t)}(x)$ denotes the output of our local polynomial interpolator for $\frac{\partial f}{\partial {\theta }_i}(·;{\theta }^{(t-1)})$ when evaluated at x (cf. (29)), and the third inequality follows from Theorem 1 and our choices of grid sizes (26) and bandwidths (27) in Section 3.2.1 (where (25) ensures that $h_t\le h^{\prime }$). Hence, we obtain

$$F({\theta }^{(T)})-F_{*}\le {\left(1-\frac{1}{\sigma }\right)}^T(F({\theta }^{(0)})-F_{*})+\frac{p}{2L_1}{\displaystyle \sum _{t=1}^T{\left(1-\frac{1}{\sigma }\right)}^{T-t}{\delta }_t^2}.$$(61)

We now select appropriate values for the constants $\{{\delta }_t\in (0,1):t\in \mathbb{N}\}$. For every $t\in \mathbb{N}$, choose

$${\delta }_t={\left(1-\frac{1}{\sigma }\right)}^{\frac{T}{2}},$$

as indicated in (33). This produces

$$\begin{array}{l}F({\theta }^{(T)})-F_{*}\le {\left(1-\frac{1}{\sigma }\right)}^T\left(F({\theta }^{(0)})-F_{*}+\frac{p}{2L_1}{\displaystyle \sum _{t=1}^T{\left(1-\frac{1}{\sigma }\right)}^{T-t}}\right)\\ \le {\left(1-\frac{1}{\sigma }\right)}^T\left(F({\theta }^{(0)})-F_{*}+\frac{p}{2L_1}{\displaystyle \sum _{t=0}^{\infty }{\left(1-\frac{1}{\sigma }\right)}^t}\right)\\ ={\left(1-\frac{1}{\sigma }\right)}^T\left(F({\theta }^{(0)})-F_{*}+\frac{p\sigma }{2L_1}\right),\end{array}$$

where the last equality computes the value of the previous geometric series.

Finally, to ensure that ${\theta }^{(T)}$ is an ϵ-approximate solution of (4), viz.

$$F({\theta }^{(T)})-F_{*}\le ϵ,$$

as defined in (13), it suffices to impose the condition

$${\left(1-\frac{1}{\sigma }\right)}^T\left(F({\theta }^{(0)})-F_{*}+\frac{p\sigma }{2L_1}\right)\le ϵ.$$

By taking logarithms of both sides, and rearranging and simplifying the resulting inequality, we obtain the equivalent inequality

$$T\ge {\left(\log \left(\frac{\sigma }{\sigma -1}\right)\right)}^{-1} \log \left(\frac{F({\theta }^{(0)})-F_{*}+\frac{p}{2\mu }}{ϵ}\right),$$

where we also use the fact that $\sigma =L_1/\mu$. Therefore, after running the LPI-GD algorithm for $T\in \mathbb{N}$ iterations with T given by (32), the updated parameter vector ${\theta }^{(T)}$ forms an ϵ-approximate solution to (4). We also note that for every $t\in \mathbb{N}$,

$${\delta }_t={\left(1-\frac{1}{\sigma }\right)}^{\frac{T}{2}}=\mathrm{\Theta }\left(\sqrt{\frac{ϵ}{p}}\right),$$

by substituting (32) into our choice of δ_t above. This yields (33). Hence, the proof is complete. □

5.2. Proof of Theorem 2

We next establish the oracle complexity bound in Theorem 2 using Proposition 2.

Proof of Theorem 2.

Because we choose the maximum allowable approximation errors $\{{\delta }_t\in (0,1):t\in [T]\}$ to be constant with respect to t as in (33), a single common uniform grid ${\mathcal{G}}_m$ is used in all iterations of the LPI-GD algorithm, where m = m_t is given by (26) (and is also constant with respect to t). The size of this grid is upper bounded by

$$\begin{array}{l}|{\mathcal{G}}_m|=m^d\\ \le {\left(\frac{110(2L_2+1)c}{b}\left(\frac{d{(3e)}^d}{\mathrm{\Lambda }{(d,l)}^2}\right){\left(\frac{\sigma }{\sigma -1}\right)}^{T/(2\eta )}+1\right)}^d\\ \le {\left(\frac{110(2L_2+1)c}{b}\left(\frac{d{(3e)}^d}{\mathrm{\Lambda }{(d,l)}^2}\right){\left(\frac{\sigma }{\sigma -1}\right)}^{1/(2\eta )}{\left(\frac{F({\theta }^{(0)})-F_{*}+\frac{p}{2\mu }}{ϵ}\right)}^{1/(2\eta )}+1\right)}^d\\ \le {\left(\frac{220(2L_2+1)c}{b}\right)}^d\left(\frac{d^d{(3e)}^{d^2}}{\mathrm{\Lambda }{(d,l)}^{2d}}\right){\left(1+\frac{\sigma }{2(L_1-\mu )}\right)}^{d/(2\eta )}{\left(\frac{p+2\mu (F({\theta }^{(0)})-F_{*})}{ϵ}\right)}^{d/(2\eta )}\\ \le \left(1+\frac{\sigma }{2(L_1-\mu )}\right){\left(\frac{220(2L_2+1)c}{b}\right)}^d\left(\frac{d^d{(3e)}^{d^2}}{\mathrm{\Lambda }{(d,l)}^{2d}}\right){\left(\frac{p+2\mu (F({\theta }^{(0)})-F_{*})}{ϵ}\right)}^{d/(2\eta )},\end{array}$$

where the second inequality follows from (26) and (33) in Proposition 2, the third inequality follows from (32) in Proposition 2, the fourth inequality follows from straightforward manipulations and the fact that $\sigma =L_1/\mu$, and the fifth inequality holds because we have assumed that $\eta >d$. The LPI-GD algorithm makes $|{\mathcal{G}}_m|$ oracle calls in each of its T iterations and produces an ϵ-approximate solution ${\theta }^{(T)}\in {\mathbb{R}}^p$ (in the sense of (13)), as argued in Proposition 2. Hence, the oracle complexity of this algorithm is upper bounded by

$$\begin{array}{l}\mathrm{\Gamma }(\text{LPI-GD})\\ \hspace{1em}\le T|{\mathcal{G}}_m|\\ \hspace{1em}\le \left(\frac{\sigma +2(L_1-\mu )}{(L_1-\mu )\log \left(\frac{\sigma }{\sigma -1}\right)}\right)\log \left(\frac{p+\mathrm{\Delta }}{2\mu ϵ}\right){\left(\frac{220(2L_2+1)c}{b}\right)}^d\left(\frac{d^d{(3e)}^{d^2}}{\mathrm{\Lambda }{(d,l)}^{2d}}\right){\left(\frac{p+\mathrm{\Delta }}{ϵ}\right)}^{d/(2\eta )}\\ \hspace{1em}=C_{\mu ,L_1,L_2,b,c}(d,l){\left(\frac{p+\mathrm{\Delta }}{ϵ}\right)}^{d/(2\eta )}\log \left(\frac{p+\mathrm{\Delta }}{2\mu ϵ}\right),\end{array}$$

where we let $\mathrm{\Delta }=2\mu (F({\theta }^{(0)})-F_{*})$ for convenience, the second inequality uses (32) and the fact that $ϵ\le \frac{(L_1-\mu )p}{2L_1\mu }$ (which is equivalent to $\frac{\sigma }{\sigma -1}\le \frac{p}{2\mu ϵ}$), and the third equality uses (34). This completes the proof. □

We remark that the proofs of Proposition 2 and Theorem 2 can be combined to determine the optimal t-dependent approximation error parameters $\{{\delta }_t\in (0,1):t\in \mathbb{N}\}$ and number of iterations $T\in \mathbb{N}$ such that ${\theta }^{(T)}$ is an ϵ-approximate solution in the sense of (13) with as few oracle calls as possible. However, this does not qualitatively improve the dominant term in the scaling law of the derived oracle complexity bound in Theorem 2.

5.3. Proof of Proposition 3

We next derive Proposition 3 from Theorem 2.

Proof of Proposition 3.

We first establish the scaling of $C_{\mu ,L_1,L_2,b,c}(d,l)$ in (34). Because $d=O(\log \log (n))$, observe that

$$\left(\frac{\sigma +2(L_1-\mu )}{(L_1-\mu ) \log \left(\frac{\sigma }{\sigma -1}\right)}\right){\left(\frac{220(2L_2+1)c}{b}\right)}^d=O(\text{polylog}(n)),$$(62)

and because $l\ge d$, observe that

$$\frac{d^d{(3e)}^{d^2}}{\mathrm{\Lambda }{(d,l)}^{2d}}=O\left({(l+d)}^{6dl{\left(\frac{e(l+d)}{d}\right)}^d}\right),$$

as d grows, where we use (20). To simplify the second scaling expression above, notice that for all sufficiently large $d\in \mathbb{N}$, we have

$$\begin{array}{l}\log \log \left(\frac{d^d{(3e)}^{d^2}}{\mathrm{\Lambda }{(d,l)}^{2d}}\right)\le \log \log (l+d)+\log (dl)+d \log \left(\frac{e(l+d)}{d}\right)+\mathrm{\Theta }(1)\\ =\mathrm{\Theta }(\log (d))+d \log \left(\frac{e(l+d)}{d}\right)\\ \le \mathrm{\Theta }(\log (d))+d \log (e(\gamma +1))\\ \le 2 \log (e(\gamma +1))d,\end{array}$$

where the second equality follows from the fact that $\eta =\mathrm{\Theta }(d)$, and the third inequality holds because $l\le \gamma d$. This implies that for all sufficiently large $d\in \mathbb{N}$,

$$\frac{d^d{(3e)}^{d^2}}{\mathrm{\Lambda }{(d,l)}^{2d}}\le \exp \left(e^{2 \log (e(\gamma +1))d}\right).$$(63)

Hence, because $d\le \frac{\log \log (n)}{4 \log (e(\gamma +1))}$, we obtain

$$\begin{array}{l}{C}_{\mu ,L_1,L_2,b,c}(d,l)\le O(\text{polylog}(n))·\exp \left(e^{\frac{1}{2}\log \log (n)}\right)\\ =O(\text{polylog}(n))·\exp \left(\sqrt{\log (n)}\right)\\ \le \exp \left(2\sqrt{\log (n)}\right),\end{array}$$

for all sufficiently large $n\in \mathbb{N}$, using (34), (62), and (63). Next, applying (36), which follows from Theorem 2, we have that for all sufficiently large $n\in \mathbb{N}$,

$$\begin{array}{l}\mathrm{\Gamma }(\text{LPI-GD})\le c_1 \exp \left(2\sqrt{\log (n)}\right){\left(\frac{p}{ϵ}\right)}^{d/(2\eta )}\log \left(\frac{p}{ϵ}\right)\\ \le c_1 \exp \left(2\sqrt{\log (n)}\right){\left(c_2{n}^{\alpha +\beta }\right)}^{d/(2\eta )}\log \left(c_2{n}^{\alpha +\beta }\right)\\ \le c_3\log (n) \exp \left(2\sqrt{\log (n)}\right)n^{1/\tau },\end{array}$$(64)

where $c_1,c_2,c_3>0$ are some positive constants, the second inequality follows from the assumptions that $ϵ=\mathrm{\Theta }(n^{-\alpha })$ and $p=\mathrm{\Theta }(n^{\beta })$, and the third inequality holds because $\eta \ge \tau (\alpha +\beta )d/2$. This proves the scaling of $\mathrm{\Gamma }(\text{LPI-GD})$ with respect to n.

Finally, from the relations (37) and (38) and the assumptions that $ϵ=\mathrm{\Theta }(n^{-\alpha })$ and $p=\mathrm{\Theta }(n^{\beta })$, we get ${\mathrm{\Gamma }}_{*}(\text{GD})=\mathrm{\Theta }(n \log (n))$ and ${\mathrm{\Gamma }}_{*}(\text{SGD})=\mathrm{\Theta }(n^{\alpha +\beta })$. Therefore, for some constants $c_4,c_5>0$, we obtain

$$\begin{array}{l}\underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{GD})}\le \underset{n\to \infty }{\lim }\frac{c_3\exp \left(2\sqrt{\log (n)}\right)n^{1/\tau }}{c_{4}n}=0,\\ \underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{SGD})}\le \underset{n\to \infty }{\lim }\frac{c_3 \log (n)\exp \left(2\sqrt{\log (n)}\right)n^{1/\tau }}{c_5{n}^{\alpha +\beta }}=0,\end{array}$$

which use the fact that $\tau >\text{max}\{1,{(\alpha +\beta )}^{-1}\}$. Furthermore, we also obtain the stronger results

$$\begin{array}{l}\underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{GD}))}\le \underset{n\to \infty }{\lim }\frac{\log (c_3)+\log \log (n)+2\sqrt{\log (n)}+{\tau }^{-1}\log (n)}{\log (c_4)+\log (n)+\log \log (n)}=\frac{1}{\tau }<1,\\ \underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{SGD}))}\le \underset{n\to \infty }{\lim }\frac{\log (c_3)+\log \log (n)+2\sqrt{\log (n)}+{\tau }^{-1}\log (n)}{\log (c_5)+(\alpha +\beta )\log (n)}=\frac{1}{\tau (\alpha +\beta )}<1,\end{array}$$

which imply the previous limits.⁸ This completes the proof. □

5.4. Proof of Proposition 4

Lastly, we derive Proposition 4 from Theorem 2.

Proof of Proposition 4.

This proof closely follows the argument for Proposition 3. Under the assumptions of the proposition statement, for some positive constant $c_3>0$ and all sufficiently large $n\in \mathbb{N}$, we have from (64) that

$$\mathrm{\Gamma }(\text{LPI-GD})\le c_3 \log (n) \exp \left(2\sqrt{\log (n)}\right)n^{1/\tau }.$$

Next, from the existing first-order oracle complexity results in Propositions A.3 and A.4 in Appendix A, we define the upper bounds:

$$\begin{array}{l}\mathrm{\Gamma }(\text{SVRG})\le {\mathrm{\Gamma }}_{*}(\text{SVRG})\triangleq \mathrm{\Theta }\left(n \log \left(\frac{p}{ϵ}\right)\right),\\ \mathrm{\Gamma }(\text{Mini-batch}\hspace{0.17em}\text{SGD})\le {\mathrm{\Gamma }}_{*}(\text{Mini-batch}\hspace{0.17em}\text{SGD})\triangleq \mathrm{\Theta }\left(\frac{p}{ϵ}\right).\end{array}$$

Using the assumptions $ϵ=\mathrm{\Theta }(n^{-\alpha })$ and $p=\mathrm{\Theta }(n^{\beta })$, we obtain that ${\mathrm{\Gamma }}_{*}(\text{SVRG})=\mathrm{\Theta }(n \log (n))$ and ${\mathrm{\Gamma }}_{*}(\text{Mini-batch}\hspace{0.17em}\text{SGD})=\mathrm{\Theta }(n^{\alpha +\beta })$. Therefore, as in the proof of Proposition 3, we have

$$\underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{SVRG})}=\underset{n\to \infty }{\lim }\frac{\mathrm{\Gamma }(\text{LPI-GD})}{{\mathrm{\Gamma }}_{*}(\text{Mini-batchSGD})}=0.$$

Furthermore, we also obtain the stronger results

$$\begin{array}{l}\underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{SVRG}))}\le \frac{1}{\tau }<1,\\ \underset{n\to \infty }{\text{lim sup}}\frac{\log (\mathrm{\Gamma }(\text{LPI-GD}))}{\log ({\mathrm{\Gamma }}_{*}(\text{Mini-batch}\hspace{0.17em}\text{SGD}))}\le \frac{1}{\tau (\alpha +\beta )}<1.\end{array}$$

This completes the proof. □

6. Conclusion

We conclude by summarizing our main contributions and reiterating some directions of future research. In order to exploit the smoothness of loss functions in data for optimization purposes, we proposed the LPI-GD algorithm (see Algorithm 1) to compute approximate solutions of the ERM problem in (3) and (4), where the loss function is strongly convex and smooth in both the parameter and the data. We then derived the iteration and oracle complexities of this algorithm and illustrated that its oracle complexity beats the oracle complexities of GD, SGD, and their variants for a broad range of dependencies between the parameter dimension p, the approximation accuracy ϵ, and the number of training samples n, when the data dimension d is sufficiently small. Finally, in the course of our convergence analysis, we also provided a careful and detailed analysis of multivariate local polynomial interpolation with supremum norm guarantees, which may be of independent interest in nonparametric statistics.

We close our discussion by stating four avenues for future work. First, as stated at the end of Section 3 and in Appendix C, it would be very valuable to improve the dependence of Theorem 2 on the data dimension d by, for example, sharpening the minimum eigenvalue lower bound in the key technical estimate in Lemma 1 (see Conjecture C.1 in Appendix C). Second, the explicit dependence of Theorem 2 on the parameter dimension p could be improved by establishing new interpolation guarantees for approximating smooth ${\mathbb{R}}^p$-valued functions using local polynomial regression or other methods. Third, as noted in Section 1, although our work focuses on the strongly convex loss function setting, similar analyses could be carried out for other conventional scenarios—for example, when $f(x;·):{\mathbb{R}}^p\to \mathbb{R}$ is nonconvex and ϵ-approximate solutions to (4) are defined via approximate first-order stationary points. Last, the LPI-GD method presented in Algorithm 1 is really a theoretical method because it depends on parameters that are often unknown in practice, such as Lipschitz and strong convexity constants. It would therefore be interesting to create a practical version of this method and empirically investigate any average-case gains it provides in running time.