# R/threshold.R In dfleis/morpca: Manifold Optimization for Robust PCA

#' Matrix Hard-Thresholding
#'
#' Computes the hard-thresholding of a given real-valued matrix \eqn{X}
#' according to thresholding percentile \eqn{\gamma}.
#'
#' @details The hard-thresholding procedure will set the \eqn{(i,j)}-th
#' entry of X to 0 if
#' \deqn{|X_{i,j}| > |X^{[\gamma]}_{i,.}|}
#' and
#' \deqn{|X_{i,j}| > |X^{[\gamma]}_{.,j}|,}
#' where \eqn{X_{i,.}} denotes the \eqn{i}-th row of \eqn{X}, and
#' \eqn{X_{.,j}} represents the \eqn{j}-th column of \eqn{X}. Additionally,
#' \eqn{X^{[\gamma]}_{i,.}} and \eqn{X^{[\gamma]}_{.,j}} represent the
#' \eqn{\gamma}-th percentile of the absolute values of the entries of the
#' \eqn{i}-th row and the \eqn{j}-th column, respectively.
#'
#' That is, the elements of \eqn{X} which simultaneously exceed the
#' \eqn{\gamma}-th most extreme fraction of entries (in terms of absolute
#' value) in the corresponding rows and columns of \eqn{X} will be removed
#' and replaced with 0.
#'
#' @param X Input real-valued data matrix.
#' @param gamma Thresholding percentile.
#' @param sparsity TO DO...
#'
#' @return Returns a real-valued matrix with \eqn{\gamma}-thresholded
#' entries.
#'
#' @export
threshold <- function(X, gamma, sparsity) {
# TO DO:
#   * Sparsity
#   * Possibly implement the apply() functionality within C++?
# n1 <- nrow(X); n2 <- ncol(X)
#
# t1 <- rep(1, n1); t2 <- rep(1, n2)
# X <- X * sparsity
#
# for (i in 1:n1) {
#   tt <- sort(abs(X[i, ]), decreasing = TRUE)
#   t1[i] <- tt[floor(gamma * sum(sparsity[i, ])) + 1]
# }
# for (j in 1:n2) {
#   tt <- sort(abs(X[, j]), decreasing = TRUE)
#   t2[j] <- tt[floor(gamma * sum(sparsity[, j])) + 1]
# }
#
# threshold1 <- abs(X) <= matrix(rep(t1, each = n2), ncol = n2, byrow = TRUE)
# threshold2 <- abs(X) <= matrix(rep(t2, each = n1), nrow = n1)
#
# X_thresholded <- X * (as.double(threshold1 + threshold2) >= 1)
#
# X_thresholded
X_abs <- abs(X)
row_pctls <- row_pctls_cpp(X_abs, 1 - gamma)
col_pctls <- col_pctls_cpp(X_abs, 1 - gamma)

percentile_threshold_cpp(X, X_abs, row_pctls, col_pctls)
}

dfleis/morpca documentation built on May 15, 2019, 3:17 p.m.