Skip to content

Matrix inverse

5 messages · Göran Broström, Douglas Bates, Brian Ripley

#
On Fri, 28 Apr 2000, Thomas Lumley wrote:

            
How about 'Ainv <- qr.solve(A)'?

I happened to read the help page for 'qr.solve' the other day, and there I
found that qr.solve(A, b) "is == but much better than solve(A) %*% b".
(I guess that 'better than' refers to numerical stability?)

Is qr.solve generally to be preferred to solve, which seems to be
indicated by the cited help page? And what is hidden behind solve?

G?ran
#
gb <gb at stat.umu.se> writes:
The "better than" refers to the fact that creating the inverse is
more-or-less equivalent to solving n systems of linear equations,
where n is the number of columns in A.  If n is large it does not make
sense to compute the solutions to n systems of equations in order to
evaluate the solution to one system of equations.
I think if you check the code for solve you will find that it usually
calls qr.solve.  qr.solve is just one method of solving a linear
system of equations.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
On 1 May 2000, Douglas Bates wrote:
I see! But then I think the help page is somewhat misleading. If it said
'...better than qr.solve(A) %*% b', there would be no doubt.
I tried
function(a, b, ...) UseMethod("solve")
[1] "solve.default" "solve.qr"     

I guess that I have to find the underlying  C  code to find out
what "solve.default" is, and when which method is used. I have
looked around in R-1.0.1/src/* with no great success. Can I
get a hint where to search?

Thanks, G?ran

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
gb <gb at stat.umu.se> writes:
solve.default is just another function in the base R library.
function(a, b, tol = 1e-7)
{
    if( !is.qr(a) )
	a <- qr(a, tol = tol)
    nc <- ncol(a$qr)
    if( a$rank != nc )
	stop("singular matrix `a' in solve")
    if( missing(b) ) {
	if( nc != nrow(a$qr) )
	    stop("only square matrices can be inverted")
	b <- diag(1,nc)
    }
    ## pre 0.63.3: b <- as.matrix(b)
    return(qr.coef(a,b))
}

As you can see, it checks if the first argument is a qr object and, if
not, it takes a qr decomposition.

It appears that solve.default, solve.qr, and qr.solve are identical.
I imagine the reason for identical functions under multiple names is 
for S-PLUS compatibility although I'm not sure.  Does anyone else
know?


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
#
On 1 May 2000, Douglas Bates wrote:

            
qr.solve does not exist in S-PLUS.  The other two do, and are different
in that solve.qr will only only work if A is already a QR decomposition.