# Variance Covariance Regularization Enforces Pairwise Independence in Self-
Supervised Representations

Grégoire Mialon
Meta AI
Paris, France
<EMAIL_ADDRESS>
Randall Balestriero
Meta AI
NYC, USA
<EMAIL_ADDRESS>
Yann LeCun
Meta AI & NYU
NYC, USA
<EMAIL_ADDRESS>

###### Abstract

Self-Supervised Learning (SSL) methods such as VICReg, Barlow Twins or W-MSE
avoid collapse of their joint embedding architectures by constraining or
regularizing the covariance matrix of their projector’s output. This study
highlights important properties of such strategy, which we coin Variance-
Covariance regularization (VCReg). More precisely, we show that VCReg enforces
pairwise independence between the features of the learned representation. This
result emerges by bridging VCReg applied on the projector’s output to kernel
independence criteria applied on the projector’s input. This provides the
first theoretical motivations and explanations of VCReg. We empirically
validate our findings where (i) we put in evidence which projector’s
characteristics favor pairwise independence, (ii) we use these findings to
obtain nontrivial performance gains for VICReg, (iii) we demonstrate that the
scope of VCReg goes beyond SSL by using it to solve Independent Component
Analysis. We hope that our findings will support the adoption of VCReg in SSL
and beyond.

##  1 Introduction

Self-Supervised Learning (SSL) via joint embedding architectures has risen to
learn visual representations outperforming their supervised counterpart. This
paradigm typically enforces similar embeddings for two augmented versions of
the same sample, thus allowing an encoder $f$ to learn a representation for a
given modality without labels. Importantly, $f$ could solve the learning task
by predicting the same embedding for every input, a failure mode known as
collapse. To avoid this, various mechanisms have been proposed hence the
diversity of SSL methods (e.g., Grill et al. (2020); Caron et al. (2020,
2021); Chen and He (2021)). Most of these are a composition $g\circ f$ of the
encoder with a projector neural network $g$. Only $f$ is retained after
training, in opposition to supervised training that never introduces $g$. The
projector was proposed by Chen et al. (2020) and significantly improved the
quality of the learned representation in terms of test accuracy on ImageNet
and other downstream tasks. Although various works (Appalaraju et al., 2020;
Cosentino et al., 2022; Bordes et al., 2022) provide empirical knowledge on
the projector, its exact role remains poorly understood.

Variance-Covariance Regularization (VCReg) was introduced in recent methods
for SSL with joint embeddings (Bardes et al., 2022; Zbontar et al., 2021;
Ermolov et al., 2021) to cope with collapse by constraining or regularizing
the covariance or cross-correlation of the projector $g$ output to be
identity. Bardes et al. (2022); Zbontar et al. (2021) observe that wider
projectors further improve the representation learned by the encoder $f$ and
the former alludes to the possibility that decorrelating a non-linear
projection of the representation improves the quality of the learned
representation. However, neither work further study these intriguing
phenomena.

This study sheds a new light on the role of the projector by demonstrating
that VC regularization of the projector’s output precisely enforces pairwise
independence for the components of the projector’s input i.e. the encoder’s
output, and connects this property to projector’s parameters such as width and
depth. Fully or partially pairwise independent representations are generally
sought for, e.g to disentangle factors of variations (Li et al., 2019; Träuble
et al., 2021). Our experimental analysis suggests that different levels of
pairwise independence of the features in the representation emerge from a
variety of SSL criteria along with mutual independence. However, as opposed to
other frameworks, VCReg allows for theoretical study and explicit control of
the learned independence amount. We prove and experimentally validate this
property for random projectors, study how it emerges in learned projectors,
and use our results to obtain new and significant performance gains on VICReg
over Bardes et al. (2022). We then ablate the SSL context and lean on our
findings to show that VCReg is sufficient to solve Independent Component
Analysis (ICA). Beyond providing a novel theoretical understanding of VCReg,
we believe that this work also leads to a better understanding of the
projector in SSL methods. The scope of VCReg is not limited to SSL: our
experiments on ICA open the way to other applications where some degree of
independence is needed.

##  2 Background

###  2.1 Measuring Statistical Independence Using Kernel Methods

Measuring the independence between two set of realizations
$\\{{\bm{X}}_{1}^{1},\dots,{\bm{X}}_{1}^{N}\\},\\{{\bm{X}}_{2}^{1},\dots,{\bm{X}}_{2}^{N}\\},{\bm{X}}_{1}^{i}\in\mathbb{R}^{M},{\bm{X}}_{2}^{i}\in\mathbb{R}^{M}$
is a fundamental task that has a long history in statistics e.g. through the
Mutual Information (MI) of the two random variables $X_{1}$ and $X_{2}$ from
which those two sets are independently drawn from (Cover, 1999). Those
variables are said independent if the realization of one does not affect the
probability distribution of the other. Computing the MI in practice is known
to be challenging (Goebel et al., 2005), which has led to considerable
interest in using alternative criteria e.g. based on functions in Reproducing
Kernel Hilbert Spaces (RKHS) (Bach and Jordan, 2002), a special case of what
is known as functional covariance or correlation (Rényi, 1959). It consists in
computing those statistics after nonlinear transformation (Leurgans et al.,
1993) as in

$\displaystyle\sup_{f_{1}\in\mathcal{F}_{1},f_{2}\in\mathcal{F}_{2}}\mathrm{Corr}\langle
f_{1}(X_{1}),f_{2}(X_{2})\rangle,$ (1)

where $f_{1},f_{2}$ are constrained to lie within some designed functional
space, and the $\mathrm{Cov}$ can be used instead of the $\mathrm{Corr}$. If
Eq. 1 is small enough, then $X_{1}$ and $X_{2}$ are independent in regard to
the functional spaces $\mathcal{F}_{1},\mathcal{F}_{2}$. For example, if
$\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ are unit balls in their respective
vector spaces, then Eq. 1 is just the norm of the usual correlation/covariance
operator (Mourier, 1953) which would be enough for independence under joint
Gaussian distributions (Melnick and Tenenbein, 1982).

#### HSIC and pairwise independence.

More recently, Gretton et al. (2005a) introduced a pairwise independence
criterion known as the Hilbert-Schmidt Independence Criterion (HSIC) which can
be estimated given empirical samples ${\bm{X}}_{1}$ and
${\bm{X}}_{2}\in\mathbb{R}^{N\times M}$ as

$\mathrm{HSIC}({\bm{X}}_{1},{\bm{X}}_{2}):=\widehat{\mathrm{HSIC}}(X_{1},X_{2})=\frac{1}{(N-1)^{2}}\operatorname{Tr}({\bm{K}}_{1}{\bm{H}}{\bm{K}}_{2}{\bm{H}}),$
(2)

with ${\bm{H}}$ the centering matrix
${\bm{I}}-\mathbf{1}\mathbf{1}^{T}\frac{1}{N}$,
$({\bm{K}}_{1})_{i,j}=k_{1}({\bm{X}}_{1}^{i},{\bm{X}}_{1}^{j})$ and
$({\bm{K}}_{2})_{i,j}=k_{2}({\bm{X}}_{2}^{i},{\bm{X}}_{2}^{j})$ the two kernel
matrices of ${\bm{X}}_{1}$ and ${\bm{X}}_{2}$ respectively, and $k_{1}$,
$k_{2}$ of $\mathcal{F}_{1},\mathcal{F}_{2}$ universal kernels such as the
Gaussian kernel (see Steinwart (2001); Micchelli et al. (2006) for other
examples). Crucially, since

$\mathrm{HSIC}(X_{1},X_{2})\geq\sup_{f_{1}\in\mathcal{F}_{1},f_{2}\in\mathcal{F}_{2}}\mathrm{Cov}\langle
f_{1}(X_{1}),g(X_{2})\rangle,$ (3)

HSIC can be used to test for (pairwise) independence as formalized below.

######  Theorem 1 (Thm. 4 from (Gretton et al., 2005a)).

$\mathrm{HSIC}(X_{1},X_{2})=0$ if and only if $X_{1}$ and $X_{2}$ are
independent.

Gretton et al. (2005a) also provide a statistical test for pairwise
independence based on HSIC. Further quantities such as upper bounds on the MI
can be found in a similar way, e.g. see Thm. 16 in Gretton et al. (2005b). In
our experiments, we will rely on HSIC under the Gaussian kernel scaled by the
median of the distribution of pairwise euclidean distances between samples.

#### dHSIC and mutual independence.

Mutual independence of a set of $D$ random variables $X_{1},\dots,X_{D}$ is a
stronger property than independence between all pairs of random variables in
the set. To evaluate it, Pfister et al. (2018) introduce dHSIC, a multivariate
extension of HSIC. In short, dHSIC measures the distance between mean
embeddings $\mu$ under a RKHS $\mathcal{F}$ (Smola et al., 2007) of the
product of distributions and the joint distribution:

$\mathrm{dHSIC}(X_{1},\dots,X_{D}):=\|\mu(\mathbb{P}^{X_{1}}\otimes\dots\otimes\mathbb{P}^{X_{D}})-\mu(\mathbb{P}^{X_{1},\dots,X_{D}})\|_{\mathcal{F}}.$
(4)

Similarly to HSIC, Pfister et al. (2018) establish the equivalence between
null dHSIC and mutual independence along with a statistical test. We provide
an estimator of dHSIC given empirical samples
${\bm{X}}_{1},\dots,{\bm{X}}_{D}$ in $\mathbb{R}^{N\times M}$ each as well as
implementations of HSIC and dHSIC in Appendix D.

#### Complexities.

In what follows, we consider scalar random variables ($M=1$). The subsequent
complexities of HSIC and dHSIC are $\mathcal{O}(N^{2})$ and
$\mathcal{O}(DN^{2})$ respectively. Testing pairwise independence between all
features in a $D$-set with HSIC is $\mathcal{O}(D^{2}N^{2})$ while dHSIC
requires $N\geq 2D$ (Pfister et al., 2018). Since competitive visual
representations typically have $D=2048$, we resort to surrogates to estimate
independence of the representations in practice. These surrogates are detailed
in Section 5.

###  2.2 Variance-Covariance Regularization in Self-Supervised Learning

SSL with joint embeddings learns visual representations by producing two
different augmented views of an input batch of images
${\bm{S}}\in\mathbb{R}^{N\times W\times H}$, denoted by ${\bm{S}}_{\rm left}$
and ${\bm{S}}_{\rm right}$. Each view is fed to an encoder $f$, typically a
ResNet50 (He et al., 2016), producing representations ${\bm{X}}_{\rm left}$
and ${\bm{X}}_{\rm right}\in\mathbb{R}^{N\times D}$ which are passed through a
projector $g$ to output embeddings ${\bm{Z}}_{\rm left}$ and ${\bm{Z}}_{\rm
right}\in\mathbb{R}^{N\times P}$. Finally, an invariance term encouraging
${\bm{Z}}_{\rm left}$ and ${\bm{Z}}_{\rm right}$ to be similar is applied.
While most SSL methods require architectural or training strategies to avoid
collapse (Grill et al., 2020; He et al., 2021), Bardes et al. (2022) and
Zbontar et al. (2021) only require to modify the loss. After training, only
$f$ is retained to be used in downstream tasks. We will denote the $(2N,P)$
matrix ${\bm{Z}}\triangleq[{\bm{Z}}_{\rm left}^{T},{\bm{Z}}_{\rm
right}^{T}]^{T}$.

#### VICReg.

In Bardes et al. (2022), an anti-collapse term $\mathcal{L}_{\rm VC}$, which
we coin VC regularization (VCReg), is added to an invariance loss to form
$\mathcal{L}_{\rm VIC}$:

$\displaystyle\mathcal{L}_{\rm VC}=$
$\displaystyle\sum_{k=1}^{P}\max\left(0,1-\sqrt{\mathrm{Cov}({\bm{Z}})_{k,k}}\right)+\alpha\sum_{j=1,j\not=k}^{P}\mathrm{Cov}({\bm{Z}})^{2}_{k,j},$
(5) $\displaystyle\mathcal{L}_{\rm VIC}=$
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\|({\bm{Z}}_{\rm
left})_{n,.}-({\bm{Z}}_{\rm right})_{n,.}\|_{2}^{2}+\mathcal{L}_{\rm VC},$ (6)

The leftmost term in $\mathcal{L}_{\rm VC}$ corresponds to regularizing the
variance of each feature in ${\bm{Z}}$ to be unit, while the second term seeks
to minimize the covariance between each pair of features in ${\bm{Z}}$. In
Bardes et al. (2022), the weight of each term in $\mathcal{L}_{\rm VIC}$ can
be tuned. Since the authors find that best results are obtained with equal
weights for the invariance and the variance terms, we only vary $\alpha$.
Although this work focus on VCReg as formulated in VICReg, we provide
background on Barlow Twins and W-MSE in the Appendix A, while next section
establishes the similarity between VICReg, Barlow Twins and W-MSE as VCReg
optimizers.

##  3 VC Regularization of SSL Projector’s Output Enforces Pairwise
Independent Features at the Encoder’s Output

In this section, we demonstrate that for random Multi-Layer Perceptrons (MLPs)
projectors, minimizing VC of the projector’s output amounts to minimizing HSIC
—a measure of pairwise dependence (see Section 2.1)— between all pairs of
variables in the learned representation, i.e the projector’s input. We then
justify how this reasoning extends to learned projectors. Our claims are
experimentally verified in Section 5.

