Is the MLE estimate of the mean of Gamma distribution same as the average? Consider the Gamma distribution:
$\Gamma(\alpha,\beta) = \frac{\beta^\alpha x^{\alpha-1} e^{-\beta x}}{\Gamma(\alpha)}$
According to many sources (e.g. wikipedia), the mean is given by $\mu=\alpha/\beta$. On the other hand, using Maximum Likelihood Estimation, we can find that (see page 11 here for example) $\hat{\beta}_{MLE}=\frac{\hat{\alpha}_{MLE}}{\bar{x}}$, where $\bar{x}$ is the average. Accordingly, $\hat{\mu}_{MLE} = \frac{\hat{\alpha}_{MLE}}{\hat{\beta}_{MLE}}=\bar{x}$.
Is this correct? If so, why do many sources go on length to do numerical estimation of the parameters? Since most of the time all what we care about is mean and variance, I don't see why we would resort to numerical estimation of individual parameters.

Indeed it's isn't necessary to solve for the shape parameter, $α$ in order to estimate $\mu$. If you reparameterize the gamma to the shape-mean parameterization (as you essentially have in a gamma GLM), then the MLE of the mean can be done without having an estimate of α and indeed will be $\bar{x}$.
You can find the shape-mean parameterization here and if you do the calculations - which are quite straightforward - you'll see the shape parameter cancel out of the MLE for $\mu$. So you can certainly estimate $\mu$ by MLE without solving for the shape parameter.
However, for the variance ($\frac{\mu^2}{\alpha}$ in this parameterization), I don't think you can avoid estimating the shape parameter (or something essentially as difficult, if you replace $\alpha$ in yet another reparameterization).
For example, I don't think it helps whether you replace $\alpha$ by the coefficient of variation or even the variance itself, you'll either end up with a function of the shape parameter or you end up with a problem of similar (or slightly greater) difficulty to estimating the shape. Might as well just be done and estimate the shape itself in that case.
