R/SURE_MSEthresh.R

In fabnavarro/gasper: Graph Signal Processing

#' Stein's Unbiased Risk Estimate with MSE
#'
#' Adaptive Threshold Selection Using Principle of SURE with the inclusion of Mean Squared Error (MSE) for comparison.
#'
#' @export SURE_MSEthresh
#' @param wcn Numeric vector of the noisy spectral graph wavelet coefficients.
#' @param wcf Numeric vector of the true spectral graph wavelet coefficients.
#' @param thresh Numeric vector of threshold values.
#' @param diagWWt Numeric vector of weights typically derived from the diagonal elements of the wavelet frame matrix.
#' @param beta A numeric value specifying the type of thresholding to be used, for example:
#'  \itemize{
#'             \item 1 for soft thresholding.
#'             \item 2 for James-Stein thresholding.
#'             }
#' @param sigma A numeric value representing the standard deviation (sd) of the noise.
#' @param hatsigma An optional numeric value providing an estimate of the noise standard deviation (default is NA).
#' @param policy A character string determining the thresholding policy. Valid options include:
#' \itemize{
#'     \item "uniform" for a global threshold applied uniformly across all coefficients.
#'     \item "dependent" for threshold values that adaptively depend on the corresponding `diagWWt` weights.
#'               }
#' @param keepwc A logical value determining if the thresholded wavelet coefficients should be returned (Default is TRUE).
#' @return A list containing:
#' \itemize{
#'         \item A dataframe with calculated MSE, SURE, and hatSURE values.
#'         \item Minima of SURE, hatSURE, and MSE, and their corresponding optimal thresholds.
#'         \item Thresholded wavelet coefficients (if \code{keepwc = TRUE}).
#'         }
#' @details
#' \code{SURE_MSEthresh} function extends the \code{SUREthresh} function by providing an MSE between the true coefficients and their thresholded versions for a given thresholding function \eqn{h}. This allows for a more comprehensive evaluation of the denoising quality in simulated scenarios where the true function is known.
#'
#' @seealso \code{\link{SUREthresh}}, \code{\link{GVN}}, \code{\link{HPFVN}}
#'
#' @examples
#' \dontrun{
#' # Compute the Laplacian matrix and its eigen-decomposition
#' L <- laplacian_mat(grid1$sA)
#' U <- eigensort(L)
#'
#' # Compute the tight frame coefficients
#' tf <- tight_frame(U$evalues, U$evectors)
#'
#' # Generate some noisy observation
#' n <- nrow(L)
#' f <- randsignal(0.01, 3, grid1$sA)
#' sigma <- 0.01
#' noise <- rnorm(n, sd = sigma)
#' tilde_f <- f + noise
#'
#' # Compute the transform coefficients
#' wcn <- analysis(tilde_f, tf)
#' wcf <- analysis(f, tf)
#'
#' # Compute the weights and use DJ trick for the SURE evaluation
#' diagWWt <- colSums(t(tf)^2)
#' thresh <- sort(abs(wcn))
#'
#' # Compute to optimal threshold
#' opt_thresh_u <- SURE_MSEthresh(wcn,
#'                                wcf,
#'                                thresh,
#'                                diagWWt,
#'                                beta=2,
#'                                sigma,
#'                                NA,
#'                                policy = "uniform",
#'                                keepwc = TRUE)
#' # Extract corresponding wavelet coefficients
#' wc_oracle_u <- opt_thresh_u$wc[, opt_thresh_u$min["xminMSE"]]
#' wc_SURE_u <- opt_thresh_u$wc[, opt_thresh_u$min["xminSURE"]]
#'
#' # Get the graph signal estimators
#' hatf_oracle_u <- synthesis(wc_oracle_u, tf)
#' hatf_SURE_u  <- synthesis(wc_SURE_u, tf)
#'
#' # Compare the perfomance according to SNR measure
#' round(SNR(f, hatf_oracle_u), 2)
#' round(SNR(f, hatf_SURE_u), 2)
#' }
#'
#' @references
#' Donoho, D. L., & Johnstone, I. M. (1995). Adapting to unknown smoothness via wavelet shrinkage. Journal of the american statistical association, 90(432), 1200-1224.
#'
#' de Loynes, B., Navarro, F., Olivier, B. (2021). Data-driven thresholding in denoising with Spectral Graph Wavelet Transform. Journal of Computational and Applied Mathematics, Vol. 389.
#'
#' Stein, C. M. (1981). Estimation of the mean of a multivariate normal distribution. The annals of Statistics, 1135-1151.