#### Notations.

In our setting, ${\bm{X}}\in\mathbb{R}^{N\times D}$ is the $D$ dimensional
output of the encoder $f$ for a batch of size $N$. ${\bm{X}}$ is then fed into
the projector $g$ to form embeddings
${\bm{Z}}=g({\bm{X}})=[g({\bm{X}})_{1},\dots,g({\bm{X}})_{P}]\in\mathbb{R}^{N\times
P}$ on which SSL criteria are generally applied. The acronym MLP will refer to
projectors typically used in SSL: neural network with three layers of same
width unless stated otherwise.

Our proof strategy consists in proving that VCReg of each MLP layer (Linear +
Batch Normalization (BN) + ReLU) output enforces pairwise independence of its
input, before composing these results. We first study nonlinear elementwise
projectors $g:\mathbb{R}\mapsto\mathbb{R}^{L}$, which belong to the wider
class of DeepSets (Zaheer et al., 2017), and of which BN followed by ReLU can
be seen as an instance. We will denote the mapping of such projectors as
${\bm{Z}}=g({\bm{X}})=[g({\bm{X}}_{:,1}),\dots,g({\bm{X}}_{:,P})]$.

######  Lemma 1 (Nonlinear elementwise projectors minimize HSIC of their
input).

Let $g:\mathbb{R}\mapsto\mathbb{R}^{L}$ be a nonlinear elementwise projector;
then, minimizing the covariance of ${\bm{Z}}$ with respect to the encoder $f$
amounts to minimizing HSIC on all feature pairs in ${\bm{X}}$ with kernels
$K_{i}=g(({\bm{X}})_{:,i})g(({\bm{X}})_{:,i})^{T}$.

###### Proof.

Let us consider the $N$ values of the $i^{\rm th}$ data feature
(${\bm{X}}_{:,i}$) as realizations of a random variable. Recalling Eq. 2,
independence of two random variables can be estimated via HSIC. Considering
the arbitrarily complicated network $g$ and
${\bm{Z}}=[g({\bm{X}}_{:,1}),\dots,g({\bm{X}}_{:,D})]\in\mathbb{R}^{N\times
DL}$, we have:

$\displaystyle\mathrm{HSIC}(({\bm{X}})_{:,i},({\bm{X}})_{:,j})=$
$\displaystyle\frac{1}{(N-1)^{2}}\operatorname{Tr}(g(({\bm{X}})_{:,i})g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,j})g(({\bm{X}})_{:,j})^{T}{\bm{H}})$
$\displaystyle=$
$\displaystyle\frac{1}{(N-1)^{2}}\left\|g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,j})\right\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\left\|\mathrm{Cov}\left(g(({\bm{X}})_{:,i}),g(({\bm{X}})_{:,j}\right)\right\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\left\|\mathrm{Cov}\left({\bm{Z}}\right)_{1+iL:1+(i+1)L,1+jl:1+(j+1)L}\right\|_{F}^{2}$
$\displaystyle\implies\sum_{i\not=j}\mathrm{HSIC}(({\bm{X}})_{:,i},({\bm{X}})_{:,j})=$
$\displaystyle\left\|\mathrm{Cov}\left({\bm{Z}}\right)\odot\left((1-{\bm{I}}_{D})\otimes\mathbf{1}_{L}\mathbf{1}_{L}^{T}\right)\right\|_{F}^{2},$

and, in the case $L=1$ we have $\mathbf{1}_{L}\mathbf{1}_{L}^{T}=1$ leading to
$\sum_{i\not=j}\mathrm{HSIC}(({\bm{X}})_{:,i},({\bm{X}})_{:,j})=\sum_{i\not=j}\mathrm{Cov}\left({\bm{Z}}\right)^{2}_{i,j}$,
concluding the proof. ∎

To rigorously obtain independence, the $K_{i}$’s must be universal kernels. To
satisfy this assumption, we may randomize Batch Normalization by drawing the
mean and variance from some distribution. Combining this operation with a
ReLU, we obtain random elementwise nonlinearities approaching random features
(Rahimi and Recht, 2007) of a universal kernel (Sun et al., 2018). Increasing
$L$ improves the approximation of such kernel (Chen and Phillips, 2017) i.e.
the larger $L$, the better approximation of HSIC will VICReg’s covariance term
be; see Appendix D for the randomized Batch Normalization implementation,
Figure 7 for its justification, and Section 5 for empirical validation.

######  Remark 2 (Necessity of variance regularization).

Although the variance regularization term does not explicitly appear when
minimizing HSIC on all pairs, it is necessary when optimizing ${\bm{X}}$ to
prevent the degenerate solution of ${\bm{X}}$ being a constant, a common
collapse mode of SSL.

######  Lemma 2 (Random linear projectors minimize HSIC of their input).

Let $g$ be a random linear projector and ${\bm{X}}$ has same variance for each
column. Then, minimizing the covariance of ${\bm{Z}}={\bm{X}}{\bm{W}}$ with
respect to the encoder $f$ amounts to minimizing HSIC with a linear kernel for
each pair of variables in ${\bm{X}}$. (Proof in Section B.1.)

Since the corresponding kernel in lemma 2 is linear, only decorrelation can be
achieved for such projector’s input. This however differs from PCA as we
optimize over ${\bm{X}}$ and not $g$’s parameters (${\bm{W}}$). Proving lemma
2 requires a projector with orthogonal weights, i.e.
${\bm{W}}^{T}{\bm{W}}={\bm{I}}$, which gets more and more accurate with random
weights $\mathcal{U}(\nicefrac{{1}}{{\sqrt{D}}},\nicefrac{{1}}{{\sqrt{D}}})$
as $P$ (the output dimension of $g$) increases. This follows since the central
limit theorem states that the dot-product between two $P$-dimensional weight
vectors tends to $0$ with rate $\mathcal{O}(\nicefrac{{1}}{{\sqrt{P}}})$.
Weight initialization in neural nets roughly follows
$\mathcal{U}(\nicefrac{{1}}{{\sqrt{D}}},\nicefrac{{1}}{{\sqrt{D}}})$, which
will be used to implement random projectors.

######  Theorem 3 (MLP projectors with random weights enforce pairwise
independence.).

Let us consider a MLP composed of alternating random linear layers and
elementwise nonlinearities. Then, minimizing the variance and covariance of
the output ${\bm{Z}}$ enforces pairwise independence between all pairs of
variables in the input ${\bm{X}}$.

###### Proof.

Let us consider the last block of such MLP, i.e. a fully-connected linear
layer followed by an elementwise nonlinearity. According to lemma 1, applying
VCReg to the MLP, i.e. at the output of the nonlinearity will enforce pairwise
independence under corresponding kernel for the input of the nonlinearity,
which is also the output of the last linear layer. If the latter is wide
enough so that it can be considered orthogonal, Theorem 11 in Comon (1994)
ensures that pairwise independence is preserved for the input of the layer. We
can then recursively extend the result backward to the whole MLP. If the last
MLP layer is a fully-connected linear, then lemma 2 applies and we go back to
the preceding elementwise nonlinearity. ∎

Fig. 5 in Section C.1 shows that each layer in the random MLP is recursively
VC-regularized from VCReg being applied only at the projector’s output.
Following Theorem 3, we can expect that wider projectors are beneficial to
enforce pairwise independence while adding layers or learning the projector is
not necessary; see Section 5 for empirical validation.

#### Extension to BarlowTwins and W-MSE, and generality of VCReg.

Our results focus on VCReg as formulated in VICReg but can in fact be extended
to methods that constrain the covariance of ${\bm{Z}}$ explicitly, namely
BarlowTwins and W-MSE (Balestriero and LeCun, 2022). Indeed, the objective in
W-MSE (Equations 3-4 in Ermolov et al. (2021)) is VICReg with explicit
constraint on the variance and covariance. Increasing the variance and
covariance hyper-parameters in VICReg produces W-MSE hence our results extend
seamlessly. The objective from Eq. 7 in BarlowTwins is also similar to VICReg.
The derivation, deferred to Section B.2, shows that minimizing the constrained
form of BarlowTwins objective from Eq. 7 is equivalent to minimizing VICReg’s
invariance term whilst explicitly constraining the variance covariance terms
as in W-MSE; hence our results also hold for BarlowTwins; see Section 5 for
empirical validation. In particular we will see that as BarlowTwins explicitly
enforces minimum VCReg, it better optimizes HSIC compared to VICReg with
standard hyper-parameters. Finally, as opposed to BarlowTwins loss and most
SSL methods, VCReg can be used and be beneficial within single branch
architectures. We provide such use case in Section 5.

#### Learned projectors.

In state-of-the-art SSL representations, the projector is learned, which is
not rigorously covered by Theorem 3. Complementing the study of Bordes et al.
(2022), we argue that learning the projector is only crucial to satisfy the
invariance criterion since random projectors are sufficient to obtain pairwise
independent features, as demonstrated in Section 5. In fact, our experiments
show that i) using VICReg, keeping the projector random generally produces
lower HSIC of their input i.e. better enforce pairwise independence than when
learning the projector’s weights, ii) using VICReg, learning the projector to
optimize VCReg more strongly than the invariance term reduces performances,
and iii) using VCReg, learning the projector e.g. in the later studied ICA
setting creates a degenerate representation that does not enforce pairwise
independence. We thus conjecture that learning the projector’s parameters to
mainly minimize VICReg’s invariance term leaves the parameters close enough to
their random initialization (Jacot et al., 2018) to maintain an accurate
estimate of HSIC. In fact, we will see that the wider the projector, the less
far away from initialization the parameters have to move, and the better
optimized will be HSIC.

##  4 Related Work

#### Feature decorrelation.

Also known as whitening, feature decorrelation ensures that correlation
between each pair of different features in a batch of feature vectors is zero,
and that each feature has unit variance. It was originally used as a data pre-
processing technique, see e.g. Hyvärinen and Oja (2000), before being extended
to deep networks (Cogswell et al., 2015) as a regularizer. In the context of
SSL, Hua et al. (2021); Ermolov et al. (2021) find that feature decorrelation
helps solving the collapse issue. The former avoids collapse via appropriate
Batch Normalization and its decorrelating variant (Huang et al., 2018) in the
projector. This variant of Batch Normalization can be seen as the hard
constraint counterpart of BarlowTwins. Practically, whitening can be
implemented by learning a fully-connected layer Husain and Bober (2019)
recovering Principal Component Analysis; as mentioned earlier (Section 3),
this differs from VCReg which keeps the layer’s parameters random and optimize
its input.

#### Independence criterion for learning features.

Enforcing mutual independence to learn a representation has been proposed for
example by Schmidhuber (1992). More recently, and in the context of supervised
learning, Chen et al. (2019) demonstrated improved training of ResNets by
reducing the pairwise dependence in the features at each layer via a
combination of Dropout and Batch Normalization. By opposition, VCReg is
applied to the projector’s output and common SSL projectors do not rely on
Dropout; we also show that Batch Normalization is not necessary to reduce
pairwise dependence.

#### HSIC in supervised and self-supervised learning.

HSIC-based losses have been employed in supervised learning e.g. see Mooij et
al. (2009); Greenfeld and Shalit (2020) to ensure independence between the
residual errors and the labels of a task at hand. Kornblith et al. (2019)
measure similarity between two neural representations via HSIC. In SSL, Tsai
et al. (2021) showed that a modified version of BarlowTwins loss maximizes
HSIC under a linear kernel between the embeddings of two augmented versions of
the same sample. In a similar fashion, Li et al. (2021) proposed a SSL
framework based on maximizing the dependence between the embeddings of
transformations of an image and the image identity via HSIC. Both these works
differ from ours, which connects VCReg to the minimization of HSIC between
pairs of features in the representation.

##  5 Experiments

We first put in evidence pairwise independence properties emerging in learned
visual representations, before validating and exploiting our theoretical
findings. Finally, we ablate the SSL context and demonstrate that the
composition of a SSL-like random projector with VCReg induces enough
independence to perform ICA. Importantly, none of the experiments with VCReg
use Dropout, and Batch Normalization is not used when the projector is random.
Hence, pairwise independence cannot be attributed to those two techniques as
opposed to Chen et al. (2019); Hua et al. (2021).

###  5.1 Pairwise independence emerges in most visual representations

#### Setup.

In these experiments, we track two metrics for statistical independence of the
components in the learned representation during ResNet50 training on ImageNet
with popular SSL frameworks, as well as a supervised baseline. The first
metric, HSIC (Gretton et al., 2005a), tracks pairwise independence. The second
metric, dHSIC (Pfister et al., 2018), tracks mutual independence. As explained
in Section 2.1, neither HSIC nor dHSIC scale to the full representation (2048
components). Therefore, we instead rely on proxys: for pairwise independence,
we compute and average HSIC on all pairs of the first $n$ components of the
representation. For mutual independence, we sample sets of $n$ components for
which we display dHSIC. We set $n=10$ and compute these statistics on the
whole ImageNet validation set. We expect both metrics to decrease during the
training, which would suggest decreasing dependence among the representations.
Although statistical tests are available for both metrics, we track bare HSIC
and dHSIC as they are continuous. However, we perform HSIC tests in this first
serie of experiments so that the bare values can be linked to concrete
portions of independent pairs. The width of the Gaussian kernel in HSIC is
scaled for each representation (see 2.1) thus allowing for comparison between
methods. See Appendix E for the detailed setup.

