Nicholas Lewin-Koh wrote:
Hi, I am writing some basic smoothers in R for cleaning some spectral data. I wanted to see if I could get close to matlab for speed, so I was trying to compare SparseM with Matrix to see which could do the choleski decomposition the fastest.
Here is the function using SparseM difsm <- function(y, lambda, d){ # Smoothing with a finite difference penalty # y: signal to be smoothed # lambda: smoothing parameter # d: order of differences in penalty (generally 2) # based on code by Paul Eilers 2002 require(SparseM) m <- length(y) E <- as(m,"matrix.diag.csr") D <- diff(E,differences=d) B <- E + (lambda * t(D)%*%D) z <- solve(B,y) z }
To do this in Matrix, I am not sure how to get an Identity matrix of dimension m. From looking at the documentation I would think that E<-new("cscMatrix", nrow=m, i=integer(m),x=rep(1,m),p=0:(m-1))) Should do what I want, but it fails?
m<-10 E<-new("cscMatrix", nrow=m, i=integer(m),x=rep(1,m),p=0:(m-1))
Error in initialize(value, ...) : Invalid names for slots of class cscMatrix: nrow
E<-new("cscMatrix", i=integer(m),x=rep(1,m),p=0:(m-1))
Error in validObject(.Object) : Invalid "cscMatrix" object: last element of slot p must match length of slots i and x
Granted I am not very fascile with the S4 classes, so I may be doing something stupid. Also to do the next step there is no diff method for matrices in Matrix, that would be nice :)
I guess the last step would be easy z <- solve((E + (lambda * crossprod(D))),y)
But I can't get the Identity matrix???
Thanks for the help, it is probably obvious, but I am missing it.
It is not really that obvious. Current versions of the Matrix package have limitations, especially with respect to constructors. Essentially that Matrix package implements what I needed at the time.
The easiest way to generate a sparse matrix is to create a tripletMatrix object and coerce it to a cscMatrix, or in your case perhaps an sscMatrix (symmetric sparse column-oriented matrix) object.
An identity matrix could be generated by
m <- 10
E <- as(new("tripletMatrix", Dim = as.integer(c(m, m)), i = 1:m,
j = 1:m, x = rep(1, m)), "cscMatrix")(I hope that is correct. I am writing this email in a system that does not do parenthesis matching for me so I may have mistakes in there.)
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