SURE_MSEthresh <- function(wcn, wcf, thresh, diagWWt, beta = 2, sigma, hatsigma = NA, policy = "uniform", keepwc=TRUE) {
  nthresh <- length(thresh)
  erisk <- dof  <- MSE <- rep(0, nthresh)
  if(keepwc){
    wcs <- matrix(0, ncol=nthresh, nrow = length(wcn))
    if (policy == "uniform") {
      for (i in 1:nthresh) {
        wc <- betathresh(wcn, thresh[i], beta)
        wcs[,i] <- wc
        MSE[i] <- sum((wcf - wc)^2)
        erisk[i] <- sum((wc - wcn)^2)
        dof[i] <- 2*sum((abs(wcn) > thresh[i])*
                          (1 + (beta - 1)*thresh[i]^beta/abs(wcn)^beta)*diagWWt,na.rm=TRUE)
      }
    }
    else if (policy == "dependent") {
      for (i in 1:nthresh) {
        wc <- betathresh(wcn, thresh[i]*sqrt(diagWWt), beta)
        wcs[,i] <- wc
        MSE[i] <- sum((wcf - wc)^2)
        erisk[i] <- sum((wc - wcn)^2)
        dof[i] <- 2*sum((abs(wcn) > thresh[i]*sqrt(diagWWt))*
                          (1 + (beta - 1)*(thresh[i]*sqrt(diagWWt))^beta/abs(wcn)^beta)*diagWWt,na.rm=TRUE)
      }
    } else {
      warning("Allowable policy are listed above")
      print("uniform")
      print("dependent")
    }
    SURE <- erisk + dof*sigma^2
    hatSURE <- erisk + dof*hatsigma^2
    minMSE <- which.min(MSE)
    minSURE <- which.min(SURE)
    minhatSURE <- which.min(hatSURE)
    opthreshMSE <- thresh[minMSE]
    opthreshSURE <- thresh[minSURE]
    opthreshhatSURE <- thresh[minhatSURE]
    res <- list("wc"=wcs,
                "res"=data.frame(MSE = MSE, SURE = SURE, hatSURE = hatSURE),
                "min"=c(xminMSE = minMSE,
                        xminSURE = minSURE,
                        xminhatSURE = minhatSURE),
                "thr"=c(opthreshMSE = opthreshMSE,
                        opthreshSURE = opthreshSURE,
                        opthreshhatSURE = opthreshhatSURE))
  }
  else{
    if (policy == "uniform") {
      for (i in 1:nthresh) {
        wc <- betathresh(wcn, thresh[i], beta)
        MSE[i] <- sum((wcf - wc)^2)
        erisk[i] <- sum((wc - wcn)^2)
        dof[i] <- 2*sum((abs(wcn) > thresh[i])*
                          (1 + (beta - 1)*thresh[i]^beta/abs(wcn)^beta)*diagWWt,na.rm=TRUE)
      }
    }
    else if (policy == "dependent") {
      for (i in 1:nthresh) {
        wc <- betathresh(wcn, thresh[i]*sqrt(diagWWt), beta)
        MSE[i] <- sum((wcf - wc)^2)
        erisk[i] <- sum((wc - wcn)^2)
        dof[i] <- 2*sum((abs(wcn) > thresh[i]*sqrt(diagWWt))*
                          (1 + (beta - 1)*(thresh[i]*sqrt(diagWWt))^beta/abs(wcn)^beta)*diagWWt,na.rm=TRUE)
      }
    } else {
      warning("Allowable policy are listed above")
      print("uniform")
      print("dependent")
    }
    SURE <- erisk + dof*sigma^2
    hatSURE <- erisk + dof*hatsigma^2
    minMSE <- which.min(MSE)
    minSURE <- which.min(SURE)
    minhatSURE <- which.min(hatSURE)
    opthreshMSE <- thresh[minMSE]
    opthreshSURE <- thresh[minSURE]
    opthreshhatSURE <- thresh[minhatSURE]
    res <- list("res"=data.frame(MSE = MSE,
                                 SURE = SURE,
                                 hatSURE = hatSURE),
                "min"=c(xminMSE = minMSE,
                        xminSURE = minSURE,
                        xminhatSURE = minhatSURE),
                "thr"=c(opthreshMSE = opthreshMSE,
                        opthreshSURE = opthreshSURE,
                        opthreshhatSURE = opthreshhatSURE))
  }
  return(res)
}