R/nHSIC.R

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

#' Normalized Hilbert-Schmidt Independence Criterion
#'
#' Computes a normalized version of the Hilbert-Schmidt Independence Criterion
#' (HSIC) between two random vectors.
#'
#' The normalization rescales HSIC to reduce the influence of marginal scale
#' effects and yields a quantity analogous in spirit to a correlation measure.
#'
#' The normalized HSIC is defined as
#'
#' \deqn{
#' \mathrm{nHSIC}(X,Y) =
#' \sqrt{
#' \frac{\mathrm{HSIC}(X,Y)}
#' {\sqrt{\mathrm{HSIC}(X,X)\mathrm{HSIC}(Y,Y)}}
#' }.
#' }
#'
#' Here, \eqn{\mathrm{HSIC}(X,Y)} denotes the empirical HSIC value computed
#' using \code{\link[dHSIC]{dhsic}}.
#'
#' Vectors are treated as one-column matrices.
#'
#' @param x A numeric vector or numeric matrix. Rows correspond to observations.
#' @param y A numeric vector or numeric matrix. Rows correspond to observations.
#'
#' @return
#' A numeric scalar giving the normalized HSIC value.
#'
#' @details
#' The function computes empirical HSIC values for \code{(x, y)},
#' \code{(x, x)}, and \code{(y, y)} and combines them according to the
#' normalization formula above.
#'
#' If the normalization denominator is numerically non-positive, the function
#' returns \code{0}.
#'
#' Small negative values caused by floating-point arithmetic are truncated to
#' zero before taking the square root.
#'
#' @references
#' Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005).
#' Measuring statistical dependence with Hilbert-Schmidt norms.
#' \emph{Algorithmic Learning Theory}, 63--77.
#' <doi:10.1007/11564089_7>.
#'
#' @examples
#' set.seed(1)
#'
#' n <- 200
#'
#' x <- matrix(rnorm(n), ncol = 1)
#' y <- x^2 + 0.3 * rnorm(n)
#'
#' # Standard HSIC
#' dHSIC::dhsic(x, y)$dHSIC
#'
#' # Normalized HSIC
#' nHSIC(x, y)
#'
#' # Independent variables
#' z <- matrix(rnorm(n), ncol = 1)
#'
#' dHSIC::dhsic(x, z)$dHSIC
#' nHSIC(x, z)
#'
#' @export
nHSIC <- function(x, y) {
  x <- as.matrix(x)
  y <- as.matrix(y)
  
  hsic_xy <- dHSIC::dhsic(x, y)$dHSIC
  hsic_xx <- dHSIC::dhsic(x, x)$dHSIC
  hsic_yy <- dHSIC::dhsic(y, y)$dHSIC
  
  denom <- sqrt(hsic_xx * hsic_yy)
  
  if (denom <= 0) {
    return(0)
  }
  
  val <- hsic_xy / denom
  
  if (val < 0) {
    val <- 0
  }
  
  sqrt(val)
}