Near-Linear Time Local Polynomial Nonparametric Estimation with Box Kernels
Yining Wang Yi Wu Simon S. Du





























Near-Linear Time Local Polynomial Nonparametric Estimation


Near-Linear Time Local Polynomial Nonparametric Estimation with Box Kernels








Yining Wang
Warrington College of Business, University of Florida, Gainesville, FL 32611, USA.
Yi Wu
Institute of Interdisciplinary Information Sciences, Tsinghua University, Beijing, 100084, China.
Simon S. Du
Paul G. Allen School of Computer Science and Engineering, University of Washington, WA 98195, USA.

 


	Local polynomial regression <cit.> is an important class of methods
	for nonparametric density estimation and regression problems.
	However, straightforward implementation of local polynomial regression has quadratic time complexity which hinders its applicability in large-scale data analysis.
	In this paper, we significantly accelerate the computation of local polynomial estimates
	by novel applications of multi-dimensional binary indexed trees <cit.>.
	Both time and space complexity of our proposed algorithm is nearly linear in the number of input data points.
	Simulation results confirm the efficiency and effectiveness of our proposed approach.







local polynomial regression, nonparametric density estimation, binary indexed trees, hashing



The overlooked role of stellar variability in LMC intermediate-age clusters

    Ricardo Salinas
    December 30, 2023
=============================================================================












§ INTRODUCTION


Big data analytics has become essential for modern operations research and operations management applications <cit.>.
Statistics methods, such as nonparametric density and function estimation <cit.>,
play important roles in the prediction of econometric trends like loan management and market profit prediction <cit.>,
as well as exploratory data analysis as preliminary steps to operations management problems, such as the understanding and estimation of customers' demand or preferences in
revenue management applications.
Naturally, big data scenarios pose unique challenges from computing aspects, as many classical statistical computational routines (especially nonparametric statistical methods)
are computationally heavy and not scalable to large amounts of data. 

The focus of this paper is on nonparametric density estimation and regression problems.
Let 𝒳⊆ℝ^d be a compact domain in ℝ^d,
which is conventionally taken to be the unit cube [0,1]^d for convenience.
In the nonparametric density estimation problem, n independent samples X_1,⋯,X_n∈𝒳
are obtained as

    X_1,⋯,X_n  i.i.d.∼  f_0,

where f_0 is the density of an unknown distribution P_0 supported on 𝒳.
In the nonparametric regression problem, on each of n “design points” X_1,⋯,X_n∈𝒳
a “response” or “measurement” Y_i is obtained according to the model

    Y_i = m_0(X_i) + ξ_i,     i=1,⋯ n,

where m_0:𝒳→ℝ is an unknown regression model of interest,
and {ξ_i} are independent sub-Gaussian noise variables satisfying 𝔼[ξ_i|X_i] = 0.
The objective in both nonparametric estimation problems is to construct estimates f̂_n or m̂_n such that 
the mean-square error (MSE) 

    ∫_𝒳|f̂_n(x)-f_0(x)|^2 x   or  ∫_𝒳|m̂_n(x)-m_0(x)|^2 x

is minimized.










The nonparametric density estimation and regression problems have a long history of study, dating back to the 1920s <cit.>.
A large family of methods have been developed and their properties analyzed,
including kernel smoothing <cit.>,
spline smoothing <cit.>,
wavelet smoothing <cit.>
and local polynomial regression <cit.>.

In this paper, we concentrate on the local polynomial regression method introduced in <cit.>.
Compared to its competitors, local polynomial regression has the benefits of being adaptive to non-uniform design densities and regression functions of unbounded values <cit.>.
In Sec. <ref> we give a succinct description of these methods and also summarize some of their basic properties.

While the statistical properties of local polynomial regression are well understood, computational efficiency has been less studied.
In particular, straightforward implementation of local polynomial regression typically has quadratic time complexity O(sn)
to calculate f̂_n(z_i) or m̂_n(z_i) on  s “testing points” z_1,⋯,z_s∈𝒳,
which is prohibitively slow for analysis of large data sets.
We give an illustrative example in Figure <ref>, which shows that the computational bottleneck of the naive implementation of local polynomial regression
significantly restricts the number of observations (training points) it processes under given time budget, leading to inaccurate function estimates.

