In my opinion the documentation or behaviour of qr.solve, qr.coef, qr.resid, and qr.fitted is not easily comprehensible and unfortunate.
We all know that a linear system Ax=b can have 0, one or infinitely many solutions. To treat all these cases uniformly we can rephrase the problem
x = argmin_u||Au-b||,
where ||.|| denotes the Euclidean norm. There is then exactly one natural and distinguished solution x, the minimizer x which is itself of minimal
Euclidean Norm. So if we want to return only one solution, I think we can agree that this should be it.
In fact this very solution can be computed from the QR-decomposition of either A (overdetermined system) or t(A) (underdetermined system).
I tried qr.solve on the underdetermined system Ax=b with
b <- c(3,5,7)
A <- rbind(
The system has infinitely many solutions. The minimal norm solution is x=c(1,0,2/3,2/3,2/3).
But qr.solve(A,b) yielded the solution x=(1,0,2,0,0) which is destinguished only by being sparse and I do not think qr.solve
tries to compute the sparsest solution. So what does qr.solve do in case of an underdetermined system?
It is not documented.
Then I tried to figure out what qr.coef, qr.resid, and qr.fitted do. I had to do actual experiments to figure out that it seems to solve the
problem x=argmin_u||Au-y|| with
qr.coef(A,y) = x but which x when there are infinitely many?
qr.fitted(A,y) = Ax
qr.resid(A,y) = y-Ax
but this is certainly not evident from the language in the documentation which conflates qr(x,...) with solving the system Ax=b then states:
"The functions qr.coef, qr.resid, and qr.fitted return the coefficients, residuals and fitted values obtained
when fitting y to the matrix with QR decomposition qr."
Since when do we call solving a system of equations "fitting the right hand side to the matrix ..." or call the solution x "the coefficients"
(which more usually are the elements of A) or introduce the "fitted values" with no definition?
Moreover the language does not fit the underdetermined case Ax=y, where we need the QR-decomposition of t(A) and not of A to compute the
minimizer x = argmin_u||Au-y|| which is itself of minimal norm.
Or, maybe this is not at all what these functions are doing.
But then, what is it and should this not be evident from the documentation?