What conditioning assumptions are hidden in the bias-variance tradeoff proof? I've read through this question and its answers quite a few times. I'm curious to know what conditioning assumptions are used that are hidden from the derivation.
Here's what I mean: let's say I have a parameter $\theta$ and an estimator of $\theta$ given by $\widehat{\theta}$. I define the bias of $\widehat{\theta}$ by
$$\text{bias}(\widehat{\theta} \mid \theta) = \mathbb{E}[\widehat{\theta} \mid \theta] - \theta\text{.}$$
Notice how in the above I am specifically indicating that $\theta$ is fixed by conditioning on it.
Additionally, one has that the mean-squared error of $\widehat{\theta}$ is given by
$$\text{MSE}(\widehat{\theta} \mid \theta) = \mathbb{E}[(\widehat{\theta} - \theta)^2 \mid \theta] = \text{bias}^2(\widehat{\theta} \mid \theta) + \text{Var}(\widehat{\theta} \mid \theta)\text{.}$$

Going back to the problem at hand, suppose we have a training set $\{(\mathbf{x}_i, y_i)\}_{i=1}^{N}$ where $\mathbf{x}_i$ is a real-valued vector so that there exists a relationship
$$y_i = f(\mathbf{x}_i) + \epsilon$$
so that $\mathbb{E}[\epsilon] = 0$ and $\text{Var}(\epsilon) = \sigma^2 > 0$.
Suppose we have have an estimator $\hat{f}$ of $f$, so that  $\hat{f}(\mathbf{x}^{\prime}) = \hat{y}^{\prime}$ attempts to estimate $y^{\prime} = f(\mathbf{x}^{\prime}) + \epsilon$ at a new point $\mathbf{x}^{\prime}$. We aim to minimize the mean-squared error of $\hat{y}^{\prime}$, which should be given by
$$\begin{align}
\mathbb{E}[(\hat{y}^{\prime} - y^{\prime})^2 \mid y^{\prime}] &= \mathbb{E}[(\hat{f}(\mathbf{x}^{\prime}) - (f(\mathbf{x}^{\prime}) + \epsilon))^2 \mid y^{\prime}] \\
&= \mathbb{E}[(\hat{f}(\mathbf{x}^{\prime}) - f(\mathbf{x}^{\prime}))^2 + \epsilon^2 + 2\epsilon(\hat{f}(\mathbf{x}^{\prime}) - f(\mathbf{x}^{\prime})) \mid y^{\prime}] \\
&= \mathbb{E}[(\hat{f}(\mathbf{x}^{\prime}) - f(\mathbf{x}^{\prime}))^2 \mid y^{\prime}]  + \mathbb{E}[\epsilon^2 \mid y^{\prime}] + 2 \cdot \mathbb{E}[\epsilon(\hat{f}(\mathbf{x}^{\prime}) - f(\mathbf{x}^{\prime})) \mid y^{\prime}]\text{.} 
\end{align}$$
Now here, I don't understand what's going on.
My understanding is that $\mathbb{E}[(\hat{f}(\mathbf{x}^{\prime}) - f(\mathbf{x}^{\prime}))^2 \mid y^{\prime}] = \text{MSE}(\hat{f}(\mathbf{x}^{\prime}) \mid f(\mathbf{x}^{\prime}))$, which decomposes into the squared bias and variance decomposition I mentioned above. But surely conditioning on $y^{\prime}$ is not the same as conditioning on $f(\mathbf{x}^{\prime})$?
Additionally, $\mathbb{E}[\epsilon^2 \mid y^{\prime}] = \text{Var}(\epsilon^2 \mid y^{\prime})$ (if we assume $\mathbb{E}[\epsilon \mid y^{\prime}] = 0$) should apparently equal $\sigma^2$. But we're told that the unconditional variance of $\epsilon$ is $\sigma^2$ - how does this apply when conditioning on $y^{\prime}$?
Additionally, $\mathbb{E}[\epsilon(\hat{f}(\mathbf{x}^{\prime}) - f(\mathbf{x}^{\prime})) \mid y^{\prime}]$ should equal $0$. The only way I can see this making sense is if conditioning on $y_i$ is equivalent to conditioning on $f(\mathbf{x}^{\prime})$, and in addition, $\epsilon$ and $\hat{f}(\mathbf{x}^{\prime})$ are conditionally independent given $f(\mathbf{x}^{\prime})$. How does this make sense?
 A: Out of all of the references I dug through, the clearest exposition of this derivation is found in Bishop's Pattern Recognition and Machine Learning (2006). I roughly follow the exposition of sections 1.5.5 (p. 46-7) and section 3.2 (pp. 147-9).
