R/sssp.R

In sssp: Semi-Symmetric Sparse Projection

SSSP_CITATION_LINE <- paste0(
  "SSSP algorithm \u2014 J. Shen & S. Damadi, Sparse projection onto ",
  "semi-symmetric sets with applications to sparse optimization, ",
  "J. Global Optimization (2026), doi:10.1007/s10898-026-01592-y"
)

#' Citation for the SSSP algorithm
#'
#' Returns the citation line for the SSSP algorithm and the Journal of Global
#' Optimization paper.
#'
#' @return A character scalar.
#' @export
sssp_citation <- function() {
  SSSP_CITATION_LINE
}

.sssp_is_integerish <- function(x) {
  is.numeric(x) && all(is.finite(x)) && all(x == as.integer(x))
}

.sssp_groups_from_sizes <- function(sizes, n) {
  if (!.sssp_is_integerish(sizes)) {
    stop("group sizes must be integers", call. = FALSE)
  }
  sizes <- as.integer(sizes)
  if (any(sizes <= 0L)) {
    stop("group sizes must be positive integers", call. = FALSE)
  }
  if (sum(sizes) != n) {
    stop(sprintf("group sizes must sum to length(v); got %d and %d", sum(sizes), n), call. = FALSE)
  }
  cuts <- c(0L, cumsum(sizes))
  indices <- vector("list", length(sizes))
  for (j in seq_along(sizes)) {
    indices[[j]] <- seq.int(cuts[[j]] + 1L, cuts[[j + 1L]])
  }
  list(indices = indices, sizes = sizes, form = "sizes")
}

.sssp_groups_from_indices <- function(groups, n) {
  if (length(groups) == 0L) {
    stop("groups must contain at least one group", call. = FALSE)
  }
  indices <- vector("list", length(groups))
  seen <- integer(0)
  for (j in seq_along(groups)) {
    idx <- groups[[j]]
    if (!.sssp_is_integerish(idx)) {
      stop(sprintf("groups[[%d]] must be an integer index vector", j), call. = FALSE)
    }
    idx <- as.integer(idx)
    if (length(idx) == 0L) {
      stop(sprintf("groups[[%d]] must be nonempty", j), call. = FALSE)
    }
    bad <- idx[idx < 1L | idx > n]
    if (length(bad) > 0L) {
      stop(
        sprintf("groups[[%d]] contains index %d; R group indices must be 1-based and within 1..%d", j, bad[[1]], n),
        call. = FALSE
      )
    }
    if (length(intersect(seen, idx)) > 0L) {
      stop("groups must form a disjoint union; an index appears more than once", call. = FALSE)
    }
    seen <- c(seen, idx)
    indices[[j]] <- idx
  }
  if (!identical(sort(seen), seq_len(n))) {
    stop("groups must form a disjoint union of all indices of v", call. = FALSE)
  }
  list(indices = indices, sizes = as.integer(vapply(indices, length, integer(1))), form = "index_lists")
}

.sssp_normalize_groups <- function(groups, n) {
  if (is.list(groups)) {
    .sssp_groups_from_indices(groups, n)
  } else {
    .sssp_groups_from_sizes(groups, n)
  }
}

.sssp_validate_levels <- function(n, p, r, s) {
  if (!.sssp_is_integerish(r) || length(r) != 1L) {
    stop("r must be an integer group sparsity level", call. = FALSE)
  }
  if (!.sssp_is_integerish(s) || length(s) != 1L) {
    stop("s must be an integer sparsity level", call. = FALSE)
  }
  r <- as.integer(r)
  s <- as.integer(s)
  if (r < 1L || r > p) {
    stop(sprintf("r must satisfy 1 <= r <= p; got r=%d, p=%d", r, p), call. = FALSE)
  }
  if (s < 1L || s > n) {
    stop(sprintf("s must satisfy 1 <= s <= n; got s=%d, n=%d", s, n), call. = FALSE)
  }
  list(r = r, s = s)
}

