R help - Jun 2004 - Matrix: Help with syntax and comparison with SparseM

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.

Nicholas

I'll defer to others on the best way to implement this in Matrix, but it
probably should be noted that for banded structure like this sort of
Whittaker problem more
specialized algorithms have a real advantage.  It  might
require very large problems to see this, however.


url:	www.econ.uiuc.edu/~roger        	Roger Koenker
email	rkoenker at uiuc.edu			Department of Economics
vox: 	217-333-4558				University of Illinois
fax:   	217-244-6678				Champaign, IL 61820

On Jun 25, 2004, at 2:30 AM, 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.
>
> Nicholas
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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

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