Suppose we have training data $(\mathbf{x}_1, y_1), \dots, (\mathbf{x}_n, y_n)$ with $\mathbf{x}_i \in A \subset \mathbb{R}^p$ and $y_i \in \mathbb{R}$. We call $A$ the predictor space, and partition it into regions $A_1, \dots, A_J$.
Consider a predictor $\hat{f}: \mathbb{R}^p \to \mathbb{R}$ given by $$\hat{f}(\mathbf{x}) = \sum_{j=1}^{J}c_j \cdot \mathbf{I}(\mathbf{x} \in A_j)\text{,}$$ where $\mathbf{I}$ denotes the indicator function.
Our aim is to minimize the loss function $$\sum_{i=1}^{n}[y_i - \hat{f}(\mathbf{x}_i)]^2\text{.}$$
Without proof, both Elements of Statistical Learning by Hastie, Tibshirani, and Freedman (2nd ed.), as well as Modern Multivariate Statistical Techniques by Izenman state that the optimal $c_j$ by the loss function above is $$\hat{c}_j = \dfrac{1}{n_j}\sum_{\{j: y_j \in A_j\}}y_j$$ i.e., take the average of all $y_i$ values in region $A_j$ (where $n_j$ is the number of $y_i$ in region $A_j$).
How is this proven? It is a standard result that for a general estimator $\hat{f}(\mathbf{X})$ of a single value $Y$, the optimal estimator under squared-error loss is $\hat{f}(\mathbf{X}) = \mathbb{E}[Y \mid \mathbf{X}]$. However, I'm not clever enough to see how this fact could be applicable.
Edit: Is this just as simple as differentiating with respect to $c_j$ and setting that equal to $0$? Funny how an obvious approach seems to make sense after asking the question...
Edit 2: Unfortunately, that approach doesn't lead to a very clean solution... I obtain \begin{equation*} \hat{c}_j = \dfrac{\sum_{i=1}^{n}y_i \cdot \mathbf{I}(\mathbf{x}_i \in A_j) - \sum_{k \neq j}\hat{c}_k n_k}{n_j}\text{.} \end{equation*} after setting the partial derivative with respect to $c_j$ equal to $0$ of the loss function.
Edit 3: I think I figured it out. I screwed up in mixing up the indices. The answer currently up is great, and I'm accepting it.
