R/getD.R

In qhengncsu/BNRPxMCMC: Bayesian l1 Trend Filtering via Proximal MCMC

#' Compute k+1-th order differencing matrix
#' @param k order of the desired piecewise polynomial; 0 constant, 1 linear, 2 quadratic
#' @param n length of signal we want to operate on
#' @param x predictor (must be strictly increasing!)
#' @export
getD <- function(k, n, x=NULL){
  if(is.null(x)){
    x <- 1:n
  }
  diags <- list(rep(-1,n),rep(1,n))
  D <- Matrix::bandSparse(n-1,n,k=c(0,1),diag=diags,symm=F)
  if(k>=1){
    for(i in 1:k){
      leftD <- Matrix::bandSparse(n-i-1,n-i,k=c(0,1),diag=diags,symm=F)
      xdiag <- Matrix::Diagonal(n-i,i/diff(x,lag=i))
      D <- leftD %*% xdiag %*% D
    }}
  return(D)
}