#### Results.

Figure 1 shows that different learning frameworks implicitly optimize pairwise
independence of the variables in the representation throughout the training
(left), and at some point mutual independence (right). As expected from
Section 3, HSIC is continuously optimized in Barlow Twins and VICReg (it
increases for DINO)111HSIC is continuously optimized for SimCLR as well: this
could be explained in light of Garrido et al. (2022), which connect VICReg and
SimCLR, and is omitted here for conciseness., and Barlow Twins better enforces
pairwise independence since the covariance matrix is constrained. Precise
values along with HSIC independence tests at level $\alpha=0.5$ and
corresponding test accuracies can be found in Table 1. Lower HSIC values do
not necessarily entail more successful tests: it is possible to allocate low
HSIC to sufficiently many pairs while having an overall larger HSIC. For
example, Supervised has less independent pairs in spite of an overall smaller
HSIC. Table 1 also suggests that pairwise independence is not predictive of
test accuracy. Although it seems to be a desirable property, it is not
sufficient to yield a good representation: for example, DINO achieves better
test accuracy than BarlowTwins while having higher HSIC. We study this trade-
off in the next subsection. Although we observe decreasing dHSIC, i.e.
improvement of mutual independence in the representations, all methods quickly
reach a floor and do not improve further although it is lower for Barlow Twins
and VICReg. Hence, we cannot conclude that this property is properly optimized
and do not study it in subsection 5.2, but conduct experiments to better
evaluate mutual independence obtained by VCReg in subsection 5.3.

Figure 1: Overall, HSIC and dHSIC of the encoder output improves for all methods during training. VCReg based methods have smoother and/or further optimization of HSIC and dHSIC. Table 1: Measure of pairwise dependence (HSIC), pairwise independence testing, and test accuracy for popular SSL or supervised representations averaged on multiple subsets of features. Method (100 ep.) |  HSIC $\downarrow$ |  % of indep. pairs $\uparrow$ |  ImageNet Top1 $\uparrow$
---|---|---|---
BarlowTwins (Zbontar et al., 2021) |  0.30 $\pm$ 0.03 |  38 $\pm$ 2 | 67.8
VICReg (Bardes et al., 2022) |  0.39 $\pm$ 0.03 |  34 $\pm$ 2 | 68.4
SimCLR (Chen et al., 2020) |  0.40 $\pm$ 0.03 |  36 $\pm$ 1 | 67.7
DINO (Caron et al., 2021) |  0.38 $\pm$ 0.02 |  35 $\pm$ 3 |  70.4222As opposed to other methods here, DINO relies on multicrop, which partly explains the performance gap.
Supervised (Wightman et al., 2021) |  0.20 $\pm$ 0.05 |  33 $\pm$ 3 | 78.1

###  5.2 Projector characteristics fostering or hurting pairwise independence

#### Setup.

These experiments validate the results from Section 3 by studying which
projector characteristics helps or hurt HSIC optimization when learning visual
representations. The setup is the same as above except that dHSIC is not used
anymore and the projector may be replaced by projectors introduced in Section
3: an elementwise nonlinear projector from lemma 1, random linear projectors
of various width from lemma 2, and a MLP taking the representation ($D=2048$)
as input, and outputting an embedding of size $P$, with two hidden layers of
width $P$ and ReLUs. The MLP can be random as in Theorem 3 or learned as in
Section 5.1. We apply an invariance loss along with VCReg on the output of the
projectors and scale the covariance coefficient with the size of the projector
for fair comparison. See Appendix E for the detailed setup.

#### Results.

Figure 2 (top left) demonstrates that nonlinear elementwise projectors, random
linear projectors and random MLPs with one layer achieve lower HSIC than
learned MLPs, in particular when resampled so as to get better HSIC estimates.
Hence, learning the projector is not necessary to obtain pairwise
independence. As expected in Section 3, for both theoretical and classical
projectors, increasing width yields lower HSIC while increasing depth hurts
HSIC, see Figure 2. Table 4 in Appendix C.1 shows that increasing the weight
of VCReg does not improve HSIC: a possible explanation is that the projector
would otherwise start to optimize too much for VCReg. This intuition will be
backed by Experiments 5.3. Modifying Batch Normalization to get closer to a
universal kernel in lemma 1 yields significant test accuracy gains for VICReg
over Bardes et al. (2022) (Figure 2, top right); the Batch Normalization
employed (see Figure 7) enables explicit control of the statistics standard
deviation and prevents gradient to flow through them. Finally, HSIC is not
correlated to test accuracy across all methods but the correlation is true
within the same method, see Figure 2 (top right).

#### Discussion.

As expected from Section 3, Batch Normalization in the projector is not
required to enforce pairwise independence: Figure 2 (top right) shows that
removing BN improves HSIC while maintaining test accuracy, which could be
explained by the fact that each activation in the projector is implicitly
regularized by VCReg (Figure 5), alleviating the need for BN. This could also
explain why adding layers in detrimental to HSIC: implicit VCReg of the
activations, an assumption of Theorem 3, will be more and more loosely
enforced. Although projectors from Section 3 learn representations with better
pairwise independence properties than classical SSL methods, most of them are
not as competitive in terms of test accuracy (Fig. 2, top right). This
complements the view of Appalaraju et al. (2020); Bordes et al. (2022), which
argue that the projector requires some learning capacity to filter out
information that is irrelevant to the invariance criterion. Our experiments
suggest that projector capacity should rather be increased via width than via
depth: adding layers can be detrimental to HSIC as seen above but also to test
accuracy. Appalaraju et al. (2020); Zbontar et al. (2021) also observe that
the performance saturates for projectors with more than three layers.

Figure 2: Top left: wider projectors and resampling both yield more pairwise
independence for the encoder output. Bottom left: increasing depth hurts HSIC
for random or learned projectors. Top right: trade-off between independence of
the representation and test accuracy (width $D=8192$). Bottom right: ICA
experiments. Resampling and increased width improve the quality of the
reconstruction.

###  5.3 Ablation: Independent Component Analysis (ICA) with VCReg

The following experiments aim to demonstrate that our findings hold outside of
SSL. To this end, we show that VCReg of a SSL-like projector’s output induces
enough pairwise independence to solve ICA. The goal of ICA (Comon, 1994) is to
find a transformation of a random vector ${\bm{Y}}$333For simplicity,
${\bm{Y}}$ directly denotes $N$ empirical $D$-dimensional observations. which
minimizes the statistical dependence of its components. Its simplest instance
is linear ICA, and is motivated by problems such as the cocktail-party, where
${\bm{Y}}\in\mathbb{R}^{N\times D}$ typically results from a linear
transformation of independent sources ${\bm{S}}$ (e.g overlapping voices) one
wants to recover. Formally, ${\bm{Y}}={\bm{S}}{\bm{A}}$, with
${\bm{A}}\in\mathbb{R}^{D\times D}$ an unknown mixing matrix that can
therefore not be inverted. Instead, the linear ICA approach searches for
${\bm{M}}\in\mathbb{R}^{D\times D}$ such that ${\bm{M}}{\bm{Y}}={\bm{S}}$ by
maximizing the statistical mutual independence between the components of
${\bm{M}}{\bm{Y}}$. Being able to recover ${\bm{S}}$ is a property known as
identifiability. In this setting, enforcing pairwise independence in the
components of ${\bm{M}}{\bm{Y}}$ is generally sufficient to recover the true
signal $s$ (Comon (1994), Theorem 11).

#### VC regularized projectors solve linear ICA.

Since pairwise independence is sufficient to solve linear ICA, VCReg of a
random projector’s output should be able to recover ${\bm{S}}$. Our model can
be seen on Figure 3: whitened batches of mixtures ${\bm{Y}}$ are fed to an
encoder $f$. Here, $f$ is a linear transformation ${\bm{M}}$ described above.
The output ${\bm{X}}=f({\bm{Y}})$ is then fed to a projector $g$, and VCReg is
applied to the output ${\bm{Z}}=g({\bm{X}})$ covariance matrix. As in
Experiments 5.2, the projector is randomly resampled at each gradient step to
get better HSIC estimates. This setting corresponds to one branch of a VICReg
network i) with a linear projection $M$ as encoder instead of a neural
network, ii) with a random projector, and iii) without invariance criterion.
We perform ICA on two datasets, a synthetic one (Brakel and Bengio, 2017) with
6 sources among which 2 are noise, and a real audio dataset (Kabal, 2002) with
3 sources among which one is noise also used in Brakel and Bengio (2017). The
evaluation metric is the maximum correlation between the true and
reconstructed sources. As this metric is not available without ground truth,
we select the model with lowest dHSIC for ${\bm{X}}$, which can be perfectly
evaluated here since the number of sources is small.

Figure 3: Linear ICA model expressed from VCReg of the projector’s output.

Our model recovers the sources, and is even competitive with methods
specifically designed for linear ICA such as Fast ICA (Hyvärinen and Oja,
2000) or Anica (Brakel and Bengio, 2017), see Table 2. Both increasing the
width of the projector and resampling it at each step (especially for smaller
projectors) improve the reconstruction as can be seen in Figure 2 (bottom
right). This is in line with our findings in Section 3 and Experiment 5.2.

