R/dacor_large.R

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

#' Distance autocovariance with optional blockwise computation
#'
#' Computes lag-\code{lag} distance autocovariance for the rows of \code{x}. 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-lag) \times (n-lag)} distance
#' matrices.
#'
#' \code{dacov_large()} returns the distance autocovariance and
#' \code{dacov2_large()} returns the corresponding squared distance
#' autocovariance.
#'
#' Vectors are treated as one-column matrices.
#'
#' @param x A numeric vector or numeric matrix. Rows are observations ordered in
#'   time.
#' @param lag A positive integer smaller than \code{nrow(x)}. The default is
#'   \code{1L}.
#' @param block_size Either \code{NULL} or one positive integer. If
#'   \code{NULL}, the standard full computation is used. If positive, the
#'   blockwise computation is used with that block size.
#'
#' @return
#' For \code{dacov_large()}, a numeric scalar giving the distance
#' autocovariance at lag \code{lag}.
#'
#' For \code{dacov2_large()}, a numeric scalar giving the squared distance
#' autocovariance at lag \code{lag}.
#'
#' @details
#' The functions compare the segment \code{x[1:(n-lag), ]} with the lagged
#' segment \code{x[(1+lag):n, ]} using distance covariance.
#'
#' If \code{block_size = NULL}, all pairwise distances for these two segments
#' are computed in the usual way.
#'
#' If \code{block_size} is a positive integer, the same quantity is computed
#' block by block. The effective sample size is \code{n - lag}, and the
#' corresponding rows are split into consecutive blocks of size at most
#' \code{block_size}. Distances are accumulated over all block pairs, so only
#' one block pair needs to be stored at a time.
#'
#' The block size does not need to divide \code{n - lag}. If \code{n - lag} is
#' not a multiple of \code{block_size}, the final block is simply smaller.
#'
#' The blockwise algorithm mainly reduces memory usage; the computational
#' complexity remains of order \eqn{O((n-lag)^2)}.
#'
#' @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>.
#'
#' Fokianos, K. and Pitsillou, M. (2017).
#' Consistent testing for pairwise dependence in time series.
#' \emph{Technometrics}, \strong{59}(2), 262--270.
#' <doi:10.1080/00401706.2016.1156024>.
#'
#' @examples
#' set.seed(1)
#' n <- 100
#' x <- matrix(rnorm(n * 2), n, 2)
#'
#' v1 <- dacov2_large(x, lag = 1L)
#' v2 <- dacov2_large(x, lag = 1L, block_size = 32L)
#'
#' v1
#' v2
#' all.equal(v1, v2, tolerance = 1e-10)
#'
#' a1 <- dacov_large(x, lag = 2L)
#' a2 <- dacov_large(x, lag = 2L, block_size = 32L)
#'
#' a1
#' a2
#' all.equal(a1, a2, tolerance = 1e-10)
#'
#' @export
dacov_large <- function(x, lag = 1L, block_size = NULL) {
  v <- dacov2_large(x, lag = lag, block_size = block_size)
  if (v < 0) {
    v <- 0
  }
  sqrt(v)
}


#' @rdname dacov_large
#' @export
dacov2_large <- function(x, lag = 1L, block_size = NULL) {
  x <- .as_numeric_matrix_dcov(x, "x")
  .check_lag_dcov(lag, nrow(x))
  
  if (is.null(block_size)) {
    return(dcov2_lag_full_cpp(x, as.integer(lag)))
  }
  
  .check_block_size_dcov(block_size)
  dcov2_lag_block_cpp(x, as.integer(lag), as.integer(block_size))
}


#' Distance autocorrelation with optional blockwise computation
#'
#' Computes lag-\code{lag} distance autocorrelation for the rows of \code{x}. 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-lag) \times (n-lag)} distance
#' matrices.
#'
#' Vectors are treated as one-column matrices.
#'
#' @param x A numeric vector or numeric matrix. Rows are observations ordered in
#'   time.
#' @param lag A positive integer smaller than \code{nrow(x)}. The default is
#'   \code{1L}.
#' @param block_size Either \code{NULL} or one positive integer. If
#'   \code{NULL}, the standard full computation is used. If positive, the
#'   blockwise computation is used with that block size.
#'
#' @return A numeric scalar giving the distance autocorrelation at lag
#'   \code{lag}.
#'
#' @details
#' Distance autocorrelation is computed as the distance correlation between
#' \code{x[1:(n-lag), ]} and \code{x[(1+lag):n, ]}.
#'
#' If \code{block_size = NULL}, the required pairwise distances are computed
#' using the standard full computation.
#'
#' If \code{block_size} is a positive integer, the effective sample of size
#' \code{n - lag} is partitioned into consecutive blocks of size at most
#' \code{block_size}, and distances are accumulated block by block. Only one
#' block pair is stored at a time, which can substantially reduce memory usage.
#'
#' The block size does not need to divide \code{n - lag}. If \code{n - lag} is
#' not a multiple of \code{block_size}, the final block is simply smaller.
#'
#' The blockwise algorithm mainly reduces memory usage; the computational
#' complexity remains of order \eqn{O((n-lag)^2)}.
#'
#' @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>.
#'
#' Fokianos, K. and Pitsillou, M. (2017).
#' Consistent testing for pairwise dependence in time series.
#' \emph{Technometrics}, \strong{59}(2), 262--270.
#' <doi:10.1080/00401706.2016.1156024>.
#'
#' @examples
#' set.seed(1)
#' n <- 100
#' x <- matrix(rnorm(n * 2), n, 2)
#'
#' r1 <- dacor_large(x, lag = 1L)
#' r2 <- dacor_large(x, lag = 1L, block_size = 32L)
#'
#' r1
#' r2
#' all.equal(r1, r2, tolerance = 1e-10)
#'
#' @export
dacor_large <- function(x, lag = 1L, block_size = NULL) {
  x <- .as_numeric_matrix_dcov(x, "x")
  .check_lag_dcov(lag, nrow(x))
  
  if (is.null(block_size)) {
    return(dcor_lag_full_cpp(x, as.integer(lag)))
  }
  
  .check_block_size_dcov(block_size)
  dcor_lag_block_cpp(x, as.integer(lag), as.integer(block_size))
}


.check_lag_dcov <- function(lag, n) {
  if (length(lag) != 1L || is.na(lag)) {
    stop("lag must be a single positive integer", call. = FALSE)
  }
  
  if (!is.numeric(lag) || lag <= 0 || lag != as.integer(lag)) {
    stop("lag must be a single positive integer", call. = FALSE)
  }
  
  if (lag >= n) {
    stop("lag must be smaller than the number of rows of x", call. = FALSE)
  }
  
  invisible(TRUE)
}