[R] Sparse matrix methods
Thomas Lumley
tlumley at u.washington.edu
Thu Mar 14 17:17:15 CET 2002
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list