The core idea behind using residuals to update the $\beta$ coefficients is to move the coefficient estimates in a direction that reduces the difference between the actual and predicted values of $y$. This is based on the principle that the best-fitting model is the one that minimizes the sum of squared residuals (RSS).
By iteratively adjusting $\beta$ based on the residuals, the algorithm seeks to find the minimum RSS in a step-wise fashion. The adjustment vector $\hat{b}$ points in the direction that most rapidly decreases the RSS, given the current estimate $\hat{\beta}^{[m-1]}$. By moving in this direction, the algorithm iteratively reduces the error between the observed and predicted values.
Under suitable conditions (e.g., $X^TX$ is non-singular), this iterative process converges to the $\hat{\beta}$ that minimizes the RSS, which is the solution to the normal equations of linear regression:
$$X^TX\hat{\beta} = X^Ty$$
This iterative method can be more numerically stable and efficient than directly solving the normal equations, especially when $X^TX$ is close to singular or when $n$ is very large. By adjusting $\hat{\beta}$ iteratively, the method can avoid the numerical difficulties associated with directly inverting $X^TX$.
In the context of linear regression, the gradient of the RSS with respect to $\beta$ is given by $-2X^T(y-X\beta)$, and the Hessian matrix (the second derivative of RSS with respect to $\beta$) is $2X^TX$. The adjustment $\hat{b}$ is derived by applying a Newton-like step, assuming the Hessian is constant (which is true for linear regression), hence moving $\beta$ towards the minimum of RSS.