The apparent difficulties of computationally efficient local polynomial regresstion motivate the following question:

	Given n training data points {(X_i,Y_i)}_i=1^n and s testing points z_1,⋯,z_s∈𝒳,
	can we compute local polynomial estimates f̂_n(z_1),⋯,f̂_n(z_s) or m̂_n(z_1),⋯,m̂_n(z_s) in 
	O((n+s)log n) time?
	

Essentially, Question <ref> concerns estimation algorithms whose time complexity is nearly linear in the total number of training and testing points n+s,
apart from possible poly-logarithmic factors.
On the other hand, a linear dependency on n+s is clearly necessary as the number of inputs is already on the order of n+s.

We give an affirmative answer to Question <ref> in this paper by accelerating local polynomial regression
using efficient data structures.
In particular, applying a multi-dimensional extension of binary indexed trees <cit.>
we are able to reduce the computation time on each testing point from O(n) to O(log^d n).
Furthermore, to avoid space complexity growing exponentially with dimension d,
we consider a lazy allocation strategy that allocates memory on the fly with the help of a Hash table,
which achieves near-linear space complexity.
Our computations of all local polynomial estimates are exact, and therefore all statistical properties of local polynomial regression
are automatically preserved.

Finally, we remark that certain existing nonparametric estimation methods can also achieve near-linear time complexity,
such as the B-spline <cit.> and/or falling factorial basis <cit.> construction in smoothing splines
and block thresholding of wavelet coefficients <cit.>.
However, both smoothing spline and wavelet thresholding (shrinkage) methods become very complicated for multi-dimensional and non-uniformly spaced
data points.
For example, the extension of smoothing splines to higher dimension requires complicated constructions of the “thin-plate” splines,
and in the case of wavelet smoothing the basis functions are complicated and difficult to compute when design points are not evenly spaced.
On the other hand, the local polynomial regression method easily handles both cases of multi-dimensionality
and unevenly spaced design points without much additional complexity.






















§ LOCAL POLYNOMIAL REGRESSION


In this section we give a succinct description of local polynomial regression <cit.> in the context of nonparametric estimation problems,
with high-level summarization of its properties.

Let k∈ℕ be the desired polynomial degree used in estimation.
Practical choices of k are small non-negative integers, such as k=0 (local mean averaging),
k=1 (local linear regression) and k=2 (local quadratic regression).
For any z∈ℝ^d, define the polynomial mapping ψ_z:ℝ^d→ℝ^D,
D=1+d+⋯+d^k as

    ψ_z(x)    := [1, x_1-z_1,⋯,x_d-z_d,(x_1-z_1)^2, (x_1-z_1)(x_2-z_2), ⋯, (x_d-z_d)^2,
       ⋯, (x_1-z_1)^k, ⋯, (x_d-z_d)^k].

At a higher level, the polynomial mapping ψ_x(·) is an expansion of basis of degree-k polynomials on ℝ^d,
appropriately offset so that ψ_z(z) = [1,0,⋯,0].
As we shall see immediately, the mapping ψ_z plays a fundamental role in local polynomial regression methods.


We first consider the regression problem Y_i=m_0(X_i)+ξ_i.
We start by introducing the concepts of kernel functions and bandwidths, which are crucial to almost every nonparametric estimation method:

	
  * A kernel function K:ℝ→ℝ that satisfies constraints ∫_-∞^∞ K(u) u = 1, ∫_-∞^∞ K(u)u u=0 and 
	∫_-∞^∞ K^2(u) u < ∞ is used to “average” neighboring data points in a certain manner.
	Common kernel functions include the box kernel K(u)=𝕀[|u|≤ 1/2], the Gaussian kernel K(u)=1/√(2π)e^-u^2/2 and the 
	Epanechnikov kernel K(u)=3/4√(5)max(0, 1-u^2/5).
	
	
  * A bandwidth parameter h>0 is used to control the size of the “neighborhood” of a certain test point z∈𝒳
	in which training points are incorporated and smoothed.
	A small bandwidth h will have fewer effective training points rendering the variance of the resulting estimate large,
	while a large bandwidth h tends to “over-smooth” the unknown function over a large domain
	and therefore leading to larger estimation bias.


