Local polynomial regression (Fan & Gijbels, 1996) 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 (Fenwick, 1994). Both time and space complexities of our proposed algorithm are nearly linear in the number of inputs. Simulation results confirm the efficiency and effectiveness of our approach.