  We study the nonparametric least squares estimator (LSE) of a multivariate
convex regression function. The LSE, given as the solution to a quadratic
program with $O(n^2)$ linear constraints ($n$ being the sample size), is
difficult to compute for large problems. Exploiting problem specific structure,
we propose a scalable algorithmic framework based on the augmented Lagrangian
method to compute the LSE. We develop a novel approach to obtain smooth convex
approximations to the fitted (piecewise affine) convex LSE and provide formal
bounds on the quality of approximation. When the number of samples is not too
large compared to the dimension of the predictor, we propose a regularization
scheme --- Lipschitz convex regression --- where we constrain the norm of the
subgradients, and study the rates of convergence of the obtained LSE. Our
algorithmic framework is simple and flexible and can be easily adapted to
handle variants: estimation of a non-decreasing/non-increasing convex/concave
(with or without a Lipschitz bound) function. We perform numerical studies
illustrating the scalability of the proposed algorithm.
