Introduction to HSIC: An Interesting Idea for Determining Correlation

Measuring Correlation

We know that the basic form of mutual information is: I(X,Y)=\iint p(x,y)\log \frac{p(x, y)}{p(x)p(y)}dxdy \label{eq:i} If I(X,Y)=0, it means p(x, y)\equiv p(x)p(y), which implies the two variables are independent; otherwise, they are correlated. However, the term \log \frac{p(x, y)}{p(x)p(y)} means we need to estimate the probability density in some way.

The role of HSIC is similar to mutual information, but unlike mutual information, it does not require estimating the probability density of the two variables. Instead, it transforms the problem into a sampling form.

Long-time readers of this blog know that "mutual information" is a frequent concept here. We can use mutual information for new word discovery (e.g., "New Word Discovery Based on Segmentation") or for unsupervised learning (e.g., "Mutual Information in Deep Learning: Unsupervised Feature Extraction"). The importance of mutual information is evident. If there is an indicator that can replace mutual information and is even more convenient, it is certainly something I must learn.

Problem Definition

Generally, we define the problem as:

Given data (x_1, y_1),(x_2, y_2),\dots,(x_n,y_n)\sim p(x, y), determine whether p(x, y) is identically equal to p(x)p(y), i.e., whether x and y are independent.

Strictly speaking, for continuous variables, "identically equal" here means "equal almost everywhere," but we will not state this strictly here.

For formal description, let x\in X, y\in Y, and f(x),g(y)\in \mathbb{R}. Note that x and y might be two variables with completely different meanings. For example, x might be "Monday" and y might be "going to work"; p(x,y) would be the probability that "today is Monday and I have to go to work." Given this, X and Y might be two completely different domains.

The basic idea is to calculate the mutual information [eq:i], but in many problems, we cannot estimate the probability or probability density well. One possible solution is to transform it into a dual problem and use an adversarial-like approach to learn mutual information (the infomax approach), but this method can be unstable and is affected by the sampling scheme. The best solution would be to have an indicator similar to a "correlation coefficient" that allows us to explicitly calculate and optimize it.

HSIC is aimed at this goal.

Basic Idea

HSIC observes that:

p(x, y)\equiv p(x)p(y) if and only if for any f, g, the expression \begin{aligned} C[f,g]=&\iint p(x,y)f(x)g(y)dxdy - \iint p(x)p(y)f(x)g(y)dxdy\\ =&\mathbb{E}_{(x,y)\sim p(x,y)}[f(x)g(y)]-\mathbb{E}_{x\sim p(x)}[f(x)]\mathbb{E}_{y\sim p(y)}[g(y)] \end{aligned} is equal to 0.

This conclusion is not difficult to understand. Interestingly, the right side of the second equals sign is in the form of an expectation (sampling). This means we transform the indicator into a sampling form, avoiding direct estimation of probability density.

Thus, we have a method to judge independence: select "enough" f, g, and then calculate L_H=\sum_{f,g} \big(C[f,g]\big)^2 \label{eq:l0} to see how close L_H is to 0. Conversely, in an optimization problem where we want features x and y to be as independent as possible, we can add L_H to the loss function.

Peeling the Onion

In fact, the form of L_H already reflects the discriminant idea of HSIC. Let us follow this path and continue to peel the onion, gradually moving towards the final form of HSIC.

First, let’s calculate \big(C[f,g]\big)^2: \begin{aligned} \big(C[f,g]\big)^2=&\big(\mathbb{E}_{(x,y)\sim p(x,y)}[f(x)g(y)]\big)^2 + \big(\mathbb{E}_{x\sim p(x)}[f(x)]\big)^2 \big(\mathbb{E}_{y\sim p(y)}[g(y)]\big)^2\\ & - 2\big(\mathbb{E}_{(x,y)\sim p(x,y)}[f(x)g(y)]\big)\big(\mathbb{E}_{x\sim p(x)}[f(x)]\big)\big(\mathbb{E}_{y\sim p(y)}[g(y)]\big) \end{aligned}

