In this work, we consider the problem of estimating a nonparametric function under a variety of smoothness and shape constraints such as monotonicity, convexity, unimodality and Lipschitz smoothness whenever some prior knowledge about the relationship between the independent and dependent variables is given. By imposing shape restrictions only on a finite set of points by discretizing the domain of the unknown function, we obtain finite-dimensional problems which are sampled versions of the infinite-dimensional function estimation problem. We solve the corresponding shape constrained finite-dimensional problems on a discretized grid using first-order methods with specialized projection operators, and then iteratively refine the solution to solve the infinite-dimensional problem to any degree of optimality. We extend the first-order algorithms developed for shape constrained univariate function estimation to solve multivariate shape constrained problems using additive modeling and provide certificates of optimality for a solution on the finite discretization with respect to the infinite-dimensional problem. Finally, using specialized projection operators for convexity constraints, we propose a fast algorithm to solve the univariate convex regression problem.