Table 2: Maximum correlation (between true and reconstructed sources (the higher the better $\uparrow$). The projector performs competitive reconstruction in the linear setting (left), where pairwise independence is sufficient, but not in PNL (right) where it is not. Linear mixture | Post Non Linear mixture
---|---
Method | Synthetic data | Audio data | Method | Synthetic data | Audio data
Whitening | 0.8074 | 0.9876 | Whitening | 0.7981 | 0.9046
FastICA | 0.9998 | 1.0 | FastICA | 0.8311 | 0.8989
Anica |  0.9987$\pm$6.5e-4 |  0.9996$\pm$4.9e-4 | Anica |  0.9794$\pm$53e-4 |  0.9929$\pm$18e-4
VCReg |  0.9986$\pm$8.2e-4 |  0.9936$\pm$64e-4 | VCReg |  0.8465$\pm$142e-4 |  0.8706$\pm$376e-4

#### VC regularized projectors do not solve nonlinear ICA.

Experiment 5.1 suggests that mutual independence also improves during
training, although not as clearly as pairwise independence. Hence, one could
ask whether VCReg also enforces the former enough to solve nonlinear ICA. To
test this hypothesis, we apply VCReg to a particular case of nonlinear ICA
which allows identifiability but does not have equivalence between pairwise
and mutual independence: the post-nonlinear mixture (PNL) (Taleb and Jutten,
1999). In PNL, the sources are linearly mixed before being fed to elementwise
nonlinear functions. For these experiments, and following Brakel and Bengio
(2017), our encoder is a MLP. During our first experiments, we observed an
informational collapse of the encoder, which produces seemingly mutually
independent variables with very poor reconstruction of the sources. To
alleviate this issue, we add a decoder taking ${\bm{X}}$ as input and
reconstructing ${\bm{Y}}$. Figure 8 in Appendix E shows our modified setup. We
compare VCReg to FastICA and Anica. Although FastICA is not meant to solve the
nonlinear case, it remains an interesting baseline. Table 2 shows that our
model fails to recover the sources, as it does only slightly better than
FastICA in the synthetic case. Hence, although mutual independence increases
during training, VCReg does not optimize it enough to solve nonlinear ICA. We
propose a simple explanation for this limitation. Indeed, each feature in
${\bm{Z}}$ is a nonlinear function of all features in ${\bm{X}}$. Hence, it is
not possible to completely decorrelate two components of ${\bm{Z}}$ as they
both contain the same set of features. It is still possible to improve mutual
independence to some extent since, in practice, only parts of the inputs are
considered at once by nonlinear mappings such as neural networks (Erhan et
al., 2009; Adebayo et al., 2018).

#### Learning the projector does not enforce pairwise independence.

ICA experiments done with a learned projector fail to recover independent
sources: this further demonstrates that learning the projector to specifically
optimize the VCReg is counter-productive, and that in the context of SSL,
learning the projector is rather useful for satisfying the invariance
criterion, which is absent in the ICA experiments. VCReg is rather optimized
by the encoder $f$.

##  6 Conclusion

Although pairwise independence in the representation emerges from many SSL
criteria, it can be precisely explained for frameworks using VCReg. This study
also shows that randomizing the nonlinearities yields significant gains on
very competitive methods. Wider and shallower projectors should also be
explored when increasing the size of the projector, a task we leave for future
work since this raises computational issues. Our experiments on ICA not only
confirm our findings but expand the scope of VCReg outside of SSL. Moreover,
VCReg seems to enforce pairwise independence only, and future work could
consist in ensuring better optimization of pairwise and/or mutual
independence.

## References

  * Adebayo et al. (2018) Julius Adebayo, Justin Gilmer, Michael Muelly, Ian Goodfellow, Moritz Hardt, and Been Kim.  Sanity checks for saliency maps.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 2018\.
  * Appalaraju et al. (2020) Srikar Appalaraju, Yi Zhu, Yusheng Xie, and István Fehérvári.  Towards good practices in self-supervised representation learning.  _arXiv preprint arXiv:2012.00868_ , 2020.
  * Bach and Jordan (2002) Francis R Bach and Michael I Jordan.  Kernel independent component analysis.  _Journal of Machine Learning Research (JMLR)_ , 2002.
  * Balestriero and LeCun (2022) Randall Balestriero and Yann LeCun.  Contrastive and non-contrastive self-supervised learning recover global and local spectral embedding methods.  _arxiv preprint arXiv:2205.11508_ , 2022.
  * Bardes et al. (2022) Adrien Bardes, Jean Ponce, and Yann LeCun.  Vicreg: Variance-invariance-covariance regularization for self-supervised learning.  _International Conference on Learning Representations (ICLR)_ , 2022\.
  * Bordes et al. (2022) Florian Bordes, Randall Balestriero, Quentin Garrido, Adrien Bardes, and Pascal Vincent.  Guillotine regularization: Improving deep networks generalization by removing their head.  _arXiv preprint arXiv:2206.13378_ , 2022.
  * Brakel and Bengio (2017) Philemon Brakel and Yoshua Bengio.  Learning independent features with adversarial nets for non-linear ica.  _arXiv preprint arXiv:1710.05050_ , 2017.
  * Caron et al. (2020) Mathilde Caron, Ishan Misra, Julien Mairal, Priya Goyal, Piotr Bojanowski, and Armand Joulin.  Unsupervised learning of visual features by contrasting cluster assignments.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 2020\.
  * Caron et al. (2021) Mathilde Caron, Hugo Touvron, Ishan Misra, Hervé Jégou, Julien Mairal, Piotr Bojanowski, and Armand Joulin.  Emerging properties in self-supervised vision transformers.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 9650–9660, 2021.
  * Chen and Phillips (2017) Di Chen and Jeff M Phillips.  Relative error embeddings of the gaussian kernel distance.  In _International Conference on Algorithmic Learning Theory_ , pages 560–576. PMLR, 2017.
  * Chen et al. (2019) Guangyong Chen, Pengfei Chen, Yujun Shi, Chang-Yu Hsieh, Benben Liao, and Shengyu Zhang.  Rethinking the usage of batch normalization and dropout in the training of deep neural networks.  _arXiv preprint arXiv:1905.05928_ , 2019.
  * Chen et al. (2020) Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In _International Conference on Machine Learning (ICML)_ , 2020.
  * Chen and He (2021) Xinlei Chen and Kaiming He.  Exploring simple siamese representation learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 15750–15758, 2021.
  * Cogswell et al. (2015) Michael Cogswell, Faruk Ahmed, Ross Girshick, Larry Zitnick, and Dhruv Batra.  Reducing overfitting in deep networks by decorrelating representations.  _arXiv preprint arXiv:1511.06068_ , 2015.
  * Comon (1994) Pierre Comon.  Independent component analysis, a new concept?  _Signal processing_ , 36(3):287–314, 1994.
  * Cosentino et al. (2022) Romain Cosentino, Anirvan Sengupta, Salman Avestimehr, Mahdi Soltanolkotabi, Antonio Ortega, Ted Willke, and Mariano Tepper.  Toward a geometrical understanding of self-supervised contrastive learning.  _arXiv preprint arXiv:2205.06926_ , 2022.
  * Cover (1999) Thomas M Cover.  _Elements of information theory_.  John Wiley & Sons, 1999.
  * Erhan et al. (2009) Dumitru Erhan, Yoshua Bengio, Aaron Courville, and Pascal Vincent.  Visualizing higher-layer features of a deep network.  _University of Montreal_ , 1341(3):1, 2009.
  * Ermolov et al. (2021) Aleksandr Ermolov, Aliaksandr Siarohin, Enver Sangineto, and Nicu Sebe.  Whitening for self-supervised representation learning.  In _International Conference on Machine Learning (ICML)_ , 2021.
  * Garrido et al. (2022) Quentin Garrido, Yubei Chen, Adrien Bardes, Laurent Najman, and Yann Lecun.  On the duality between contrastive and non-contrastive self-supervised learning.  _arXiv preprint arXiv:2206.02574_ , 2022.
  * Goebel et al. (2005) Bernhard Goebel, Zaher Dawy, Joachim Hagenauer, and Jakob C Mueller.  An approximation to the distribution of finite sample size mutual information estimates.  In _IEEE International Conference on Communications, 2005. ICC 2005\. 2005_ , volume 2, pages 1102–1106. IEEE, 2005.
  * Greenfeld and Shalit (2020) Daniel Greenfeld and Uri Shalit.  Robust learning with the hilbert-schmidt independence criterion.  In _International Conference on Machine Learning (ICML)_ , 2020.
  * Gretton et al. (2005a) Arthur Gretton, Olivier Bousquet, Alex Smola, and Bernhard Schölkopf.  Measuring statistical dependence with hilbert-schmidt norms.  In _International conference on algorithmic learning theory_ , pages 63–77. Springer, 2005a.
  * Gretton et al. (2005b) Arthur Gretton, Ralf Herbrich, Alexander Smola, Olivier Bousquet, Bernhard Schölkopf, et al.  Kernel methods for measuring independence.  _Journal of Machine Learning Research (JMLR)_ , 2005b.
  * Grill et al. (2020) Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin Tallec, Pierre Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Avila Pires, Zhaohan Guo, Mohammad Gheshlaghi Azar, et al.  Bootstrap your own latent-a new approach to self-supervised learning.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 2020\.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * He et al. (2021) Kaiming He, Xinlei Chen, Saining Xie, Yanghao Li, Piotr Dollár, and Ross Girshick.  Masked autoencoders are scalable vision learners.  _arXiv preprint arXiv:2111.06377_ , 2021.
  * Hua et al. (2021) Tianyu Hua, Wenxiao Wang, Zihui Xue, Sucheng Ren, Yue Wang, and Hang Zhao.  On feature decorrelation in self-supervised learning.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 9598–9608, 2021.
  * Huang et al. (2018) Lei Huang, Dawei Yang, Bo Lang, and Jia Deng.  Decorrelated batch normalization.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 791–800, 2018.
  * Husain and Bober (2019) Syed Sameed Husain and Miroslaw Bober.  Remap: Multi-layer entropy-guided pooling of dense cnn features for image retrieval.  _IEEE Transactions on Image Processing_ , 28(10):5201–5213, 2019.
  * Hyvärinen and Oja (2000) Aapo Hyvärinen and Erkki Oja.  Independent component analysis: algorithms and applications.  _Neural networks_ , 13(4-5):411–430, 2000.
  * Jacot et al. (2018) Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Kabal (2002) Peter Kabal.  Tsp speech database, 2002.  McGill University, Database Version, 1(0):09-02.
  * Kornblith et al. (2019) Simon Kornblith, Mohammad Norouzi, Honglak Lee, and Geoffrey Hinton.  Similarity of neural network representations revisited.  In _International Conference on Machine Learning (ICML)_ , 2019.
  * Leurgans et al. (1993) Sue E Leurgans, Rana A Moyeed, and Bernard W Silverman.  Canonical correlation analysis when the data are curves.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 55(3):725–740, 1993.
  * Li et al. (2021) Yazhe Li, Roman Pogodin, Danica J Sutherland, and Arthur Gretton.  Self-supervised learning with kernel dependence maximization.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 2021\.
  * Li et al. (2019) Zejian Li, Yongchuan Tang, Wei Li, and Yongxing He.  Learning disentangled representation with pairwise independence.  In _Proceedings of the AAAI conference on artificial intelligence_ , 2019.
  * Melnick and Tenenbein (1982) Edward L Melnick and Aaron Tenenbein.  Misspecifications of the normal distribution.  _The American Statistician_ , 36(4):372–373, 1982\.
  * Micchelli et al. (2006) Charles A Micchelli, Yuesheng Xu, and Haizhang Zhang.  Universal kernels.  _Journal of Machine Learning Research (JMLR)_ , 2006.
  * Mooij et al. (2009) Joris Mooij, Dominik Janzing, Jonas Peters, and Bernhard Schölkopf.  Regression by dependence minimization and its application to causal inference in additive noise models.  In _International Conference on Machine Learning (ICML)_ , 2009.
  * Mourier (1953) Edith Mourier.  Eléments aléatoires dans un espace de banach.  In _Annales de l’institut Henri Poincaré_ , volume 13, pages 161–244, 1953.
  * Pfister et al. (2018) Niklas Pfister, Peter Bühlmann, Bernhard Schölkopf, and Jonas Peters.  Kernel-based tests for joint independence.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 80(1):5–31, 2018.
  * Rahimi and Recht (2007) Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 2007\.
  * Rényi (1959) Alfréd Rényi.  On measures of dependence.  _Acta mathematica hungarica_ , 10(3-4):441–451, 1959.
  * Schmidhuber (1992) Jürgen Schmidhuber.  Learning factorial codes by predictability minimization.  _Neural computation_ , 4(6):863–879, 1992.
  * Smola et al. (2007) Alex Smola, Arthur Gretton, Le Song, and Bernhard Schölkopf.  A hilbert space embedding for distributions.  In _International Conference on Algorithmic Learning Theory_ , pages 13–31. Springer, 2007.
  * Steinwart (2001) Ingo Steinwart.  On the influence of the kernel on the consistency of support vector machines.  _Journal of Machine Learning Research (JMLR)_ , 2001.
  * Sun et al. (2018) Yitong Sun, Anna Gilbert, and Ambuj Tewari.  On the approximation properties of random relu features.  _arXiv preprint arXiv:1810.04374_ , 2018.
  * Taleb and Jutten (1999) Anisse Taleb and Christian Jutten.  Source separation in post-nonlinear mixtures.  _IEEE Transactions on signal Processing_ , 47(10):2807–2820, 1999.
  * Träuble et al. (2021) Frederik Träuble, Elliot Creager, Niki Kilbertus, Francesco Locatello, Andrea Dittadi, Anirudh Goyal, Bernhard Schölkopf, and Stefan Bauer.  On disentangled representations learned from correlated data.  In _International Conference on Machine Learning (ICML)_ , 2021.
  * Tsai et al. (2021) Yao-Hung Hubert Tsai, Shaojie Bai, Louis-Philippe Morency, and Ruslan Salakhutdinov.  A note on connecting barlow twins with negative-sample-free contrastive learning, 2021.
  * Wightman et al. (2021) Ross Wightman, Hugo Touvron, and Hervé Jégou.  Resnet strikes back: An improved training procedure in timm.  _arXiv preprint arXiv:2110.00476_ , 2021.
  * Zaheer et al. (2017) Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola.  Deep sets.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 2017\.
  * Zbontar et al. (2021) Jure Zbontar, Li Jing, Ishan Misra, Yann LeCun, and Stéphane Deny.  Barlow twins: Self-supervised learning via redundancy reduction.  _arXiv preprint arXiv:2103.03230_ , 2021.

##  Appendix A Further Background on VCReg methods

#### BarlowTwins.

Zbontar et al. (2021) propose a slightly different approach based on
regularizing ${\bm{C}}$, the $P\times P$ cross-correlation matrix between
${\bm{Z}}_{\rm left}$ and ${\bm{Z}}_{\rm right}$ by optimizing

$\displaystyle\mathcal{L}_{\rm
BT}=\sum_{k=1}^{K}(({\bm{C}})_{k,k}-1)^{2}+\alpha\sum_{k^{\prime}\not=k}({\bm{C}})_{k,k^{\prime}}^{2}\;.$
(7)

Here, the leftmost term corresponds to regularizing the cross-correlation of
the same feature in the two views to be unit, while the rightmost term
regularizes the cross-correlation between pairs of different features in the
two views. Importantly, $({\bm{C}})_{i,j}$ falls back to measuring the cosine
similarity between the $i^{\rm th}$ column of ${\bm{Z}}_{\rm left}$ and the
$j^{\rm th}$ column of ${\bm{Z}}_{\rm right}$ i.e.
$({\bm{C}})_{i,j}=\frac{\langle({\bm{Z}}_{\rm left})_{.,i},({\bm{Z}}_{\rm
right})_{.,j}\rangle}{\|({\bm{Z}}_{\rm left})_{.,i}\|_{2}\|({\bm{Z}}_{\rm
right})_{.,j}\|_{2}}$. Hence, the leftmost term is also an invariance term:
all features must be similar for both views.

#### W-MSE.

Ermolov et al. (2021) use the following loss

$\displaystyle\min\;$ $\displaystyle 2-2\frac{\langle{\bm{Z}}_{\rm
left},{\bm{Z}}_{\rm right}\rangle}{\|{\bm{Z}}_{\rm
left}\|_{2}\dot{\|}{\bm{Z}}_{\rm right}\|_{2}}\;$ $\displaystyle\text{s.t.
}\mathrm{Cov}({\bm{Z}}_{\rm left})=Id,\mathrm{Cov}({\bm{Z}}_{\rm right})=Id.$

where the cosine similarity can be replaced by the Euclidean distance.

##  Appendix B Proofs

###  B.1 Proof of lemma 2: Random Linear Projectors maximize pairwise
independence

###### Proof.

Let’s consider the linear regime for which
$g(({\bm{X}})_{:,i})=({\bm{X}})_{:,i}{\bm{w}}^{T}_{i}$ with orthogonal weights
i.e. $\langle{\bm{w}}_{i},{\bm{w}}_{j}\rangle=1_{\\{i=j\\}}$. This assumption
is realistic since we assume that ${\bm{w}}\in\mathbb{R}^{K}$ with $K$ very
large, and that those are randomly sampled. From that, we then have

$\displaystyle\sum_{i\not=j}\mathrm{HSIC}(({\bm{X}})_{:,i},({\bm{X}})_{:,j})=$
$\displaystyle\frac{1}{(N-1)^{2}}\sum_{i\not=j}\operatorname{Tr}(g(({\bm{X}})_{:,i})g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,j})g(({\bm{X}})_{:,j})^{T}{\bm{H}})$
$\displaystyle=$
$\displaystyle\frac{1}{(N-1)^{2}}\sum_{i\not=j}\|g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,j})\|_{F}^{2},$

