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 [EMAIL PROTECTED] 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
______________________________________________
[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
______________________________________________
[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
