What's wrong with this derivation of the GLS of a linear model? Question:

Find the GLS estimator $\hat{\beta}$ for the linear model
$$Y=\beta X+\epsilon$$ where $\mathbb{E}\epsilon=0$ and
$\text{Cov}(\epsilon)=\sigma^2\Sigma$, where $\Sigma$ is positive
definite.

From Wikipedia, the correct answer is $$\hat{\beta}=(X^T\Sigma^{-1}X)^{-1}X^T\Sigma^{-1}Y,$$ but what is wrong with this derivation?
Set $B=Y-X\beta$, then we wish to minimise $f(B)=B^T\Sigma^{-1}B$. Now $\nabla f=2\Sigma^{-1}B$ (I'll skip the working here,  but I can insert it if this might be the source of error), so we have $$\nabla f=0\iff B=0\iff Y=X\beta\iff X^TY=(X^TX)\beta\iff \beta=(X^TX)^{-1}X^TY$$
This doesn't seem correct, since it is a significant simplification of the above result (it's just the OLS estimator). Where have I gone wrong?
Note: If I don't cancel the $\Sigma^{-1}$ in $\Sigma^{-1}B=0$ and carry out the same working I get the right answer, but I don't see why I can't cancel this.
 A: The vector $B$, which may be thought of as a vector of residuals, has expectation 0. In fact, the cross of $B$ with any column of the design matrix $X$ is 0 because regression is a projection. But as @whuber points out, saying $B=0$ is probably interpreted to mean $B^T = [0, 0, \ldots, 0]$ as a $n \times 1$ vector. This means you have perfect fit. In this case it doesn't matter what $\Sigma$ is, and that's why you're able to prove that OLS is "optimal" (any arbitrary weighting would also be "optimal" in this case, so the WLS optimal solution is not unique).
The proof that the WLS is a best linear unbiased estimator is by contradiction. We do not directly optimize the loss function. Your step $\nabla f \implies B=0$ needs serious consideration. If $f = B^T \Sigma^{-1} B$ can you prove with a lemma or reference that actually $\nabla f = 2 \Sigma{-1} B $? Your proposed $f$ is a $p \times p$ matrix. If you derive with respect to a vector, shouldn't it be a $p \times p \times p$ manifold? It's far too complicated to construct the WLS by optimization when the answer is fairly evident from the design.
Seber & Lee is a good reference.
A: You're writing your objective function as a function of the error, not of the parameter $\beta$, so that's why it's behaving weirdly.
I'll use $b$ as a candidate parameter value and $e_b$ as the corresponding residual vector. The loss then is
$$
f(e_b) = e_b^T\Sigma^{-1}e_b = (y - Xb)^T\Sigma^{-1}(y - Xb).
$$
Differentiating w.r.t. $e_b$ we do indeed get
$$
\nabla_e f = 2 \Sigma^{-1} e_b
$$
but the issue is that $e_b$ is not free to vary over all of $\mathbb R^n$. This means that in general we can't solve $\nabla_e f = \mathbf 0$.
If we instead differentiate $f(e_b)$ w.r.t. $b$ we get
$$
\nabla_b f = -2 y^T\Sigma^{-1}X + 2 X^T\Sigma^{-1}Xb
$$
and now $b$ is free to vary over all of $\mathbb R^p$ so this can be solved for zero (uniquely, when $X$ is full column rank)
