Please allow me to expand on Prof. Venables' reply. Here is a comparison of various approaches to solving (or inverting) a linear system:
> set.seed(123) > n <- 1000 > > amat <- matrix(rnorm(n*n), n, n) > > amat <- amat %*% t(amat) # a symmetric positive-definite matrix > > x <- 1:n # correct solution > > b <- c(amat %*% x) # R.H.S. > > # Default `solve' > system.time(ans1 <- solve(amat, b))[1] user.self 0.48 > > max(abs(x - ans1)) #maximum error [1] 4.657716e-08 > > # Using LAPACK's QR decomposition > system.time(ans2 <- solve(qr(amat, LAPACK=TRUE), b))[1] user.self 1.34 > > max(abs(x - ans2)) #maximum error [1] 1.047159e-07 > > # Using Cholesky decomposition > # Note: Cholesky decomp works only for symmetric positive-definite matrices > system.time({ + ch <- chol(amat) + ans3 <- backsolve(ch, forwardsolve(t(ch), b)) + })[1] user.self 0.36 > > max(abs(x - ans3)) #maximum error [1] 4.907562e-08 > > # Using SVD from MASS > system.time(ans4 <- c(ginv(amat, tol=1.e-14) %*% b))[1] user.self 10.04 > max(abs(x - ans4)) #maximum error [1] 6.043251e-08 > > # Using SVD by hand > system.time({ + sv <- svd(amat) + ans5 <- c(sv$v %*% diag(1/sv$d) %*% crossprod(sv$u, b)) + })[1] user.self 8.18 > > max(abs(x - ans5)) #maximum error [1] 5.131403e-08 > >From this we can learn a few things: 1. SVD, as Prof. Venables said, is very slow, although this really matters only in larger matrices of order 10^3 or greater. 2. Cholesky decomposition approach is fastest, but is limited to SPD matrices. 3. All methods are reasonably accurate for SPD (and non-singular) matrices. But when the system is ill-conditioned or singular, only `ginv' approach would work. It provides the "minimum-norm" solution. Hope this helps, Ravi. ____________________________________________________________________ Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvarad...@jhmi.edu ----- Original Message ----- From: bill.venab...@csiro.au Date: Sunday, July 4, 2010 8:50 pm Subject: Re: [R] if using ginv function, does it mean there is no need to use solve function any more? To: rprojecth...@gmail.com, r-help@r-project.org > ginv() is slower than solve(). This is the price you pay for more generality. > > -----Original Message----- > From: r-help-boun...@r-project.org [ On Behalf Of song song > Sent: Monday, 5 July 2010 10:21 AM > To: r-help@r-project.org > Subject: [R] if using ginv function, does it mean there is no need to > use solve function any more? > > since ginv can deal with both singular and non-singular conditions, > is there > any other difference between them? > > if I use ginv only, will be any problem? > > thanks > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. > > ______________________________________________ > R-help@r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.