.sssp_build_H <- function(v, group_data, s) {
  p <- length(group_data$indices)
  H_flat <- numeric((s + 1L) * p)
  sorted_indices <- vector("list", p)
  for (j in seq_len(p)) {
    idx <- group_data$indices[[j]]
    ordered <- idx[order(abs(v[idx]), decreasing = TRUE)]
    sorted_indices[[j]] <- ordered
    n_elem <- min(s, length(ordered))
    if (n_elem > 0L) {
      csum <- cumsum(v[ordered[seq_len(n_elem)]]^2)
      for (k in seq_len(n_elem)) {
        H_flat[k * p + j] <- csum[[k]]
      }
    }
  }
  list(H_flat = H_flat, sorted_indices = sorted_indices)
}

.sssp_project_normalized <- function(v, group_data, r, s, verbose = FALSE) {
  start <- proc.time()[["elapsed"]]
  p <- length(group_data$indices)
  H_data <- .sssp_build_H(v, group_data, s)
  core <- .Call(C_sssp_solve_h, H_data$H_flat, group_data$sizes, as.integer(r), as.integer(s), as.integer(p))
  counts <- as.integer(core$counts)
  x <- numeric(length(v))
  selected_indices <- integer(0)
  for (j in seq_along(counts)) {
    q <- counts[[j]]
    if (q > 0L) {
      chosen <- H_data$sorted_indices[[j]][seq_len(q)]
      x[chosen] <- v[chosen]
      selected_indices <- c(selected_indices, chosen)
    }
  }
  wall_time <- proc.time()[["elapsed"]] - start
  info <- list(
    objective = sum((x - v)^2),
    tuple = counts,
    selected_groups = which(counts > 0L),
    selected_indices = sort(selected_indices),
    n_tuples = core$n_tuples,
    wall_time = unname(wall_time),
    selected_energy = core$selected_energy,
    citation = SSSP_CITATION_LINE
  )
  res <- list(x = x, info = info)
  class(res) <- "sssp_projection"
  if (isTRUE(verbose)) {
    print(res)
  }
  res
}

#' Exact joint sparsity and group sparsity projection
#'
#' Computes an exact element of \eqn{P_{G_r \cap C_s}(v)}.
#'
#' The argument `v` is \eqn{v \in \mathbb{R}^n}. The argument `groups` encodes
#' \eqn{\mathcal{L}_1,\ldots,\mathcal{L}_p}; pass either a vector of group
#' sizes \eqn{(|\mathcal{L}_1|,\ldots,|\mathcal{L}_p|)} for contiguous groups
#' or a list of 1-based index vectors that form a disjoint union of all
#' indices of `v`. The argument `r` is the group sparsity level \eqn{r}. The
#' argument `s` is the sparsity level \eqn{s}.
#'
#' @param v Numeric vector \eqn{v \in \mathbb{R}^n}.
#' @param groups Group sizes or a list of 1-based index vectors.
#' @param r Group sparsity level \eqn{r}.
#' @param s Sparsity level \eqn{s}.
#' @param verbose If `TRUE`, print a citation-bearing result.
#' @return A list with `x`, the projection, and `info`, including
#'   `objective`, `tuple`, `selected_groups`, `n_tuples`, and `wall_time`.
#' @export
sssp_project <- function(v, groups, r, s, verbose = FALSE) {
  v <- as.numeric(v)
  if (length(v) == 0L) {
    stop("v must be nonempty", call. = FALSE)
  }
  group_data <- .sssp_normalize_groups(groups, length(v))
  levels <- .sssp_validate_levels(length(v), length(group_data$indices), r, s)
  .sssp_project_normalized(v, group_data, levels$r, levels$s, verbose = verbose)
}

