R/dcov_large.R

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

#' Distance covariance with optional blockwise computation
#'
#' Computes distance covariance between \code{x} and \code{y}. If
#' \code{block_size = NULL}, the standard full pairwise-distance computation is
#' used. If \code{block_size} is a positive integer, a blockwise algorithm is
#' used that avoids allocating the full \eqn{n \times n} distance matrices.
#'
#' \code{dcov_large()} returns the distance covariance and thus behaves like
#' \code{\link[energy]{dcov}}. \code{dcov2_large()} returns the corresponding squared
#' distance covariance.
#'
#' Vectors are treated as one-column matrices.
#'
#' @param x A numeric vector or numeric matrix. Rows are observations.
#' @param y A numeric vector or numeric matrix. Rows are observations.
#' @param block_size Either \code{NULL} or one positive integer. If \code{NULL},
#'   the standard full computation is used, except in the univariate case where
#'   \href{https://CRAN.R-project.org/package=dccpp}{dccpp} is used. If positive, 
#'   the blockwise computation is used with that block size.
#'
#' @return
#' For \code{dcov_large()}, a numeric scalar giving the distance covariance.
#'
#' For \code{dcov2_large()}, a numeric scalar giving the squared distance
#' covariance.
#'
#' @details
#' The functions compute the biased V-statistic version of distance covariance
#' based on pairwise Euclidean distances between the rows of \code{x} and
#' \code{y}.
#'
#' If \code{block_size = NULL}, all pairwise distances are computed in the
#' usual way. In the special case where both \code{x} and \code{y} are
#' univariate, the optimized \href{https://CRAN.R-project.org/package=dccpp}{dccpp} implementation is used instead. This
#' implementation has computational complexity of order \eqn{O(n \log n)},
#' making it preferable to the blockwise approach in the univariate setting.
#'
#' If \code{block_size} is a positive integer, the same quantity is computed
#' block by block. The observations are split into consecutive blocks of size
#' at most \code{block_size}, and pairwise distances are accumulated over all
#' block pairs. This reduces memory usage because only distances for one block
#' pair are stored at a time rather than the full distance matrices.
#'
#' The block size does not need to divide the sample size. If \code{n} is not
#' a multiple of \code{block_size}, the final block simply contains the
#' remaining observations and is processed in the same way as the other blocks.
#'
#' The blockwise result is intended to agree with the standard computation up
#' to numerical differences due to floating-point arithmetic. It mainly reduces memory usage; the computational
#' complexity remains of order \eqn{O(n^2)}.
#'
#' Smaller block sizes reduce peak memory usage but may increase computation
#' time due to additional block overhead. Larger block sizes are closer to the
#' standard full computation in memory behavior.
#'
#' @references
#' Székely, G. J., Rizzo, M. L., and Bakirov, N. K. (2007).
#' Measuring and testing dependence by correlation of distances.
#' \emph{Annals of Statistics}, \strong{35}(6), 2769--2794.
#' <doi:10.1214/009053607000000505>.
#'
#' Székely, G. J. and Rizzo, M. L. (2009).
#' Brownian distance covariance.
#' \emph{Annals of Applied Statistics}, \strong{3}(4), 1236--1265.
#' <doi:10.1214/09-AOAS312>.
#'
#' Chaudhuri, A. and Hu, W. (2019).
#' A fast algorithm for computing distance correlation.
#' \emph{Computational Statistics & Data Analysis}, \strong{135}, 15--24.
#' <doi:10.1016/j.csda.2019.01.016>.
#'
#' @examples
#' n <- 100
#' x <- matrix(rnorm(n * 3), n, 3)
#' y <- matrix(rnorm(n * 2), n, 2)
#'
#' v1 <- dcov2_large(x, y)
#' v2 <- dcov2_large(x, y, block_size = 32L)
#'
#' v1
#' v2
#' all.equal(v1, v2, tolerance = 1e-10)
#'
#' d1 <- dcov_large(x, y)
#' d2 <- dcov_large(x, y, block_size = 32L)
#'
#' d1
#' d2
#' all.equal(d1, d2, tolerance = 1e-10)
#' 
#' # Compare with energy::dcov()
#' if (requireNamespace("energy", quietly = TRUE)) {
#'  
#'  d_energy <- energy::dcov(x, y)
#'  d_large  <- dcov_large(x, y)
#'  
#'  d_energy
#'  d_large
#'  
#'  all.equal(d_energy, d_large, tolerance = 1e-10)
#' }
#'
#' @export
dcov_large <- function(x, y, block_size = NULL) {
  v <- dcov2_large(x, y, block_size = block_size)
  
  if (v < 0) {
    v <- 0
  }
  
  sqrt(v)
}


#' @rdname dcov_large
#' @export
dcov2_large <- function(x, y, block_size = NULL) {
  x <- .as_numeric_matrix_dcov(x, "x")
  y <- .as_numeric_matrix_dcov(y, "y")
  
  if (nrow(x) != nrow(y)) {
    stop("x and y must have the same number of rows", call. = FALSE)
  }
  
  if (is.null(block_size)) {
    
    if (.is_univariate_dcov(x, y)) {
      return(dccpp::dcov(as.numeric(x), as.numeric(y))^2)
    }
    
    return(dcov2_full_cpp(x, y))
  }
  
  .check_block_size_dcov(block_size)
  
  dcov2_block_cpp(x, y, as.integer(block_size))
}