To estimate the value of the unknown regression function m_0 at a specific testing point z, first solve
a weighted least-squares problem
[Although in the literature the ℓ_2 norm z-X_i_2=√(∑_j=1^d(z_j-X_ij)^2) is more commonly used,
	replacing the ℓ_2 norm with ℓ_∞ does not affect the statistical properties of θ̂_n,
	because ℓ_2 and ℓ_∞ are equivalent norms when dimension d is fixed.]

    θ̂_n = min_θ∈ℝ^D∑_i=1^n(Y_i-θ^⊤ψ_z(X_i))^2 K(z-X_i_∞/h)

where z-X_i_∞ = max_1≤ j≤ d|z_j-X_ij|.
Afterwards, m̂_n(z) is taken to be θ̂_n^⊤ψ_z(z), which corresponds to the first component of θ̂_n.


	When K(·) is taken to be the box kernel K(u)=𝕀[|u|≤ 1/2], 
	Eq. (<ref>) can be re-formulated as
	
    ĝ_n = min_g∈𝒫_k∑_X_i∈ B_h(z) (Y_i-g(X_i))^2

	where 𝒫_k is the class of all degree-k polynomials on ℝ^d and B_h(z)={u∈𝒳: u-z_∞≤ h/2}
	denotes the neighborhood of z with radius h.
	The re-formulation of ĝ_n is a clear interpretation of the polynomial fitting (in 𝒫_k) and the local properties (in B_h(z)).


The density estimation problem is slightly more involved than local polynomial regression.
Here we adopt the approach in <cit.> by re-casting the density estimation problem as a regression problem.
For every X_i∈𝒳 define empirical cumulative density function (CDF)

    F̂(X_i) := 1/n∑_j=1^n𝕀[X_j1≤ X_i1, ⋯, X_jd≤ X_id]

that approximates the true CDF F(X_i) = ∫_𝒳𝕀[u_1≤ X_i1,⋯,u_d≤ X_id]f_0(u) u.
Afterwards, a local polynomial fit of F̂ is computed as

    β̂_n = min_β∈ℝ^D∑_i=1^n(F̂(X_i)-θ^⊤ψ_z(X_i))^2 K(z-X_i_∞/h).

Because f_0 = ∂^d F/∂ x_1⋯∂ x_d,
an estimate f̂_n(z) can be obtained by reading off the component in β̂_n corresponding to (z_1-x_1)(z_2-x_2)⋯(z_d-x_d)
multiplied by d!.
Note that k≥ d is required for this density estimation method.

Local polynomial regression enjoys a wide range of desired properties.
When bandwidths h are properly tuned (for example by cross-validation), 
local polynomial regression is minimax optimal when the underlying model f_0 or m_0
belongs to Hölder or Sobolev smoothness classes.
In addition, when bandwidths are selected locally, the resulting estimation is optimal even for spatially heterogeneous functions <cit.>.
The local polynomial regression estimator also adapts to non-uniform and non-smooth design densities and boundaries <cit.>,
properties that do not hold for kernel smoothing estimators.




§ OUR METHOD




We describe our proposed method for local polynomial regression with the box kernel that achieves
O((n+s)log^d n) time complexity and O(nlog^d n) space complexity,
which is the main contribution of this paper.

We start by formulating sufficient statistics required to compute local polynomial estimates,
and giving a simple nearly linear-time algorithm for the univariate case d=1 as a warm-up exercise.
We then proceed with a formal discretization argument for higher dimensions, and explain how binary indexed trees
can be applied to enable fast calculations.
Finally, a lazy memory allocation scheme is employed to avoid memory explosion in high-dimensional binary indexed trees.



 §.§ Sufficient statistics of local polynomial regression


When the box kernel K(u)=𝕀[|u|≤ 1/2] is used, 
the local polynomial regression formulation in Eq. (<ref>) is an ordinary least squares (OLS)
problem on ℬ(z,h) := {(X_i,Y_i): X_i-z_∞≤ h/2}, where z∈𝒳 is the testing point.
The solution θ̂_n admits the following closed form:

    θ̂_n = [∑_ℬ(z,h)ψ_z(X_i)ψ_z(X_i)^⊤]^-1[∑_ℬ(z,h)Y_iψ_z(X_i)].