we will now push the sum inside the norm by considering the following
equality:

$\displaystyle\|\sum_{i\not=j}f(i,j)\|_{F}^{2}=$
$\displaystyle\sum_{i\not=j}\|f(i,j)\|_{F}^{2}+\sum_{i\not=j}\sum_{k\not=\ell,(i,j)\not=(k,\ell)}\operatorname{Tr}\left(f(i,j)^{T}f(k,\ell)\right)=\sum_{i\not=j}\|f(i,j)\|_{F}^{2},$

since in our case

$\displaystyle f(i,j)$
$\displaystyle=g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,j})={\bm{w}}_{i}({\bm{X}})_{:,i}^{T}{\bm{H}}({\bm{X}})_{:,j}{\bm{w}}_{j}^{T},$
$\displaystyle f(i,j)^{T}f(k,\ell)$
$\displaystyle={\bm{w}}_{j}({\bm{X}})_{:,j}^{T}{\bm{H}}({\bm{X}})_{:,i}{\bm{w}}_{i}^{T}{\bm{w}}_{k}({\bm{X}})_{:,k}^{T}{\bm{H}}({\bm{X}})_{:,\ell}{\bm{w}}_{\ell}^{T},$
$\displaystyle\operatorname{Tr}(f(i,j)^{T}f(k,\ell))$
$\displaystyle=1_{\\{i=k\wedge
j=\ell\\}}({\bm{X}})_{:,j}^{T}{\bm{H}}({\bm{X}})_{:,i}({\bm{X}})_{:,k}^{T}{\bm{H}}({\bm{X}})_{:,\ell},$

leading to

$\displaystyle\sum_{i\not=j}\mathrm{HSIC}(({\bm{X}})_{:,i},({\bm{X}})_{:,j})=$
$\displaystyle\frac{1}{(N-1)^{2}}\|\sum_{i\not=j}g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,j})\|_{F}^{2}$
$\displaystyle=\frac{1}{(N-1)^{2}}\|(\sum_{i}g(({\bm{X}})_{:,i}))^{T}{\bm{H}}(\sum_{j}g(({\bm{X}})_{:,j}))-\sum_{i}g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,i})\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\|\frac{1}{N-1}(\sum_{i}g(({\bm{X}})_{:,i}))^{T}{\bm{H}}(\sum_{j}g(({\bm{X}})_{:,j}))-\operatorname{diag}(\mathrm{Cov}({\bm{X}}{\bm{W}}))\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\|\mathrm{Cov}({\bm{X}}{\bm{W}})-\operatorname{diag}(\mathrm{Cov}({\bm{X}}{\bm{W}}))\|_{F}^{2}$
$\displaystyle=$
$\displaystyle\sum_{i\not=j}\mathrm{Cov}({\bm{X}}{\bm{W}})_{i,j}^{2},$

since we assume that all columns of ${\bm{X}}$ have same variance and that
${\bm{W}}$ is orthogonal we have for the pre-last equality

$\displaystyle\frac{1}{N-1}\sum_{i}g(({\bm{X}})_{:,i})^{T}{\bm{H}}g(({\bm{X}})_{:,i})=\mathrm{Var}(({\bm{X}})_{:,1})\sum_{i}{\bm{w}}_{i}{\bm{w}}_{i}^{T}=\mathrm{Var}(({\bm{X}})_{:,1}),$

and for the last equality we use the fact that since $g$ is linear,
$\sum_{i}g(({\bm{X}})_{:,i})=\sum_{i}({\bm{X}})_{:,i}{\bm{w}}_{i}^{T}={\bm{X}}{\bm{W}}$
with ${\bm{W}}=[{\bm{w}}_{1},\dots,{\bm{w}}_{K}]^{T}$.

∎

###  B.2 On the equivalence between VICReg and BarlowTwins objectives

We can express Barlow Twins objective as:

$\displaystyle\min$ $\displaystyle\sum_{k=1}^{K}(\mathrm{Cov}({\bm{Z}}_{\rm
left},{\bm{Z}}_{\rm
right})_{k,k}-1)^{2}+\alpha\sum_{k^{\prime}\not=k}\mathrm{Cov}({\bm{Z}}_{\rm
left},{\bm{Z}}_{\rm right})_{k,k^{\prime}}^{2}\;$ $\displaystyle\text{s.t.
}\mathrm{Cov}({\bm{Z}}_{\rm left})={\bm{I}},\mathrm{Cov}({\bm{Z}}_{\rm
right})={\bm{I}}.$

Assuming ${\bm{Z}}_{\rm left}^{T}{\bm{Z}}_{\rm left}={\bm{I}},{\bm{Z}}_{\rm
right}^{T}{\bm{Z}}_{\rm right}={\bm{I}}$ i.e. perfect minimization of the
variance and covariance terms, we have

