R/dcorICA.R

In dcorBSS: Distance-Correlation Based Methods for Blind Source Separation and Dependence Analysis

#' Independent Component Analysis via Distance Correlation
#'
#' Estimate an unmixing matrix by sequentially minimizing a distance-correlation
#' criterion on whitened data. The method can be run either with the ordinary
#' distance correlation criterion or with a robustified criterion based on the
#' bowl transformation.
#'
#' The returned object has class \code{"bss"} so that common extractor methods
#' from package \pkg{JADE}, such as [JADE::bss.components()] and [stats::coef()]
#' can be used directly.
#'
#' @param X A numeric data matrix or an object coercible to a matrix. Rows are
#'   observations and columns are variables.
#' @param mu Optional location vector used for whitening and centering. If
#'   \code{NULL}, the column means of \code{X} are used.
#' @param scatter Optional positive definite scatter matrix used for whitening.
#'   If \code{NULL}, the sample covariance matrix is used.
#' @param transform Character string specifying the dependence criterion.
#'   Possible values are \code{"none"} for ordinary distance correlation and
#'   \code{"bowl"} for the bowl-transformed robust criterion.
#' @param alpha Probability level passed to [bowl_transform()] when
#'   \code{transform = "bowl"}.
#' @param sweeps Number of sweeps through the sequential optimization scheme.
#' @param seed Optional integer random seed for the random sweep permutations.
#' @param block_size Either \code{NULL} or one positive integer. If
#'   \code{NULL}, the standard full distance-correlation computation is used in
#'   the objective. If positive, the blockwise computation is used with that
#'   block size.
#' @param optimizer Character string specifying the optimizer to use. Either
#'   \code{"cobyla"} for [nloptr::cobyla()] or \code{"bobyqa"} for 
#'   [minqa::bobyqa()]. Defaults to \code{"cobyla"}.
#' @param control A named list of control parameters passed to
#'   [nloptr::cobyla()] or [minqa::bobyqa()].
#'
#' @return
#' An object of class \code{"bss"}, implemented as a list with components:
#' \describe{
#'   \item{W}{Estimated unmixing matrix in the convention
#'   \eqn{S = (X - \mu) W^T}.}
#'   \item{S}{Estimated independent components.}
#'   \item{mu}{Location vector used for centering.}
#'   \item{transform}{Selected dependence criterion.}
#'   \item{alpha}{Bowl-transform tuning parameter.}
#'   \item{block_size}{Block size used in the distance-correlation criterion,
#'   or \code{NULL} if the standard full computation was used.}
#'   \item{convergence}{Optimization diagnostics from the sequential optimization
#'   steps.}
#' }
#'
#' @details
#' The algorithm first whitens the observations using the supplied location and
#' scatter estimates. It then searches for an orthogonal rotation that minimizes
#' the distance correlation between one component and the remaining components,
#' proceeding sequentially until all components have been extracted.
#'
#' The dependence criterion is evaluated using [dcor_large()]. If
#' \code{block_size = NULL}, the standard full pairwise-distance computation is
#' used. If \code{block_size} is a positive integer, the same criterion is
#' computed blockwise in order to reduce memory usage. The block size does not
#' need to divide the sample size; if it does not, the final block is simply
#' smaller.
#'
#' When the sample mean and sample covariance matrix are used for whitening and
#' the ordinary, non-transformed distance-correlation criterion is minimized,
#' the resulting procedure is closely related in spirit to the distance
#' covariance ICA method of Matteson and Tsay (2017). That paper develops
#' ICA based on distance covariance and distance correlation under classical
#' whitening. When robust location and scatter estimates such as MCD are used
#' for whitening and the bowl-transformed distance-correlation criterion is
#' minimized, the procedure corresponds to the robust approach referred to as
#' PICARD (Leyder et al. 2026).
#'
#' @references
#' Matteson, D. S. and Tsay, R. S. (2017). Independent component analysis via
#' distance covariance. \emph{Journal of the American Statistical Association},
#' \strong{112}(518), 623--637. <doi:10.1080/01621459.2016.1150851>
#'
#' Leyder, S., Raymaekers, J., Rousseeuw, P. J., Van Deuren, T., and Verdonck,
#' T. (2026). Independent component analysis by robust distance correlation.
#' \emph{Advances in Data Analysis and Classification}, \strong{20}, 727--765.
#' <doi:10.1007/s11634-026-00674-9>
#'
#' @examples
#' \donttest{
#' n <- 400
#' S <- cbind(runif(n), rnorm(n), rchisq(n, df = 3))
#'
#' A <- matrix(rnorm(ncol(S)^2), ncol = ncol(S))
#' X <- tcrossprod(S, A)
#'
#' fit_default <- dcorICA(X, seed = 1)
#' fit_block <- dcorICA(X, seed = 1, block_size = 128L)
#' fit_bobyqa <- dcorICA(X, seed = 1, optimizer = "bobyqa")
#'
#' fit_default$W
#' fit_block$W
#' fit_bobyqa$W
#' head(fit_default$S)
#'
#' if (requireNamespace("JADE", quietly = TRUE)) {
#'   coef(fit_default)
#'   plot(fit_default)
#' }
#' 
#' # PICARD: Robust whitening + bowl criterion
#' if (requireNamespace("robustbase", quietly = TRUE)) {
#'   mcd <- robustbase::covMcd(X)
#'   fit_picard <- dcorICA(
#'     X,
#'     mu = mcd$center,
#'     scatter = mcd$cov,
#'     transform = "bowl",
#'     seed = 1
#'   )
#'   if (requireNamespace("JADE", quietly = TRUE)) {
#'     JADE::MD(fit_picard$W, A)
#'   }
#' }
#' }
#'
#' @importFrom stats cov
#' @export
#' @seealso [bowl_transform()], [dcor_large()], [JADE::bss.components()],
#'   [stats::coef()]
dcorICA <- function(X,
                    mu = NULL,
                    scatter = NULL,
                    transform = c("none", "bowl"),
                    alpha = 0.9975,
                    sweeps = ncol(X)+1,
                    seed = NULL,
                    block_size = NULL,
                    optimizer = c("cobyla", "bobyqa"),
                    control = list()) {
  X <- validate_matrix_input(X)
  
  if (is.null(mu)) {
    mu <- colMeans(X)
  }
  if (length(mu) != ncol(X)) {
    stop("`mu` must have length ncol(X).", call. = FALSE)
  }
  
  if (is.null(scatter)) {
    scatter <- stats::cov(X)
  }
  
  transform <- match.arg(transform)
  optimizer <- match.arg(optimizer)
  
  if (!is.numeric(alpha) || length(alpha) != 1L || is.na(alpha) ||
      alpha <= 0 || alpha >= 1) {
    stop("`alpha` must be a single number satisfying 0 < alpha < 1.",
         call. = FALSE)
  }
  
  if (!is.numeric(sweeps) || length(sweeps) != 1L || is.na(sweeps) ||
      sweeps < 1) {
    stop("`sweeps` must be a positive integer.", call. = FALSE)
  }
  sweeps <- as.integer(sweeps)
  
  .check_block_size_or_null_dcov(block_size)
  
  white <- my_whitener(X, center = mu, scatter = scatter)
  K <- white$K
  Z <- white$Z
  
  fit <- optimize_sequential_dcorICA(
    Z = Z,
    transform = transform,
    alpha = alpha,
    sweeps = sweeps,
    seed = seed,
    block_size = block_size,
    optimizer = optimizer,
    control = control
  )
  
  W <- fit$W_rot %*% K
  X_centered <- sweep(X, 2L, mu, FUN = "-")
  S <- tcrossprod(X_centered, W)
  colnames(S) <- paste0("IC.", seq_len(ncol(S)))
  rownames(S) <- rownames(X)
  
  all_steps <- unlist(
    lapply(fit$convergence, function(sw) sw$steps),
    recursive = FALSE
  )
  ierr_vals <- vapply(all_steps, function(x) x$convergence, numeric(1L))
  
  out <- list(
    W = W,
    S = S,
    mu = mu,
    transform = transform,
    alpha = alpha,
    block_size = block_size,
    convergence = list(
      converged = if (optimizer == "cobyla") all(ierr_vals %in% 1:4) else all(ierr_vals == 0),
      codes = ierr_vals,
      details = fit$convergence
    )
  )
  
  class(out) <- c("bss")
  out
}