I've been led to believe that it is in fact possible to convert the least squares approximation from their partial derivatives = 0 sum form to the matrix form $\overline{b}=(M^TM)^{-1}M^T\overline{y}$. Apparently one does not get the matrix M, but combined $M^TM$ and $M^T\overline{y}$.
How might one go about the job?
For example, let the model be: $a+b\cdot \sin x$ $$S(a,b) = \sum_1^n{(y_i-a-b\cdot \sin x_i)^2}$$ $$\frac{\partial S(a,b)}{\partial a}=\sum_1^n{2(y_i-a-b\cdot \sin x_i)(-1)}$$ $$\frac{\partial S(a,b)}{\partial b}=\sum_1^n{2(y_i-a-b\cdot \sin x_i)(-\sin x_i)}$$
$$\sum y = \sum a + \sum b\sin x$$ $$\sum y\sin x = \sum a\sin x + \sum b\sin x\sin x$$
Then what?