Recall that ψ_z(X_i) is a multivariate polynomial function in X_i. 
The sufficient statistics required to compute Eq. (<ref>) can then be organized as follows:

	
  * ∑_(X_i,Y_i)∈ℬ(z,h) X_i1^j_1X_i2^j_2⋯ X_id^j_d for j_1+⋯+j_d≤ 2k;
	
  * ∑_(X_i,Y_i)∈ℬ(z,h) Y_i X_i1^j_1⋯ X_id^j_d for j_1+⋯+j_d≤ k.


Let M be the total number of terms summarized above.
For each ℓ∈{1,2,⋯,M}, define T_ℓ(X_i,T_i) to be either X_i1^j_1X_i2^j_2⋯ X_id^j_d
or Y_i X_i1^j_1⋯ X_id^j_d, where j_1+⋯+j_d≤ 2k or j_1+⋯+j_d≤ k are indices corresponding to term ℓ≤ M.
The terms T_ℓ(X_i,Y_i) for ℓ≤ M are then sufficient statistics required to compute the local polynomial estimate θ̂_n
in Eq. (<ref>).


Once the sufficient statistics 𝒯_ℓ(X_i,Y_i), (X_i,Y_i)∈ℬ(z,h) are accumulated for a specific testing point z∈𝒳 and bandwidth h>0,
the estimate θ̂_n and subsequently m̂_n(z) or f̂_n(z) can be computed in O(D^3) time and O(D^2) space
via Eq. (<ref>), both treated as constants under fixed-dimension (d) settings.

Thus, the computational question reduces to efficient computation of sufficient statistics for every testing point z and bandwidth h.
[In this paper we consider constant bandwidths h for every testing point. Nevertheless, our algorithm framework can be 
easily modified to handle heterogeneous bandwidths as well.]
To simplify our presentation, we abstract the sufficient statistics as

    𝒯_ℓ^z,h = ∑_ℬ(z,h) T_ℓ(X_i,Y_i),       ℓ=1,2,⋯,M,

where M is the total number of sufficient statistics required.
A simple brute-force algorithm loops over all {(X_i,Y_i)}_i=1^n and takes O(n) time to calculate 𝒯_ℓ^z,h for each testing point z,
which results in an overall O(sn) quadratic time complexity.
Further optimization of computational procedures of 𝒯_ℓ^z,h is the focus of the rest of this section.



 §.§ The univariate case: a warm-up exercise


Before giving the full description of our algorithm, we first consider the basic univariate case (d=1) and show
a simple algorithm that computes 𝒯_ℓ^z,h for each z∈𝒳 and h>0 in O(log n) time.
Though simple, the univariate algorithm serves as a conceptual starting point for our follow-up development of general high-dimensional algorithms.

Let {(X_i,Y_i)}_i=1^n be the training data points where X_i,Y_i∈ℝ.
As pre-processing steps, (X_i,Y_i) are sorted in ascending order with respect to X_i, such that X_1≤ X_2≤⋯≤ X_n,
and for each ℓ=1,⋯,M the cumulative statistics

    𝒜_ℓ(i) = ∑_j≤ iT_ℓ(X_j,Y_j)       i=1,2,⋯,n

are computed by a sweeping algorithm from i=1 to i=n.
It is clear that the pre-processing can be done in O(nlog n+Mn) time with O(Mn) space.

Afterwards, for each testing point z∈ℝ and h>0, the left and right “boundaries”
L = max{i: X_i<z-h/2} and U=max{i: X_i≤ z+h/2} are found by binary searches.
The sufficient statistic 𝒯_ℓ^z,h=∑_z-h/2≤ X_i≤ z+h/2T_ℓ(X_i,Y_i) can then be computed as

    𝒯_ℓ^z,h = 𝒜_ℓ(U) - 𝒜_ℓ(L).

The time required to compute m̂_n(z) can then be upper bounded by O(log n + M).
The overall time complexity on s testing points is then O((n+s)(log n+M)),
which is simply O((n+s)log n) because M only depends on d and k, both treated as constants in this paper.



 §.§ Discretization in multiple dimensions



