[R] Sparse matrix methods

[Thomas Lumley]

On Tue, 12 Mar 2002, Doug Nychka wrote:

>
> Does anyone know of contributions to R for solving sparse linear systems?
> In particular for spatial stats I am interested in solving large
> positive definite symmetric systems.

I have a matrix-free implementation of the conjugate gradient algorithm
(below). That is, to solve Ax=b you supply b and a function that computes
Ax.

I wrote this while working on my dissertation, exactly to solve large
symmetric systems from spatial statistics, but ended up using the Aztec
sparse matrix library from Sandia Labs, which has a wider set of methods.

	-thomas

Thomas Lumley			Asst. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle


cgsparse <-function(Atimes,b,guess,epsilon=0.0001,maxit=50,verbose=F,...){
## Conjugate gradient iterative solver.
## Atimes is a function that returns A%*%x for argument x
## Finds x such that Ax=b with relative error epsilon using
## at most maxit iterations starting from x=guess
##
## Algorithm from Axelsson "Iterative Solution Methods" p470
##
## This algorithm works much better on correlation matrices than
## covariance matrices -- it's worth rescaling first.
##
  x<-guess
  goodenough<-norm(b)*epsilon
  r<-b
  r<-Atimes(x,...)-r
  delta0<-norm(r)
  if (delta0<goodenough) return(cbind(x,r))
  d<- -r
  i<-0
  while(i<maxit){
    i<-i+1
    h<-Atimes(d,...)
    tau<-delta0/crossprod(d,h)
    x<-x+tau*d
    r<-r+tau*h
    delta1<-norm(r)
    if ((verbose) & (delta1>2*delta0))
      warning("?SPD matrix")
    if (delta1<goodenough) break
    beta<-delta1/delta0
    delta0<-delta1
    d<- -r+beta*d
  }
        return(cbind(x,r))

}


-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._