blockband() -> sparsify_block_diag(), sparsify_block_banded()

Author: Carol-seven

DESCRIPTION

@@ -30,7 +30,6 @@ Encoding: UTF-8
 LazyData: true
 Imports: 
     igraph,
-    MASS,
     Rcpp,
     Rdpack
 LinkingTo:

NAMESPACE

@@ -1,9 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
-export(blockband)
 export(gen_prec_sbm)
 export(grasps)
-importFrom(MASS,mvrnorm)
+export(sparsify_block_banded)
+export(sparsify_block_diag)
 importFrom(Rcpp,evalCpp)
 importFrom(Rdpack,reprompt)
 importFrom(igraph,as_adjacency_matrix)

R/blockband.R

@@ -0,0 +1,61 @@
+#' Groupwise Block-Banded Sparsifier
+#'
+#' @description
+#' Make a precision-like matrix block-banded according to group membership,
+#' keeping only entries within specified group neighborhoods.
+#'
+#' @param mat A p-by-p precision-like matrix specifying the base matrix to be
+#' masked.
+#'
+#' @param membership An integer vector specifying the group membership.
+#' The length of \code{membership} must be consistent with the dimension p.
+#'
+#' @param neighbor.range An integer (default = 1) specifying the neighbor range,
+#' where groups whose labels differ by at most \code{neighbor.range} are
+#' considered neighbors and kept in the mask.
+#'
+#' @return A list containing:
+#' \describe{
+#' \item{Omega}{The masked precision matrix.}
+#' \item{Sigma}{The covariance matrix, i.e., the inverse of \code{Omega}.}
+#' \item{sparsity}{Proportion of zero entries in \code{Omega}.}
+#' \item{membership}{An integer vector specifying the group membership.}
+#' }
+#'
+#' @examples
+#' ## reproducibility for everything
+#' set.seed(1234)
+#'
+#' ## precision matrix estimation
+#' X <- matrix(rnorm(200), 10, 20)
+#' membership <- c(rep(1,5), rep(2,5), rep(3,4), rep(4,6))
+#' est <- grasps(X, membership = membership, penalty = "lasso", crit = "BIC")
+#'
+#' ## default: keep blocks within ±1 of each group
+#' res1 <- sparsify_block_banded(est$hatOmega, membership, neighbor.range = 1)
+#'
+#' ## wider band: keep blocks within ±2 of each group
+#' res2 <- sparsify_block_banded(est$hatOmega, membership, neighbor.range = 2)
+#'
+#' @export
+
+sparsify_block_banded <- function(mat, membership, neighbor.range = 1) {
+
+  p <- ncol(mat)
+  if (length(membership) != p) {
+    stop(sprintf("Length of 'membership' (%d) must equal the matrix dimension p (%d).",
+                 length(membership), p))
+  }
+
+  ## determine which entries to keep: membership within 'neighbor.range'
+  mask <- abs(outer(membership, membership, `-`)) <= neighbor.range
+
+  ## apply mask to matrix
+  Omega <- mat * mask
+  ## compute covariance
+  Sigma <- solve(Omega)
+
+  return(list(Omega = Omega, Sigma = Sigma,
+              sparsity = sum(Omega == 0) / length(Omega),
+              membership = membership))
+}

R/sparsify_block_banded.R

@@ -0,0 +1,73 @@
+#' Groupwise Block-Diagonal Sparsifier
+#'
+#' @description
+#' Make a precision-like matrix block-diagonal according to group membership,
+#' keeping only within-group blocks. Optionally, further reduce selected
+#' within-group blocks to diagonal-only (identity-like) structure.
+#'
+#' @param mat A p-by-p precision-like matrix specifying the base matrix to be
+#' masked.
+#'
+#' @param membership An integer vector specifying the group membership.
+#' The length of \code{membership} must be consistent with the dimension p.
+#'
+#' @param group.diag (default = \code{NULL})
+#' \enumerate{
+#' \item \code{NULL}: Make \code{mat} block-diagonal, i.e., only keep all
+#' within-group blocks.
+#' \item An integer vector specifying which within-group blocks are further
+#' reduced to diagonal-only structure (based on the block-diagonal form).
+#' }
+#'
+#' @return A list containing:
+#' \describe{
+#' \item{Omega}{The masked precision matrix.}
+#' \item{Sigma}{The covariance matrix, i.e., the inverse of \code{Omega}.}
+#' \item{sparsity}{Proportion of zero entries in \code{Omega}.}
+#' \item{membership}{An integer vector specifying the group membership.}
+#' }
+#'
+#' @examples
+#' ## reproducibility for everything
+#' set.seed(1234)
+#'
+#' ## precision matrix estimation
+#' X <- matrix(rnorm(200), 10, 20)
+#' membership <- c(rep(1,5), rep(2,5), rep(3,4), rep(4,6))
+#' est <- grasps(X, membership = membership, penalty = "lasso", crit = "BIC")
+#'
+#' ## block-diagonalizel; within-group block for group 3 made diagonal-only.
+#' res <- sparsify_block_diag(est$hatOmega, membership, group.diag = 3)
+#'
+#' @export
+
+sparsify_block_diag <- function(mat, membership, group.diag = NULL) {
+
+  p <- ncol(mat)
+  if (length(membership) != p) {
+    stop(sprintf("Length of 'membership' (%d) must equal the matrix dimension p (%d).",
+                 length(membership), p))
+  }
+
+  ## block-diagonal mask: TRUE only when in same group
+  mask <- outer(membership, membership, `==`)
+
+  ## adjust within-group blocks
+  if (!is.null(group.diag)) {
+    block_idx <- split(seq_along(membership), membership)
+    for (g in group.diag) {
+      idx <- block_idx[[as.character(g)]]
+      mask[idx, idx] <- FALSE
+      diag(mask[idx, idx]) <- TRUE
+    }
+  }
+
+  ## apply mask to matrix
+  Omega <- mat * mask
+  ## compute covariance
+  Sigma <- solve(Omega)
+
+  return(list(Omega = Omega, Sigma = Sigma,
+              sparsity = sum(Omega == 0) / length(Omega),
+              membership = membership))
+}