Convex regression (CR) problem deals with fitting a convex function to a finite number of observations. It has many applications in various disciplines, such as statistics, economics, operations research, and electrical engineering. Computing the least squares (LS) estimator via solving a quadratic program (QP) is the most common technique to fit a piecewise-linear convex function to the observed data. Since the number of constraints in the QP formulation increases quadratically in N, the number of observed data points, computing the LS estimator is not practical using interior point methods when N is very large. The first-order method proposed in this paper carefully manages the memory usage through parallelization, and efficiently solves large-scale instances of CR.