Bias-variance decomposition and independence of $X$ and $\epsilon$ I took a look at a couple of derivations of this decomposition and I think they all require that $E[\epsilon \hat{f}] = 0$.
The most transparent one I found was from Wikipedia here.
I reproduce the steps below for your convenience.
The data generating process is $y = f(x) + \epsilon$. As a shorthand, let $f = f(x)$ and let $\hat{f} = \hat{f}(x)$ be the fitted model.
$
\begin{align}
\mathrm{E}\big[(y - \hat{f})^2\big]
 & = \mathrm{E}[y^2 + \hat{f}^2 - 2 y\hat{f}] \\
 & = \mathrm{E}[y^2] + \mathrm{E}[\hat{f}^2] - \mathrm{E}[2y\hat{f}] \\
 & = \mathrm{Var}[y] + \mathrm{E}[y]^2 + \mathrm{Var}[\hat{f}] + \mathrm{E}[\hat{f}]^2 - 2f\mathrm{E}[\hat{f}] \\
 & = \mathrm{Var}[y] + \mathrm{Var}[\hat{f}] + (f - \mathrm{E}[\hat{f}])^2 \\
 & = \mathrm{Var}[y] + \mathrm{Var}[\hat{f}] + \mathrm{E}[f - \hat{f}]^2 \\
 & = \sigma^2 + \mathrm{Var}[\hat{f}] + \mathrm{Bias}[\hat{f}]^2
\end{align}$
My question is about the third step, where they use
$ E[2y\hat{f}] = 2fE[\hat{f}].$
I think what they are assuming is that $E[\epsilon \hat{f}] = 0$.
Is it because we are assuming that $X$ is independent of $\epsilon$? Is that necessary to do this derivation? It seems such a big assumption that I feel I must be missing something. Is there another decomposition that doesn't require this?
 A: I find the most difficult parts of this subject to be in keeping track of what is is being integrated over (the expectation with respect to what?).  Many expositions take this lightly, and the end result can be highly confusing.  Here is The Elements of Statistical Learning on this point

Discussions of error rate estimation can be confusing, because we have to make clear which quantities are fixed and which are random

I think you're seeing a bit of that in the wikipedia exposition.
My Preferred Derivation
Here's a detailed break down.  Working out the exposition in this way is what finally cleared up the various concepts at play for me.  I hope it helps you too.
The first step is to condition on $x$, and take the expectation with respect to $y$.  There are no assumptions for this part, $g$ is a general function
$$\begin{align*}
    ESE(f; x) &= E_Y \left[ \left( y - g(x) \right)^2 \mid x \right] \\
      & = E_Y \left[ \left( y - E[y \mid x] + E[y \mid x]
          - g(x) \right)^2 \mid x \right] \\
      & = E_Y \left[ \left( y - E[y \mid x] \right)^2 \mid x \right] 
          + E_Y \left[ \left(E[y \mid x] - g(x) \right)^2 \mid x \right] \\      
      & \quad + 2 E_Y \left[ 
          \left( y - E[y \mid x] \right)
          \left( E[y \mid x] - g(x) \right)
        \mid x \right] \\
\end{align*}$$
The trick here is introducing the $E[ y \mid x ]$ into the middle of the sum.  This expectation is over $y$ (the only thing that makes sense in this context), so $E[ y \mid x ]$ is just a function of $x$, same as $g$.
Now we may observe that in
$$ 2 E_Y \left[ 
          \left( y - E[y \mid x] \right)
          \left( E[y \mid x] - g(x) \right)
        \mid x \right] $$
the second factor has no dependence on $y$, and the first factor is zero in expectation, so this term dies, and we are left with
$$\begin{align*}
    ESE(f; x) = E_Y \left[ \left( y - E[y \mid x] \right)^2 \mid x \right] 
                + E_Y \left[ \left(E[y \mid x] - g(x) \right)^2 \mid x \right] \\      
