Greetings,
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
as
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)
and
A <- rbind(
c(1,1,1,1,1),
c(1,2,2,2,2),
c(1,2,3,3,3)
)
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?
Sincerely,
Michael Meyer
