[R] Matrix: Help with syntax and comparison with SparseM

[Douglas Bates]

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