We consider the problem of estimating rectangular piecewise constant functions in multiple dimensions using the method of least squares. We introduce and study two least squares estimators over: (i) the class of all completely monotone functions (in [0,1]^d), and (ii) the class of functions (in [0,1]^d) with finite Hardy and Krause variation (anchored at 0). The above two classes can be thought of as extensions of the class of all monotone functions and the class of all functions with finite bounded variation on the real line, respectively.
We show that the two estimators can be computed by solving a nonnegative least squares problem and a LASSO problem respectively, and furthermore the design matrices in these two finite-dimensional problems are the same. Our main results show that the worst case risk of these estimators do not depend on the dimension d (up to logarithmic factors). Further, these estimators exhibit a parametric (k/n) rate of convergence (up to logarithmic factors) for estimating piecewise constant functions in multiple dimensions, where k is the number of constant "pieces" and n is the sample size.