Then we use a trick: we know that \mathbb{E}_{x\sim p(x)}[f(x)]=\mathbb{E}_{x'\sim p(x')}[f(x')], which means the result of the expectation has nothing to do with the notation of the random variable. So we have: \begin{aligned} \big(\mathbb{E}_{x\sim p(x)}[f(x)]\big)^2=&\big(\mathbb{E}_{x_1\sim p(x)}[f(x_1)]\big)\big(\mathbb{E}_{x_2\sim p(x)}[f(x_2)]\big)\\ =&\mathbb{E}_{x_1\sim p(x),x_2\sim p(x)}[f(x_1)f(x_2)] \end{aligned}

By transforming the other terms in the same way, we eventually get: \begin{aligned} \big(C[f,g]\big)^2=&\mathbb{E}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}[f(x_1)f(x_2)g(y_1)g(y_2)] \\ & + \mathbb{E}_{x_1\sim p(x),x_2\sim p(x),y_1\sim p(y),y_2\sim p(y)}[f(x_1)f(x_2)g(y_1)g(y_2)]\\ & - 2 \mathbb{E}_{(x_1,y_1)\sim p(x,y),x_2\sim p(x),y_2\sim p(y)}[f(x_1)f(x_2)g(y_1)g(y_2)] \end{aligned} \label{eq:c} In this way, every term is an expectation of f(x_1)f(x_2)g(y_1)g(y_2), just with different sampling distributions for the variables.

Eigenfunctions

The question now is: Which f, g should we choose? How many are "enough"?

Analogous to the knowledge of vector spaces, all possible f(x) form a vector space \mathcal{F}, and all g(y) form a vector space \mathcal{G}. If we can traverse all the "bases" of these two spaces, that would certainly be enough. So the question is: how to find all the bases?

This is where the "kernel function" comes in. A kernel function is—well, it’s actually quite complicated, and I don’t fully understand it either. Simply put, a kernel function is similar to a "positive definite matrix" in linear algebra. It is a binary symmetric function K(x_1, x_2)=K(x_2, x_1) defined on X\times X. If we analogize a unary function f(x) to a vector, then \int K(x_1,x_2) f(x_2)dx_2 is equivalent to a matrix-vector multiplication. Like the eigenvalues and eigenvectors of a matrix, a kernel function can also define eigenvalues and eigenfunctions. A unary function \psi that satisfies the following identity is called an eigenfunction of this kernel function: \int K(x_1,x_2) \psi(x_2)dx_2=\alpha \psi(x_1)

The above content is preparatory; the strict definition belongs to the category of "Reproducing Kernel Hilbert Spaces." What we will use later are actually two properties:

1. All eigenfunctions \psi_1,\psi_2,\dots of a kernel function constitute an orthogonal basis for that space;

2. All eigenvalues \alpha_1,\alpha_2,\dots of a kernel function are positive and satisfy K(x_1,x_2)=\sum_i \alpha_i \psi_i(x_1)\psi_i(x_2) \label{eq:k}

HSIC Appears

After the above preparation, HSIC is basically ready to appear!

First, suppose we have a kernel function K_X(x_1,x_2) defined on X\times X. We can then calculate the eigenvalues \alpha_1,\alpha_2,\dots and eigenfunctions \psi_1,\psi_2,\dots corresponding to K_X(x_1,x_2). Similarly, with a kernel function K_Y(y_1,y_2) defined on Y\times Y, we can calculate the eigenvalues \beta_1,\beta_2,\dots and eigenfunctions \phi_1,\phi_2,\dots corresponding to K_Y(y_1,y_2).

Then, because the eigenfunctions constitute the basis, in [eq:l0], we can replace f, g with the corresponding eigenfunctions \psi_i, \phi_j: L_H=\sum_{i,j}\big(C[\psi_i, \phi_j]\big)^2 Since all eigenvalues are positive, we can also use the eigenvalues as weights for a weighted sum without changing the purpose of L_H: L_H=\sum_{i,j}\alpha_i \beta_j\cdot\big(C[\psi_i, \phi_j]\big)^2 Now we substitute [eq:c] into the above to get: \begin{aligned} L_H=&\mathbb{E}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}\left[\sum_{i,j}\alpha_i \beta_j\psi_i(x_1)\psi_i(x_2)\phi_j(y_1)\phi_j(y_2)\right] \\ & + \mathbb{E}_{x_1\sim p(x),x_2\sim p(x),y_1\sim p(y),y_2\sim p(y)}\left[\sum_{i,j}\alpha_i \beta_j\psi_i(x_1)\psi_i(x_2)\phi_j(y_1)\phi_j(y_2)\right]\\ & - 2 \mathbb{E}_{(x_1,y_1)\sim p(x,y),x_2\sim p(x),y_2\sim p(y)}\left[\sum_{i,j}\alpha_i \beta_j\psi_i(x_1)\psi_i(x_2)\phi_j(y_1)\phi_j(y_2)\right] \end{aligned}

