R/leftMultiplyByXt.R

In jointseg: Joint Segmentation of Multivariate (Copy Number) Signals

#' leftMultiplyByXt
#' 
#' Compute X'*Y where X is the n*(n-1) design matrix for the weighted group
#' fused Lasso, with weights defined by the vector w, and Y is any n*p matrix.
#' The computation is done in O(np).
#' 
#' This implementation is derived from the MATLAB code of Vert and Bleakley:
#' \url{http://cbio.ensmp.fr/GFLseg}.\
#' 
#' Contrary to \code{\link{getUnivStat}} it does not handle missing values.
#' 
#' @param Y A n*p matrix
#' @param w (n-1)*1 vector of weights
#' @param verbose A \code{logical} value: should extra information be output ?
#' Defaults to \code{FALSE}.
#' @return The (n-1)*p matrix equal to X'*Y
#' @author Morgane Pierre-Jean and Pierre Neuvial
#' @references Bleakley, K., & Vert, J. P. (2011). The group fused lasso for
#' multiple change-point detection. arXiv preprint arXiv:1106.4199.\ Vert, J.
#' P., & Bleakley, K. (2010). Fast detection of multiple change-points shared
#' by many signals using group LARS. Advances in Neural Information Processing
#' Systems, 23, 2343-2351.
#' @keywords internal
#' @examples
#' 
#' Y <- matrix(rnorm(20), ncol=2)
#' C <- leftMultiplyByXt(Y)
#' 
#' @export leftMultiplyByXt
leftMultiplyByXt <- function(Y,
    w=defaultWeights(nrow(Y)),
    verbose=FALSE){
    n <- as.numeric(dim(Y)[1])
    p <- ncol(Y)
    u <- apply(Y, 2, cumsum)
    if(length(w)!= (n-1)) {
        stop("w needs to be of length nrow(Y)-1")
    }
    C <- apply(u, 2, function(x){
        w*((1:(n-1))*x[n]/n-x[1:(n-1)])
    })
    dim(C) <- c(n-1,p)  ## so that the code also works with n=2...
    return(C)
}

############################################################################
## HISTORY:
## 2013-01-30
## o Now an internal function (not exported anymore).
## o Added 'jointSeg:::' to example because function is not exported anymore.
## 2012-12-27
## o Some code and doc cleanups.
## 2012-12-06
## o replaced the function defaultWeigths by a (n-1)*1 vector of weigths 
## o BUG FIX: when n=2 the result would loose its 'dim' attribute.
## 2012-09-13
## o Some code cleanups.
## 2012-08-13
## o Created.
############################################################################