\end{align*}$$
The first term is the irreducible error of the classifier.  It's your $\sigma^2$ above.
Its worth noting that only the integrand of the first term in this result has any dependence on $y$, so the $E_Y$ may be dropped from the second term.
To get to the bias-variance breakdown, we introduce an expectation over the sampling distribution of $x$.  That is, we consider $g$ itself as a random variable; training a learning algorithm takes us from a data set $D$ to a function $g$.  Thus, it makes sense to take expectations over this sampling distribution with $g$ as a random variate.  I'll write $g(x, D)$ to emphasize this dependence of $g$ on the training dataset $D$.  With this notation, the second term in our final equation above becomes
$$ \left[ \left(E[y \mid x] - g(x,D) \right)^2 \mid x \right] $$
Taking the expectation of everything over $D$, and writing $Eg(x)$ for the expectation of $g(x, D)$ with respect to this sampling distribution, we can break this down further
$$\begin{align*}
    & E_{D} \left[ \left(E[y \mid x] - g(x,D) \right)^2 \mid x \right] \\
      & = E_{D} \left[ \left(E[y \mid x] - Eg(x) + Eg(x) - g(x,D) \right)^2 \mid x \right] \\
      & = E_{D} \left[ \left(E[y \mid x] - Eg(x) \right)^2 
            \mid x \right] 
           + E_{D} \left[ \left(Eg(x) - g(x, D) \right)^2 \mid x \right] \\
      & \quad + 2 E_{D} \left[ 
          \left(E[y \mid x] - Eg(x) \right)
         \left( Eg(x) - g(x,D) \right)
        \mid x \right] \\
  \end{align*}
$$
The same trick applies.  In this breakdown the first factor has no dependence on $D$, and the second is zero in expectation, so the cross term dies.  Therefore
$$\begin{align*}
    E_{D} \left[ \left(E[y \mid x] - g(x,D) \right)^2 \mid x \right] = E_{D} \left[ \left(E[y \mid x] - Eg(x) \right)^2 
            \mid x \right] 
           + E_{D} \left[ \left(Ef(x) - g(x, D) \right)^2 \mid x \right]
\end{align*}$$
The first term here is the bias, and the second is the variance.
Of course, I kept the conditioning on $x$ the entire time, so what we have really arrived at is the pointwise decomposition.  To get the usual full decomposition of the total error, all that is needed is to integrate out $x$ as well.
How This Relates to the Wikipedia Breakdown
Wikipedia writes $f(x)$ for my $E[ y \mid x ]$.  I think it makes the exposition easier to follow to make this explicit.  Wikipedia writes
$$ y = f(x) + \epsilon $$
with the assumption that $E(\epsilon) = 0$.  What is really meant is $E(\epsilon \mid x) = 0$.  This is implicit in my exposition, I simply use $E[y \mid x]$ throughout.
Note that no independence assumptions are necessary for the bias variance decomposition.  What wikipedia says

Thus, since $\epsilon$ and $\hat{f}$ are independent

what they really mean is that $\hat{f}$ can be considered a constant when conditioning on $x$, so it can be taken out front of the expectation. 
 
