HSIC: An Interesting Approach to Judging Correlation
A few days ago, I came across a post on Synced (机器之心): "A New Method Solves the Gradient Explosion Problem Once and for All—Training ResNets Without Backpropagation". As usual, let's ignore the sensationalist framing that media outlets tend to favor, and just focus on the content. What this Synced article introduces is the work from the paper The HSIC Bottleneck: Deep Learning without Back-Propagation, which proposes an algorithm for training neural networks via an HSIC Bottleneck.
To be honest, I haven't fully understood this paper myself, since it introduces quite a few new concepts that are unfamiliar to me. However, the concept of "HSIC" in the paper caught my interest. After some study, I've finally gained a basic understanding of what HSIC means and where it comes from, and that's how this post came about—an attempt to give as accessible (though perhaps not fully rigorous) an explanation of HSIC as possible.
Background
HSIC stands for "Hilbert-Schmidt Independence Criterion." Like mutual information, it can be used to measure the independence between two variables.more
On measuring correlation
As we know, the basic form of mutual information is
$$\begin{equation}I(X,Y)=\iint p(x,y)\log \frac{p(x, y)}{p(x)p(y)}dxdy\label{eq:i}\end{equation}$$
If $I(X,Y)=0$, then this tells us that $p(x, y)\equiv p(x)p(y)$, meaning the two variables are independent of each other; otherwise they are correlated. But the term $\log \frac{p(x, y)}{p(x)p(y)}$ implies that we need some way of estimating a probability density.
HSIC serves a purpose similar to mutual information, but unlike mutual information, it doesn't require estimating the probability densities of the two variables—it can be directly reformulated in terms of sampling.
Long-time readers of this blog know that "mutual information" is a concept that comes up frequently here. We've used mutual information for new-word discovery (e.g., New Word Discovery Based on Segmentation), and for unsupervised learning (e.g., Mutual Information in Deep Learning: Unsupervised Feature Extraction)—its importance is evident. So if there's a metric that could replace mutual information, or be even more convenient to use, that's certainly something I need to study.
Problem definition
Generally, we can define the problem as follows:
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 always equal to $p(x), p(y)$, i.e., whether $x,y$ are independent.
Strictly speaking, for continuous variables, "always equal to" here should mean "equal almost everywhere," but we won't dwell on this technicality.
For the sake of clear notation, let $x\in X, y\in Y$ and $f(x),g(y)\in \mathbb{R}$. Note that $x,y$ could be two variables with completely different meanings—for example, $x$ might be "Monday," $y$ might be "going to work," and $p(x,y)$ is the probability of "today is Monday, and today I go to work." Given this, $X,Y$ could be two entirely different domains.
The basic idea is to compute the mutual information $\eqref{eq:i}$, but in many problems we cannot estimate the probability or probability density well. One possible approach is to convert this into a dual problem and learn the mutual information via something like an adversarial approach (the infomax idea), but this method can be unstable and is sensitive to the sampling scheme. The best solution would be to have something like a "correlation coefficient"—an indicator that we can explicitly compute and optimize.
HSIC is aimed precisely at this goal~
HSIC
Here we'll try to introduce the concept of HSIC as clearly as possible. However, "as clearly as possible" doesn't mean "as short as possible"—in fact, the following section will still be fairly long and contain quite a few mathematical formulas. Still, compared to standard tutorials that dive straight into Hilbert spaces, reproducing kernels, various operators, and so on, the introduction here should be friendlier to readers unfamiliar with these concepts.
The basic idea
HSIC is built on the following observation:
$p(x, y)\equiv p(x)p(y)$ if and only if, for any $f,g$, the expression
$$\begin{equation}\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}\end{equation}$$
is equal to 0.
This conclusion is not hard to understand. What's interesting is that the right-hand side of the second equality takes the form of a sampling expression—that is, we've converted this criterion into a sampling-based form, avoiding the need to directly estimate probability densities.
This gives us a way to test independence: pick "enough" $f,g$, then compute
$$\begin{equation}L_H=\sum_{f,g} \big(C[f,g]\big)^2\label{eq:l0}\end{equation}$$
and see how close $L_H$ is to 0. Conversely, if in an optimization problem we want a set of features $x,y$ to be as mutually independent as possible, we can add $L_H$ to the loss function.
Peeling back the layers
Actually, the form of $L_H$ already captures the essence of HSIC's approach to independence testing. Below, we'll continue along this line of thought, peeling back the layers step by step until we arrive at the final form of HSIC.
First, let's work out $\big(C[f,g]\big)^2$:
$$\begin{equation}\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}\end{equation}$$
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 tells us that the value of this expectation has nothing to do with the notation we use for the random variable. So we have
$$\begin{equation}\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}\end{equation}$$
Transforming all the remaining terms in the same way, we ultimately obtain
$$\begin{equation}\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}\end{equation}\label{eq:c}$$
In this way, each term is an expectation of $f(x_1)f(x_2)g(x_1)g(x_2)$, just with the variables sampled from different distributions.
Eigenfunctions
The remaining question is: which $f,g$ should we choose? And what counts as "enough"?
By analogy with vector space theory, all possible $f(x)$ form a vector space $\mathcal{F}$, and likewise all possible $g(y)$ form a vector space $\mathcal{G}$. If we could enumerate all the "basis vectors" of both spaces, that would certainly be enough. So the question becomes: how do we find all the basis vectors?
This is where the "kernel function" comes into play. What exactly is a kernel function—well, that's actually quite complicated, and I don't fully understand it myself. Roughly speaking, a kernel function is analogous to a "positive definite matrix" in linear algebra: it's a symmetric bivariate function $K(x_1, x_2)=K(x_2, x_1)$ defined on $X\times X$. If we think of a univariate function $f(x)$ as analogous to a vector, then
$$\begin{equation}\int K(x_1,x_2) f(x_2)dx_2\end{equation}$$
is akin to a matrix-vector multiplication. Just as with the eigenvalues and eigenvectors of a matrix, a kernel function also has eigenvalues and eigenfunctions: a univariate function $\psi$ satisfying the identity below is called an eigenfunction of this kernel:
$$\begin{equation}\int K(x_1,x_2) \psi(x_2)dx_2=\alpha \psi(x_1)\end{equation}$$
Everything above is just groundwork; the rigorous definition falls under the theory of "reproducing kernel Hilbert spaces." What we'll actually need going forward are two properties:
1. All the eigenfunctions $\psi_1,\psi_2,\dots$ of a kernel function form an orthogonal basis of the space;
2. All the eigenvalues $\alpha_1,\alpha_2,\dots$ of a kernel function are positive, and satisfy
$$\begin{equation}K(x_1,x_2)=\sum_i \alpha_i \psi_i(x_1)\psi_i(x_2)\end{equation}\label{eq:k}$$
Enter HSIC
With this groundwork laid, HSIC is now essentially ready to make its entrance~
First, suppose we already have a kernel function $K_X(x_1,x_2)$ defined on $X\times X$; then we can compute the corresponding eigenvalues $\alpha_1,\alpha_2,\dots$ and eigenfunctions $\psi_1,\psi_2,\dots$ of $K_X(x_1,x_2)$. Likewise, once we have a kernel function $K_Y(y_1,y_2)$ defined on $Y\times Y$, we can compute the corresponding eigenvalues $\beta_1,\beta_2,\dots$ and eigenfunctions $\phi_1,\phi_2,\dots$ of $K_Y(y_1,y_2)$.
Then, since the eigenfunctions form a basis, in $\eqref{eq:l0}$ we can replace $f,g$ with the corresponding eigenfunctions $\psi_i,\phi_j$:
$$\begin{equation}L_H=\sum_{i,j}\big(C[\psi_i, \phi_j]\big)^2\end{equation}$$
Since all the eigenvalues are positive, we can also take a weighted sum using the eigenvalues as weights without changing the role of $L_H$:
$$\begin{equation}L_H=\sum_{i,j}\alpha_i \beta_j\cdot\big(C[\psi_i, \phi_j]\big)^2\end{equation}$$
Now substituting $\eqref{eq:c}$ into the above, we get
$$\begin{equation}\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}\end{equation}$$
Finally, using the identity $\eqref{eq:k}$, the expression inside the brackets is actually just $K_X(x_1,x_2)K_Y(y_1,y_2)$, and so HSIC makes its entrance:
$$\begin{equation}\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}\end{equation}\label{eq:hsic}$$
This is the metric of correlation we've been looking for—it's purely expressed in terms of sampling, and since $K_X,K_Y$ are all given in advance and typically differentiable, this is an indicator that can be explicitly computed via sampling and directly optimized!
In practical computation, there are many candidate kernel functions we could use; a commonly used one is
$$\begin{equation}K(x_1, x_2) = \exp\left(-\frac{\Vert x_1 - x_2\Vert_2^2}{\sigma^2}\right)\end{equation}\label{eq:gk}$$
where $\sigma > 0$ is a constant. The paper mentioned at the beginning of this post, The HSIC Bottleneck: Deep Learning without Back-Propagation, also uses this kernel function. Different kernel functions produce somewhat different results, but all guarantee $HSIC(X,Y)=0 \Leftrightarrow p(x,y)\equiv p(x)p(y)$.
Matrix form
Finally, let's derive the matrix form of $\eqref{eq:hsic}$ under a finite sample.
Following the idea of estimating an expectation via sampling, $\mathbb{E}_{(x_1,y_1)\sim p(x,y)}$ is actually just the average of the results over 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 essentially this averaging operation performed twice. So
$$\begin{equation}\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]\end{equation}$$
Here $K_X(x_i,x_j),K_Y(y_i,y_j)$ are actually symmetric $n\times n$ matrices, which we denote respectively as $K_X,_y$; then the above operation can be written as the matrix product $\frac{1}{n^2}\text{Tr}(K_X K_Y)$, where $\text{Tr}$ denotes the trace of the matrix. Following the same idea, the second term is essentially "the average of all elements of $K_X$ multiplied by the average of all elements of $K_Y$." If we insist on writing this in matrix form, it becomes $\frac{1}{n^4}\text{Tr}(K_X \boldsymbol{1}K_Y \boldsymbol{1})$, where the bold $\boldsymbol{1}$ denotes the all-ones matrix of size $n\times n$. Correspondingly, the last term is "twice the product of the averages of all elements of $K_X K_Y$ and $1/n$," i.e. $\frac{2}{n^3}\text{Tr}(K_X K_Y \boldsymbol{1})$.
So, if we express HSIC in matrix form, we get
$$\begin{equation}\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}\end{equation}$$
Here $J = \boldsymbol{I}-\boldsymbol{1}/n$, and $\boldsymbol{I}$ is the $n$-order identity matrix. As discussed in A Brief Discussion of Biased and Unbiased Estimation, this is actually a biased estimate; replacing the leading $1/n$ with $1/(n-1)$ gives the unbiased estimate:
$$\begin{equation}HSIC(X,Y)=\frac{1}{(n-1)^2}\text{Tr}(K_X J K_Y J)\end{equation}\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 replaced with $1/(n-1)$).
Miscellaneous
Here I'll first give a reference implementation and run a simple experiment to verify HSIC's effectiveness, and then in the next section, we'll reflect on some potential issues with HSIC.
Reference implementation
Given the kernel matrix $K_X,K_Y$, a reference implementation for computing HSIC 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 follows the meaning of each term in $\eqref{eq:hsic}$ directly, rather than the matrix form $\eqref{eq:hsic-m}$. In fact, the matrix form $\eqref{eq:hsic-m}$ is not very efficient (it involves three matrix multiplications).
Let's now run a simple experiment to verify HSIC's effectiveness:
# 产生两组独立无关的随机变量
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) # 计算核矩阵
Ky = np.expand_dims(y, 0) - np.expand_dims(y, 1)
Ky = np.exp(- Ky**2) # 计算核矩阵
print(hsic(Kx, Ky)) # 计算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 implies that $x,y$ have a fairly strong correlation, and HSIC's result reflects this too, coming out to about 0.096—more than two orders of magnitude larger than 0.0002. This confirms that HSIC is indeed effective. (Note that HSIC's output value is generally only meaningful in relative comparisons; its absolute value has no clear interpretation on its own.)
Some personal reflections
Clearly, the HSIC value computed via $\eqref{eq:hsic}$ depends on the choice of kernel function. Regardless of which kernel function is used, in theory we can always guarantee
$$\begin{equation}HSIC(X,Y)=0 \Leftrightarrow p(x,y)\equiv p(x)p(y)\end{equation}$$
But the question is: when $HSIC(X,Y) > 0$, exactly how correlated are $X,Y$?
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 $\eqref{eq:gk}$, we can roughly sense that a kernel function amounts to a similarity measure between two samples. The issue is: what kind of similarity definition truly matches the background of the problem at hand? There's no standard answer to this, and it's usually a difficult question.
For example, suppose $x_1,x_2,x_3$ each represent an image. We know that $\Vert x_1 - x_2\Vert_2 = 0$ means the two images $x_1,x_2$ are exactly identical, but when $\Vert x_1 - x_2\Vert_2,\Vert x_1 - x_3\Vert_2$ are all nonzero, we cannot conclude, simply because $\Vert x_1 - x_2\Vert_2 < \Vert x_1 - x_3\Vert_2$, that image $x_2$ necessarily "looks" more like $x_3$ than image $x_1$ does—because the norm $\Vert\cdot\Vert_2$ is not a perfect measure of visual similarity.
I think this is actually a common shortcoming of all kernel methods: kernel methods can only guarantee that when a certain indicator equals 0, our ideal target has been achieved—but when the indicator is nonzero, it cannot properly measure how far we are from that ideal. A good measure needs to be carefully designed for the specific problem at hand, or automatically learned from the dataset via some GAN-like approach.
That said, this doesn't mean HSIC has no value. HSIC's value lies in its usefulness as an auxiliary objective for optimization—much like when training an image autoencoder, even if we adopt a GAN-based approach, we would still typically use the MSE between the original and reconstructed images as an auxiliary loss.
Closing remarks
All in all, this post has introduced the concept of HSIC in a relatively accessible and straightforward way. The explanation involves some mathematical content, but I've omitted rigorous mathematical definitions and derivations, trying to retain only the core ideas. I believe this approach makes the material easier for readers to digest. For readers who value rigor, please bear with me~
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.