$\displaystyle C_{i,j}=$ $\displaystyle\frac{\langle({\bm{Z}}_{\rm
left})_{:,i},({\bm{Z}}_{\rm left})_{:,j}\rangle}{\|({\bm{Z}}_{\rm
left})_{:,i}\|_{2}\|({\bm{Z}}_{\rm
left})_{:,j}\|_{2}}=\frac{1}{2}\langle({\bm{Z}}_{\rm
left})_{:,i},({\bm{Z}}_{\rm right})_{:,j}\rangle=-\|({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2}-1,$

and thus

$\displaystyle\sum_{i}(C_{i,i}-1)^{2}=\sum_{i}\|({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm right})_{:,i}\|_{2}^{4}=\|{\bm{Z}}_{\rm
left}-{\bm{Z}}_{\rm right}\|_{F}^{4}\propto I({\bm{Z}}_{\rm
left},{\bm{Z}}_{\rm right}),$

so we see that we recover the invariance loss exactly with the diagonal terms
of BarlowTwins. Now the only thing left to show is that this is actually
enough to minimize to also minimize the off-diagonal terms. To see that, let’s
follow the derivations below

$\displaystyle\sum_{i\not=j}C_{i,j}^{2}=$
$\displaystyle\sum_{i\not=j}(1-\frac{1}{2}\|({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2})^{2}$ $\displaystyle=$
$\displaystyle\sum_{i\not=j}(1-\frac{1}{2}\|({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm left})_{:,j}+({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2})^{2}$ $\displaystyle=$
$\displaystyle\sum_{i\not=j}\Big{(}1-\frac{1}{2}\|({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm left})_{:,j}\|_{2}^{2}-\frac{1}{2}\|({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2}$
$\displaystyle-\langle({\bm{Z}}_{\rm left})_{:,i}-({\bm{Z}}_{\rm
left})_{:,j},({\bm{Z}}_{\rm left})_{:,j}-({\bm{Z}}_{\rm
right})_{:,j}\rangle\Big{)}^{2}$ $\displaystyle=$
$\displaystyle\sum_{i\not=j}\Big{(}-\frac{1}{2}\|({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2}-\langle({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm left})_{:,j},({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\rangle\Big{)}^{2}$ $\displaystyle=$
$\displaystyle\sum_{i\not=j}\Big{(}\frac{1}{2}\|({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2}+\langle({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm left})_{:,j},({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\rangle\Big{)}^{2}$
$\displaystyle\leq$
$\displaystyle\sum_{i\not=j}\Big{(}\frac{1}{2}\|({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2}+\|({\bm{Z}}_{\rm
left})_{:,i}-({\bm{Z}}_{\rm left})_{:,j}\|_{2}^{2}\|({\bm{Z}}_{\rm
left})_{:,j}-({\bm{Z}}_{\rm right})_{:,j}\|_{2}^{2}\Big{)}^{2}$
$\displaystyle=$ $\displaystyle\frac{25(K-1)}{4}\|{\bm{Z}}_{\rm
left}-{\bm{Z}}_{\rm right}\|_{2}^{4},$

hence minimizing the BarlowTwins loss with the explicit whitening constraint
is equivalent to minimizing VICReg with explicit whitening constraint.

##  Appendix C Additional experimental results

###  C.1 HSIC of representations learned from ImageNet

Table 3: Measure of pairwise dependence (the lower the better) and test accuracy for the projectors studied in Section 3. $L$ is the number of nonlinear projections in Lemma 1, $D$ is the width of the projector and $l$ is the number of layers. | Resampling | No resampling
---|---|---
Projector type |  Mean HSIC $\downarrow$ |  ImageNet Top1 $\uparrow$ |  Mean HSIC $\downarrow$ |  ImageNet Top1 $\uparrow$
No projector | N.A. | N.A. | 0.48 | 48.6
Elementwise nonlinear ($L=8$) | N.A. | N.A. | 0.34 | 51.0
Random linear ($D=16384$) | 0.16 | 54.3 | 0.13 | 54.2
Random linear ($D=8192$) | 0.14 | 53.6 | 0.30 | 52.9
Random linear ($D=4096$) | 0.16 | 52.3 | 0.24 | 52.9
Random linear ($D=2048$) | 0.18 | 51.7 | 0.38 | 52.3
Random linear ($D=1024$) | 0.27 | 51.1 | 0.39 | 52.0
Random linear ($D=512$) | 0.49 | 50.0 | 0.54 | 51.3
Random ($D=8192$, $l=1$) | 0.06 | 50.9 | 0.48 | 53.6
Random ($D=8192$, $l=2$) | 0.63 | 53.3 | 0.52 | 53.4
Random ($D=8192$, $l=3$) | 0.56 | 52.8 | 0.62 | 54.9
Learned ($D=8192$, $l=1$) | 0.17 | 58.4 | N.A. | N.A.
Learned ($D=8192$, $l=2$) | 0.36 | 65.1 | N.A. | N.A.
Learned ($D=8192$, $l=3$) | 0.30 | 68.4 | N.A. | N.A.
Same without BN | 0.22 | 68.1 | N.A. | N.A.
Same with Random BN | 0.31 | 69.1 | N.A. | N.A.
Figure 4: The projectors proposed in Section 3 better optimize HSIC than their classical counterpart. Table 4: VICReg with 2048-2048-2048 projector and varying covariance coefficient. Increasing this coefficient does not necessarily decreases HSIC. Indeed, the projector may start to optimize too much for the VC criterion thus getting trapped in bad solutions. Covariance coeff |  Mean HSIC $\downarrow$ |  ImageNet Top1 $\uparrow$
---|---|---
2 | 0.33 | 66.6
4 | 0.42 | 66.7
6 | 0.37 | 66.6
8 | 0.29 | 66.3
10 | 0.28 | 66.5
12 | 0.37 | 66.2
Table 5: VICReg with learned MLP projector and varying sizes. The wider the projector, the lower HSIC, which conforms to our theoretical results in Section 3. Covariance coeff |  Mean HSIC $\downarrow$ |  ImageNet Top1 $\uparrow$
---|---|---
512 | 0.42 | 63.5
1024 | 0.37 | 65.4
2048 | 0.33 | 66.6
8192 | 0.31 | 68.4
Figure 5: The covariance matrices for some features of each activation in a
random VICReg projector (same as in Bardes et al. (2022) with width $8192$)
before corresponding hidden layer are close to diagonal, suggesting that each
hidden layer output is implicitly VC-regularized.  Figure 6: The covariance
matrices for some features of each activation in a learned VICReg projector
(same as in Bardes et al. (2022) with width $8192$) before corresponding
hidden layer are close to diagonal, suggesting that each hidden layer output
is implicitly VC-regularized.

##  Appendix D Implementations

###  D.1 HSIC Eq. 2

⬇

1def GaussianKernelMatrix(X, sigma):

2 pairwise_distances = torch.cdist(X, X)

3 return torch.exp( -pairwise_distances / (2 * sigma**2))

4

5def HSIC(X_1, X_2, sigma_1, sigma_2):

6 N = x.size(0) # batch size, should be the same for X_1 and X_2

7 K_1 = GaussianKernelMatrix(X_1, sigma_1) # Gaussian kernel matrix of X_1
with bandwidth sigma_1

8 K_2 = GaussianKernelMatrix(X_2, sigma_2) # Gaussian kernel matrix of X_2
with bandwidth sigma_2

9 H = torch.eye(N) - 1.0 / N # centering matrix

10 HSIC = torch.trace(K_1 @ H @ K_2 @ H)/((N-1)**2)

11 return HSIC

In our study, the bandwidth $\sigma$ of the Gaussian kernel is determined by
the median of the distribution of pairwise euclidean distances between
samples.

###  D.2 dHSIC Eq. 4

dHSIC can be estimated given empirical samples
${\bm{X}}_{1},\dots,{\bm{X}}_{d}$ with respective kernel matrices
${\bm{K}}_{1},\dots,{\bm{K}}_{d}$ as:

$\displaystyle\mathrm{dHSIC}({\bm{X}}_{1},\dots,{\bm{X}}_{d})=$
$\displaystyle\frac{1}{N^{2}}\sum_{i,j}(\bigodot_{k=1}^{d}{\bm{K}}_{k})_{i,j}+\frac{1}{N^{2d}}\prod_{k=1}^{d}\sum_{i,j}({\bm{K}}_{k})_{i,j}-\frac{2}{N^{d+1}}\sum_{i}\bigodot_{k=1}^{d}\sum_{j}({\bm{K}}_{k})_{i,j}.$

When $d=2$, the first term corresponds to a biased $\mathrm{HSIC}$ estimator.
The implementation of dHSIC is given by:

⬇

1import GaussianKernelMatrix # defined above

2

3def dHSIC(X, sigma):

4 length = X.shape(0)

5 term_1 = 1.0

6 term_2 = 1.0

7 term_3 = 2.0 / length

8 for j in range(D):

9 K_j = GaussianKernelMatrix(X[j])

10 term_1 = torch.mul(term_1, K_j)

11 term_2 = 1.0 / length / length / term_2 * torch.sum(K_j)

12 term_3 = 1.0 / length * term_3 * K_j.sum(axis=0)

13

14 term_1 = (1.0 / length) ** 2 * torch.sum(term_1)

15 term_3 = torch.sum(term_3)

16

17 return term_1 + term_2 - term_3 # the three terms of the estimator

As these evaluations already take a few minutes, scaling it to significantly
larger portions of components would be impractical.

###  D.3 Randomized Batch Normalization

⬇

1class RandBatchNorm1d(torch.nn.BatchNorm1d):

2 def forward(input):

3 # the mean and variance term in classical batch norm are

4 # randomized

5 m = torch.randn_like(m) * self.estimator_mean.std +

6 self.estimator_mean.mean

7 std = torch.distributions.Gamma(

8 self.estimator_std.mean.square() / self.estimator_std.var,

9 self.estimator_std.mean / self.estimator_std.var,

10 ).sample()

11 return self.weight * (input - m) / torch.sqrt(std.square_() + self.eps) +
self.bias

###  D.4 Linear ICA model 3

⬇

1for p in projector.parameters(): # freeze the projector

2 p.requires_grad = False

3

4for y in loader:

5 x = encoder(y) # matrix multiplication with M

6

7 z = projector(x) # compute embedding

8

9 C = torch.cov(z) # covariance matrix

10 loss = torch.MSE(C, torch.identity()) # VCReg loss

11 loss.backward()

12 optimizer.step()

13

14 projector.__init__() # to resample the projector

15 for p in projector.parameters(): # freeze the projector

16 p.requires_grad = False

###  D.5 Nonlinear ICA model 8

The implementation of our nonlinear ICA model 8 differs from 3 by the addition
of a reconstruction constraint:

⬇

1for p in projector.parameters(): # freeze the projector

2 p.requires_grad = False

3

4for y in loader:

5 x = encoder(y) # MLP encoder

6

7 z = projector(x) # compute embedding

8

9 y_rec = reconstructor(x) # MLP reconstructor

10

11 C = torch.cov(z) # covariance matrix

12 loss = torch.MSE(C, torch.identity()) + lmda * torch.MSE(y, y_rec) # VCReg
loss with reconstruction

13 loss.backward()

14 optimizer.step()

15

16 projector.__init__() # to resample the projector

17 for p in projector.parameters(): # freeze the projector

18 p.requires_grad = False

##  Appendix E Experimental details

Figure 7: Distributions of $\mu^{\ell,k}$ (top row) and $\sigma^{\ell,k}$
(bottom row) which are the BN centering and scaling statistics computed at
each mini-batch for three units in a DN layer with
${\bm{W}}^{\ell}{\bm{z}}^{\ell-1}\sim\mathcal{N}([1,0,-1],\operatorname{diag}([1,3,0.1]))$
($\ell$ is not relevant in this context). The empirical distributions in black
are obtained by repeatedly sampling mini-batches of size $64$ from a training
set of size $1000$. The analytical distributions in green can be obtained
analysis, and its empirical distributions in red closely match the training
one. However by doing the amount of noise can be controlled and is now
independent of the mini-batch size.

###  E.1 Details on ImageNet experiments

#### Architectures and hyper-parameters.

We use a Resnet50 bottleneck, optimized during $100$ epochs with LARS and an
initial learning rate of $0.3$ for all methods. The batch size is $1024$ for
all methods except SimCLR. All methods have projector $8192-8192-8192$ except
DINO and Supervised. Throughout the training, we evaluate the learned
representation using an online linear classifier. Details specific to each
method:

  * •

Barlow Twins: for Experiment 5.1, we have off-diagonal coefficient $0.0051$.

  * •

VICReg: for Experiment 5.1, we have projector size $8192-8192-8192$,
Invariance coeff $25$, Variance coeff $25$ and Covariance coeff $1$. For the
rest of the experiments, and following Bardes et al. (2022), we only move the
covariance coefficient when changing the size of the projector. We scale it in
the square root of the output size of the projector in order to have similar
magnitude for the covariance term for all projector sizes.

  * •

SimCLR: for Experiment 5.1, we choose a higher batch size of $2048$ as SimCLR
is sensitive to this parameter (Chen et al., 2020). The temperature is $0.15$.

  * •

DINO: we use $8$ crops. The head has $4$ layers of widths
$2048-2048-2048-256-65536$.

#### Augmentations.

At train time, the resolution is $160$ and we use

  * •

RandomHorizontalFlip().

  * •

ColorJitter(0.8, 0.4, 0.4, 0.2, 0.1).

  * •

Greyscale(0.2).

  * •

NormalizeImage with ImageNet mean and standard deviation.

  * •

GaussianBlur() with kernel size $(5,9)$ and sigma $(0.1,2)$.

At validation time, the resolution is $224$ and the images are cropped and
normalized with ImageNet mean and standard deviation.

###  E.2 ICA setup

#### Detailed setup.

In these experiments, we keep the same optimizer, train on $100$ epochs and
choose a batch size of $64$ for both datasets. We tune the learning rate
according to a logarithmic grid $[1.0,10.0,100.0]$ (then rescaled according to
$\frac{lr\times batchsize}{256}$). The MLP can be either a SSL-like projector
or simply a fully-connected layer followed by a ReLU. The wider the MLP, the
better the result hence it does not require selection. For the synthetic
dataset, the MLP has width $1024$, and $8192$ for the audio one.

  * •

Linear ICA: We tune the standard and covariance coefficients according to a
logarithmic grid $[1,10,100]$.

  * •

PNL ICA: Our architecture for the nonlinear ICA experiments is presented in
Figure 8. The encoder $f$ is a MLP with $3$ layers and width $128$. The
decoder $h$ is an learnable nonlinearity (a MLP with $3$ hidden layers and
depth $16$) followed by a fully-connected layer. We tune the standard,
covariance and reconstruction coefficients according to a larger logarithmic
grid $[1,3,10,30,100]$.

Figure 8: Nonlinear ICA model.

#### Reconstructed sources for ICA.

In this paragraph, we provide ground truth sources for the synthetic data
along with examples of reconstructed sources.

Figure 9: Data before mixing. Figure 10: Data after mixing, before feeding to
the various ICA models. Figure 11: This level of max correlation is typically
achieved by FastICA, VCReg or Anica. Figure 12: This level of max correlation
is typically achieved by Whitening.


# Data-based optimal bandwidth for kernel density estimation of statistical
samples

Zhen-Wei Li  Center for Theoretical Physics and College of Physics, Jilin
University, Changchun, 130012, China  Changchun Observatory, National
Astronomical Observatories, CAS, Changchun 130117, China  Ping He
<EMAIL_ADDRESS>Center for Theoretical Physics and College of Physics, Jilin
University, Changchun, 130012, China  Center for High Energy Physics, Peking
University, Beijing 100871, China

###### Abstract

Abstract It’s a common practice to evaluate probability density function or
matter spatial density function from statistical samples. Kernel density
estimation is a frequently used method, but to select an optimal bandwidth of
kernel estimation, which is completely based on data samples, is a long-term
issue that has not been well settled so far. There exist analytic formulae of
optimal kernel bandwidth, but they cannot be applied directly to data samples,
since they depend on the unknown underlying density functions from which the
samples are drawn. In this work, we devise an approach to pick out the totally
data-based optimal bandwidth. First, we derive correction formulae for the
analytic formulae of optimal bandwidth to compute the roughness of the
sample’s density function. Then substitute the correction formulae into the
analytic formulae for optimal bandwidth, and through iteration we obtain the
sample’s optimal bandwidth. Compared with analytic formulae, our approach
gives very good results, with relative differences from the analytic formulae
being only $2\sim 3\%$ for sample size larger than $10^{4}$. This approach can
also be generalized easily to cases of variable kernel estimations.

Key words: numerical methods, kernel density estimation, optimal bandwidth,
large-scale structure of Universe

##  I Introduction

In statistics and many other scientific disciplines such as modern cosmology,
it is a common practice to estimate probability density function from
statistical data, or matter density function in a given spatial region from
observational or simulation samples. In particular, $N$-body simulations are
extensively used to investigate the distribution, formation and evolution of
large-scale nonlinear structures of the Universe, such as galaxies, galaxy
clusters, filaments and inter-galactic gas, and so forth efs85 ; bert98 ;
Springel01 . For example, density profiles of dark matter halos are frequently
investigated, and these density profiles of dark matter halos are usually
considered as universal nfw95 ; nfw96 ; nfw97 ; moore99 ; einasto65 ;
navarro10  and widely used in various cosmological studies. Hence, it is of
great significance to evaluate matter densities of nonlinear structures from
simulation data samples.

There are many excellent methods to evaluate matter density from cosmological
simulation samples, such as Delaunay Tessellation Field Estimator Schaap2000 ,
or the dark-matter sheet method Abel2012 ; Shandarin2012 . In statistics or
cosmology, however, the widely used ones are still histogram method and kernel
estimation method silverman86 . They are parametric methods, i.e. they depend
on the binwidth for histogram or on the bandwidth for kernel method. So how to
make the optimal binwidth or bandwidth is crucial for density estimation from
statistical samples, e.g. Refs. scott92 ; gentle09 .

Kernel estimation is believed to be superior to histogram method jones96 ,
since if the kernel is a derivable function, then the estimated density
function is also derivable. Hence, in this work, we only consider the kernel
estimation method.

A non-negative real-valued function $K_{1}(u)$ is called a kernel function, if
satisfying the following two conditions:

$\int_{-\infty}^{+\infty}K_{1}(u){\rm d}u=1;{\hskip
20.0pt}K_{1}(u)=K_{1}(-u).$ (1)

The density function $f(x)$ of a univariate statistical sample, $x_{i}$, with
$i$ running from 1 to the sample size $N_{p}$, can be estimated by using the
kernel function as:

$\hat{f}_{h}(x)=\frac{1}{N_{p}h}\sum\limits_{i=1}^{N_{p}}K_{1}(\frac{x-x_{i}}{h}),$
(2)

in which $h$ is the bandwidth of the kernel.

As aforementioned, for a proper density estimation, it is crucial to select an
appropriate bandwidth for the kernel function. A small bandwidth $h$ makes the
estimated $\hat{f}_{h}$ too rough and produces spurious features, while a
large bandwidth over-smoothes the data sample, so that useful features are
lost jones96 . Hence, the most important task for kernel estimation is to
choose the optimal bandwidth $h$ for the data sample.

It is possible to select an optimal bandwidth of kernel method by computing
the entropy $S[\hat{f}_{h}]=-\int\hat{f}_{h}(x)\ln\hat{f}_{h}(x){\rm d}x$ with
the estimated density $\hat{f}_{h}(x)$ wolpert13 ; hep14 . Yet, it is very
time-consuming to pick out the optimal bandwidth by evaluating the entropy,
and hence it is not a proper approach to select the optimal bandwidth.

There are indeed analytical formulae of optimal bandwidth for kernel
estimation, which are obtained by minimizing the so-called asymptotic mean
integrated squared error (AMISE) silverman86 ; jones96 . The optimal bandwidth
formulae, however, cannot be applied directly to data samples, since they
depend on the underlying density function $f(x)$, from which the samples are
drawn, but $f$ is usually not known. So some data-based methods of selecting
the optimal bandwidth have been developed, such as the ‘cross validation’
methods Rudemo1982 ; Bowman1984 ; Hall1992 , and ‘plug-in’ methods
Sheather1991 ; Botev2010 . For a comprehensive review about selecting the
optimal bandwidth of kernel estimation, see Ref. jones96 .

In this work, based on the analytic formulae of optimal bandwidth, we devise a
practical method on how to select an optimal bandwidth from data samples
without the knowledge of the underlying density function. In Section II, we
provide the analytic optimal bandwidth of kernel estimation for univariate
samples. We introduce the method in Section III, and present the results, the
generalization to three-dimensional (3D) cases, and a cosmological application
in Section IV. In Section V, we give summary and conclusions of our results.

Figure 1: The one-dimensional mass assignment function of NGP, CIC and TSC
scheme. In this work, we do not use them to assign particle mass to mesh
points, but use them as kernel functions for density estimation of univariate
samples.

##  II Analytic Optimal Bandwidth of Kernel Estimation for Univariate Samples

For a univariate statistical sample, a common way of measuring the error in
the estimation process is the mean integrated squared error (MISE),

${\rm MISE}(h)=E\int(\hat{f}_{h}(x)-f(x))^{2}{\rm d}x,$ (3)

where $E$ denotes the expected value with respect to that sample. MISE is
asymptotically (as $N_{p}\rightarrow\infty$) approximated by the AMISE
silverman86 ; jones96 ,

$\displaystyle{\rm
AMISE}(h)=\frac{R_{1}(K_{1})}{hN_{p}}+h^{4}R_{1}(f^{\prime\prime})\left(\int
x^{2}\frac{K_{1}(x)}{2}{\rm d}x\right)^{2},$ (4)

in which $f^{\prime\prime}$ is the second derivative of $f(x)$ with respect to
$x$, and the functional $R_{1}(g)$ is defined as

$R_{1}(g)\equiv\int_{-\infty}^{+\infty}g(x)^{2}{\rm d}x.$ (5)

By minimizing the AMISE, i.e. ${\rm d}({\rm AMISE})/{\rm d}h=0$, one obtains
the optimal bandwidth $h_{\rm opt}$ for univariate statistical samples as

$h_{\rm opt}=\left[\frac{R_{1}(K_{1})}{R_{1}(f^{\prime\prime})(\int
x^{2}K_{1}(x){\rm d}x)^{2}}\right]^{1/5}N_{p}^{-1/5}.$ (6)

Eq. (6) cannot be directly applied to statistical samples, since it depends on
the unknown density function $f$ and its second derivative $f^{\prime\prime}$.

Figure 2: The optimal bandwidth $h$ vs. the sample size $N_{p}$ for univariate
samples. The red solid lines indicate the analytic bandwidth of Eq. (6). While
the blue solid lines indicate results of our approach, which is purely data-
based. The percentage differences between the two are explicitly shown in the
figure. The density function is the 1D Gaussian function of Eq. (11). Panels
(a), (b) and (c) are for kernel functions of 1D NGP, 1D CIC and 1D TSC,
respectively.

Figure 3: The same as Fig. 2, but the density function for the experimental
sample is the 1D TSC function of Eq. (9). The kernel function we adopted is
also 1D TSC function.

##  III Method

In $N$-body simulation codes such as those based on particle-mesh scheme
hockney81 , particle masses are assigned to a set of pre-defined mesh points
according to some mass assignment scheme. The usual mass assignment approaches
are Nearest Grid Point (NGP), Cloud In Cell (CIC) and Triangular Shaped Cloud
(TSC) schemes. The NGP function is

$\displaystyle W_{\rm NGP}(x)=\left\\{\begin{matrix}&1,&|x|\leq 1/2,\\\
&0,&{\rm otherwise},\end{matrix}\right.$ (7)

the CIC function is

$\displaystyle W_{\rm CIC}(x)=\left\\{\begin{matrix}&1-|x|,&|x|\leq 1,\\\
&0,&{\rm otherwise},\end{matrix}\right.$ (8)

and the TSC function is

$\displaystyle W_{\rm
TSC}(x)=\left\\{\begin{matrix}&\frac{3}{4}-x^{2},&|x|\leq 1/2,\\\
&\frac{1}{2}(\frac{3}{2}-|x|)^{2},&1/2\leq|x|\leq 3/2,\\\ &0,&{\rm
otherwise}.\end{matrix}\right.$ (9)

These functions are shown in Fig. 1. It is interesting to notice that these
mass assignment functions can also be used as kernel functions, since they
satisfy the two conditions of Eq. (1).

In order to demonstrate our approach of picking out the optimal bandwidth $h$
for univariate samples, we generate some experimental samples with known
density functions, from which we draw $N_{p}$ random data points with Monte
Carlo technique.

For a bandwidth $h$, the roughness $R_{1}(f^{\prime\prime})$ in Eq. (6) is
approximated by

$R_{1}(f^{\prime\prime})\approx
R_{1}(\hat{f}^{\prime\prime}_{h})-\frac{6}{wh^{5}N_{p}},$ (10)

in which $w$ is the width of the kernel function. For a kernel function of
finite extension, We define the non-zero range of a kernel function as the
width of the kernel, and hence $w=1$, $2$ and $3$ for NGP, CIC and TSC,
respectively. We give the derivation of Eq. (10) in Appendix-A.

Our approach to pick out the optimal bandwidth $h$ is briefly described as
follows.

  1. 1.

With $W_{\rm NGP}(x)$, $W_{\rm CIC}(x)$ or $W_{\rm TSC}(x)$ as one-dimensional
(1D) kernel functions, we construct the density estimator $\hat{f}_{h}(x)$ of
univariate samples by Eq. (2).

  2. 2.

With this estimated $\hat{f}_{h}(x)$, we calculate $R_{1}(f^{\prime\prime})$
by Eq. (10), and put it into Eq. (6) to update the bandwidth $h$.

  3. 3.

Iterate step 1 to 2 until the relative difference of $h$ between two adjacent
iterations is smaller than a tolerance parameter, say $\leq 0.1\%$.

  4. 4.

Since an extremely small $h$ will make $R_{1}(f^{\prime\prime})$ in Eq. (10)
negative and hence meaningless, we choose a slightly larger $h$ as the initial
value of the bandwidth.

Since the iteration is convergent, in this way, we obtain the optimal
bandwidth $h$ that is completely data-based.

In the next section, we check whether or not our approach can work properly.

Figure 4: The estimated density functions are compared with analytic ones.
Panels (a) and (b) are for 1D Gaussian function, and (c) and (d) are for 1D
TSC function, respectively. The sample size $N_{p}$, the optimal bandwidth $h$
derived from our approach, are shown in the four panels. In all cases, we use
1D TSC as kernel function for the density estimation.

##  IV Results and Application

###  IV.1 Results for univariate samples with unimodal density

As aforementioned, we need some experimental samples to verify our approach.
For this purpose, we consider two sets of experimental data: (1) The first set
of experimental data are drawn from the 1D standard normal distribution,

$f(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}.$ (11)

(2) The 1D TSC function can also be treated as a density function, and we draw
the second set of samples from it by the usual acceptance-rejection method
press92 .

With our approach, it takes just several iterations to obtain the optimal
bandwidth $h$. In Figs. 2 and 3, we compare our results, which are purely
data-based, with the analytic bandwidth $h$ of Eq. (6). We can see that
results obtained with 1D NGP as kernel function are slightly worse for small
sample size $N_{p}$, but those based on 1D TSC functions give the best
results, with the relative differences being only $2\sim 3\%$. In Fig. 4, we
show the density functions estimated by using our optimal bandwidth $h$, to
compare with the analytic density functions. We can see that the agreements
between the estimated and the analytic densities are satisfactory.

Figure 5: The same as Figs. 2 and 3, but the density function for the
experimental sample is the multimodal density of Eq. (12). We use 1D TSC as
kernel function for the density estimation.

Figure 6: The estimated density functions are compared with analytic ones for
multimodal density as Eq. (12). The sample size $N_{p}$, the optimal bandwidth
$h$ derived from our approach, are shown in the figure. We use 1D TSC as
kernel function for the density estimation.

###  IV.2 Results for univariate samples with multimodal density

The approach is not restricted to just unimodal density, but can be also
applied to cases of multimodal density functions. To check its validity for
multimodal densities, we construct the test density function as the following,

$f(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}+\frac{1}{\sqrt{2\pi}\sigma_{1}}e^{-\frac{(x-x_{1})^{2}}{2\sigma^{2}_{1}}}+\frac{1}{\sqrt{2\pi}\sigma_{2}}e^{-\frac{(x-x_{2})^{2}}{2\sigma^{2}_{2}}},$
(12)

with $\sigma_{1}=2$, $x_{1}=-4$, $\sigma_{2}=0.5$, and $x_{2}=4$,
respectively. In Fig. 5, we compare our data-based bandwidth with the analytic
results, and in Fig. 6, we compare the density functions estimated by our
optimal bandwidth $h$ with the analytic density of Eq. (12). It can be seen
that the method is still be valid for estimation of multimodal densities.

Figure 7: The same as Figs. 2 and 3, but the density function for the
experimental sample is the 3D Gaussian function of Eq. (19). Panels (a), (b)
and (c) are for kernel functions of 3D NGP, 3D CIC and 3D TSC, respectively.

###  IV.3 Results for three-dimensional samples

With $|x|$ replaced by $r=(x^{2}_{1}+x^{2}_{2}+x^{2}_{3})^{1/2}$, and with
proper normalization, the 3D generalizations of the 1D kernel functions
$W_{\rm NGP}(x)$, $W_{\rm CIC}(x)$, $W_{\rm TSC}(x)$ of Eqs. (7)-(9) are as
follows

$\displaystyle W_{\rm NGP3}(\textbf{\emph{x}})$ $\displaystyle=$
$\displaystyle\frac{6}{\pi}W_{\rm NGP}(r),$ $\displaystyle W_{\rm
CIC3}(\textbf{\emph{x}})$ $\displaystyle=$ $\displaystyle\frac{3}{\pi}W_{\rm
CIC}(r),$ $\displaystyle W_{\rm TSC3}(\textbf{\emph{x}})$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}W_{\rm TSC}(r),$ (13)

in which $\textbf{\emph{x}}=(x_{1},x_{2},x_{3})$. All these 3D kernel
functions should satisfy

$\int K_{3}(\textbf{\emph{x}}){\rm d}\textbf{\emph{x}}=1;{\hskip
20.0pt}K_{3}(\textbf{\emph{x}})=K_{3}(-\textbf{\emph{x}}).$ (14)

Hence the estimated $\hat{f}_{h}(\textbf{\emph{x}})$ is

$\hat{f}_{h}(\textbf{\emph{x}})=\frac{1}{N_{p}h^{3}}\sum\limits_{i=1}^{N_{p}}K_{3}(\frac{\textbf{\emph{x}}-\textbf{\emph{x}}_{i}}{h}).$
(15)

Eq. (6) is just the optimal bandwidth for univariate samples. Its 3D
generalization is

$h_{\rm opt}=\left[\frac{3R_{3}(K_{3})}{R_{3}(\nabla^{2}f)(\int
x_{1}^{2}K_{3}(\textbf{\emph{x}}){\rm
d}\textbf{\emph{x}})^{2}}\right]^{1/7}N_{p}^{-1/7},$ (16)

in which $K_{3}=K_{3}(\textbf{\emph{x}})$, $f=f(\textbf{\emph{x}})$,
$\nabla^{2}=\partial^{2}_{x_{1}}+\partial^{2}_{x_{2}}+\partial^{2}_{x_{3}}$,
and $R_{3}$ is defined as

$R_{3}(g)\equiv\int g(\textbf{\emph{x}})^{2}{\rm d}\textbf{\emph{x}}.$ (17)

About Eq. (16), we refer the reader to Section-4.3.1 of Ref. silverman86  for
further details.

For a bandwidth $h$, the roughness $R_{3}(\nabla^{2}f)$ in Eq. (16) is
approximated by

$R_{3}(\nabla^{2}f)\approx
R_{3}(\nabla^{2}\hat{f}_{h})-\frac{42}{w^{3}h^{7}N_{p}}.$ (18)

Please see the explanation about Eq. (18) in Appendix-A.

Proceed in the similar way as in Section III, we use 3D standard normal
distribution

$f(\textbf{\emph{x}})=\frac{1}{(2\pi)^{3/2}}e^{-\frac{\textbf{\emph{x}}^{2}}{2}}$
(19)

as the experimental density function, from which we draw a random sample with
size $N_{p}$, and with Eqs. (16) and (18), we iterate to pick out the optimal
bandwidth $h$.

In Fig. 7, we show our results to compare with the analytic bandwidth $h$ of
Eq. (16). It can be seen that results obtained with 3D NGP as kernel function
are obviously worse for small sample size $N_{p}$, with relative errors being
17% and 12% for $N_{p}=10^{3}$ and $10^{4}$ respectively. While computations
based on 3D TSC functions give satisfactory results, with the relative
differences being only $2\sim 3\%$ for $N_{p}>10^{4}$.

###  IV.4 A Cosmological Application

We apply our approach to a realistic model, the Hernquist model, to further
verify its applicability. Hernquist hernquist90  proposed an analytic model of
density profile for spherical galaxies and bulges, as

$\rho(r)=\frac{M_{\rm T}}{2\pi}\frac{r_{c}}{r}\frac{1}{(r_{c}+r)^{3}},$ (20)

where $M_{\rm T}$ is the total mass, and $r_{c}$ is a scale length. With this
density profile, we can generate the spatial distribution of a set of
experimental sample. The code we used for generating this distribution is the
initial-condition generator, taken from the $N$-body simulation code of Ref.
londrillo03 . In a previous work, we have shown that the code works very well
when $r/r_{c}>0.05$ hep17 .

Fig. 8 compares the estimated density profile by the 1D kernel estimation
approach with the analytic model. For Hernquist model, the analytic optimal
bandwidth given by Eq. (6) is $h=0.1712$, while our approach gives $h=0.1678$,
with only -1.9% relative error from the analytic one.

It can be seen that our approach gives very satisfactory result.

Figure 8: The estimated density profile is compared with the analytic one,
i.e. the Hernquist hernquist90  model. The optimal bandwidth $h$ derived from
our approach, is shown in the figure. The total mass in Eq. (20) of Hernquist
model is set to be $M_{\rm T}=1$. The sample is generated with the initial-
condition generator, taken from the $N$-body simulation code of londrillo03 ,
with the sample size $N_{p}=1.05\times 10^{6}$. We use 1D TSC as kernel
function for the density estimation.

##  V Summary and Conclusions

In many scientific disciplines such as statistics and modern cosmology, it is
a common practice to estimate the probability density function or matter
density function in a given spatial region from statistical, observational or
simulation samples. Kernel density estimation is widely used among all the
approaches, but to select an optimal bandwidth of kernel estimation, which is
completely based on data samples, is a long-term issue that has not been well
settled so far. The formulae of Eqs. (6) and (16) do provide the analytic
optimal bandwidth, but they cannot be applied directly to data samples, since
they depend on the unknown underlying density functions from which the samples
are drawn.

In this work, we devise an approach to pick out the totally data-based optimal
bandwidth. First, we derive the correction formulae, Eq. (10) for 1D or Eq.
(18) for 3D samples, to compute the roughness of the sample’s density
function. Then substitute the correction formulae into the analytic formulae
Eqs. (6) or (16) for optimal bandwidth, and through iteration we obtain the
sample’s optimal bandwidth. The whole process terminates just after several
iterations, and compared with the analytic formulae, our approach gives very
good results, whose relative differences from the analytic formulae are only
$2\sim 3\%$ for sample size larger than $10^{4}$.

Technically, we use NGP, CIC and TSC mass assignment schemes as kernel
functions for density estimation, and use the Gaussian and TSC functions as
well as Ref. hernquist90 ’s model as experimental density functions, but we
emphasis that the results should not depend on the adopted kernel functions or
density functions.

Incidentally, if we divide the whole sample into several parts, and apply our
approach of selecting optimal bandwidth to sub-samples part by part, then our
approach is not hard to be generalized to cases of variable kernel
estimations.

###### Acknowledgements.

This work is supported by the National Science Foundation of China (no.
11273013), and also supported by the Natural Science Foundation of Jilin
Province, China (no. 20180101228JC).

##  Appendix A Correction Formula for the Roughness of Density Functions

We give the derivation of the correction formula in Eq. (10) for the roughness
of a 1D density function.

The estimated density $\hat{f}_{h}(x)$ of the sample is a random variable for
any $x$, and should obey Poisson distribution. $\hat{f}_{h}(x)$ can be written
as

$\hat{f}_{h}(x)=f(x)+\Delta\hat{f}_{h}(x),$ (21)

where $f$ is the unknown underlying density function, from which the sample is
drawn, and $\Delta\hat{f}_{h}$ is the difference between the two, whose
ensemble average is of course $\overline{\Delta\hat{f}_{h}}=0$.

We require both $\hat{f}_{h}$ and $f$ satisfy the normalization
$\int\hat{f}_{h}(x)dx=\int f(x)dx=1$. So we have $\hat{f}_{h}(x_{i})\Delta
x={n_{i}}/{N_{p}}$, in which $n_{i}$ is the counted particle number at $x_{i}$
within the interval $\Delta x$. We roughly estimate $\Delta x\approx wh$, in
which $w$ is the width of the kernel function. From Fig. 1, we see that $w=1$,
$2$ and $3$ for NGP, CIC and TSC, respectively. Hence, $\hat{f}_{h}(x_{i})$ is
estimated as

$\hat{f}_{h}(x_{i})\approx\frac{n_{i}}{whN_{p}},$ (22)

and hence

$f(x_{i})=\overline{\hat{f}_{h}(x_{i})}\approx\frac{\overline{n_{i}}}{whN_{p}}.$
(23)

Calculating the second derivative and taking the square of both sides of Eq.
(21), we have

$\hat{f}^{\prime\prime}_{h}(x)^{2}=f^{\prime\prime}(x)^{2}+2f^{\prime\prime}(x)\Delta\hat{f}^{\prime\prime}_{h}(x)+\Delta\hat{f}^{\prime\prime}_{h}(x)^{2},$
(24)

in which the second derivative of the estimated density
$\hat{f}^{\prime\prime}_{h}$ is evaluated by

$\hat{f}^{\prime\prime}_{h}(x_{i})=\frac{\hat{f}_{h}(x_{i+1})+\hat{f}_{h}(x_{i-1})-2\hat{f}_{h}(x_{i})}{h^{2}},$
(25)

in which $x_{i}=x_{i+1}-h=x_{i-1}+h$. Making the ensemble average of
$\hat{f}^{\prime\prime}_{h}(x)^{2}$, we obtain

$\displaystyle\overline{\hat{f}^{\prime\prime}_{h}(x)^{2}}$ $\displaystyle=$
$\displaystyle\overline{f^{\prime\prime}(x)^{2}}+2\overline{f^{\prime\prime}(x)\Delta\hat{f}^{\prime\prime}_{h}(x)}+\overline{\Delta\hat{f}^{\prime\prime}_{h}(x)^{2}}$
(26) $\displaystyle=$ $\displaystyle
f^{\prime\prime}(x)^{2}+\overline{\Delta\hat{f}^{\prime\prime}_{h}(x)^{2}},$

in which $f^{\prime\prime}(x)$ is not random variable, and
$\overline{\Delta\hat{f}^{\prime\prime}_{h}(x)}=0$. Further, from Eq. (25), we
have

$\displaystyle\overline{\Delta\hat{f}^{\prime\prime}_{h}(x_{i})^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{h^{4}}(\overline{\Delta\hat{f}_{h}(x_{i+1})^{2}}+\overline{\Delta\hat{f}_{h}(x_{i-1})^{2}}+4\overline{\Delta\hat{f}_{h}(x_{i})^{2}}$
(27)
$\displaystyle+2\overline{\Delta\hat{f}_{h}(x_{i+1})\Delta\hat{f}_{h}(x_{i-1})}-4\overline{\Delta\hat{f}_{h}(x_{i+1})\Delta\hat{f}_{h}(x_{i})}$
$\displaystyle-4\overline{\Delta\hat{f}_{h}(x_{i})\Delta\hat{f}_{h}(x_{i-1})})=\frac{1}{h^{4}}(\overline{\Delta\hat{f}_{h}(x_{i+1})^{2}}$
$\displaystyle+\overline{\Delta\hat{f}_{h}(x_{i-1})^{2}}+4\overline{\Delta\hat{f}_{h}(x_{i})^{2}}).$

Since the deviation $\Delta\hat{f}_{h}$’s are independent at different $x$, to
arrive at the last equality in Eq. (27), we use

$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i+1})\Delta\hat{f}_{h}(x_{i-1})}$
$\displaystyle\approx$
$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i+1})}\cdot\overline{\Delta\hat{f}_{h}(x_{i-1})}=0,$
$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i+1})\Delta\hat{f}_{h}(x_{i})}$
$\displaystyle\approx$
$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i+1})}\cdot\overline{\Delta\hat{f}_{h}(x_{i})}=0,$
$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i})\Delta\hat{f}_{h}(x_{i-1})}$
$\displaystyle\approx$
$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i})}\cdot\overline{\Delta\hat{f}_{h}(x_{i-1})}=0,$

