[R] if using ginv function, does it mean there is no need to use solve function any more?

[Ravi Varadhan]

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: rvaradhan at jhmi.edu


----- Original Message -----
From: Bill.Venables at 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: rprojecthelp at gmail.com, r-help at r-project.org


> ginv() is slower than solve().  This is the price you pay for more generality.
>  
>  -----Original Message-----
>  From: r-help-bounces at r-project.org [ On Behalf Of song song
>  Sent: Monday, 5 July 2010 10:21 AM
>  To: r-help at 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 at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.