What is the expected value of the logarithm of Gamma distribution? If the expected value of $\mathsf{Gamma}(\alpha, \beta)$ is $\frac{\alpha}{\beta}$, what is the expected value of $\log(\mathsf{Gamma}(\alpha, \beta))$? Can it be calculated analytically? 
The parametrisation I am using is shape-rate.
 A: This one (maybe surprisingly) can be done with easy elementary operations (employing Richard Feynman's favorite trick of differentiating under the integral sign with respect to a parameter).

We are supposing $X$ has a $\Gamma(\alpha,\beta)$ distribution and we wish to find the expectation of $Y=\log(X).$  First, because $\beta$ is a scale parameter, its effect will be to shift the logarithm by $\log\beta.$  (If you use $\beta$ as a rate parameter, as in the question, it will shift the logarithm by $-\log\beta.$) This permits us to work with the case $\beta=1.$
After this simplification, the probability element of $X$ is
$$f_X(x) = \frac{1}{\Gamma(\alpha)} x^\alpha e^{-x} \frac{\mathrm{d}x}{x}$$
where $\Gamma(\alpha)$ is the normalizing constant
$$\Gamma(\alpha) = \int_0^\infty x^\alpha e^{-x} \frac{\mathrm{d}x}{x}.$$
Substituting $x=e^y,$ which entails $\mathrm{d}x/x = \mathrm{d}y,$ gives the probability element of $Y$,
$$f_Y(y) = \frac{1}{\Gamma(\alpha)} e^{\alpha y - e^y} \mathrm{d}y.$$
The possible values of $Y$ now range over all the real numbers $\mathbb{R}.$
Because $f_Y$ must integrate to unity, we obtain (trivially)
$$\Gamma(\alpha) = \int_\mathbb{R} e^{\alpha y - e^y} \mathrm{d}y.\tag{1}$$
Notice $f_Y(y)$ is a differentiable function of $\alpha.$ An easy calculation gives
$$\frac{\mathrm{d}}{\mathrm{d}\alpha}e^{\alpha y - e^y} \mathrm{d}y =  y\, e^{\alpha y - e^y} \mathrm{d}y = \Gamma(\alpha) y\,f_Y(y).$$
The next step exploits the relation obtained by dividing both sides of this identity by $\Gamma(\alpha),$ thereby exposing the very object we need to integrate to find the expectation; namely, $y f_Y(y):$
$$\eqalign{
\mathbb{E}(Y) &= \int_\mathbb{R} y\, f_Y(y) = \frac{1}{\Gamma(\alpha)} \int_\mathbb{R} \frac{\mathrm{d}}{\mathrm{d}\alpha}e^{\alpha y - e^y} \mathrm{d}y \\
&= \frac{1}{\Gamma(\alpha)} \frac{\mathrm{d}}{\mathrm{d}\alpha}\int_\mathbb{R} e^{\alpha y - e^y} \mathrm{d}y\\
&= \frac{1}{\Gamma(\alpha)} \frac{\mathrm{d}}{\mathrm{d}\alpha}\Gamma(\alpha)\\
&= \frac{\mathrm{d}}{\mathrm{d}\alpha}\log\Gamma(\alpha)\\
&=\psi(\alpha),
}$$
the logarithmic derivative of the gamma function (aka "polygamma").  The integral was computed using identity $(1).$
Re-introducing the factor $\beta$ shows the general result is
$$\mathbb{E}(\log(X)) = \log\beta + \psi(\alpha)$$
for a scale parameterization (where the density function depends on $x/\beta$) or
$$\mathbb{E}(\log(X)) = -\log\beta + \psi(\alpha)$$
for a rate parameterization (where the density function depends on $x\beta$).
A: The answer by @whuber is quite nice; I will essentially restate his answer in a more general form which connects (in my opinion) better with statistical theory, and which makes clear the power of the overall technique. 
Consider a family of distributions $\{F_\theta : \theta \in \Theta\}$ which consitute an exponential family, meaning they admit a density 
$$
f_\theta(x) = \exp\left\{s(x)\theta - A(\theta) + h(x)\right\}
$$
with respect to some common dominating measure (usually, Lebesgue or counting measure). Differentiating both sides of
$$
\int f_\theta(x) \ dx = 1
$$
with respect to $\theta$ we arrive at the score equation
$$
\int f'_\theta(x) = 
\int \frac{f'_\theta(x)}{f_\theta(x)} f_\theta(x)
=
\int u_\theta(x) \, f_\theta(x) \ dx = 0 \tag{$\dagger$}
$$
where $u_\theta(x) = \frac d {d\theta} \log f_\theta(x)$ is the score function and we have defined $f'_\theta(x) = \frac{d}{d\theta} f_\theta(x)$. In the case of an exponential family, we have 
$$
u_\theta(x) = s(x) - A'(\theta) 
$$
where $A'(\theta) = \frac d {d\theta} A(\theta)$; this is sometimes called the cumulant function, as it is evidently very closely related to the cumulant-generating function. It follows now from $(\dagger)$ that $E_\theta[s(X)] = A'(\theta)$. 
We now show this helps us compute the require expectation. We can write the gamma density with fixed $\beta$ as an exponential family
$$
f_\theta(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}
= 
\exp\left\{\log(x) \alpha + \alpha \log \beta - \log \Gamma(\alpha) - \beta x \right\}.
$$
This is an exponential family in $\alpha$ alone with $s(x) = \log x$ and $A(\alpha) = \log \Gamma(\alpha) - \alpha \log \beta$. It now follows immediately by computing $\frac d {d\alpha} A(\alpha)$ that 
$$
E[\log X] = \psi(\alpha) - \log \beta. 
$$
