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?
Based on Wikipedia provided in the comments I got it as it was actually pretty close. So by taking the constants out of the summations, one can easily bring the equation system to first a vectorized form and finally to a matrix form.
$$\sum y = a\sum 1 + b\sum \sin x$$ $$\sum y\sin x = a\sum \sin x + b\sum \sin^2 x$$ And so: $$\begin{bmatrix}\sum y \\ \sum y\sin x\end{bmatrix} = a\begin{bmatrix}\sum 1 \\ \sum \sin x\end{bmatrix} + b\begin{bmatrix}\sum\sin x \\ \sum\sin^2x\end{bmatrix}$$ $$\begin{bmatrix}\sum y \\ \sum y\sin x\end{bmatrix} = \begin{bmatrix}\sum 1 & \sum\sin x \\ \sum\sin x & \sum \sin^2 x\end{bmatrix}\begin{bmatrix}a\\b\end{bmatrix}$$