A: Here is a derivation of the bias-variance decomposition, in which I make use of the independence of $X$ and $\epsilon$. 
True model
Suppose that a target variable $Y$ and a feature variable $X$ are related via $Y = f(X) + \epsilon$, where $X$ and $\epsilon$ are independent random variables and the expected value of $\epsilon$ is zero, $E[\epsilon] = 0$. 
We can  use this mathematical relationship to generate a data set $\cal D$. Because data sets are always of finite size, we may think of $\cal D$ as a random variable, the realizations of which take the form
$d = \{ (x_1,y_1), \ldots , (x_m,y_m) \}$, where $x_i$ and $y_i$ are realizations of $X$ and $Y$. 
Estimated model
Machine learning uses a particular realization $d$ of $\cal D$ to train an estimate of the function $f(x)$, called the hypothesis $h_d(x)$. The subscript $d$ reminds us that the hypothesis is a random function that varies over training data sets. 
Test error of the estimated model
Having learned an hypothesis for a particular training set $d$, we next evaluate the error made in predicting the value of $y$ on an unseen test value $x$. In linear regression, that test error is quantified by taking a test data set (also drawn from the distribution of $\cal D$) and computing the average of $(Y - h_d)^2$ over the data set. If the size of the test data set is large enough, this average is approximated by $E_{X,\epsilon} [ (Y(X,\epsilon) - h_{d}(X))^2 ]$. As the training data set $d$ varies, so does the test error; in other words, test error is a random variable, the average of which over all training sets is given by 
\begin{equation*}
\text{expected test error} = E_{\cal D} \left[ E_{X,\epsilon} \left[ (Y(X,\epsilon) - h_{\cal D}(X))^2 \right] \right].
\end{equation*}
In the following sections, I will show how this error arises from three sources: a bias that quantifies how much the average of the hypothesis deviates from $f$; a variance term that quantifies how much the hypothesis varies among training data sets; and an irreducible error that describes the fact that one's ability to predict is always limited by the noise $\epsilon$.  
Establishing a useful order of integration
To compute the expected test error analytically, we rewrite the expectation operators in two steps. The first step is to recognize that 
$ E_{X,\epsilon} [\ldots] = E_X \left[ E_\epsilon [ \ldots ] \right],$
since $X$ and $E$ are independent. The second step is to use Fubini's theorem to reverse the order in which $X$ and $D$ are integrated out. The final result is that the expected test error is given by 
$$ 
\text{expected test error} = 
E_X \left[ E_{\cal D} \left[ E_\epsilon \left[
(Y - h)^2 
\right] \right] \right], $$
where I have dropped the dependence of $Y$ and $h$ on $X$, $\epsilon$ and $\cal D$ in the interests of clarity.
Reducible and irreducible error
We fix values of $X$ and $\cal D$ (and therefore $f$ and $h$) and compute the inner-most integral in the expected test error:
\begin{align*}
E_\epsilon \left[ (Y - h)^2 \right] 
& = E_\epsilon \left[ (f + \epsilon - h)^2 \right]\\
& = E_\epsilon \left[ (f-h)^2 + \epsilon^2 + 2\epsilon (f-h) \right]\\
& = (f-h)^2 + E_\epsilon\left[ \epsilon^2 \right] + 0 \\
& = (f-h)^2 + Var_\epsilon \left[ \epsilon \right].
\end{align*}
The last term remains unaltered by subsequent averaging over $X$ and $D$. It represents the irreducible error contribution to the expected test error.
The average of the first term, 
$E_X \left[ E_{\cal D} \left[ \left( f-h\right)^2 \right] \right]$,
is sometimes called the reducible error.
Decomposing the reducible error into 'bias' and 'variance'
We relax our constraint that $\cal D$ is fixed (but keep the constraint that $X$ is fixed) and compute the innermost integral in the reducible error:
\begin{align*}
E_{\cal D} \left[ (f-h)^2 \right] 
& = E_{\cal D} \left[ f^2 + h^2 - 2fh \right]\\
& = f^2 + E_{\cal D} \left[ h^2 \right] - 2f E_{\cal D} \left[h\right]\\
\end{align*}
Adding and subtracting $E_{\cal D} \left[ h^2 \right]$, and rearranging terms, we may write the right-hand side above as 
$$
\left( f - E_{\cal D} \left[ h \right] \right)^2 + Var_{\cal D} \left[ h \right].
$$
Averaging over $X$, and restoring the irreducible error, yields finally: 
$$
\boxed{
\text{expected test error} = 
E_X \left[ \left( f - E_{\cal D} \left[ h \right] \right)^2 \right]
+ E_X \left[ Var_{\cal D} \left[ h \right] \right] 
+ Var_\epsilon \left[ \epsilon \right].
}
$$
The first term is called the bias and the second term is called the variance. 
The variance component of the expected test error is a consequence of the finite size of the training data sets. In the limit that training sets contain an infinite number of data points, there are no fluctuations in $h$ among the training sets and the variance term vanishes. Put another way, when the size of the training set is large, the expected test error is expected to be solely due to bias (assuming the irreducible error is negligible). 
More info
An excellent exposition of these concepts and more can be found here.