The sorting and partial sum idea in Eq. (<ref>) can be formally generalized to multiple dimensions d>1
via the argument of discretization. (See the leftmost plot of Figure <ref> for a graphical illustration.)
Specifically, for each dimension j∈{1,⋯,d}, 
the values X_1j, X_2j,⋯,X_nj∈ℝ are sorted in ascending order.
We also denote λ_j(i) as the ith smallest value in {X_1j,⋯, X_nj} (repetitions counted as multiple values).
By performing such “discretization” on all dimensions, 
every training point X_i∈𝒳 can be mapped to a unique k-tuple χ_i = (χ_i1,⋯,χ_id) ∈ [n]^d,
where [n]={1,2,⋯,n}.




Similar to the definition of 𝒜_ℓ(i) in the d=1 setting, for the general multivariate case define

    𝒜_ℓ(i_1,⋯,i_d) = ∑_t=1^n𝕀[χ_t1≤ i_1,⋯,χ_td≤ i_d] T_ℓ(X_t,Y_t)

for i_1,⋯,i_d∈[n].
Suppose for now that 𝒜_ℓ(i_1,⋯,i_d) are available for all ℓ and i_1,⋯,i_d by some pre-processing procedure.
To compute the estimate m̂_n(z), first find L_1,⋯,L_d and U_1,⋯,U_d by binary searches, where (for j=1,⋯,d)

    L_j    = max{i: λ_j(i) < z_j-h/2};
    
    	U_j    = max{i:λ_j(i)≤ z_j+h/2}.

The sufficient statistics 𝒯_ℓ^z,h can then be computed as 

    𝒯_ℓ^z,h = ∑_ν_1,⋯,ν_d∈{0,1}^d (-1)^ν_1+⋯+ν_d𝒜_ℓ(i_1^(ν_1),⋯,i_d^(ν_d)),

where i_j^(0)=U_j and i_j^(1)=L_j.
The validity of the computational formula for 𝒯_ℓ^z,h is implied by the inclusion-exclusion principle,
and involves 2^d evaluations of 𝒜_ℓ(i_1,⋯,i_d).
An illustrative example of the d=2 case is given in the middle plot of Figure <ref>.

The question remains to design fast algorithms that compute 𝒜_ℓ(i_1,⋯,i_d)
for arbitrary input tuple (i_1,⋯,i_d)∈ [n]^d.
In the basic univariate case of d=1, this is accomplished by a sweeping pre-processing that records 𝒜_ℓ(i) for all i∈[n].
Such an approach, however, no longer works for d>1 because a sweeping algorithm that records all 𝒜_ℓ(i_1,⋯,i_d)
takes O(n^d) space and time complexity and is clearly impractical.
To overcome this difficulty, we consider a data structure named binary indexed trees and show how it can 
compute every 𝒜_ℓ(i_1,⋯,i_d) evaluation in O(log^d n) time.






 §.§ Binary indexed trees



The binary indexed tree or Fenwick's tree, invented and first documented by <cit.>,
is a simple yet powerful data structure that supports fast queries and updates of partial sums,
corresponding to 𝒜_ℓ(i_1,⋯,i_d) defined in this paper.

In the rightmost panel of Figure <ref> we give a toy example with n=8 data points and d=1.
The V_ℓ(i) entries for i=1,⋯,n are the main data structure kept by the binary indexed tree,
each corresponding to the cumulative statistics of data points in a specific interval.
To query or interrogate a particular partial sum 𝒜_ℓ(i), one starts with i and repeatedly rips off the least significant bit (LSB) in the binary expansion of i
and accumulates the entries of V_ℓ on the path.
For example, for i=5 the query/interrogation path would be 5→ 4→ 0.
Formally, the following update rule is applied when interrogating partial sum 𝒜_ℓ(i):

    Interrogation:    i  i - LSB(i).

Here LSB(i) denotes the least significant bit of i in its binary expansion.

To add (update) a data point 
[Here T_ℓ(i) denotes T_ℓ(X,Y) for the data tuple (X,Y) that corresponds to the ith position in the sorted data locations.
	T_ℓ(i_1,⋯,i_d) has a similar meaning in higher dimensions.]
T_ℓ(i) to the binary indexed tree,
one has to undo the interrogation path in Eq. (<ref>).
It is easy to verify that the following update procedure is an exact mirroring of the interrogation rule in Eq. (<ref>):

    Update:    i  i + LSB(i+1).

