The math behind minimizing the loss function for regression trees 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.

Intuitively, if you were to choose a constant as your prediction, the best choice would be $E[y|x]$ as you quoted. If you can divide the predictor space into $J$ regions (pre-defined) and have the freedom to choose $J$ constants, then region specific conditional averages would your best choice. This is like taking the mean of the samples in the leaves of a decision tree.
$$E[Y|\mathbf X]=\sum_{j=1}^J \mathbb P(\mathbf X \in A_j)\mathbb E[Y|\mathbf X, \mathbf X\in A_j]=c_m$$
because only one of the terms in this sum is non-zero and that is the region of $\mathbf X$, say $m$.