Also, since $n(x)=N_{p}\hat{f}_{h}(x)wh$ satisfies the Poisson distribution at
$x$, so

$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i+1})^{2}}$ $\displaystyle=$
$\displaystyle\frac{\overline{n_{i+1}}}{w^{2}h^{2}N^{2}_{p}}\approx\frac{f(x_{i+1})}{whN_{p}},$
$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i})^{2}}$ $\displaystyle=$
$\displaystyle\frac{\overline{n_{i}}}{w^{2}h^{2}N^{2}_{p}}\approx\frac{f(x_{i})}{whN_{p}},$
$\displaystyle\overline{\Delta\hat{f}_{h}(x_{i-1})^{2}}$ $\displaystyle=$
$\displaystyle\frac{\overline{n_{i-1}}}{w^{2}h^{2}N^{2}_{p}}\approx\frac{f(x_{i-1})}{whN_{p}},$

in which we considered Eqs. (22) and (23). So Eq. (27) is simplified as

$\overline{\Delta\hat{f}^{\prime\prime}_{h}(x_{i})^{2}}\approx\frac{1}{wh^{5}N_{p}}(f(x_{i+1})+4f(x_{i})+f(x_{i-1})).$
(28)

Substitute Eq. (28) into Eq. (26), and integrate both sides, we obtain