For example, starting with i=5 the updating procedure would take us to 5→ 6→ 8,
meaning that V_ℓ(5), V_ℓ(6) and V_ℓ(8) are updated with the data entry T_ℓ(5) added.

To simplify notations, we write ^I(i) and ^U(i) as the interrogation and update paths starting from i, respectively.
For example, ^I(5) = {5,4} and ^U(5) = {5,6,8} if n=8.
It is a simple observation that for any i∈[n], |^I(i)|, |^U(i)| = O(log n).
Thus, both interrogation and updating operations (for d=1) can be completed in O(log n) time.

To extend the univariate binary indexed tree to higher dimensions, consider tuple i_1,⋯,i_d∈[n]^d
at which interrogation or update is needed.
For evaluation of 𝒜_ℓ(i_1,⋯,i_d), 
accumulate all V_ℓ(i_1',⋯,i_d') such that i_1'∈^I(i_1), i_2'∈^I(i_2),⋯, i_d'∈^I(i_d)
and add them together.
Similarly, when adding a data point T_ℓ(i_1,⋯,i_d), consider all i_1'∈^U(i_1),i_2'∈^U(i_2),⋯,i_d'∈^U(i_d)
and update each V_ℓ(i_1',⋯,i_d') with T_ℓ(i_1,⋯,i_d).
The time complexity for both interrogation and updating is O(log^d n).



 §.§ Lazy memory allocation via Hashing


Although the time complexity of d-dimensional binary indexed trees is O(log^d n) and is small,
the space complexity is another story.
A naive implementation of the multi-dimensional binary indexed tree would require O(n^d) memory to 
store all V_ℓ(i_1,⋯,i_d) values.
It is imperative to design more intelligent space allocation algorithms that significantly reduce the O(n^d) space complexity
for practical usages.

An important observation of the multi-dimensional binary indexed tree is that the V_ℓ values are extremely sparse,
with most of the entries equal to zero.
More specifically, with n points in the training data set, the number of non-zero entries of V_ℓ is at most O(nlog^d n)
because each data point is associated with a update path of length at most O(log^d n),
which is significantly smaller than O(n^d).
This motivates us to consider a lazy memory allocation scheme based on Hash tables,
which would subsequently reduce the space complexity from O(n^d) to nearly O(nlog^d n).

In particular, we construct M hash tables {ℋ_ℓ}_ℓ=1^M each corresponding to a sufficient statistics 𝒯_ℓ to be maintained.
When adding a new data point T_ℓ(i_1,⋯,i_d), we insert all entries i_1'∈^U(i_1),⋯,i_d'∈^U(i_d) into the corresponding Hash table ℋ_ℓ;
similarly when interrogating a point (i_1,⋯,i_d) we read values off all relevant entries i_1'∈^I(i_1),⋯,i_d'∈^I(i_d) from the Hash table ℋ_ℓ and add them together.

Let b∈ℕ be a pre-determined capacity parameter of hash tables {ℋ_ℓ}_ℓ=1^M to be constructed, which satisfies b=Ω(nlog^d n). 
Let τ∈ℕ be another pre-defined integer.
Suppose ħ_1,⋯,ħ_d:[n]→[b] are τ d-wise independent hash functions, meaning that for all j=1,⋯,d and any i_1,i_2,⋯,i_τ d∈[n], a_1,⋯,a_τ d∈[b], 

    [ħ_j(i_1)=a_1∧⋯∧ħ_j(i_τ d)=a_τ d] = b^-τ d.

A “composite” hash function H:[n]^d→[b] can then be defined as

    H(i_1,⋯,i_d) := (ħ_1(i_1) + ⋯ + ħ_d(i_d))  b.

It is a well-known result that H is a τ-wise independent hash function on [n]^d (see for example <cit.>).
The corresponding entry of (i_1,⋯,i_d)∈[n]^d in the Hash table ℋ_ℓ can then be located at

    ℋ_ℓ(H(i_1,⋯,i_d)).

When collision occurs, standard methods such as hash chains or linear probing can be employed to resolve collision.
In particular, when τ=5 for linear probing, the expected number of collision would be at a constant level <cit.>.




 §.§ Summary of time and space complexity


The capacity of the Hash table constructed in the previous section, b, can be set as b≍ nlog^d n,
leading to the expected number of collision to be upper bounded at a constant level.
Hence, the space complexity of the proposed algorithm is O(nlog^d n).

