What is this bias-variance tradeoff for regression coefficients and how to derive it? In this paper, (Bayesian Inference for Variance Components Using Only Error Contrasts, Harville, 1974), the author claims
$$(y-X\beta)'H^{-1}(y-X\beta)=(y-X\hat\beta)'H^{-1}(y-X\hat\beta)+(\beta-\hat\beta)'(X'H^{-1}X)(\beta-\hat\beta)$$
to be a "well-known relationship", for a linear regression
$$y=X\beta+\epsilon,$$
where
$$\epsilon\sim\mathcal{N}(0, H).$$
How is this well-known? What is the simplest way of proving this?
 A: They arrive at this identity by a technique called completing the square.  The left hand side is in a quadratic form, so start by multiplying it out
$$ (y-X\beta)'H^{-1}(y-X\beta)= y'H^{-1}y - 2y'H^{-1}X\beta + \beta'X'H^{-1} X\beta $$
continue on and then rewrite in terms of $\hat{\beta} = (X'H^{-1}X)^{-1}X'H^{-1}y$. The algebra is kind of long but googling completing the square in Bayesian regression and you can find plenty of hints.  For example, see the wikipedia on Bayesian linear regression, and other CrossValided answers regarding completing the square, like here.  
A: The last term in the equation can be written as
$$
(X\beta - X\hat{\beta})'H^{-1}(X\beta - X\hat{\beta}).
$$
In this form the equation is saying something interesting.  Assuming $H$ is positive definite and symmetric, so is its inverse.  Therefore,we can define an inner product $<x, y>_{H^{-1}} = x'H^{-1}y$, giving us geometry.  Then the above equality is essentially saying that,
$$
(X\beta - X\hat{\beta}) \perp (y - X\hat{\beta}).
$$
I wanted to give you this bit of intuition since a commenter has already left a link to the derivation.
Edit: For Posterity
LHS:
\begin{eqnarray}
(y-X \beta)'H^{-1}(y-X \beta) &=& y'H^{-1}y &-& 2y'H^{-1}X \beta &+& \beta'X'H^{-1}X\beta \\
&=& (A) &-& (B) &+& (C)
\end{eqnarray}
RHS:
$$
(y-X\hat\beta)'H^{-1}(y-X\hat\beta)+(\beta-\hat\beta)'(X'H^{-1}X)(\beta-\hat\beta)
$$
\begin{eqnarray}
&=& y'H^{-1}y &- 2y'H^{-1}X\hat{\beta} &+ \hat{\beta}'X'H^{-1}X\hat{\beta} &+ \beta X'H^{-1}X\beta &- 2\hat{\beta}X'H^{-1}X\beta &+ \hat{\beta}'X'H^{-1}X\hat{\beta} \\
&=& (A) &- (D) &+ (E) &+ (C) &- (F) &+ (E)
\end{eqnarray}
Relation:
$$
\hat{\beta} = (X'H^{-1}X)^{-1}X'H^{-1}y
$$
By plugging in the relation you can show that (B) = (F), and that 2(E) = (D).  All done.
A: If you know your matrix algebra, then this should be doable by multiplying everything out and verifying that you indeed have the same on both sides. This is what jlimahaverford has demonstrated.
To be able to do this you need the formula for the estimate of $\hat{\beta}$. We can derive the formula in a similar manner as for linear regression when we have uncorrelated error terms. The trick is to standardize.
Here is some information on how to standardize a RV that comes from a multivariate normal distribution. Let's assume that you have 
$$
\mathbf{X}\sim \mathcal{N}(\mu,\Sigma).
$$
$\Sigma$ is positive definite, so you can factorize it as $\Sigma = PP^T$. Now the random variable 
$$
\mathbf{Y}=P^{-1}(\mathbf{X}-\mu)
$$
comes from the distribution $\mathcal{N}(0,I)$. Now we can use this trick for our problem to find $\hat{\beta}$. Let's factorize $H=PP^T$. We have
$$
\begin{align}
y&=X\beta+\epsilon\\
P^{-1}y &= P^{-1}X\beta + P^{-1}\epsilon
\end{align} 
$$
Now $\epsilon$ has been standardized, such that $\text{cov}(P^{-1}\epsilon)=I$, so we can now treat this as a simple multiple linear regression model where:
$$
\tilde{X}=P^{-1}X,\qquad \tilde{y}=P^{-1}y\quad\text{and}\quad \tilde{\epsilon}=P^{-1}\epsilon.
$$
So we have the regression problem:
$$
\tilde{y}=\tilde{X}\beta+\tilde{\epsilon}
$$
The formula for $\hat{\beta}$ is
$$
\begin{align}
\hat{\beta} &= (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^T\tilde{y}\\
&=((P^{-1}X)^TP^{-1}X)^{-1}(P^{-1}X)^TP^{-1}y\\
&=(X^T(PP^T)^{-1}X)^{-1}X(PP^T)^{-1}y\\
&=(X^TH^{-1}X)^{-1}XH^{-1}y
\end{align}
$$
This is the key to do this, the rest is the algebraic manipulation demonstrated in the solution by jlimahaverford.