$\displaystyle\int\overline{\hat{f}^{\prime\prime}_{h}(x)^{2}}{\rm d}x$
$\displaystyle=$ $\displaystyle\int f^{\prime\prime}(x)^{2}{\rm
d}x+\int\overline{\Delta\hat{f}^{\prime\prime}_{h}(x)^{2}}{\rm d}x,$ (29)
$\displaystyle\approx$ $\displaystyle\int f^{\prime\prime}(x)^{2}{\rm
d}x+\frac{6}{wh^{5}N_{p}}.$

Since the ensemble average $\overline{\hat{f}^{\prime\prime}_{h}(x)^{2}}$ is
not easy to derive, in practice we can simply approximate it as
$\hat{f}^{\prime\prime}_{h}(x)^{2}$, and finally we obtain Eq. (10).

It’s not hard to generalize the one-dimensional formula Eq. (10) to 3D case.
For 3D samples, Eq. (22) should be adjusted as

$\hat{f}_{h}(\textbf{\emph{x}}_{i})\approx\frac{n_{i}}{w^{3}h^{3}N_{p}}.$ (30)

In the expression of $R_{3}(\nabla^{2}\hat{f}_{h})$, there are total 3 squared
terms of the second derivative $(\partial^{2}\hat{f}_{h}/\partial
x^{2}_{\alpha})^{2}$ ($\alpha=1,2,3$), and each contributes a correction
$6/w^{3}h^{7}N_{p}$. Also, there are total 6 mixed products of the second
derivative as $(\partial^{2}\hat{f}_{h}/\partial
x^{2}_{\alpha})(\partial^{2}\hat{f}_{h}/\partial x^{2}_{\beta})$
($\alpha,\beta=1,2,3,\alpha\neq\beta$), and each contributes a correction
$4/w^{3}h^{7}N_{p}$. So finally, these corrections are added up to
$42/w^{3}h^{7}N_{p}$, and we obtain Eq. (18).

## References

  * (1) G. Efstathiou, M. Davis, C. S. Frenk, and S. D. M. White, Astrophys. J. 57 (1985) 241.
  * (2) E. Bertschinger, ARA&A 36 (1998) 599.
  * (3) V. Springel, N. Yoshida, and S. D. M. White, New Astronomy 6 (2001) 79.
  * (4) J. F. Navarro, C. S. Frenk, and S. D. M. White, MNRAS 275 (1995) 720.
  * (5) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462 (1996) 563.
  * (6) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490 (1997) 493.
  * (7) B. Moore, T. Quinn, F. Governato, J. Stadel, and G. Lake, MNRAS 310 (1999) 1147.
  * (8) J. Einasto, Trudy Inst. Astrofizicheskogo Alma-Ata 51 (1965) 87.
  * (9) J. F. Navarro, A. D. Ludlow, V. Springel, J. Wang, M. Vogelsberger, S. D. M. White, A. Jenkins, C. S. Frenk, and A. Helmi, MNRAS 402 (2010) 21.
  * (10) W. E. Schaap and R. van de Weygaert, A&A 363 (2000) L29.
  * (11) T. Abel, O. Hahn, and R. Kaehler, MNRAS 427 (2012) 61.
  * (12) S. Shandarin, S. Habib, and K. Heitmann, Phys. Rev. D85 (2012) 083005.
  * (13) B. W. Silverman, Density Estimation for Statistics and Data Analysis, Chapman and Hall, London (1986).
  * (14) D. W. Scott, Multivariate Density Estimation: Theory, Practice, and Visualization, Wiley, New York (1992).
  * (15) J. E. Gentle, Computational Statistics. Springer Science+Business Media, New York (2009).
  * (16) M. C. Jones, J. S. Marron, and S. J. Sheather, J. Am. Stat. Assoc. 91 (1996) 401.
  * (17) D. H. Wolpert and S. DeDeo, Entropy 15 (2013) 4668.
  * (18) N. Sui, M. Li, and P. He, MNRAS 445 (2014) 4211.
  * (19) M. Rudemo, Scand. J. Statist. 9 (1982) 65.
  * (20) A. W. Bowman, Biometrika 71 (1984) 353.
  * (21) P. Hall, J. S. Marron, and B. U. Park, Probab. Theory Relat. Fields 92 (1992) 1.
  * (22) S. J. Sheather and M. C. Jones, J. R. Statist. Soc. B 53 (1991) 683.
  * (23) Z. I. Botev, J. F. Grotowski, and D. P. Kroese, Annals of Statistics 38 (2010) 2916.
  * (24) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran - The Art of Scientific Computing, 2nd edt., Cambridge University Press (1992).
  * (25) R. W. Hockney and J. W. Eastwood, Computer Simulations Using Particles, McGraw-Hill, New York (1981).
  * (26) L. Hernquist, Astrophys. J. 356 (1990) 359.
  * (27) P. Londrillo, C. Nipoti, and L. Ciotti, Mem. Soc. Astron. Ital. 1 (2003) 18.
  * (28) N. Sui and P. He, Int. J. Mod. Phys. D 26 (2017) 1750130.