For the time complexity, note that each insertion of a data point x_i,i∈[n] requires M update steps along the update paths in Eq. (<ref>),
and each estimate requires 2^d M interrogation steps along the interrogation paths in Eq. (<ref>).
Because M≤ (2k)^d+k^d and the length of each update/interrogation path is upper bounded by O(log^d n),
the total time complexity of Algorithm <ref> is upper bounded by O((2k)^d (n+s)log^d n)=O((n+s)log^d n), where both k and d
are treated as constants.









§ NUMERICAL RESULTS






























































































































In this section we use simulations to verify the effectiveness of our proposed algorithm.
All source codes of our implementation are provided in the supplementary materials accompanying this paper.
The raw data and results in the plotted figures and tables are also provided in a separate online supplement.



We implement algorithm <ref> using C++.

All simulations are done on a computer installed with uBuntu 14.04 LTS 64bit OS, 31.4 Gib Memory, Intel Core i7-4970 CPU @3.60GHz x 8,
except for the running time reported for our own implementation of the brute-force algorithm
which is experimented on a MacBook Pro laptop with macOS Mojave 10.14.6, 16GB memory, 2.3 GHz Intel Core i9
because the former server has been de-serviced at the time this part of the experiment is carried out.



Two baselines are considered and compared against our proposed method.
The first baseline is a widely used R library [https://cran.r-project.org/web/packages/locpol], which is implemented in C.
A detailed description of the  package can be found in <cit.>.
The relevant routines in  include , ,  and ,
corresponding to local polynomial estimation of k=0, k=1, k=2 and k>2, respectively.
Note that for d≥ 2 only local constant regression (i.e., k=0) is supported in .
While the detailed implementation strategy of  is not available, judging from the high computational costs in numerical experiments
we conjecture that the naive brute-force implementation is used which computes θ̂_n separately for each data point.

The second baseline is our own implementation of the brute-force computation of local polynomial estimates,
which serves as a second baseline for comparison.
We implement the brute-force approach using C++, which supports general data dimension d and polynomial degree k.

















We first compare the running times of our fast local polynomial regression implementation (Algorithm <ref>, denoted as )
with the baseline R package  under various simulation settings in Figure <ref> (d=1,k=0),
Figure <ref> (d=1,k=1) and Figure <ref> (d=2,k=0).

The ground-truth functions f:ℝ^d→ℝ are defined as f(x_1,⋯,x_d) = ∑_i=1^d sin(x_i).
Note that for d=2 we only experiment with local mean averaging (k=0) because the  package only supports k=0 for d>1.
For all experimental setups, the running times are reported under varying number of training and testing points,
all uniformly sampled from the unit cube [0,1]^d.




From Figures <ref>, <ref> and <ref>, we observe that
as the number of testing points (s) or the training points become larger, the advantages of our proposed algorithm are more significant.
In particular, when 1,024,000 training and testing points are present under the setting of k=1 and d=1, Algorithm <ref> achieves 40,000× speed-up
over the naive implementation in .
The reason is that the naive implementation has Θ(ns) time complexity whereas ours is O((n+s)log^d n).


We compare our algorithm with a baseline implementation of local polynomial regression in Figures <ref>, <ref>, <ref>
on 3-dimensional data (i.e., d=3) with polynomial degrees ranging from k=0 to k=2.
Because the  package in R only supports k=0 for d>1, we compare our proposed algorithm with our own brute-force implementation
of local polynomial regression in C++.
As we can observe from the figures, our proposed algorithm still consistently outperforms the baseline implementation for all polynomial degrees
in the case of d=3.

We also compare memory consumptions between our method and the naive method in Table <ref>.
For all experiments we fix d=3 and s = 4000 (we found different s do not affect the memory consumption of both methods by much).
As we can observe from Table <ref>, our implementation consumes significantly more memory compared to baseline
because the multi-dimensional binary indexed trees is quite a memory-intensive data structure,
even with lazy memory allocation schemes in place.
Nevertheless, the highest amount of memory consumption (approximately 7GB for d=3, k=2 and n=64,000)
is still well within the limit of modern main memories on desktop or workspace computers.
Such memory consumption is worthy because of the significant improvement in running time our proposed algorithm brings.







Lastly, we report in Figure <ref> the function reconstruction errors of both FastLPR and the baseline method in  under fixed running time budgets.
We fix number of testing points to be 4000 and the bandwidth parameter h takes values of n^-1/3, n^-1/4 and n^-1/5.
Figure <ref> shows under given time budgets, our proposed FastLPR algorithm processes much more observations (training samples)
and therefore achieves significantly smaller reconstruction error.



§ DISCUSSION




  
Kernel functions.
In our accelerated local polynomial estimation algorithm, the kernel K(·) is assumed to be the box kernel K(u)=𝕀[|u|≤ 1/2]
and for multiple dimensions (d≥ 2) the “composite kernel” is defined in terms of ℓ_∞ norm, meaning that the kernel evaluation between x,z∈ℝ^d
is K(x-z_∞/h).
The box kernel ensures that the sufficient statistics can be expressed as a linear combination of unweighted partial sums,
and the ℓ_∞ norm in multi-variate kernel evaluations lead to rectangle regions whose edges are parallel to the standard basis.
Both properties are crucial for our algorithmic development that enables fast sufficient statistics computation via binary indexed trees.

In general, box kernels with ℓ_∞ distance measures are sufficient as they achieve the same minimax statistical efficiency
as other kernel functions and/or equivalent distance measures.
However, for certain applications it might be desirable to consider kernel functions smoother than the box kernel (e.g., the Gaussian kernel K(u)=e^-u^2/2/√(2π))
to obtain smoother function fits.
We believe significantly different techniques are required to handle such kernels that are non-uniform and not restricted to parallel rectangular neighborhoods of testing points.




  
Incremental training and testing sets.
While in the description of our algorithm the sizes of both training and testing sets (n and s) are fixed a priori,
we remark that our algorithm can perfectly handle training and testing data with increasing sizes.
In particular, for every new training or testing point, the corresponding update or interrogation paths are calculated in O(log^d n) time
and the Hash tables can be updated/queried using a similar amount of running time.
This property is particularly useful in applications where training data constantly grow/change, 
and estimates on incoming test points have to be made in a real-time fashion.




  
Further acceleration with approximate computation. 
In cases where exact computations of local polynomial estimates are not mandatory and small error in the estimates can be tolerated,
it is possible to further reduce the time and space complexity beyond O(nlog^d n).

One approach is to consider significantly smaller (shorter) Hash tables with capacity b≪ nlog^d n.
Because the capacity of Hash tables is significantly smaller than the number of entries created by binary indexed trees, collisions are unavoidable.
Instead of resolving collisions completely, one can borrow the ideas of CountSketch <cit.> that multiplies an additional Rademacher 
variable
[A Rademacher random variable takes on values of ± 1 with equal probability.]
σ(i_1,⋯,i_d) in Hash table updates:

    ℋ_ℓ(H(i_1,⋯,i_d))    ℋ_ℓ(H(i_1,⋯,i_d)) 
            + σ(i_1,⋯,i_d)T_ℓ(i_1,⋯,i_d).





  
Dependency on data dimension d.
The time complexity of our proposed algorithm has a leading O((2k)^d) term which scales exponentially with data dimension d,
which is dropped in the main O((n+s)log^d n) time complexity bound because both polynomial degree k and data dimension d are treated
as constants in this paper.
On the other hand, the naive approach has an O(n^2) time complexity and does not explicitly depend exponentially on d.
Nevertheless, in nonparametric regression/estimation it is typical that an exponential number of sample points are required 
(e.g., to estimate a Lipschitz continuous function up to precision ε, Ω(ε^-(2+d)) samples are needed <cit.>).
Hence, our O((n+s)log^d n) time complexity is generally much more preferable than the O(n^2) time complexity of naive approaches.





§ CONCLUSION


In this paper we propose nearly linear time algorithms that compute local polynomial estimates for nonparametric density and regression function estimates.
Our algorithm is based on the novel application of multivariate binary indexed trees together with discretization and hashing techniques.
Simulation results demonstrate an up to 40000× speed-up over state-of-the-art R implementation of local polynomial regression.
Future directions including general kernel functions and further algorithmic acceleration via sketching approximation are discussed.

informs2014