#' Projected gradient descent with the SSSP projection
#'
#' Applies \eqn{x^{k+1} \in P_{G_r \cap C_s}(x^k - \gamma \nabla f(x^k))}
#' to \eqn{f(x)=\frac{1}{2}\|Ax-b\|_2^2}.
#'
#' @param A Numeric matrix.
#' @param b Numeric response vector.
#' @param groups Group sizes or a list of 1-based index vectors.
#' @param r Group sparsity level \eqn{r}.
#' @param s Sparsity level \eqn{s}.
#' @param gamma Step length \eqn{\gamma}; if `NULL`, use \eqn{1/\|A\|_2^2}.
#' @param max_iter Maximum iterations.
#' @param tol Relative step stopping tolerance.
#' @param x0 Optional initial vector.
#' @param verbose If `TRUE`, print a citation-bearing result.
#' @return A list with `x` and `info`.
#' @export
sssp_solve <- function(A, b, groups, r, s, gamma = NULL, max_iter = 200L,
                       tol = 1e-6, x0 = NULL, verbose = FALSE) {
  A <- as.matrix(A)
  b <- as.numeric(b)
  if (length(b) != nrow(A)) {
    stop("length(b) must equal nrow(A)", call. = FALSE)
  }
  n <- ncol(A)
  group_data <- .sssp_normalize_groups(groups, n)
  levels <- .sssp_validate_levels(n, length(group_data$indices), r, s)
  if (is.null(x0)) {
    x <- numeric(n)
  } else {
    x <- as.numeric(x0)
    if (length(x) != n) {
      stop("x0 must have length ncol(A)", call. = FALSE)
    }
  }
  if (is.null(gamma)) {
    d <- svd(A, nu = 0L, nv = 0L)$d
    gamma <- if (length(d) == 0L || d[[1]] == 0) 1 else 1 / (d[[1]]^2)
  }
  if (!is.numeric(gamma) || length(gamma) != 1L || gamma <= 0) {
    stop("gamma must be positive", call. = FALSE)
  }
  start <- proc.time()[["elapsed"]]
  last_projection <- NULL
  iterations <- 0L
  for (k in seq_len(as.integer(max_iter))) {
    residual <- as.vector(A %*% x - b)
    gradient <- as.vector(t(A) %*% residual)
    projected <- .sssp_project_normalized(x - gamma * gradient, group_data, levels$r, levels$s)
    x_next <- projected$x
    last_projection <- projected$info
    step_norm <- sqrt(sum((x_next - x)^2))
    x_norm <- max(1, sqrt(sum(x^2)))
    x <- x_next
    iterations <- k
    if (step_norm <= tol * x_norm) {
      break
    }
  }
  residual <- as.vector(A %*% x - b)
  info <- list(
    objective = 0.5 * sum(residual^2),
    iterations = iterations,
    gamma = gamma,
    wall_time = unname(proc.time()[["elapsed"]] - start),
    projection_info = last_projection,
    citation = SSSP_CITATION_LINE
  )
  res <- list(x = x, info = info)
  class(res) <- "sssp_solve"
  if (isTRUE(verbose)) {
    print(res)
  }
  res
}

#' @export
print.sssp_projection <- function(x, ...) {
  cat(SSSP_CITATION_LINE, "\n", sep = "")
  cat("objective ||x - v||_2^2: ", format(x$info$objective, digits = 12), "\n", sep = "")
  cat("optimal tuple r*: ", paste(x$info$tuple, collapse = ", "), "\n", sep = "")
  cat("selected groups: ", paste(x$info$selected_groups, collapse = ", "), "\n", sep = "")
  cat("tuples enumerated: ", format(x$info$n_tuples, scientific = FALSE), "\n", sep = "")
  cat("wall time: ", format(x$info$wall_time, digits = 6), " s\n", sep = "")
  invisible(x)
}

#' @export
print.sssp_solve <- function(x, ...) {
  cat(SSSP_CITATION_LINE, "\n", sep = "")
  cat("objective 0.5||Ax - b||_2^2: ", format(x$info$objective, digits = 12), "\n", sep = "")
  cat("iterations: ", x$info$iterations, "\n", sep = "")
  cat("gamma: ", format(x$info$gamma, digits = 12), "\n", sep = "")
  cat("wall time: ", format(x$info$wall_time, digits = 6), " s\n", sep = "")
  invisible(x)
}