[Bioc-devel] (no subject)

[Gordon K Smyth]

Hi Markus,

If you're worried about efficiency, you might try the loessFit() function 
in the limma package.  I suspect it might be faster, as well as more 
accurate, because it would avoid the need to subset and interpolate.

Cheers
Gordon

> Date: Tue, 24 Jun 2008 16:40:00 +0200 (CEST)
> From: lgautier at altern.org
> Subject: Re: [Bioc-devel] normalize.loess effizient method
> To: schmidb at ibe.med.uni-muenchen.de
> Cc: bioc-devel at stat.math.ethz.ch
> Message-ID: <57511.90.27.82.241.1214318400.squirrel at www.altern.org>
> Content-Type: text/plain;charset=iso-8859-1
>
>
> Thanks a lot Markus.
>
> To facilitate our work (review and inclusion), could you directly post a
> patch for your changes ?
>
> (Something like svn diff > fixnormalizeloess.patch should make it).
>
>
> Thanks,
>
>
> L.
>
>
>
>
>
>> Hello,
>>
>> I parallelized the normalize loess function (see affyPara package).
>> Thereby I found a memory inefficient implementation of normalize.loess.
>> The matrices fs and newdata are of the size of the complete intensity
>> matrix and not necessary.  I removed this two matrices and the oldfs
>> matrix. Attached the changed code.
>>
>> Best
>> Markus
>>
>>
>> normalize.loess <- function(mat, subset=sample(1:(dim(mat)[1]),
>> min(c(5000, nrow(mat)))),
>>         epsilon=10^-2, maxit=1, log.it=TRUE, verbose=TRUE, span=2/3,
>>         family.loess="symmetric"){
>>
>>     J <- dim(mat)[2]
>>     II <- dim(mat)[1]
>>     if(log.it)
>>         mat <- log2(mat)
>>
>>     change <- epsilon +1
>>     iter <- 0
>>     w <- c(0, rep(1,length(subset)), 0) ##this way we give 0 weight to the
>>     ##extremes added so that we can interpolate
>>
>>     while(iter < maxit){
>>         iter <- iter + 1
>>         means <- matrix(0,II,J) ##contains temp of what we substract
>>
>>         for (j in 1:(J-1)){
>>             for (k in (j+1):J){
>>                 y <- mat[,j] - mat[,k]
>>                 x <- (mat[,j] + mat[,k]) / 2
>>                 index <- c(order(x)[1], subset, order(-x)[1])
>>                 ##put endpoints in so we can interpolate
>>                 xx <- x[index]
>>                 yy <- y[index]
>>                 aux <-loess(yy~xx, span=span, degree=1, weights=w,
>> family=family.loess)
>>                 aux <- predict(aux, data.frame(xx=x)) / J
>>                 means[, j] <- means[, j] + aux
>>                 means[, k] <- means[, k] - aux
>>                 if (verbose)
>>                     cat("Done with",j,"vs",k," in iteration ",iter,"\n")
>>             }
>>         }
>>         mat <- mat - means
>>         change <- max(colMeans((means[subset,])^2))
>>
>>         if(verbose)
>>             cat(iter, change,"\n")
>>
>>     }
>>
>>     if ((change > epsilon) & (maxit > 1))
>>         warning(paste("No convergence after", maxit, "iterations.\n"))
>>
>>     if(log.it) {
>>         return(2^mat)
>>     } else
>>         return(mat)
>> }
>>
>> --
>> Dipl.-Tech. Math. Markus Schmidberger
>>
>> Ludwig-Maximilians-Universit?t M?nchen
>> IBE - Institut f?r medizinische Informationsverarbeitung,
>> Biometrie und Epidemiologie
>> Marchioninistr. 15, D-81377 Muenchen
>> URL: http://ibe.web.med.uni-muenchen.de
>> Mail: Markus.Schmidberger [at] ibe.med.uni-muenchen.de