In OLS, $X_{m\times n}$ matrix is the form $\begin{bmatrix}\mathbf x_1 & \dots & \mathbf x_{n-1} & \mathbf 1\end{bmatrix}$ where each element represents an $m\times 1$ vector. First $n-1$ are the feature vectors and the last one is all-1 vector that is to be multiplied with the bias. So, the formula for OLS solves for the bias already using the model:

$$\hat y = X\begin{bmatrix}\mathbf w\\ b\end{bmatrix}+\epsilon$$

In OLS, $X_{m\times n}$ matrix is the form $\begin{bmatrix}\mathbf x_1 & \dots & \mathbf x_{n-1} & \mathbf 1\end{bmatrix}$ where each element represents an $m\times 1$ vector. First $n-1$ are the feature vectors and the last one is all-1 vector that is to be multiplied with the bias. So, the formula for OLS solves for the bias already using the model:

$$y = X\begin{bmatrix}\mathbf w\\ b\end{bmatrix}+\epsilon$$