Bishop assumes $\mathbf{X}$, $Y$ are continuous random vectors/variables. Let $\mathbf{X} \in \mathbb{R}^k$ be the new input vector. Squared-error loss between a value $Y := Y(\mathbf{X})$ and a proposed estimator $\hat{Y}(\mathbf{X})$ is given by
$$\int_{\mathbb{R}}\int_{\mathbb{R}^k}\left[\hat{Y}(\mathbf{x}) - y \right]^2f_{\mathbf{X}, Y}(\mathbf{x}, y)\text{ d}\mathbf{x}\text{ d}y = \mathbb{E}\left[\left(\hat{Y}(\mathbf{X})-Y\right)^2 \right]$$
where $f_{\mathbf{X}, Y}(\mathbf{x}, y)$ denotes the joint density of $(\mathbf{X}, Y)$.
We know from standard results that the above quantity  leads to the decomposition
$$\begin{align}
\mathbb{E}\left[\left(\hat{Y}(\mathbf{X})-Y\right)^2 \right] &= \mathbb{E}\left\{\left[\hat{Y}(\mathbf{X}) - \mathbb{E}\left[Y \mid \mathbf{X} \right]\right]^2\right\} + \mathbb{E}\left\{\left[\mathbb{E}\left[Y \mid \mathbf{X} \right] - Y\right]^2\right\} \tag{1}
\end{align}$$
and thus we obtain, letting $g(\mathbf{x}) = \mathbb{E}[Y \mid \mathbf{X} = \mathbf{x}]$,
$$
\mathbb{E}\left[\left(\hat{Y}(\mathbf{X})-Y\right)^2 \right] = \int_{\mathbb{R}^k}\left[\hat{Y}(\mathbf{x}) - g(\mathbf{x}) \right]^2f_{\mathbf{X}}(\mathbf{x}) \text{ d}\mathbf{x}+ \int_{\mathbb{R}}\int_{\mathbb{R}^k}\left[g(\mathbf{x})  - y  \right]^2f_{\mathbf{X}, Y}(\mathbf{x}, y)\text{ d}\mathbf{x}\text{ d}y\text{.}$$
We next observe that
$$\left(\mathbb{E}\left[Y \mid \mathbf{X} \right] - Y\right)^2 = \left[f(\mathbf{X})-f(\mathbf{X}) -\epsilon\right]^2 = \epsilon^2$$
with the assumption that $\mathbb{E}[\epsilon \mid \mathbf{X}] = 0$.
Thus, in $(1)$, the term
$$\mathbb{E}\left\{\left[\mathbb{E}\left[Y \mid \mathbf{X} \right] - Y\right]^2\right\} = \mathbb{E}[\epsilon^2] = \sigma^2$$
if we assume either one of the following are true: either $\mathbb{E}[\epsilon^2] = \sigma^2$ or $\mathbb{E}[\epsilon^2 \mid \mathbf{X}] = \sigma^2$.
What we haven't taken into account here at all yet is the use of the training data. Let $\mathcal{T} = \{(\mathbf{X}_i, Y_i)\}_{i=1}^{N} = \{\mathbf{P}_i\}_{i=1}^{N}$ denote the training data, which we consider continuous random vectors/variables. In this case, our estimator $\hat{Y}(\mathbf{X})$ will not only depend on $\mathbf{X}$, but the given training data $\mathcal{T}$ as well, so we denote $\hat{Y}(\mathbf{X} \mid \mathcal{T}) := \hat{Y}(\mathbf{X})$.
Following notation similar to Bishop's, let
$$\begin{align}
&\mathbb{E}_{\mathcal{T}}\left[\hat{Y}(\mathbf{X} \mid \mathcal{T})\right] =\\
&\int_{\mathbb{R}^{k+1}} \cdots \int_{\mathbb{R}^{k+1}} \hat{Y}(\mathbf{X} \mid \mathbf{P}_1 = \mathbf{p}_1, \dots, \mathbf{P}_N = \mathbf{p}_N)f_{\mathbf{P}_1, \dots, \mathbf{P}_N}(\mathbf{p}_1, \dots, \mathbf{p}_N)\text{ d}\mathbf{p}_1 \cdots \text{ d}\mathbf{p}_N
\end{align}$$
denote the expectation of $\hat{Y}(\mathbf{X} \mid \mathcal{T})$ over all possible training sets.
With some algebra laid out in (3.39) and (3.40) of Bishop, one can build from $(1)$ and apply double expectation to obtain
$$\mathbb{E}\left\{\left[\hat{Y}(\mathbf{X} \mid \mathcal{T}) - \mathbb{E}\left[Y \mid \mathbf{X} \right]\right]^2\right\} =\mathbb{E}\left[\mathbb{E}\left(\mathbb{E}\left\{\left[\hat{Y}(\mathbf{X} \mid \mathcal{T}) - \mathbb{E}\left[Y \mid \mathbf{X} \right]\right]^2 \mid \mathbf{X}, \mathcal{T}\right\} \mid \mathbf{X}\right)\right] $$
add and subtract $\mathbb{E}_{\mathcal{T}}\left[\hat{Y}(\mathbf{X} \mid \mathcal{T})\right]$ inside the squared term (similar to what is done in the link above) and use the work laid out in Bishop to show that the above may be decomposed into, if $h(\mathbf{x}) = \mathbb{E}_{\mathcal{T}}\left[\hat{Y}(\mathbf{X} \mid \mathcal{T})\right](\mathbf{x})$,
$$\int_{\mathbb{R}^k}\left\{ h(\mathbf{x}) - \mathbb{E}[Y \mid \mathbf{X} = \mathbf{x} ]\right\}^2f_{\mathbf{X}}(\mathbf{x})\text{ d}\mathbf{x} + \int_{\mathbb{R}^k}\mathbb{E}_{\mathcal{T}}\left[\left\{\hat{Y}(\mathbf{x}) - h(\mathbf{x}) \right\}^2 \right]f_{\mathbf{X}}(\mathbf{x})\text{ d}\mathbf{x}$$
The first term above is the squared bias term, and the second term above is the variance term. If $\mathbf{X}$ is deterministic, the PDFs can be ignored, and the integrals are just sums over one point.
