hsic: Pairwise Hilbert-Schmidt Independence Criterion

Description

Computes pairwise kernel dependence measures for the numeric columns of a matrix or data frame using the Hilbert-Schmidt Independence Criterion (HSIC). By default, hsic() returns a normalised kernel independence correlation in [0, 1], analogous to a distance-correlation matrix.

Usage

hsic(
  data,
  kernel = c("gaussian", "linear", "laplace", "polynomial"),
  bandwidth = c("median", "silverman", "scott"),
  normalise = TRUE,
  estimator = c("biased", "unbiased"),
  na_method = c("error", "pairwise", "complete"),
  p_value = FALSE,
  B = 499L,
  seed = NULL,
  n_threads = getOption("matrixCorr.threads", 1L),
  output = c("matrix", "sparse", "edge_list"),
  threshold = 0,
  diag = TRUE
)
# S3 method for hsic
print(
  x,
  digits = 4,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)
# S3 method for hsic
summary(
  object,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)
# S3 method for summary.hsic
print(
  x,
  digits = NULL,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)
# S3 method for hsic
plot(
  x,
  title = NULL,
  low_color = "white",
  high_color = "steelblue1",
  value_text_size = 4,
  show_value = TRUE,
  ...
)

Value

A matrix-like correlation result. With normalise = TRUE, entries are normalised kernel independence correlations and the diagonal is 1

when the variable has positive kernel self-dependence. With normalise = FALSE, entries are raw HSIC estimates. Attributes record the kernel, bandwidth rule, estimator, raw HSIC matrix, diagnostics, and optional permutation inference.

Arguments

data: A numeric matrix or data frame with at least two numeric columns. Non-numeric columns are dropped.
kernel: Kernel used to build univariate Gram matrices. Supported options are "gaussian" (RBF), "linear", "laplace", and "polynomial".
bandwidth: Bandwidth rule for Gaussian and Laplace kernels. "median" uses the upper median of finite off-diagonal pairwise distances divided by $\sqrt{2}$, with fallback 0.001 for degenerate columns. "silverman" uses Silverman's univariate rule, and "scott" uses Scott's univariate rule. Linear and polynomial kernels ignore this argument.
normalise: Logical. If TRUE, return the normalised kernel independence correlation. If FALSE, return raw HSIC.
estimator: HSIC estimator. "biased" uses the stable V-statistic $n^{-2}\mathrm{tr}(KHLH)$. "unbiased" uses the U-statistic estimator with zero diagonal Gram matrices and requires at least four complete observations per pair.
na_method: Missing-data handling. "error" rejects missing, NaN, and infinite values; "complete" performs one listwise finite-row deletion; "pairwise" recomputes each pair on its pairwise complete rows.
p_value: Logical. If TRUE, attach permutation p-values for the HSIC independence test.
B: Positive integer. Number of permutations used when p_value = TRUE.
seed: Optional positive integer seed for reproducible permutation inference. The seed is passed to the C++ permutation engine and does not mutate the user's global R RNG state.
n_threads: Integer $\geq 1$. Number of OpenMP threads. Defaults to getOption("matrixCorr.threads", 1L).
output: Output representation: "matrix", "sparse", or "edge_list".
threshold: Non-negative absolute-value filter for non-matrix outputs. Must be 0 when output = "matrix".
diag: Logical; whether to include diagonal entries in sparse and edge-list outputs.
x: An object of class hsic.
digits: Integer; number of decimal places to print.
n: Optional row threshold for compact preview output.
topn: Optional number of leading/trailing rows to show when truncated.
max_vars: Optional maximum number of visible columns.
width: Optional display width.
show_ci: One of "yes" or "no"; accepted for compatibility with shared matrix-printing methods.
...: Additional arguments passed to print or plot methods.
object: An object of class hsic.
title: Optional plot title.
low_color: Colour for low HSIC values.
high_color: Colour for high HSIC values.
value_text_size: Font size for tile labels.
show_value: Logical; whether to overlay numeric values on the heatmap.

Author

Thiago de Paula Oliveira

Details

For paired observations $(x_i, y_i)$, HSIC measures the squared Hilbert-Schmidt norm of the RKHS cross-covariance operator. With centred Gram matrices $K_c = HKH$ and $L_c = HLH$, the biased empirical estimator is $$ \widehat{\mathrm{HSIC}}_b = \frac{1}{n^2}\mathrm{tr}(KHLH) = \frac{1}{n^2}\sum_{ij} (K_c)_{ij}(L_c)_{ij}. $$

With normalise = TRUE, the returned matrix contains the kernel independence correlation $$ \widehat{\mathrm{kCor}}(X,Y) = \frac{\widehat{\mathrm{HSIC}}(X,Y)} {\sqrt{\widehat{\mathrm{HSIC}}(X,X) \widehat{\mathrm{HSIC}}(Y,Y)}}. $$ The raw HSIC matrix is stored in attr(x, "hsic_raw"). When estimator = "unbiased", raw HSIC can be slightly negative in finite samples; these signed raw values are retained in hsic_raw, while the displayed normalised coefficient is clipped only for the user-facing matrix.

Kernel formulas are: Gaussian $\exp(-|x_i - x_j|^2 / (2\sigma^2))$, Laplace $\exp(-|x_i - x_j| / \sigma)$, linear $x_i x_j$, and polynomial $(x_i x_j + 1)^2$. For bandwidth = "median", $\sigma$ is the upper median of $\{|x_i - x_j|: i < j\}$ divided by $\sqrt{2}$; if this is non-finite or non-positive, $\sigma = 0.001$.

For p_value = TRUE, the test statistic is the raw HSIC estimate. The null distribution is generated by permuting one variable within each pair, and the p-value is $(1 + \#\{T_b \geq T_{\mathrm{obs}}\})/(B + 1)$. This is an $O(B n^2)$ pairwise procedure and can be expensive for large matrices or large B.

The polynomial kernel currently uses the fixed form $(x_i x_j + 1)^2$. More polynomial controls can be added later without changing the main HSIC object contract.

References

Gretton, A., Bousquet, O., Smola, A., & Schoelkopf, B. (2005). Measuring statistical dependence with Hilbert-Schmidt norms. Algorithmic Learning Theory.

Gretton, A., Fukumizu, K., Teo, C. H., Song, L., Schoelkopf, B., & Smola, A. J. (2007). A kernel statistical test of independence. Advances in Neural Information Processing Systems.

Gretton, A., Bousquet, O., Smola, A., & Schoelkopf, B. (2005). Kernel methods for measuring independence. Journal of Machine Learning Research, 6, 2075-2129.

Pfister, N., Buhlmann, P., Scholkopf, B., & Peters, J. (2018). Kernel-based tests for joint independence. Journal of the Royal Statistical Society: Series B, 80(1), 5-31.

Examples

Run this code

set.seed(1)
x <- rnorm(200)
y <- rnorm(200)
H0 <- hsic(cbind(x = x, y = y))
H0

set.seed(2)
x <- rnorm(200)
y <- x^2 + rnorm(200, sd = 0.1)
H1 <- hsic(cbind(x = x, y = y))
H1

set.seed(3)
X <- cbind(a = rnorm(80), b = rnorm(80), c = rnorm(80))
H_perm <- hsic(X, p_value = TRUE, B = 19, seed = 1)
summary(H_perm)
estimate(H_perm)
tidy(H_perm)

Run the code above in your browser using DataLab