Why is error of OLS not zero? Consider ordinary least squares (OLS): We have $n$ real datapoints $\mathbf{x}\in\mathbb{R}^d$ ($d$ features) organized in an ($n \times d$) matrix $X = (\mathbf{x}_1, \dots, \mathbf{x}_n)^T$ and we want to perform linear regression on $\mathbf{y} \in \mathbb{R}^n$ using a parameter vector $\boldsymbol{\theta} \in \mathbb{R}^d$. We write the loss function
\begin{equation}
L(\boldsymbol{\theta}) = \lVert \mathbf{y} - X\boldsymbol{\theta}\rVert^2 \tag{1}
\end{equation}
and set $\partial_\boldsymbol{\theta}L = 0$ to get the normal equations yielding the optimal parameter vector
\begin{equation}
\boldsymbol{\theta}^* = (X^T X)^{-1}X^T \mathbf{y} \tag{2}
\end{equation}
However, it seems to me that we can rearrange this using SVD: Write $X = U \Sigma V^T$ with orthogonal $n\times n$ matrix $U$, orthogonal $d\times d$ matrix $V$, and diagonal positive matrix $\Sigma$, then
\begin{align}
(X^T X)^{-1}X^T &= (V \Sigma^T U^T U \Sigma V^T)^{-1} V \Sigma^T U^T \tag{3a-e}
\\&= (V \Sigma^T \Sigma V^T)^{-1} V \Sigma^T U^T
\\&= V \Sigma^{-1} U^T
\\&= V \left(\Sigma^T U^T U (\Sigma^T)^{-1}\right)\Sigma^{-1}U^T
\\&= V \Sigma^T U^T (U \Sigma \Sigma^T U^T)^{-1}
\\&= X^T (X X^T)^{-1}
\end{align}
plugging (3e) into (2) gives $\boldsymbol{\theta}^* = X^T(X X^T)^{-1} \mathbf{y}$, and plugging this into (1) gives
\begin{align}
L(\boldsymbol{\theta}^*) = \lVert \mathbf{y} - X\left( X^T(X X^T)^{-1}\mathbf{y}\right)\rVert^2 = 0 \tag{4}
\end{align}
But this is nonsense, we can't possibly fit arbitrary $\mathbf{y}$ with a line. The only thing I can think of is that I'm somehow implicitly assuming $n < d$, but I can't see why this won't work for $n > d$.

My question is, why is this reasoning incorrect: what assumptions for OLS have I violated and what specific mathematical steps become invalid because of those errors?
 A: The major issue in your proof is $ (X^TX)^{-1}X^T = X^T(XX^T)^{-1} $
One obvious proof of the statement does not hold is by choosing a random matrix and then compute value on both sides. You would find they are different.
(I try X= np.array([[1,2,3],[4,5,6]])).
The reason of your proof goes wrong is because when you performing SVD, the dimension of matrix $\sum$ is $ n \times d$, which is not square if $ d<n $. This means that the firse line should be $ (X^TX)^{-1}X^T = (V \sum^T U^T...) $ instead of $ (V \sum U^T...)$
Also, the quantity $ V \sum^{-1} U^T$ does not exists at all if $d<n$, since there is not inverse for non-square matrix. This lead to false result.
The only case your proof holds is when the inverse exist, which is the case of $ n=d$ and we would expect a perfect OLS fit too.
A: How do you go from 3b to 3c?
You seem to be using the following $$(V\Sigma\Sigma V^T)^{-1}V\Sigma = V\Sigma^{-1}$$ But, that is only allowed when $\Sigma$ is square (ie when $n=p$). In that case $$(V\Sigma\Sigma V^T)^{-1} = (V^T)^{-1}\Sigma^{-1}\Sigma^{-1} V^{-1}$$
For $p=n$ it is correct that your cost function is zero.