Finally, using the identity [eq:k], the term inside the square brackets is exactly K_X(x_1,x_2)K_Y(y_1,y_2). Thus, HSIC appears: \begin{aligned} HSIC(X,Y)=&\mathbb{E}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right] \\ & + \mathbb{E}_{x_1\sim p(x),x_2\sim p(x),y_1\sim p(y),y_2\sim p(y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right]\\ & - 2 \mathbb{E}_{(x_1,y_1)\sim p(x,y),x_2\sim p(x),y_2\sim p(y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right] \end{aligned} \label{eq:hsic} This is the indicator we were looking for to measure correlation. It is purely in sampling form, and K_X, K_Y are pre-given and usually differentiable. Therefore, this is an indicator that can be explicitly calculated via sampling and directly optimized!

In practical calculations, there are many kernel functions to choose from. A commonly used one is: K(x_1, x_2) = \exp\left(-\frac{\Vert x_1 - x_2\Vert_2^2}{\sigma^2}\right) \label{eq:gk} where \sigma > 0 is a constant. The paper mentioned at the beginning, "The HSIC Bottleneck: Deep Learning without Back-Propagation", also uses this kernel function. Different kernel functions have slightly different effects, but they all guarantee that HSIC(X,Y)=0 \Leftrightarrow p(x,y)\equiv p(x)p(y).

Matrix Form

Finally, let’s derive the matrix form of [eq:hsic] for a finite sample.

According to the idea of calculating expectations through sampling, \mathbb{E}_{(x_1,y_1)\sim p(x,y)} is actually the average of the results for all sample pairs (x_i,y_i), and \mathbb{E}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)} is just performing this averaging operation twice. So: \mathbb{E}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right]=\frac{1}{n^2}\sum_{i=1}^n \sum_{j=1}^n \left[K_X(x_i,x_j)K_Y(y_i,y_j)\right] Here K_X(x_i,x_j) and K_Y(y_i,y_j) are actually n\times n symmetric matrices, denoted as K_X and K_Y. Then the above operation can be written as the matrix multiplication \frac{1}{n^2}\text{Tr}(K_X K_Y), where \text{Tr} denotes the trace of the matrix. Based on the same idea, the second term is actually "the average of all elements in K_X multiplied by the average of all elements in K_Y." If written in matrix form, it is \frac{1}{n^4}\text{Tr}(K_X \boldsymbol{1}K_Y \boldsymbol{1}), where the bold \boldsymbol{1} denotes an n\times n matrix of all ones. Correspondingly, the last term is "twice the average value of all elements in K_X K_Y divided by n," i.e., \frac{2}{n^3}\text{Tr}(K_X K_Y \boldsymbol{1}).

