Note: I've posted an expanded version of this answer on my website.
Would you kindly consider posting a similar answer with the actual R engine exposed?
Sure! Down the rabbit hole we go.
The first layer is
lm, the interface exposed to the R programmer. You can look at the source for this by just typing
lm at the R console. The majority of it (like the majority of most production level code) is busy checking of inputs, setting of object attributes, and throwing of errors; but this line sticks out
lm.fit(x, y, offset = offset, singular.ok = singular.ok,
lm.fit is another R function, you can call it yourself. While
lm conveniently works with formulas and data frame,
lm.fit wants matrices, so that's one level of abstraction removed. Checking the source for
lm.fit, more busywork, and the following really interesting line
z <- .Call(C_Cdqrls, x, y, tol, FALSE)
Now we are getting somewhere.
.Call is R's way of calling into C code. There is a C function, C_Cdqrls in the R source somewhere, and we need to find it. Here it is.
Looking at the C function, again, we find mostly bounds checking, error cleanup, and busy work. But this line is different
F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
REAL(coefficients), REAL(residuals), REAL(effects),
&rank, INTEGER(pivot), REAL(qraux), work);
So now we are on our third language, R has called C which is calling into fortran. Here's the fortran code.
The first comment tells it all
c dqrfit is a subroutine to compute least squares solutions
c to the system
c (1) x * b = y
(interestingly, looks like the name of this routine was changed at some point, but someone forgot to update the comment). So we're finally at the point where we can do some linear algebra, and actually solve the system of equations. This is the sort of thing that fortran is really good at, which explains why we passed through so many layers to get here.
The comment also explains what the code is going to do
c on return
c x contains the output array from dqrdc2.
c namely the qr decomposition of x stored in
c compact form.
So fortran is going to solve the system by finding the $QR$ decomposition.
The first thing that happens, and by far the most important, is
This calls the fortran function
dqrdc2 on our input matrix
x. Whats this?
c dqrfit uses the linpack routines dqrdc and dqrsl.
So we've finally made it to linpack. Linpack is a fortran linear algebra library that has been around since the 70s. Most serious linear algebra eventualy finds its way to linpack. In our case, we are using the function dqrdc2
c dqrdc2 uses householder transformations to compute the qr
c factorization of an n by p matrix x.
This is where the actual work is done. It would take a good full day for me to figure out what this code is doing, it is as low level as they come. But generically, we have a matrix $X$ and we want to factor it into a product $X = QR$ where $Q$ is an orthogonal matrix and $R$ is an upper triangular matrix. This is a smart thing to do, because once you have $Q$ and $R$ you can solve the linear equations for regression
$$ X^t X \beta = X^t Y $$
very easily. Indeed
$$ X^t X = R^t Q^t Q R = R^t R $$
so the whole system becomes
$$ R^t R \beta = R^t Q^t y $$
but $R$ is upper triangular and has the same rank as $X^t X$, so as long as our problem is well posed, it is full rank, and we may as well just solve the reduced system
$$ R \beta = Q^t y $$
But here's the awesome thing. $R$ is upper triangular, so the last linear equation here is just
constant * beta_n = constant, so solving for $\beta_n$ is trivial. You can then go up the rows, one by one, and substitute in the $\beta$s you already know, each time getting a simple one variable linear equation to solve. So, once you have $Q$ and $R$, the whole thing collapses to what is called backwards substitution, which is easy. You can read about this in more detail here, where an explicit small example is fully worked out.