#' Distance correlation with optional blockwise computation
#'
#' Computes the distance correlation between \code{x} and \code{y}. If
#' \code{block_size = NULL}, the standard full pairwise-distance computation is
#' used. If \code{block_size} is a positive integer, a blockwise algorithm is
#' used that avoids allocating the full \eqn{n \times n} distance matrices.
#'
#' If \code{block_size = NULL} and both inputs are univariate, the computation
#' is delegated to the optimized \href{https://CRAN.R-project.org/package=dccpp}{dccpp} implementation. In this special
#' case, the computational complexity is of order \eqn{O(n \log n)}, whereas
#' the standard and blockwise multivariate implementations both have quadratic
#' complexity.
#'
#' Vectors are treated as one-column matrices.
#'
#' @param x A numeric vector or numeric matrix. Rows are observations.
#' @param y A numeric vector or numeric matrix. Rows are observations.
#' @param block_size Either \code{NULL} or one positive integer. If \code{NULL},
#'   the standard full computation is used, except in the univariate case where
#'   \href{https://CRAN.R-project.org/package=dccpp}{dccpp} is used. If positive, the blockwise computation is used with
#'   that block size.
#'
#' @return A numeric scalar, the distance correlation.
#'
#' @details
#' Distance correlation is obtained by combining squared distance covariance
#' terms for \code{(x, y)}, \code{(x, x)}, and \code{(y, y)}.
#'
#' If \code{block_size} is a positive integer, the same quantities are computed
#' block by block. The rows of \code{x} and \code{y} are partitioned into
#' consecutive blocks of size at most \code{block_size}. Distances are then
#' computed only for one block pair at a time, and the contributions needed for
#' distance covariance are accumulated across all block pairs. This can greatly
#' reduce memory requirements for large \code{n}.
#'
#' The block size does not need to divide the sample size. If the sample size
#' is not a multiple of \code{block_size}, the final block is simply smaller.
#'
#' The blockwise computation is intended for larger data sets where forming the
#' full pairwise distance matrices may be memory intensive. It computes the same
#' sample quantity as the standard version, up to possible small differences
#' from floating-point arithmetic. It mainly reduces memory usage; the computational
#' complexity remains of order \eqn{O(n^2)}.
#'
#' Smaller block sizes use less memory but may be somewhat slower. Larger block
#' sizes use more memory and approach the standard full computation.
#'
#' @references
#' Székely, G. J., Rizzo, M. L., and Bakirov, N. K. (2007).
#' Measuring and testing dependence by correlation of distances.
#' \emph{Annals of Statistics}, \strong{35}(6), 2769--2794.
#' <doi:10.1214/009053607000000505>.
#'
#' Székely, G. J. and Rizzo, M. L. (2009).
#' Brownian distance covariance.
#' \emph{Annals of Applied Statistics}, \strong{3}(4), 1236--1265.
#' <doi:10.1214/09-AOAS312>.
#'
#' Chaudhuri, A. and Hu, W. (2019).
#' A fast algorithm for computing distance correlation.
#' \emph{Computational Statistics & Data Analysis}, \strong{135}, 15--24.
#' <doi:10.1016/j.csda.2019.01.016>.
#'
#' @examples
#' n <- 100
#' x <- matrix(rnorm(n * 3), n, 3)
#' y <- matrix(rnorm(n * 2), n, 2)
#'
#' r1 <- dcor_large(x, y)
#' r2 <- dcor_large(x, y, block_size = 32L)
#'
#' r1
#' r2
#' all.equal(r1, r2, tolerance = 1e-10)
#'
#' @export
dcor_large <- function(x, y, block_size = NULL) {
  x <- .as_numeric_matrix_dcov(x, "x")
  y <- .as_numeric_matrix_dcov(y, "y")
  
  if (nrow(x) != nrow(y)) {
    stop("x and y must have the same number of rows", call. = FALSE)
  }
  
  if (is.null(block_size)) {
    
    if (.is_univariate_dcov(x, y)) {
      r <- dccpp::dcor(as.numeric(x), as.numeric(y))
      if (is.nan(r) && all(is.finite(x)) && all(is.finite(y))) {
        r <- 0   # constant x or y, as in the multivariate case
      }
      return(r)
    }
    
    return(dcor_full_cpp(x, y))
  }
  
  .check_block_size_dcov(block_size)
  
  dcor_block_cpp(x, y, as.integer(block_size))
}


.as_numeric_matrix_dcov <- function(z, name) {
  
  if (is.data.frame(z)) {
    z <- as.matrix(z)
  } else if (is.vector(z)) {
    z <- matrix(z, ncol = 1L)
  } else {
    z <- as.matrix(z)
  }
  
  if (!is.numeric(z)) {
    stop(sprintf("%s must be numeric", name), call. = FALSE)
  }

  if (anyNA(z)) {
    stop(sprintf("%s must not contain missing values", name), call. = FALSE)
  }
  
  z
}


.is_univariate_dcov <- function(x, y) {
  ncol(x) == 1L && ncol(y) == 1L
}


.check_block_size_dcov <- function(block_size) {
  
  if (length(block_size) != 1L || is.na(block_size)) {
    stop("block_size must be NULL or a single positive integer",
         call. = FALSE)
  }
  
  if (!is.numeric(block_size) ||
      block_size <= 0 ||
      block_size != as.integer(block_size)) {
    
    stop("block_size must be NULL or a single positive integer",
         call. = FALSE)
  }
  
  invisible(TRUE)
}