So, if HSIC is expressed in matrix form, it is: \begin{aligned} HSIC(X,Y)=&\frac{1}{n^2}\text{Tr}(K_X K_Y)+\frac{1}{n^4}\text{Tr}(K_X \boldsymbol{1}K_Y \boldsymbol{1})-\frac{2}{n^3}\text{Tr}(K_X K_Y \boldsymbol{1})\\ =&\frac{1}{n^2}\text{Tr}(K_X J K_Y J) \end{aligned} where J = \boldsymbol{I}-\boldsymbol{1}/n, and \boldsymbol{I} is the n-th order identity matrix. Similar to the discussion in "A Brief Description of Unbiased and Biased Estimation", this is actually a biased estimate. By replacing the leading 1/n with 1/(n-1), we get the unbiased estimate: HSIC(X,Y)=\frac{1}{(n-1)^2}\text{Tr}(K_X J K_Y J) \label{eq:hsic-m} This is the final matrix form of the HSIC formula (note that the 1/n inside J does not need to be changed to 1/(n-1)).

Reference Implementation

If the kernel matrices K_X, K_Y are already known, the implementation of HSIC calculation is as follows:

import numpy as np

def hsic(Kx, Ky):
    Kxy = np.dot(Kx, Ky)
    n = Kxy.shape[0]
    h = np.trace(Kxy) / n**2 + np.mean(Kx) * np.mean(Ky) - 2 * np.mean(Kxy) / n
    return h * n**2 / (n - 1)**2

Note that this implementation is based on the meaning of each term in [eq:hsic], rather than the matrix form [eq:hsic-m]. In fact, the matrix form [eq:hsic-m] is not very efficient (it involves three matrix multiplications).

Below is a simple experiment to verify the effectiveness of HSIC:

# Generate two sets of independent random variables
x = np.random.randn(1000)
y = np.random.randn(1000)

Kx = np.expand_dims(x, 0) - np.expand_dims(x, 1)
Kx = np.exp(- Kx**2) # Calculate kernel matrix

Ky = np.expand_dims(y, 0) - np.expand_dims(y, 1)
Ky = np.exp(- Ky**2) # Calculate kernel matrix

print(hsic(Kx, Ky)) # Calculate HSIC

The output is roughly around 0.0002. If we change x, y to:

x = np.random.randn(1000)
y = x + 0.1 * np.random.randn(1000)

This means x and y have a strong correlation. The HSIC result also reflects this, being approximately 0.096, which is more than two orders of magnitude larger than 0.0002. This indicates that HSIC is indeed effective. (Note: The output value of HSIC generally only has meaning for relative comparison; its absolute value has no obvious meaning.)

Personal Reflections

Clearly, the HSIC calculation result given by [eq:hsic] depends on the choice of the kernel function. Regardless of which kernel function is used, it can theoretically be guaranteed that: HSIC(X,Y)=0 \Leftrightarrow p(x,y)\equiv p(x)p(y) But the question is, when HSIC(X,Y) > 0, how correlated are X and Y exactly?

This depends heavily on the choice of kernel function and the background of the original problem. From the form of the commonly used kernel function [eq:gk], we can roughly perceive that the kernel function is equivalent to the similarity between two samples. The question is what kind of similarity definition truly fits the problem background; there is no standard answer, and it is usually a very difficult problem.

For example, if x_1, x_2, x_3 represent three images, we know that \Vert x_1 - x_2\Vert_2 = 0 means images x_1 and x_2 are identical. However, when \Vert x_1 - x_2\Vert_2 and \Vert x_1 - x_3\Vert_2 are both non-zero, we cannot say that image x_2 "looks" more like x_1 than image x_3 just because \Vert x_1 - x_2\Vert_2 < \Vert x_1 - x_3\Vert_2, because the \Vert\cdot\Vert_2 norm is not a perfect measure of our vision.

In fact, I believe this is a common problem with all kernel methods: kernel methods can only guarantee that when a certain indicator is zero, it meets our ideal pursuit, but when the indicator is non-zero, it cannot well measure the gap between us and the ideal. A good metric must be carefully designed according to the specific problem, or automatically learned from the dataset using methods like GANs.

Of course, this is not to say that HSIC has no value. The value of HSIC lies in the fact that it can be used as an auxiliary objective for optimization. Just as we train an image autoencoder, even if we use the GAN idea, we usually use the MSE between the original image and the reconstructed image as an auxiliary loss.