On Discriminative Probabilistic Modeling for Self-Supervised Representation Learning

## Abstract

Abstract We study the discriminative probabilistic modeling problem on a continuous domain for (multimodal) self-supervised representation learning. To address the challenge of computing the integral in the partition function for each anchor data, we leverage the multiple importance sampling (MIS) technique for robust Monte Carlo integration, which can recover InfoNCE-based contrastive loss as a special case. Within this probabilistic modeling framework, we conduct generalization error analysis to reveal the limitation of current InfoNCE-based contrastive loss for self-supervised representation learning and derive insights for developing better approaches by reducing the error of Monte Carlo integration. To this end, we propose a novel non-parametric method for approximating the sum of conditional probability densities required by MIS through convex optimization, yielding a new contrastive objective for self-supervised representation learning.
Moreover, we design an efficient algorithm for solving the proposed objective.
We empirically compare our algorithm to representative baselines on the contrastive image-language pretraining task. Experimental results on the CC3M and CC12M datasets demonstrate the superior overall performance of our algorithm.

## 1 Introduction

Recently, self-supervised learning (SSL) of large models has emerged as a prominent paradigm for building artificial intelligence (AI) systems (Bommasani et al., ; Ozbulak et al., ; Zhou et al., ; Zong et al., ). Unlike traditional supervised learning that relies on human-labeled datasets, SSL can fully exploit the vast multimodal data readily available on the internet via self-supervision (pretext tasks such as instance discrimination, masked modeling, and inpainting) to learn large foundation models that are useful for many downstream tasks. Although self-supervision differs from human supervision, self-supervised and supervised learning share similarities. For instance, many successful self-supervised learning models, e.g., CLIP (Radford et al., ), still use the softmax function and the cross-entropy loss to define their objective functions, similar to traditional multi-class classification in supervised learning. The key difference is that self-supervised learning focuses on predicting relevant data instead of relevant labels.

Figure: Figure 1: DPM for supervised learning and self-supervised representation learning.
Refer to caption: x1.png

Discriminative probabilistic modeling (DPM) uses a parameterized model to capture the conditional probability $\Pr(\mathbf{y}|\mathbf{x})$ of a target $\mathbf{y}\in\mathcal{Y}$ given an input data point $\mathbf{x}$. It is a fundamental approach for supervised learning (see, e.g., Chapter 4.3 in Bishop ). For example, logistic regression for multi-class classification (MCC) uses $\Pr(\mathbf{y}|\mathbf{x})$ to define the probability of a label $\mathbf{y}$ given a data point $\mathbf{x}$, whose maximum likelihood estimation (MLE) yields the cross-entropy (CE) loss. Similarly, probabilistic modeling approaches such as ListNet (Cao et al., ) have been used for learning to rank (L2R) to model the probability of a candidate $\mathbf{y}$ in a list given a query $\mathbf{x}$. In these supervised learning tasks, the target $\mathbf{y}$ is from a discrete and finite set $\mathcal{Y}$ (e.g. class labels or a list of candidates).

What if the target $\mathbf{y}$ in DPM is from a continuous domain $\mathcal{Y}$? This is particularly useful for modeling the prediction task of self-supervised representation learning. Considering that each underlying object in the real world generates various forms of observational data, such as images, texts, and audio, DPM is a natural choice to model the probability of observing a data point from a continuous domain (e.g., the space of natural images, audio, or the continuous input embedding space of texts) given an “anchor” data point. The anchor data may come from a different modality. However, solving DPM over a continuous domain is deemed as a challenging task (c.f. Section 1.3 in LeCun et al. ).
Compared to the probabilistic modeling over discrete and finite sets, such as in traditional supervised learning tasks like MCC and L2R, the DPM problem over a continuous domain (real vector space) necessitates computing the partition function (i.e., the normalizing constant) for each anchor. This involves an integration over an underlying continuous space, rather than a finite summation. Previous works tackle or circumvent this issue using Markov Chain Monte Carlo (MCMC) sampling (Song and Kingma, ; Du and Mordatch, ; Ta et al., ), noise-contrastive estimation (NCE) (Gutmann and Hyvärinen, ), and energy-based models (EBMs) (LeCun et al., ; Assran et al., ). However, these approaches have drawbacks that hinder their wide adoption in self-supervised representation learning. Specifically, MCMC approaches are computationally expensive and prone to divergence during training (Yang and Ji, ; Kim and Ye, ), NCE suffers from slow convergence for high-dimensional data when the noise distribution is not well-configured (Liu et al., ), and EBMs neglects the partition functions that may throw away some discriminative power of DPM. The state-of-the-art discriminative self-supervised representation learning algorithms rely on InfoNCE-based contrastive losses (Oord et al., ; Chen et al., ; Radford et al., ). However, the connection between these losses and DPMs is not fully revealed, limiting our understanding of their strengths and weaknesses.

In this work, we study DPM over a continuous domain for self-supervised representation learning by investigating a computational framework of robust Monte Carlo integration of the partition functions based on the multiple importance sampling (MIS) approach (Veach and Guibas, ; Veach, ). Our contributions can be summarized as follows:

$\bullet$ We construct an empirical risk for maximum likelihood estimation (MLE) based on MIS. Unlike the setting in MIS for traditional stochastic simulation, the sums of conditional densities on training data are not directly accessible, requiring approximations of these sums. We show that an InfoNCE-based contrastive loss, referred to as the global contrastive loss (GCL) (Yuan et al., ), can be derived as a special case of this empirical risk when a simple uniform approximation is used.

$\bullet$ Our finite-sample generalization analysis reveals that GCL incurs a non-diminishing error due to the rather crude uniform approximation. To reduce this error, we introduce a novel non-parametric method to approximate the sum of conditional densities required by the MIS approach through optimization. In a toy experiment, we demonstrate that the proposed non-parametric method results in better approximations and significantly smaller generalization errors than GCL.

$\bullet$ The MLE framework based on MIS and our non-parametric method inspires a new contrastive objective for self-supervised representation learning. We design an efficient algorithm, NUCLR, to optimize the proposed objective, which does not require the costly MCMC steps or the elusive noise distribution selection. Furthermore, NUCLR can be interpreted as a contrastive representation learning algorithm with *nonuniform* margins that are dynamically learned for negative pairs.

$\bullet$ Experimental results on language-image pretraining using CC3M and CC12M datasets show the superior overall performance of NUCLR compared to baselines on downstream retrieval tasks.

### 1.1 Related Work

Probabilistic Models for Self-Supervised Representation Learning: Discriminative probabilistic models learn the *conditional* probability mass/density function $p(\mathbf{y}\mid\mathbf{x})$ of $\mathbf{y}$ given data $\mathbf{x}$. Recently, some works have focused on modeling the conditional probability density function $p(\mathbf{y}\mid\mathbf{x})$ for the unsupervised representation learning task, where both $\mathbf{x}$ and $\mathbf{y}$ may belong to uncountable spaces. Khemakhem et al. () studied the identifiability (i.e., the learned representations
are unique up to a linear transformation) of DPM and showed its connection to nonlinear ICA models. Ta et al. () improved the Langevin MCMC method to handle the partition function in DPM for learning implicit representations of behavior-cloned policies in robotics. Discarding the partition function, Assran et al. (); Bardes et al. () proposed the energy-based models I-JEPA and V-JEPA to learn visual representations by predicting the relevance between data representations. Although the high-level concept of JEPA is similar to our work in that both aim to predict the relevance between data representations, our approach is grounded in discriminative probabilistic modeling, whereas JEPA is an energy-based model that omits the partition function. Consequently, JEPA lacks some statistical guarantees of probabilistic models, such as the convergence of the maximum likelihood estimator, which have implications for performance on downstream tasks (See Remark 2). Furthermore, JEPA is designed specifically for the visual modality whereas our algorithm applies to multimodal data.

By modeling the joint distribution $p(\mathbf{x},\mathbf{y})$, hybrid models (Grathwohl et al., ; Wang et al., ; Kim and Ye, ; Bizeul et al., ) simultaneously perform discriminative and generative probabilistic modeling.
Although the generative component in hybrid models might offer some benefits for representation learning, such as achieving reasonably good performance with small batch size, Kim and Ye () pointed out that current hybrid models significantly increase the computational burden and cannot be applied to large-scale datasets such as ImageNet1k due to the expensive inner loops of stochastic gradient Langevin dynamics or adversarial training. Furthermore, Bizeul et al. () mentioned(^1^11[https://openreview.net/forum?id=QEwz7447tR](https://openreview.net/forum?id=QEwz7447tR)) that generative representation learning approaches face difficulties scaling to large-scale, complex datasets, as “learning representations for complex data distributions under a generative regime remains a challenge compared to discriminative approaches.”

InfoNCE-based Losses and Theoretical Guarantees: The InfoNCE loss is arguably the most widely used loss function in contrastive learning (Chopra et al., ; Oord et al., ; Chen et al., ; Radford et al., ). Given a dataset $\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{n}$ from two views or modalities, the minibatch-based InfoNCE loss contrasts each positive pair with $k$ negative pairs in the sampled batch. Both empirical observations (Chen et al., ; Radford et al., ; Yuan et al., ) and theoretical analysis (Yuan et al., ) demonstrate that algorithms based on InfoNCE perform well when the batch size is sufficiently large (e.g. 32,768 for CLIP training), which demands a lot of
computational resources. Besides, several works analyze the generalization error of the minibatch-based InfoNCE loss (Arora et al., ; Lei et al., ). However, these analyses have a critical limitation: the generalization error increases with $k$, which contradicts practical observations.

To address the dependency on large batch sizes, Yuan et al. () studied the global contrastive loss (GCL),
which can be viewed as a variant of the InfoNCE loss that contrasts each positive pair with *all* other negative pairs in the whole dataset. Based on the exponential moving average (EMA) technique, they developed the SogCLR algorithm that converges to a neighborhood of a stationary point of GCL even with small batch sizes (e.g., 256). Recently, Yau et al. () introduced the MCMC-based EMC^2 algorithm that converges to the stationary point with a small batch size. However, EMC^2 appears to empirically perform worse than SogCLR on large datasets such as ImageNet1k. Besides, Waida et al. () established the generalization bound of the kernel contrastive loss (KCL), which is a lower bound of GCL when the kernel is bilinear.

## 2 Preliminaries

Notations: We use $[n]$ to denote the set $\{1,\dotsc,n\}$. For a vector $\mathbf{v}\in\mathbb{R}^{n}$, $v^{(i)}$ represents its $i$-th coordinate and we define $\exp(\mathbf{v})\coloneqq(\exp(v^{(1)}),\dotsc,\exp(v^{(n)}))^{\top}$. Let $\mathcal{X}$ and $\mathcal{Y}$ be Lebesgue-measurable data spaces of different views or modalities. When $\mathcal{X}$ or $\mathcal{Y}$ is finite, the Lebesgue measure $\mu$ is replaced by the counting measure. For a sequence of random variables $\{\mathbf{x}_{n}\}_{n\in\mathbb{N}}$, we write $\mathbf{x}_{n}\stackrel{{\scriptstyle p}}{{\rightarrow}}\mathbf{x}$ to denote that $\{\mathbf{x}_{n}\}_{n\in\mathbb{N}}$ converges in probability to $\mathbf{x}$.

Multiple Importance Sampling (MIS): Next, we briefly review the multiple importance sampling (MIS) approach (Veach and Guibas, ; Veach, ) for robust Monte Carlo integration. MIS was originally introduced to address the glossy highlights problem for image rendering in computer graphics, which involves computing several integrals of the form $g(\mathbf{a},r,s)=\int_{\mathcal{X}}f(\mathbf{x};r,s)d\mu(\mathbf{x})$, which represents the light exiting a point $\mathbf{a}$ on a glossy surface under variations in light size $s$ and surface glossiness $r$. For Monte Carlo integration of $g(\mathbf{a},r,s)$, importance sampling based on a sample from a single distribution may lead to a large variance under some light size/surface glossiness. To address this issue, the MIS approach constructs an unbiased estimator $\hat{g}(\mathbf{a},r,s)=\sum_{j=1}^{n}\frac{1}{m}\sum_{l=1}^{m}\omega^{(j)}(
\mathbf{x}_{j,l})\frac{f(\mathbf{x}_{j,l};\mathbf{a},r,s)}{p_{j}(\mathbf{x}_{j
,l})}$ by combining samples $\hat{\mathbf{X}}_{1}\,\dotsc,\hat{\mathbf{X}}_{n}$ from distributions $p_{1},\dotsc,p_{n}$, where $\hat{\mathbf{X}}_{j}\coloneqq\{\mathbf{x}_{j,1},\dotsc,\mathbf{x}_{j,m}\}$ is sampled from the $j$-th distribution, and $\boldsymbol{\omega}=(\omega^{(1)},\dotsc,\omega^{(n)})$ is a weighting function satisfies that $\sum_{j=1}^{n}\omega^{(j)}(\mathbf{x})=1$ when $f(\mathbf{x};\mathbf{a},r,s)\neq 0$ and $\omega^{(j)}(\mathbf{x})=0$ when $p_{j}(\mathbf{x})=0$. Moreover, Veach and Guibas () proposed the “balance heuristic” $\omega^{(j)}(\mathbf{x})=\frac{p_{j}(\mathbf{x})}{\sum_{j^{\prime}=1}^{n}p_{j^
{\prime}}(\mathbf{x})}$ and proved that this choice of $\boldsymbol{\omega}$ is near-optimal in terms of variance among all possible weighting functions and the resulting MIS estimator is $\hat{g}(\mathbf{a},r,s)=\sum_{j=1}^{n}\frac{1}{m}\sum_{l=1}^{m}\frac{f(\mathbf
{x}_{j,l};\mathbf{a},r,s)}{\sum_{j^{\prime}=1}^{n}p_{j^{\prime}}(\mathbf{x}_{j
,l})}$. Empirically, they showed that MIS with the balance heuristic leads to better rendering performance than importance sampling with a single distribution.

## 3 DPM over a Continuous Domain

When choosing $\mathcal{X}$ as the anchor space, we model the probability density $p(\mathbf{y}\mid\mathbf{x})$ of an object $\mathbf{y}\in\mathcal{Y}$ given an anchor object $\mathbf{x}\in\mathcal{X}$ by the following DPM parameterized by $\mathbf{w}$:

$$ $\displaystyle p_{\mathbf{w}}(\mathbf{y}\mid\mathbf{x})=\frac{\exp(E_{\mathbf{w }}(\mathbf{x},\mathbf{y})/\tau)}{\int_{\mathcal{Y}}\exp(E_{\mathbf{w}}(\mathbf {x},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^{\prime})},$ (1) $$

where $\tau>0$ is a temperature parameter for flexibility, $E_{\mathbf{w}}:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}$ is a parameterized prediction function, e.g., $E_{\mathbf{w}}(\mathbf{x},\mathbf{y})=E_{1}(\mathbf{w}_{1},\mathbf{x})^{\top}E
_{2}(\mathbf{w}_{2},\mathbf{y})$, where $E_{1}$ and $E_{2}$ are encoder networks. We assume that $\exp(E_{\mathbf{w}}(\mathbf{x},\mathbf{y})/\tau)$ is Lebesgue-integrable for $\mathbf{w}\in\mathcal{W}$, $\mathcal{W}\subset\mathbb{R}^{d}$. Here $p_{\mathbf{w}}(\mathbf{y}\mid\mathbf{x})$ is a valid probability density function because $\int_{\mathcal{Y}}p_{\mathbf{w}}(\mathbf{y}\mid\mathbf{x})d\mu(\mathbf{y})=1$.
Given $\{(\mathbf{x}_{1},\mathbf{y}_{1}),\dotsc,(\mathbf{x}_{n},\mathbf{y}_{n})\}$ sampled from the joint distribution $p_{\mathbf{x},\mathbf{y}}$, the maximum likelihood estimation (MLE) aims to solve

$$ $\displaystyle\min_{\mathbf{w}}\left\{-\frac{1}{n}\sum_{i=1}^{n}\tau\log\frac{ \exp(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i})/\tau)}{\int_{\mathcal{Y}} \exp(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^{ \prime})}\right\}.$ (2) $$

Corresponding to the empirical objective above, the true (expected) risk can be defined as

$$ $\displaystyle\mathcal{L}(\mathbf{w})\coloneqq\mathbb{E}_{\mathbf{x},\mathbf{y} }\left[-\tau\log\frac{\exp(E_{\mathbf{w}}(\mathbf{x},\mathbf{y})/\tau)}{\int_{ \mathcal{Y}}\exp(E_{\mathbf{w}}(\mathbf{x},\mathbf{y}^{\prime})/\tau)d\mu( \mathbf{y}^{\prime})}\right].$ (3) $$

###### Remark 1 .

When $\mathcal{Y}$ is a finite set, (3) reduces to the cross-entropy loss commonly used in supervised learning. Besides, the true risk in (3) is subtly different from the asymptotic of the InfoNCE ($\infty$-InfoNCE) loss (Wang and Isola, ), which can be expressed as $\mathcal{L}_{\infty\text{-InfoNCE}}(\mathbf{w})=\mathbb{E}_{\mathbf{x},\mathbf
{y}}\left[-\tau\log\frac{\exp(E_{\mathbf{w}}(\mathbf{x},\mathbf{y})/\tau)}{
\mathbb{E}_{\mathbf{y}^{\prime}}[\exp(E_{\mathbf{w}}(\mathbf{x},\mathbf{y}^{
\prime})/\tau)]}\right]$. The difference lies in the denominator within the logarithm. It is not clear how to interpret the $\infty$-InfoNCE loss from the probabilistic modeling viewpoint, since $p^{\prime}_{\mathbf{w}}(\mathbf{y}\mid\mathbf{x})=\frac{\exp(E_{\mathbf{w}}(
\mathbf{x},\mathbf{y})/\tau)}{\mathbb{E}_{\mathbf{y}^{\prime}}[\exp(E_{\mathbf
{w}}(\mathbf{x},\mathbf{y}^{\prime})/\tau)]}$ is not a valid probability density function.

###### Remark 2 .

Learning the DPM $p_{\hat{\mathbf{w}}_{*}}$ via MLE for self-supervised pretraining naturally provides some performance guarantees for downstream discriminative tasks. Suppose that the true conditional density function is parameterized by some $\mathbf{w}_{*}\in\mathcal{W}$, i.e., $p(\mathbf{y}\mid\mathbf{x})=p_{\mathbf{w}_{*}}(\mathbf{y}\mid\mathbf{x})=\frac
{\exp(E_{\mathbf{w}_{*}}(\mathbf{x},\mathbf{y})/\tau)}{\int_{\mathcal{Y}}\exp(
E_{\mathbf{w}_{*}}(\mathbf{x},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^{
\prime})}$ for any $\mathbf{x}\in\mathcal{X},\mathbf{y}\in\mathcal{Y}$. Then, the maximum likelihood estimator $\hat{\mathbf{w}}_{*}=\operatorname*{arg\,max}_{\mathbf{w}\in\mathcal{W}}\frac{
1}{n}\sum_{i=1}^{n}\log p_{\mathbf{w}}(\mathbf{y}_{i}\mid\mathbf{x}_{i})$ with the sample $\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{n}$ of size $n$ converges in probability to $\mathbf{w}_{*}$ under some mild assumptions (see Theorem 2.1 in Newey and McFadden ). Due to the continuous mapping theorem, the learned model satisfies $E_{\hat{\mathbf{w}}_{*}}(\mathbf{x},\mathbf{y})\stackrel{{\scriptstyle p}}{{
\rightarrow}}E_{\mathbf{w}_{*}}(\mathbf{x},\mathbf{y})$ and $p_{\hat{\mathbf{w}}_{*}}(\mathbf{y}|\mathbf{x})\stackrel{{\scriptstyle p}}{{
\rightarrow}}p_{\mathbf{w}_{*}}(\mathbf{y}|\mathbf{x})$ if the parameterized models $E_{\mathbf{w}}$ and $p_{\mathbf{w}}$ have measure-zero discontinuity points on $\mathcal{W}$, which naturally provides a statistical guarantee for cross-modality retrieval. In App. C, we also discuss the performance of DPM on downstream classification tasks.

When choosing $\mathcal{Y}$ as the anchor space, we can also model the probability density of an object $\mathbf{x}\in\mathcal{X}$ given an anchor $\mathbf{y}\in\mathcal{Y}$ by the parameterized model $p_{\mathbf{w}}(\mathbf{x}\mid\mathbf{y})=\frac{\exp(E_{\mathbf{w}}(\mathbf{x},
\mathbf{y})/\tau)}{\int_{\mathcal{X}}\exp(E_{\mathbf{w}}(\mathbf{x}^{\prime},
\mathbf{y})/\tau)d\mu(\mathbf{x}^{\prime})}$ similar to (1). Based on a sample $\hat{\mathbf{S}}=\hat{\mathbf{X}}\times\hat{\mathbf{Y}}=\{(\mathbf{x}_{1},
\mathbf{y}_{1}),\dotsc,(\mathbf{x}_{n},\mathbf{y}_{n})\}$ of size $n$ from the joint distribution $p_{\mathbf{x},\mathbf{y}}$, we can simultaneously model $p_{\mathbf{w}}(\mathbf{y}\mid\mathbf{x})$ and $p_{\mathbf{w}}(\mathbf{x}\mid\mathbf{y})$ via the objective below, which resembles the symmetric InfoNCE loss in Radford et al. ():

$$ $\displaystyle\min_{\mathbf{w}}-\frac{1}{n}\sum_{i=1}^{n}\left(\tau\log\frac{ \exp(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i})/\tau)}{\int_{\mathcal{Y}} \exp(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^{ \prime})}+\tau\log\frac{\exp(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i})/ \tau)}{\int_{\mathcal{X}}\exp(E_{\mathbf{w}}(\mathbf{x}^{\prime},\mathbf{y}_{i })/\tau)d\mu(\mathbf{x}^{\prime})}\right).$ $$

###### Remark 1 .

###### Remark 2 .

### 3.1 An MIS-based Empirical Risk for Maximum Likelihood Estimation

For simplicity, we focus on the case where $\mathcal{X}$ is the anchor space and aim to model $p(\mathbf{y}\mid\mathbf{x})$. The main challenge of MLE in (2) lies in computing the integral $g(\mathbf{w};\mathbf{x}_{i},\mathcal{Y})\coloneqq\int_{\mathcal{Y}}\exp\left(E
_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}^{\prime})/\tau\right)d\mu(\mathbf{y}^{
\prime})$ for each $i\in[n]$, which is infeasible unless $\mathcal{Y}$ is finite. For Monte Carlo integration based on the importance sampling, it is difficult, if not impossible, to select a single distribution that works well for all integrals $g(\mathbf{w};\mathbf{x}_{i},\mathcal{Y})$, $i\in[n]$. Assume that each data $\mathbf{y}_{j}$ is sampled from $p_{j}=p(\cdot\mid\mathbf{x}_{j})$ given data $\mathbf{x}_{j}$ sampled from the marginal, $j=1,2,\dotsc,n$.
Thus, we employ the MIS method
with balance heuristic (Veach and Guibas, ) to construct the following estimator of $g(\mathbf{w};\mathbf{x}_{i},\mathcal{Y})$ by combining sampled data $\mathbf{y}_{1},\dotsc,\mathbf{y}_{n}$ from $n$ distributions $p_{1},\dotsc,p_{n}$:

$$ $\displaystyle\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}})$ $\displaystyle=\sum_{j=1}^{n}\frac{1}{\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{j} \mid\mathbf{x}_{j^{\prime}})}\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y }_{j})/\tau\right),\quad\hat{\mathbf{Y}}=\{\mathbf{y}_{1},\dotsc,\mathbf{y}_{n }\}.$ (4) $$

In App. B, we show that the estimator $\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}})$ in (4) is an unbiased estimator of $g(\mathbf{w};\mathbf{x}_{i},\mathcal{Y})$ and explain why we choose the balance heuristic over other possible weighting functions for MIS.

###### Remark 3 .

The variance of MIS-based estimator $\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}})$ can be further reduced if we sample multiple data $\{\mathbf{y}_{j,l}\}_{l=1}^{m}$ from each $p_{j}$, $j\in[n]$. The MIS-based estimator is then constructed as

$$ $\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}})=\sum_{j=1}^{n}\frac{1}{m} \sum_{l=1}^{m}\frac{1}{\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{j,l}\mid\mathbf{x} _{j^{\prime}})}\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j,l})/\tau \right),\quad\hat{\mathbf{Y}}=\bigcup_{j=1}^{n}\{\mathbf{y}_{j,1},\dotsc, \mathbf{y}_{j,m}\}.$ $$

In practice, $\mathbf{y}_{j,1},\dotsc,\mathbf{y}_{j,m}$ can be approximated by random augmentations of $\mathbf{y}_{j}$, assuming that the neighborhood of a high-density point also has relatively high density. Using a sample of $\{\mathbf{y}_{j,1},\dotsc,\mathbf{y}_{j,m}\}$ of size $m>1$, as opposed to a single $\mathbf{y}_{j}$, has been empirically shown to improve performance in bimodal contrastive learning (Fan et al., ; Li et al., ).

However, a remaining issue prevents us from using the MIS-based estimator in (4). Unlike the rendering problem considered in Veach and Guibas (), we do not have access to the conditional probability density $p(\mathbf{y}_{j}\mid\mathbf{x}_{j^{\prime}})$, $j,j^{\prime}\in[n]$. Thus, there is a need for a cheap approximation $\tilde{q}^{(j)}$ of the sum of conditional probability densities $q^{(j)}\coloneqq\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{j}\mid\mathbf{x}_{j^{
\prime}})$, $\forall j\in[n]$, where $q^{(j)}$ can be viewed as a measure of popularity of $\mathbf{y}_{j}$ if we consider that all data in $\hat{\mathbf{S}}=\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{n}$ form a graph and the weight on edge from node $\mathbf{x}$ to node $\mathbf{y}$ is $p(\mathbf{y}|\mathbf{x})$. With a general approximation $\tilde{\mathbf{q}}=(\tilde{q}^{(1)},\dotsc,\tilde{q}^{(n)})^{\top}$ of $\mathbf{q}=(q^{(1)},\dotsc,q^{(n)})^{\top}$, the MLE objective in (2) with MIS can be written as

$$ $\displaystyle\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}},\hat{\mathbf{S}}) =-\frac{1}{n}\sum_{i=1}^{n}\tau\log\frac{\exp(E_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y}_{i})/\tau)}{\tilde{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}})}, \quad\tilde{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}})=\sum_{j=1}^{n}\frac {\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j})/\tau\right)}{\tilde{q }^{(j)}}.$ (5) $$

###### Remark 4 .

If we simply choose the uniform approximation $\tilde{\mathbf{q}}=nc\mathbf{1}_{n}$ with some $c>0$, minimizing $\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}},\hat{\mathbf{S}})$ in (5) is equivalent to minimizing the global contrastive loss (GCL) (Yuan et al., ).

###### Remark 3 .

###### Remark 4 .

### 3.2 Finite-Sample Generalization Analysis

Next, we analyze the error between the empirical risk $\hat{\mathcal{L}}(\mathbf{w},\tilde{\mathbf{q}};\hat{\mathbf{S}})$ in (5) with a *general* approximation $\tilde{\mathbf{q}}$ and the true risk $\mathcal{L}(\mathbf{w})$ in (3) for DPM. This analysis provides (i) insights into the statistical error of GCL (Yuan et al., ), and (ii) guidance on finding an approximation $\tilde{\mathbf{q}}$ better than the uniform one used by GCL as discussed in Remark 4. First, we state the necessary assumptions.

###### Assumption 1 .

There exist $c_{1},c_{2}>0$ such that $\left\|\mathbf{x}\right\|_{2}\leq c_{1}$, $\left\|\mathbf{y}\right\|_{2}\leq c_{2}$ for any $\mathbf{x}\in\mathcal{X},\mathbf{y}\in\mathcal{Y}$.

We focus on representation learning such that $E_{\mathbf{w}}(\mathbf{x},\mathbf{y})=E_{1}(\mathbf{w}_{1};\mathbf{x})^{\top}E
_{2}(\mathbf{w}_{2};\mathbf{y})$, where $\mathbf{w}_{1}$ and $\mathbf{w}_{2}$ are the encoders+projection heads of the first and second views/modalities, respectively. In our theory, we consider the case that both $E_{1}$ and $E_{2}$ are $L$-layer neural networks(^2^22Our results could potentially be extended to other neural networks, such as ConvNets, using the corresponding Rademacher complexity bounds (See e.g., Truong, ).) with positive-homogeneous and 1-Lipschitz continuous activation function $\sigma(\cdot)$ (e.g. ReLU).

###### Assumption 2 .

Suppose that $E_{1}(\mathbf{w}_{1};\mathbf{x})\in\mathbb{R}^{d_{L}}$, $E_{2}(\mathbf{w}_{2};\mathbf{y})\in\mathbb{R}^{d_{L}}$ for some $d_{L}\geq 1$. Moreover, we have $\left\|E_{1}(\mathbf{w}_{1};\mathbf{x})\right\|_{2}\leq 1$, $\left\|E_{2}(\mathbf{w}_{2};\mathbf{y})\right\|_{2}\leq 1$ such that $E_{\mathbf{w}}(\mathbf{x},\mathbf{y})\in[-1,1]$.

Based on the assumptions above, we provide a finite-sample generalization error bound between the empirical risk $\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}},\hat{\mathbf{S}})$ in (5) and the true risk $\mathcal{L}(\mathbf{w})$ in (3).

###### Theorem 1 .

Suppose that Assumptions (1) and (2) hold. Consider the prediction function $E_{\mathbf{w}}$ parameterized by $L$-layer deep neural networks $\mathbf{w}_{1},\mathbf{w}_{2}$ and an approximation $\tilde{\mathbf{q}}$ of $\mathbf{q}$, where $q^{(j)}=\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{j}\mid\mathbf{x}_{j^{\prime}})
\geq\Omega(n)$ almost surely, $\forall j\in[n]$. With probability at least $1-\delta$, $\delta\in(0,1)$,

$$ $\displaystyle\left|\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}},\hat{ \mathbf{S}})-\mathcal{L}(\mathbf{w})\right|\leq O\left(\frac{1}{n}+\sqrt{\frac {d_{L}}{n}}+\sqrt{\frac{\log(1/\delta)}{n}}+\mathcal{E}_{\mathbf{w}}(\tilde{ \mathbf{q}},\mathbf{q};\hat{\mathbf{S}})\right),$ (6) $$

where $\mathcal{E}_{\mathbf{w}}(\tilde{\mathbf{q}},\mathbf{q};\hat{\mathbf{S}})
\coloneqq\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\left|\frac{1}{\tilde{q}^{(j)}
}-\frac{1}{q^{(j)}}\right|\exp((E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i})-
1)/\tau)$ is an error term.

###### Remark 5 .

(i) The uniform approximation $\tilde{q}^{(j)}=nc$ for all $j\in[n]$ used by the GCL leads to a non-diminishing error term $\mathcal{E}(\tilde{\mathbf{q}},\mathbf{q};\hat{\mathbf{S}})$ for solving DPM over a continuous domain;
(ii) Moreover, the error term $\mathcal{E}(\tilde{\mathbf{q}},\mathbf{q};\hat{\mathbf{S}})$ vanishes when $\mathcal{Y}$ is a finite set. Then, the result reproduces the classical result in the literature for supervised learning (Boucheron et al., ).

###### Assumption 1 .

###### Assumption 2 .

###### Theorem 1 .

###### Remark 5 .

### 3.3 Non-parametric Method for Approximating the Measure of Popularity

In Section 3.2, we show that simply choosing a uniform $\tilde{\mathbf{q}}$ as in GCL to approximate the measure of popularities $\mathbf{q}$ (i.e., the sum of conditional probability densities) leads to a non-diminishing term in generalization error. Next, we present a new way to approximate the measure of popularities $\mathbf{q}$. For brevity, we denote by $E(\cdot,\cdot)=E_{\mathbf{w}_{*}}(\cdot,\cdot)$ that corresponds to the real conditional density $p(\mathbf{y}\mid\mathbf{x})=p_{\mathbf{w}_{*}}(\mathbf{y}\mid\mathbf{x})=\frac
{\exp(E_{\mathbf{w}_{*}}(\mathbf{x},\mathbf{y})/\tau)}{\int_{\mathcal{Y}}\exp(
E_{\mathbf{w}_{*}}(\mathbf{x},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^{
\prime})}$. Thus, for any $j\in[n]$ we have

$$ $\displaystyle q^{(j)}\!=\!\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{j}\mid\mathbf{x }_{j^{\prime}})\!=\!\sum_{j^{\prime}=1}^{n}\frac{\exp(E(\mathbf{x}_{j^{\prime} },\mathbf{y}_{j})/\tau)}{\int_{\mathcal{Y}}\exp(E(\mathbf{x}_{j^{\prime}}, \mathbf{y})/\tau)d\mu(\mathbf{y})}\!\stackrel{{\scriptstyle\diamondsuit}}{{ \approx}}\!\sum_{j^{\prime}=1}^{n}\frac{\exp(E(\mathbf{x}_{j^{\prime}},\mathbf {y}_{j})/\tau)}{\sum_{i^{\prime}=1}^{n}\frac{1}{q^{(i^{\prime})}}\exp(E( \mathbf{x}_{j^{\prime}},\mathbf{y}_{i^{\prime}})/\tau)},$ (7) $$

where the last step $\diamondsuit$ is due to the MIS-based Monte Carlo integration and its error decreases when $n$ increases (See Prop. 1 in App. B). Since the expression in (7) is implicit, we propose a non-parametric method to approximate $\mathbf{q}$ by solving the following convex optimization problem:

$$ $\displaystyle\min_{\boldsymbol{\zeta}\in\mathbb{R}^{n}}\left\{-\frac{1}{n}\sum _{i=1}^{n}\tau\log\left(\frac{\exp(E(\mathbf{x}_{i},\mathbf{y}_{i})/\tau)}{ \sum_{j=1}^{n}\exp((E(\mathbf{x}_{i},\mathbf{y}_{j})-\zeta^{(j)})/\tau)}\right )+\frac{1}{n}\sum_{j=1}^{n}\zeta^{(j)}\right\}.$ (8) $$

The following theorem characterizes the set of optimal solutions to (8) and its relationship to $\mathbf{q}$.

###### Theorem 2 .

The optimal solution to (8) is unique up to an additive scalar. In other words, if $\boldsymbol{\zeta}_{*}$ belongs to the optimal solutions set $\mathcal{Z}_{*}$ of (8), then $\boldsymbol{\zeta}_{*}+c\mathbf{1}_{n}$ is also in $\mathcal{Z}_{*}$ for any $c\in\mathbb{R}$. Moreover, if we define the set $\mathcal{Q}_{*}\coloneqq\{\exp(\boldsymbol{\zeta}_{*}/\tau)\mid\boldsymbol{
\zeta}_{*}\in\mathcal{Z}_{*}\}$, then any $\bar{\mathbf{q}}\in\mathcal{Q}_{*}$ satisfies the following equation

$$ $\displaystyle\bar{q}^{(j)}=\sum_{j^{\prime}=1}^{n}\frac{\exp(E(\mathbf{x}_{j^{ \prime}},\mathbf{y}_{j})/\tau)}{\sum_{i^{\prime}=1}^{n}\frac{1}{\bar{q}^{(i^{ \prime})}}\exp(E(\mathbf{x}_{j^{\prime}},\mathbf{y}_{i^{\prime}})/\tau)},\quad \forall j\in[n].$ (9) $$

and the set $\mathcal{Q}_{*}$ is closed under scalar multiplication, i.e., any $\bar{\mathbf{q}}\in\mathcal{Q}_{*}$ and $C>0$ implies $C\bar{\mathbf{q}}\in\mathcal{Q}_{*}$.

###### Remark 6 .

Due to Theorem 2 and the resemblance between (7) and (9), there exists a specific $\tilde{\mathbf{q}}\in\mathcal{Q}_{*}$ that approximates the real $\mathbf{q}$ in (7). Thus, we can approximate the real $\mathbf{q}$ (up to a scaling factor) by solving the optimization problem in (8). Specifically, for $\tilde{\mathbf{q}}^{\prime}=\exp(\boldsymbol{\zeta}_{*}/\tau)$ computed from any solution $\boldsymbol{\zeta}_{*}\in\mathcal{Z}_{*}$ to (8) with finite $\ell_{\infty}$ norm, we can conclude that $\tilde{\mathbf{q}}^{\prime}$ and $\tilde{\mathbf{q}}$ differ only by a scaling factor, i.e., $\tilde{\mathbf{q}}^{\prime}=Z\tilde{\mathbf{q}}$ for some constant $0<Z<\infty$. There is no need to worry about the scaling factor $Z$ since minimizing the empirical risk $\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}},\hat{\mathbf{S}})$ and the corresponding true risk $\mathcal{L}(\mathbf{w})$ is equivalent to minimizing $\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}}^{\prime},\hat{\mathbf{S}})=
\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}},\hat{\mathbf{S}})+\tau\log(1/Z)$ and $\mathcal{L}(\mathbf{w})+\tau\log(1/Z)$.

We design a synthetic experiment to verify the effectiveness of our non-parametric method. Consider anchor data space and $\mathcal{X}=\left\{(x_{1},x_{2})\mid x_{1}^{2}+y_{1}^{2}\leq 1,x_{1}\in\left[-
1,1\right],x_{2}\in\left[0,1\right]\right\}$ and contrast data space $\mathcal{Y}=\left\{(y_{1},y_{2})\mid y_{1}\in\left[0,1\right],y_{2}\in\left[0,
1\right]\right\}$.
Let $\mathbf{x}$ be uniformly distributed on $\mathcal{X}$ and the conditional density of an $\mathbf{y}\in\mathcal{Y}$ given $\mathbf{x}$ be $p(\mathbf{y}\mid\mathbf{x})=\frac{\exp(E(\mathbf{x},\mathbf{y})/\tau)}{\int_{
\mathcal{Y}}\exp(E(\mathbf{x},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^{
\prime})}$, where $\tau=0.2$ and $E(\mathbf{x},\mathbf{y})\coloneqq\mathbf{x}^{\top}\mathbf{y}$, where the partition function $\int_{\mathcal{Y}}\exp(E(\mathbf{x},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^
{\prime})$ can be exactly computed.

Figure: Figure 2: Left: Illustration of spaces $\mathcal{X}$ and $\mathcal{Y}$; Middle: RBF interpolated heatmaps of the true $\mathbf{q}$ and approximation $\tilde{\mathbf{q}}$ when $n=100$; Right: Comparing the generalization error $|\hat{\mathcal{L}}(\tilde{\mathbf{q}},\hat{\mathbf{S}})-\mathcal{L}|$ of our method and GCL across various $n$. “MLE” refers to the MLE objective in (2) with the exact partition function.
Refer to caption: x2.png

We construct a dataset $\hat{\mathbf{S}}=\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{n}$ by uniformly sampling $\mathbf{x}_{1},\dotsc,\mathbf{x}_{n}$ from $\mathcal{X}$ and then sampling each $\mathbf{y}_{j}$ from $p_{j}=p(\cdot\mid\mathbf{x}_{j})$ by rejection sampling. The ground-truth $\mathbf{q}$ can be computed by the analytic expression of $p(\mathbf{y}\mid\mathbf{x})$. To solve the optimization problem in (8), we initialize $\boldsymbol{\zeta}_{0}=\mathbf{0}_{n}$ and obtain $\boldsymbol{\zeta}_{*}$ and $\tilde{\mathbf{q}}^{\prime}=\exp(\boldsymbol{\zeta}_{*}/\tau)$ by running gradient descent until the gradient norm $\leq 10^{-15}$. Besides, we estimate the true risk $\mathcal{L}=-\mathbb{E}_{\mathbf{x},\mathbf{y}}[\tau\log p(\mathbf{y}\mid
\mathbf{x})]$ by $-\frac{1}{N}\sum_{i=1}^{N}\tau\log p(\mathbf{y}_{i}\mid\mathbf{x}_{i})$ on $N=50,000$ sampled pairs and estimate $Z>0$ by $\frac{\|\tilde{\mathbf{q}}^{\prime}\|_{\infty}}{\|\mathbf{q}\|_{\infty}}$ to obtain $\tilde{\mathbf{q}}=\frac{\tilde{\mathbf{q}}^{\prime}}{Z}$. Then we plug $\tilde{\mathbf{q}}$ into (5) to calculate our empirical risk $\hat{\mathcal{L}}$. It is worth noting that estimating the true risk $\mathcal{L}$ and the constant $Z$ is only for the plots in Figure 2, which is not necessary for real-world empirical risk minimization problems.

As shown in the middle of Figure 2, our method effectively approximates the true $\mathbf{q}$ up to a constant $Z$. Moreover, we plot the generalization error of different methods in the right column of Figure 2, which confirms the result in Theorem 1 and Remark 4 that the uniform approximation of $\mathbf{q}$ in GCL results in a non-diminishing term in generalization error as $n$ increases. In contrast, our method achieves a significantly smaller generalization error, almost matching the MLE objective in (2) with the exact partition function when $n$ increases.

###### Theorem 2 .

###### Remark 6 .

### 3.4 The NUCLR Algorithm for Self-Supervised Representation Learning

By substituting the $\tilde{\mathbf{q}}$ from the non-parametric method described in Section 3.3, the problem of minimizing the empirical risk of DPM in (5) becomes

$$ $\displaystyle\min_{\mathbf{w}\in\mathcal{W}}\hat{\mathcal{L}}(\mathbf{w};\hat{ \mathbf{S}}),\quad\hat{\mathcal{L}}(\mathbf{w};\hat{\mathbf{S}})\coloneqq- \frac{1}{n}\sum_{i=1}^{n}\tau\log\left(\sum_{j=1}^{n}\frac{\exp(E_{\mathbf{w}} (\mathbf{x}_{i},\mathbf{y}_{i})/\tau)}{\exp((E_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y}_{j})-\zeta_{*}^{(j)})/\tau)}\right),$ $$

where $\boldsymbol{\zeta}_{*}$ is solved from (8). Unfortunately, the true similarity function $E:\mathcal{X}\times\mathcal{Y}\rightarrow[-1,1]$ in (8) is unknown. To address this, we can adopt a Jacobi-type alternating minimization approach. In the $t$-th iteration, we replace the $E(\cdot,\cdot)$ in (8) by a fixed model $E_{\mathbf{w}_{t}}(\cdot,\cdot)$ and solve $\boldsymbol{\zeta}_{t+1}$ from

$$ $\displaystyle\min_{\boldsymbol{\zeta}\in\mathbb{R}^{n}}\Phi_{t}(\boldsymbol{ \zeta}),\quad\Phi_{t}(\boldsymbol{\zeta})\coloneqq\left\{-\frac{1}{n}\sum_{i=1 }^{n}\tau\log\left(\frac{\exp(E_{\mathbf{w}_{t}}(\mathbf{x}_{i},\mathbf{y}_{i} )/\tau)}{\sum_{j=1}^{n}\exp((E_{\mathbf{w}_{t}}(\mathbf{x}_{i},\mathbf{y}_{j}) -\zeta^{(j)})/\tau)}\right)+\frac{1}{n}\sum_{j=1}^{n}\zeta^{(j)}\right\}.$ (10) $$

Then, we fix the auxiliary variable $\boldsymbol{\zeta}$ to $\boldsymbol{\zeta}_{t}$ and solve $\mathbf{w}_{t+1}$ from the following problem:

$$ $\displaystyle\min_{\mathbf{w}\in\mathcal{W}}\Psi_{t}(\mathbf{w}),\quad\Psi_{t} (\mathbf{w})\coloneqq\left\{-\frac{1}{n}\sum_{i=1}^{n}\tau\log\left(\frac{\exp (E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i})/\tau)}{\sum_{j=1}^{n}\exp((E_{ \mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j})-\zeta_{t}^{(j)})/\tau)}\right) \right\}.$ (11) $$

However, solving the subproblems (10) and (11) exactly is computationally infeasible. Instead, we obtain the new iterate $(\mathbf{w}_{t+1},\boldsymbol{\zeta}_{t+1})$ by one step of gradient-based update(^3^33E.g., the gradient descent or momentum-based update. It could be extended to multiple steps.) from the previous iterate $(\mathbf{w}_{t},\boldsymbol{\zeta}_{t})$. Defining $\phi_{i}(\mathbf{w},\boldsymbol{\zeta})\coloneqq\sum_{j\neq i}\exp((E_{\mathbf
{w}}(\mathbf{x}_{i},\mathbf{y}_{j})-E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{
i})-\zeta^{(j)})/\tau)$ and $\varepsilon(\zeta^{(i)})\coloneqq\exp(-\zeta^{(i)}/\tau)$, the objectives $\Phi_{t}(\boldsymbol{\zeta})$ and $\Psi_{t}(\mathbf{w})$ in (10) and (11) can be rewritten as follows:

$$ $\displaystyle\Phi_{t}(\boldsymbol{\zeta})=\frac{1}{n}\sum_{i=1}^{n}\tau\log( \varepsilon(\zeta^{(i)})+\phi_{i}(\mathbf{w}_{t},\boldsymbol{\zeta}))+\frac{1} {n}\sum_{j=1}^{n}\zeta^{(j)},~{}\Psi_{t}(\mathbf{w})=\frac{1}{n}\sum_{i=1}^{n} \tau\log(\varepsilon(\zeta_{t}^{(i)})+\phi_{i}(\mathbf{w},\boldsymbol{\zeta}_{ t})).$ $$

Given a size-$B$ minibatch $\mathcal{B}_{t}\subset[n]$, the stochastic gradients of $\Phi_{t}(\boldsymbol{\zeta}_{t})$, $\Psi_{t}(\mathbf{w}_{t})$ can be computed as

$$ $\displaystyle\frac{\partial}{\partial\zeta^{(j)}}\Phi_{t}(\boldsymbol{\zeta}_{ t};\mathcal{B}_{t})=\frac{1}{B}\sum_{i\in\mathcal{B}_{t}}\frac{\tau}{ \varepsilon_{t}^{(i)}+\hbox{\pagecolor{gray!30}$\phi_{i}(\mathbf{w}_{t}, \boldsymbol{\zeta}_{t};\mathcal{B}_{t})$}}\frac{\partial}{\partial\zeta^{(j)}} (\varepsilon_{t}^{(i)}+\phi_{i}(\mathbf{w}_{t};\boldsymbol{\zeta}_{t};\mathcal {B}_{t}))+\frac{1}{n},\quad j\in\mathcal{B}_{t},$ (12) $\displaystyle\nabla_{\mathbf{w}}\Psi_{t}(\mathbf{w}_{t};\mathcal{B}_{t})=\frac {1}{B}\sum_{i\in\mathcal{B}_{t}}\frac{\tau}{\varepsilon_{t}^{(i)}+\hbox{ \pagecolor{gray!30}$\phi_{i}(\mathbf{w}_{t},\boldsymbol{\zeta}_{t};\mathcal{B} _{t})$}}\nabla_{\mathbf{w}}\phi_{i}(\mathbf{w}_{t},\boldsymbol{\zeta}_{t}; \mathcal{B}_{t}),$ (13) $$

where $\varepsilon(\zeta_{t}^{(i)})$ is shortened to $\varepsilon_{t}^{(i)}$, $\phi_{i}(\mathbf{w}_{t};\boldsymbol{\zeta}_{t};\mathcal{B}_{t})\!=\!\frac{n-1}
{B-1}\sum_{j\in\mathcal{B}_{t}\backslash\{i\}}\exp(\frac{E_{\mathbf{w}_{t}}(
\mathbf{x}_{i},\mathbf{y}_{j})-E_{\mathbf{w}_{t}}(\mathbf{x}_{i},\mathbf{y}_{i
})-\zeta_{t}^{(j)})}{\tau})$ is a stochastic estimator of $\phi_{i}(\mathbf{w}_{t},\boldsymbol{\zeta}_{t})$. However, $\frac{\partial}{\partial\zeta^{(j)}}\Phi_{t}(\boldsymbol{\zeta}_{t};\mathcal{B
}_{t})$ and $\nabla_{\mathbf{w}}\Psi_{t}(\mathbf{w}_{t};\mathcal{B}_{t})$ are still *biased* estimators of $\frac{\partial}{\partial\zeta^{(j)}}\Phi_{t}(\boldsymbol{\zeta}_{t})$, $\nabla_{\mathbf{w}}\Psi_{t}(\mathbf{w}_{t})$ because the gray parts are in the denominators, where the bias is small only when the batch size $B$ is large. To resolve this issue, we adopt the moving average technique from the SogCLR algorithm (Yuan et al., ) to keep track of $\phi_{i}(\mathbf{w},\boldsymbol{\zeta})$ for each $i\in[n]$. To be specific, we maintain a scalar $u^{(i)}$ for each $i$ and update $u_{t+1}^{(i)}=(1-\gamma)u_{t}^{(i)}+\gamma\phi_{i}(\mathbf{w}_{t},\boldsymbol{
\zeta}_{t};\mathcal{B}_{t})$ for those $i\in\mathcal{B}_{t}$
while fixing $u_{t+1}^{(i)}=u_{t}^{(i)}$ for $i\notin\mathcal{B}_{t}$. Then, we obtain the following stochastic gradient estimators $G(\zeta_{t}^{(j)})$ and $G(\mathbf{w}_{t})$ by replacing the gray parts in (12) and (13) with $u_{t+1}^{(i)}$:

$$ $\displaystyle G(\zeta_{t}^{(j)})\coloneqq\frac{1}{B}\sum_{i\in\mathcal{B}_{t}} \frac{\tau}{\varepsilon_{t}^{(i)}+u_{t+1}^{(i)}}\frac{\partial}{\partial\zeta^ {(j)}}(\varepsilon_{t}^{(i)}+\phi_{i}(\mathbf{w}_{t};\boldsymbol{\zeta}; \mathcal{B}_{t}))+\frac{1}{n},\quad j\in\mathcal{B}_{t},$ (14) $\displaystyle G(\mathbf{w}_{t})\coloneqq\frac{1}{B}\sum_{i\in\mathcal{B}_{t}} \frac{\tau}{\varepsilon_{t}^{(i)}+u_{t+1}^{(i)}}\nabla_{\mathbf{w}}\phi_{i}( \mathbf{w}_{t},\boldsymbol{\zeta}_{t};\mathcal{B}_{t}),$ (15) $$

We also adopt an additional trick to improve the quality of the model $\mathbf{w}$, especially during the early training epochs.
The denominator $\sum_{j=1}^{n}\exp((E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j})-\zeta_{t}^{(
j)})/\tau)$ in the objective $\Psi_{t}(\mathbf{w})$ of  (11) can be viewed as the weighted version $\sum_{j=1}^{n}\varepsilon_{t}^{(j)}\exp((E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf
{y}_{j}))/\tau)$ of the corresponding term in GCL, where $\varepsilon_{t}^{(j)}=\exp(-\zeta_{t}^{(j)}/\tau)$ can be viewed as the “strength” of pushing $\mathbf{y}_{j}$ away from $\mathbf{x}_{i}$. For less parameter tuning, we prefer initializing $\boldsymbol{\zeta}$ from the same value $\zeta_{0}\in\mathbb{R}$. Consequently, nearly equal weights are assigned to both the positive pair $(\mathbf{x}_{i},\mathbf{y}_{i})$ and negative pairs $\{(\mathbf{x}_{i},\mathbf{y}_{j})\}_{j\neq i}$ during early training epochs, which may slow down the learning process of model $\mathbf{w}$. To address this issue, we introduce a scalar $\xi_{t}=\max\{\xi_{t-1},\left\|\boldsymbol{\zeta}_{t}\right\|_{\infty}\}$ and reduce $\varepsilon_{t}^{(i)}$ in (15) of positive pair $(\mathbf{x}_{i},\mathbf{y}_{i}))$ from $\exp(-\zeta_{t}^{(j)}/\tau)$ to $\exp(-\xi_{t}/\tau)$ to avoid aggresively pushing the positive data $\mathbf{y}_{i}$ away from $\mathbf{x}_{i}$.

Figure: Algorithm 1 NUCLR Algorithm for Self-Supervised Representation Learning

Then, we can update $\boldsymbol{\zeta}$ and $\mathbf{w}$ based on $G(\zeta_{t}^{(j)})$ and $G(\mathbf{w}_{t})$. The full update procedure is in Algorithm 1, referred to as NUCLR. The novelty of NUCLR lies in lines 5 and 10 of Algorithm 1. If we fix $\boldsymbol{\zeta}_{t}=\mathbf{0}_{n}$ and $\xi_{t}=0$, NUCLR becomes the SogCLR algorithm (Yuan et al., ). As discussed in App. F.1, NUCLR incurs only minor computational and memory overheads compared to SogCLR. Moreover, we offer an intuitive margin-based interpretation of NUCLR in App. F.2.

## 4 Experiments on Bimodal Representation Learning

Settings: In our experiments, we apply our algorithm to bimodal self-supervised representation learning on the Conceptual Captions (CC3M) (Sharma et al., ) and Conceptual 12M (CC12M) (Changpinyo et al., ) datasets. Because some data links have expired, our downloaded training set of CC3M contains $n=2,723,200$ image-text pairs, while that of CC12M contains $n=9,184,256$ image-text pairs. We evaluate the performance of trained models on downstream zero-shot image-text retrieval and image classification tasks. Retrieval performance is evaluated on the test splits of the Flickr30k (Plummer et al., ) and MSCOCO (Lin et al., ) datasets, in terms of the average Recall@1 score of image-to-text and text-to-image retrievals. The top-1 classification accuracy is evaluated on the ImageNet1k (Russakovsky et al., ) and CIFAR100 (Krizhevsky et al., ) datasets. We compare our proposed NUCLR algorithm with representative baselines CLIP (Radford et al., ), SigLIP (Zhai et al., ), DCL (Chuang et al., ), CyCLIP (Goel et al., ), and SogCLR (Yuan et al., )(^4^44For CLIP and SigLIP, we adapt implementations from the OpenCLIP repository. For DCL and CyCLIP, we use the authors’ official implementations. For SogCLR, we use the implementation in the updated codebase of Qiu et al. (), which yields better results of SogCLR than that reported in Qiu et al. ().)
. All experiments utilize distributed data-parallel (DDP) training on two NVIDIA A100 GPUs with 40GB memory and the total batch size $B$ in each iteration is 512. Besides, we use ResNet-50 as the vision encoder and DistilBert as the text encoder. The output embedding of each encoder is projected by a linear layer into a 256-dimensional feature representation for computing the losses. We run each algorithm 3 times with different random seeds and each run contains 30 epochs. Hyperparameters of all algorithms are tuned based on the validation performance. The optimizer for the model parameter $\mathbf{w}$ is AdamW (Loshchilov and Hutter, ) with a weight decay of 0.02 and a cosine learning rate schedule (Loshchilov and Hutter, ). For all algorithms, we choose a fixed temperature parameter $\tau$ tuned within $\{0.01,0.03,0.05,0.07\}$.
For SogCLR and NUCLR, we set $\gamma=0.8$ as in the SogCLR paper  (Yuan et al., ). For our NUCLR, we select $\zeta_{0}=-0.05$ on the CC3M dataset and $\zeta_{0}=0$ on the CC12M dataset. Besides, we freeze $\boldsymbol{\zeta}$ in the first 5 epochs and update $\boldsymbol{\zeta}$ by the SGDm optimizer with a cosine learning rate schedule.

**Table 1: A comparison of test performance. The best result in each column is highlighted in black.**
| Dataset | Algorithm | MSCOCO | Flickr30k | CIFAR100 | ImageNet1k | Mean |
| --- | --- | --- | --- | --- | --- | --- |
| CC3M | CLIP | 24.23 $\pm$ 0.14 | 46.33 $\pm$ 0.76 | 33.94 $\pm$ 0.87 | 35.91 $\pm$ 0.33 | 35.10 $\pm$ 0.22 |
| DCL | 24.44 $\pm$ 0.20 | 46.03 $\pm$ 0.75 | 32.78 $\pm$ 0.46 | 35.90 $\pm$ 0.20 | 34.79 $\pm$ 0.29 |  |
| SigLIP | 23.21 $\pm$ 0.14 | 44.95 $\pm$ 0.45 | 35.70 $\pm$ 0.84 | 37.53 $\pm$ 0.09 | 35.35 $\pm$ 0.31 |  |
| CyCLIP | 24.47 $\pm$ 0.25 | 47.10 $\pm$ 0.83 | 37.27 $\pm$ 0.61 | 36.63 $\pm$ 0.04 | 36.37 $\pm$ 0.42 |  |
| SogCLR | 28.54 $\pm$ 0.25 | 52.20 $\pm$ 0.64 | 35.50 $\pm$ 1.71 | 40.40 $\pm$ 0.12 | 39.16 $\pm$ 0.33 |  |
| NUCLR (Ours) | 29.55 $\pm$ 0.26 | 53.55 $\pm$ 0.22 | 37.45 $\pm$ 0.45 | 40.49 $\pm$ 0.30 | 40.26 $\pm$ 0.19 |  |
| CC12M | CLIP | 30.30 $\pm$ 0.15 | 55.21 $\pm$ 0.45 | 25.35 $\pm$ 0.64 | 44.28 $\pm$ 0.22 | 38.79 $\pm$ 0.30 |
| DCL | 30.23 $\pm$ 0.21 | 54.63 $\pm$ 0.50 | 25.55 $\pm$ 0.61 | 44.32 $\pm$ 0.07 | 38.68 $\pm$ 0.14 |  |
| SigLIP | 30.13 $\pm$ 0.45 | 55.40 $\pm$ 0.32 | 26.60 $\pm$ 1.89 | 46.12 $\pm$ 0.12 | 39.56 $\pm$ 0.68 |  |
| CyCLIP | 30.35 $\pm$ 0.24 | 54.63 $\pm$ 0.20 | 26.71 $\pm$ 2.09 | 44.94 $\pm$ 0.02 | 39.15 $\pm$ 0.50 |  |
| SogCLR | 33.91 $\pm$ 0.26 | 59.28 $\pm$ 0.07 | 26.10 $\pm$ 0.88 | 49.82 $\pm$ 0.14 | 42.28 $\pm$ 0.27 |  |
| NUCLR (Ours) | 34.36 $\pm$ 0.13 | 60.45 $\pm$ 0.03 | 28.16 $\pm$ 1.35 | 49.82 $\pm$ 0.23 | 43.20 $\pm$ 0.39 |  |

Figure: Figure 3: Validation Recall@1 performance of our algorithm and baseline SogCLR during training on the CC3M (left two columns) and CC12M datasets (right two columns).
Refer to caption: x5.png

Comparison with Baselines: We compare the testing performance of our algorithm on downstream retrieval and classification tasks with various baselines in Table 1. Moreover, we compare the Recall@1 of our NUCLR with that of SogCLR across the training epochs in Figure 3 on the validation sets of MSCOCO and Flickr30k. Compared to baselines, our NUCLR achieves superior overall performance. In particular, the improvement of our method is more significant in retrieval downstream tasks and zero-shot classification on CIFAR100, while its performance in zero-shot classification on ImageNet1k is only marginally better or comparable to SogCLR.

Ablation Study: For further understanding the key components of Algorithm 1, we compare the performance of NUCLR against three variants: (i) NUCLR-$\dagger$, which fixes $\boldsymbol{\zeta}$ to $\zeta_{0}\mathbf{1}_{n}$, i.e., removing step 5 in Algorithm 1; (ii) NUCLR-$\diamondsuit$, which does not reduce the strength $\varepsilon_{t}^{(i)}$ of pushing away the $i$-th positive pair to $\exp(-\xi_{t}/\tau)$ as described in the penultimate paragraph of Section 3.4; and (iii) NUCLR-$\clubsuit$, which does not freeze $\boldsymbol{\zeta}$ in the first 5 epochs; to answer the following questions.

(i) What is the advantage of updating $\boldsymbol{\zeta}$ using gradient-based updates compared to fixing it at $\zeta_{0}$?

(ii) Is it beneficial to reduce the strength of pushing away each positive pair by introducing $\xi$?

(iii) What is the benefit of freezing $\boldsymbol{\zeta}$ during the initial training epochs?

The results are shown in Figure 4. First, we can observe that NUCLR with the learned $\boldsymbol{\zeta}$ as in Algorithm 1 outperforms NUCLR-$\dagger$ with a fixed $\boldsymbol{\zeta}=\zeta_{0}\mathbf{1}_{n}$. Moreover, we present some examples of images from CC3M with large or small $\tilde{q}^{\prime}=\exp(\zeta/\tau)$ in Figure 5, showing that the learned $\tilde{q}^{\prime}$ effectively captures data popularity: Images depicting human life, such as portraits, pets, and landscapes, tend to be more popular than abstract symbols and shapes (see App. G.2.4 for more examples). Second, introducing $\xi$ to avoid pushing away positive pairs improves performance in both tasks compared to NUCLR-$\diamondsuit$. Lastly, freezing $\boldsymbol{\zeta}$ during the first 5 epochs yields better performance. This improvement may arise because model $\mathbf{w}_{t}$ is far from $\mathbf{w}_{*}$ in the early stage of training, making the $\boldsymbol{\zeta}$ solved from (10) less accurate. We provide more experimental results in App. G.2.

Figure: Figure 4: Compare the NUCLR algorithm with variants NUCLR-$\dagger$, NUCLR-$\diamondsuit$, and NUCLR-$\clubsuit$. “Downstream Retrieval” refers to the average test recall@1 on MSCOCO and Flickr30k datasets; “Downstream Classification” refers to the average test top-1 accuracy on CIFAR100 and ImageNet1k datasets.
Refer to caption: x9.png

Figure: Figure 5: Examples of CC3M images with large (in red) and small learned popularities $\tilde{q}^{\prime}$ (in blue).
Refer to caption: x13.png

## 5 Conclusion

In this paper, we tackle the discriminative probabilistic modeling problem on a continuous space by leveraging the multiple importance sampling (MIS) method and proposing a novel nonparametric method for approximating the sum of conditional probability densities required by MIS, which yields a new contrastive loss
for self-supervised representation learning. Then, we design an efficient algorithm called NUCLR to optimize our new loss. Experimental results on bimodal pretraining demonstrate the improved overall performance of our method compared to baseline approaches on downstream tasks. An open question remains whether we can learn a generative model in this discriminative framework that can conditionally generate high-quality new data. Additionally, our NUCLR algorithm requires storing additional $2n$ floating-point numbers for language-image pretraining. Another open question is how to learn a neural network to predict the popularities.

## Appendix A Other Related Work

In this section, we discuss previous research on discriminative and generative probabilistic modeling and theoretical analyses of self-supervised learning that are relevant to our work.

### A.1 Discriminative and Generative Probabilistic Models

Over the past decades, numerous probabilistic models have been developed for machine learning and pattern recognition problems [Grenander, , Wu et al., ]. Our work is closely related to previous works on discriminative probabilistic models for supervised and unsupervised learning, generative probabilistic models, and the recently proposed hybrid probabilistic models.

#### A.1.1 Discriminative Probabilistic Models for Supervised Learning

For supervised learning, discriminative probabilistic models learn the *conditional* probability mass function $p(\mathbf{y}\mid\mathbf{x})$ of $\mathbf{y}$ given data $\mathbf{x}$. Since the 2010s, discriminative probabilistic models based on deep neural networks, such as ConvNets, have achieved tremendous success in supervised learning tasks [LeCun et al., ], where $\mathbf{y}$ is a class label from a finite set of categories. Despite the availability of million-scale datasets with label annotations [Deng et al., , Kemelmacher-Shlizerman et al., ], scaling of these models is hindered by the enormous cost of collecting even larger curated datasets with label annotations.

#### A.1.2 Generative Probabilistic Models

Generative probabilistic models typically focus on modeling the marginal distribution $p(\mathbf{x})$ of real data. In particular, an energy-based generative probabilistic model can be defined as $p_{\theta}(\mathbf{x})=\frac{\exp(-E_{\theta}(\mathbf{x}))}{Z_{\theta}}$, where $E_{\theta}(\mathbf{x})$ is called the energy parameterized by $\theta$ and $Z_{\theta}$ is the partition function [LeCun et al., ]. The intractable partition function is the key challenge in training an energy-based generative probabilistic model. Markov Chain Monte Carlo (MCMC) sampling is a popular approach for calculating the gradient of the model parameter of the log-likelihood [Younes, , Hinton, , Gao et al., , Du and Mordatch, , Du et al., ]. However, MCMC sampling is computationally expensive. MCMC-free approaches such as score-matching [Hyvärinen and Dayan, , Song et al., , Song and Ermon, , Song and Kingma, , Song et al., ] have been proposed to learn a model that can be used to generate the data through MCMC sampling in the inference phase. On one hand, generative modeling tasks are generally more challenging than discriminative ones. For discriminative tasks, such as classification and retrieval, it suffices to have the predicted conditional density $p(\cdot\mid\mathbf{x})$ of a relevant $\mathbf{y}$ be higher than that of an irrelevant
$\mathbf{y}$. In contrast, generative tasks require a much more accurate high-dimensional density estimation to produce high-quality new data. On the other hand, training and inference on generative models typically demand more time and computational resources compared to discriminative models with a similar number of parameters.

#### A.1.3 Hybrid Probabilistic Models

Grathwohl et al. [] propose the Joint Energy-based Model (JEM) to simultaneously perform discriminative and generative learning by maximizing the likelihood of the joint distribution $p(\mathbf{x},\mathbf{y})$ for each image-label pair $(\mathbf{x},\mathbf{y})$, which can be factorized into an MCC problem solved by minimizing the CE loss and a generative modeling problem solved by the stochastic gradient Langevin dynamic (SGLD). When label annotations are absent, Wang et al. [] and Kim and Ye [] model the joint distribution $p(\mathbf{x},\mathbf{x}^{\prime})$ of a data point $\mathbf{x}$ and its augmented copy $\mathbf{x}^{\prime}$. Specifically, Wang et al. [] directly maximize the log-likelihood of the joint distribution $p(\mathbf{x},\mathbf{x}^{\prime})$ via adversarial learning, while Kim and Ye [] follow the idea of JEM to decompose the log-likelihood $\log p(\mathbf{x},\mathbf{x}^{\prime})$ into the sum of a generative term $\log p(\mathbf{x})$ and a discriminative term $\log p(\mathbf{x}^{\prime}\mid\mathbf{x})$. In Kim and Ye [], the gradient of the generative term is computed by multi-stage SGLD, and the discriminative term is approximated by the InfoNCE loss with a finite number of negative data. More recently, Bizeul et al. [] proposed the SimVAE method to model the *joint* distribution of semantically related data (e.g., different augmented copies of an image). The SimVAE maximizes the evidence lower bound (ELBO) of the log-likelihood by introducing the implicit latent structure of SimCLR into the variational autoencoder.

### A.2 Theory of Self-Supervised Representation Learning

In this section, we review additional papers on the theoretical properties of self-supervised representation learning from broader perspectives. The pioneering work by Arora et al. [] established a generalization bound for the minibatch-based InfoNCE loss that contrasts each positive pair with $k$ negative pairs in the sampled batch. They formally proved that the unsupervised error can bound the classification error of a constructed mean classifier. Following the settings in Arora et al. [], Lei et al. [] refined the analysis and improved the generalization bound by a factor of $k$. However, both results above depend on the assumption that the negative data points are independent from the anchor data points, which does not hold for contrastive learning algorithms used in practice(^5^55In practice, the negative data points for each anchor are from positive pairs of other anchors.). Moreover, their results suffer from a limitation that the generalization error increases as the number of negatives $k$ increases, which contradicts the practice. HaoChen et al. [] proposed a spectral contrastive loss and established generalization bounds for both representation learning and the downstream classification. Based on a graph-dependent McDiarmid inequality and a U-statistics reformulation, Waida et al. [] proved the generalization bound of a kernel contrastive loss (KCL) without assuming the independence between anchors and negative data. They also provided theoretical guarantees on downstream zero-shot classification. Recently, Chen et al. [] established a generalization bound of the symmetric minibatch-based InfoNCE loss used in CLIP, under the assumption that an image and a text are conditionally independent given their common shared feature. However, the convergence rate of the generalization error is also worse when the minibatch size increases. They also proved that an approximate minimizer of the population loss leads to infinitesimal zero-shot classification error under the assumption that “good”
encoding functions exist.

Besides, Wang and Isola [] proved that the asymptotics of the InfoNCE loss can be decomposed into two terms that optimize alignment and uniformity, respectively. Tian et al. [] showed that minimizing the InfoNCE loss can maximize the mutual information between two views/modalities since the loss is a lower bound of the mutual information. Tschannen et al. [] demonstrated that the success of InfoNCE loss cannot be fully attributed to information maximization because maximizing a tighter bound of mutual
information leads to worse empirical performance. Nakada et al. [] show that the gradient descent steps of minimizing the symmetric InfoNCE loss of CLIP can be viewed as performing singular value decomposition (SVD) on a contrastive cross-covariance matrix when the representation is linear. Shwartz Ziv and LeCun [] proposed a unified framework for self-supervised learning based on information theory. They surveyed numerous existing approaches and demonstrated how their framework can encompass these approaches.

## Appendix B MIS with A General Weight Function for DPM

We consider the following MIS-based estimator with a size-$m$ sample from each distribution $p(\cdot\mid\mathbf{x}_{j})$ and a general weight function $\boldsymbol{\omega}$ for the integral $g(\mathbf{w};\mathbf{x}_{i},\mathcal{Y})=\int_{\mathcal{Y}}\exp(E_{\mathbf{w}}
(\mathbf{x}_{i},\mathbf{y})/\tau)d\mu(\mathbf{y})$:

$$ $\displaystyle\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}},\boldsymbol{ \omega})=\sum_{j=1}^{n}\frac{1}{m}\sum_{l=1}^{m}\frac{\omega^{(j)}(\mathbf{y}_ {j,l})}{p(\mathbf{y}_{j,l}\mid\mathbf{x}_{j})}\exp\left(E_{\mathbf{w}}(\mathbf {x}_{i},\mathbf{y}_{j,l})/\tau\right),\quad\hat{\mathbf{Y}}=\bigcup_{j=1}^{n} \{\mathbf{y}_{j,1},\dotsc,\mathbf{y}_{j,m}\},$ (16) $$

where $\boldsymbol{\omega}$ is a weighting function such that $\omega(\mathbf{y})$ is on a probability simplex, $\forall\mathbf{y}\in\mathcal{Y}$. We denote $\hat{\mathbf{X}}\coloneqq\{\mathbf{x}_{1},\dotsc,\mathbf{x}_{n}\}$, $\Xi_{i,j}(\boldsymbol{\omega},\mathbf{y}_{j,l})\coloneqq\frac{\omega^{(j)}(
\mathbf{y}_{j,l})}{p(\mathbf{y}_{j,l}\mid\mathbf{x}_{j})}\exp\left(E_{\mathbf{
w}}(\mathbf{x}_{i},\mathbf{y}_{j,l})/\tau\right)$. We consider the “balance heuristic” $\omega_{\text{bl}}^{(j)}(\mathbf{y})=\frac{p(\mathbf{y}\mid\mathbf{x}_{j})}{
\sum_{j^{\prime}=1}^{n}p(\mathbf{y}\mid\mathbf{x}_{j^{\prime}})}$, $\forall\mathbf{y}\in\mathcal{Y}$ and $\forall j\in[n]$ proposed in Veach and Guibas []. Proposition 1 shows the unbiasedness of the estimator in (16) and justifies why we choose the balance heuristic.

###### Proposition 1 .

For each $\boldsymbol{\omega}$, we have that $\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}},\boldsymbol{\omega})$ is an unbiased estimator of the integral $g(\mathbf{w};\mathbf{x}_{i},\mathcal{Y})$; (ii) The balance heuristic $\boldsymbol{\omega}_{\text{bl}}$ minimizes $\mathbb{E}[\sum_{j=1}^{n}\frac{1}{m}\sum_{l=1}^{m}\Xi_{i,j}(\boldsymbol{\omega
},\mathbf{y}_{j,l})^{2}\mid\hat{\mathbf{X}}]$ among all possible weighting functions for any $i$, where $\mathbb{E}[\sum_{j=1}^{n}\frac{1}{m}\sum_{l=1}^{m}\Xi_{i,j}(\boldsymbol{\omega
},\mathbf{y}_{j,l})^{2}\mid\hat{\mathbf{X}}]$
is an upper bound of the variance $\mathrm{Var}[\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}},\boldsymbol{
\omega})\mid\hat{\mathbf{X}}]$; (iii) If $\sum_{j^{\prime}=1}^{n}p(\mathbf{y}\mid\mathbf{x}_{j^{\prime}})\geq\Omega(n)$ almost surely for any $\mathbf{y}\in\mathcal{Y}$ and Assumptions 2 holds, the variance goes to zero when $n\rightarrow\infty$ or $m\rightarrow\infty$.

###### Proof.

Since for any $j\in[n]$ $\mathbf{y}_{j,1},\dotsc,\mathbf{y}_{j,m}$ are i.i.d. distributed, we have

$$ $\displaystyle\mathbb{E}\left[\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y} },\boldsymbol{\omega})\mid\hat{\mathbf{X}}\right]=\sum_{j=1}^{n}\mathbb{E} \left[\frac{\omega^{(j)}(\mathbf{y}_{j,1})}{p(\mathbf{y}_{j,1}\mid\hat{\mathbf {X}})}\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j,1})/\tau\right) \mid\hat{\mathbf{X}}\right]$ $\displaystyle=\sum_{j=1}^{n}\int_{\mathcal{Y}}\frac{\omega^{(j)}(\mathbf{y})}{ p(\mathbf{y}\mid\mathbf{x}_{j})}p(\mathbf{y}\mid\mathbf{x}_{j})\exp\left(E_{ \mathbf{w}}(\mathbf{x}_{i},\mathbf{y})/\tau\right)d\mu(\mathbf{y})\stackrel{{ \scriptstyle\star}}{{=}}\int_{\mathcal{Y}}\sum_{j=1}^{n}\omega^{(j)}(\mathbf{y })\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y})/\tau\right)d\mu(\mathbf{ y})$ $\displaystyle=\int_{\mathcal{Y}}\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y})/\tau\right)d\mu(\mathbf{y}),$ (17) $$

where $\star$ is due to Tonelli’s theorem. Since $\{\mathbf{y}_{j,l}\}_{j\in[n],l\in[m]}$ are mutually independent and for a specific $j$, $\mathbf{y}_{j,1},\dotsc,\mathbf{y}_{j,l}$ are i.i.d., the variance of the estimator in (16) can be upper bounded as

$$ $\displaystyle\text{Var}[\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}}, \boldsymbol{\omega})\mid\hat{\mathbf{X}}]=\frac{1}{m}\sum_{j=1}^{n}\mathbb{E}[ \Xi_{i,j}(\boldsymbol{\omega},\mathbf{y}_{j,1})^{2}\mid\hat{\mathbf{X}}]-\frac {1}{m}\sum_{j=1}^{n}\mathbb{E}[\Xi_{i,j}(\omega^{(j)},\mathbf{y}_{j,1})\mid \hat{\mathbf{X}}]^{2}$ (18) $\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\leq\frac{1}{m}\sum_{j=1} ^{n}\mathbb{E}[\Xi_{i,j}(\boldsymbol{\omega},\mathbf{y}_{j,1})^{2}\mid\hat{ \mathbf{X}}]=\frac{1}{m}\sum_{j=1}^{n}\int_{\mathcal{Y}}\frac{\omega^{(j)}( \mathbf{y})^{2}\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y})/\tau\right) ^{2}}{p(\mathbf{y}\mid\mathbf{x}_{j})}d\mu(\mathbf{y}).$ $$

Due to Tonelli’s theorem, we have

$$ $\sum_{j=1}^{n}\int_{\mathcal{Y}}\frac{\omega^{(j)}(\mathbf{y})^{2}\exp\left(E_ {\mathbf{w}}(\mathbf{x}_{i},\mathbf{y})/\tau\right)^{2}}{p(\mathbf{y}\mid \mathbf{x}_{j})}d\mu(\mathbf{y})=\int_{\mathcal{Y}}\sum_{j=1}^{n}\frac{\omega^ {(j)}(\mathbf{y})^{2}\exp\left(E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y})/\tau \right)^{2}}{p(\mathbf{y}\mid\mathbf{x}_{j})}d\mu(\mathbf{y}).$ $$

We can instead minimize the variance upper bound at each $\mathbf{y}$ pointwise. Then, minimizing $\sum_{j=1}^{n}\frac{\omega^{(j)}(\mathbf{y})^{2}\exp\left(E_{\mathbf{w}}(
\mathbf{x}_{i},\mathbf{y})/\tau\right)^{2}}{p(\mathbf{y}\mid\mathbf{x}_{j})}$ subject to the simplex constraint leads to $\omega_{\text{bl}}^{(j)}(\mathbf{y})=\frac{p(\mathbf{y}\mid\mathbf{x}_{j})}{
\sum_{j^{\prime}=1}^{n}p(\mathbf{y}\mid\mathbf{x}_{j^{\prime}})}$. Plugging this into (18) and using Assumption 2 and $\sum_{j^{\prime}=1}^{n}p(\mathbf{y}\mid\mathbf{x}_{j^{\prime}})\geq\Omega(n)$ a.s., we have

$$ $\displaystyle\text{Var}[\hat{g}(\mathbf{w};\mathbf{x}_{i},\hat{\mathbf{Y}}, \boldsymbol{\omega}_{\text{bl}})\mid\hat{\mathbf{X}}]\leq\frac{1}{m}\sum_{j=1} ^{n}\int_{\mathcal{Y}}\frac{p(\mathbf{y}\mid O_{j})\exp\left(E_{\mathbf{w}}( \mathbf{x}_{i},\mathbf{y})/\tau\right)^{2}}{(\sum_{j^{\prime}=1}^{n}p(\mathbf{ y}\mid\mathbf{x}_{j^{\prime}}))^{2}}d\mu(\mathbf{y})=O\left(\frac{1}{mn}\right).$ $$

∎

Interestingly, the minimizer $\boldsymbol{\omega}_{\text{bl}}$ of $\frac{1}{m}\mathbb{E}[\sum_{j=1}^{n}\sum_{j=1}^{m}\Xi_{i,j}(\boldsymbol{\omega
},\mathbf{y}_{j,l})^{2}\mid\hat{\mathbf{X}}]$ does not depend on $\mathbf{x}_{i}$. We can obtain the estimator in (4) by plugging the balance heuristic $\boldsymbol{\omega}_{\text{bl}}$ into (16) and letting $m=1$.

###### Proposition 1 .

###### Proof.

## Appendix C Performance of DPM on Downstream Zero-Shot Classification

Suppose that the true conditional density function $p(\mathbf{y}\mid\mathbf{x})$ is generated by some $\mathbf{w}_{*}\in\mathcal{W}$, i.e., $p(\mathbf{y}\mid\mathbf{x})=p_{\mathbf{w}_{*}}(\mathbf{y}\mid\mathbf{x})=\frac
{\exp(E_{\mathbf{w}_{*}}(\mathbf{x},\mathbf{y})/\tau)}{\int_{\mathcal{Y}}\exp(
E_{\mathbf{w}_{*}}(\mathbf{x},\mathbf{y}^{\prime})/\tau)d\mu(\mathbf{y}^{
\prime})}$. Then, the maximum likelihood estimator $\hat{\mathbf{w}}_{*}=\operatorname*{arg\,max}_{\mathbf{w}\in\mathcal{W}}\frac{
1}{n}\sum_{i=1}^{n}\log p_{\mathbf{w}}(\mathbf{y}_{i}\mid\mathbf{x}_{i})$ with the random sample $\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{n}$ converges in probability to $\mathbf{w}_{*}$ under some mild assumptions (see Theorem 2.1 in Newey and McFadden []).

Consider the downstream multi-class classification with $K>1$ distinct classes. The task is to predict the ground-truth label $c(\mathbf{x})\in\{1,\dotsc,K\}$ of a data $\mathbf{x}\in\mathcal{X}$. Suppose that there are $K$ subsets $\mathcal{Y}_{1},\dotsc,\mathcal{Y}_{K}$ of $\mathcal{Y}$ and any $\mathbf{y}\in\mathcal{Y}_{k}$ belongs to the $k$-th class. Moreover, the ground-truth label $c(\mathbf{x})$ of data $\mathbf{x}$ is $c(\mathbf{x})=\operatorname*{arg\,max}_{k\in[K]}\mathrm{Pr}(k\mid\mathbf{x})$. Given the model $\hat{\mathbf{w}}_{*}$ trained via MLE, the predicted label $c_{\hat{\mathbf{w}}_{*}}(\mathbf{x})$ of a data $\mathbf{x}\in\mathcal{X}$ can be obtained by the following 1-nearest neighbor (1-NN) classifier:

$$ $c_{\hat{\mathbf{w}}_{*}}(\mathbf{x})=\operatorname*{arg\,max}_{k\in[K]}E_{\hat {\mathbf{w}}_{*}}(\mathbf{x},\mathbf{y}_{k}),$ $$

where $\mathbf{y}_{k}\in\mathcal{Y}$ is an example of the $k$-th class. For instance, the example $\mathbf{y}_{k}$ of the $k$-th class of the downstream image classification could be “a photo of {class_k}” when $\mathcal{X}$ is the image domain and $\mathcal{Y}$ is the text domain [Radford et al., ]. Due to the monotonicity of the function $\exp(\cdot/\tau)$ and the expression of $p_{\mathbf{w}}$ in (1), we have $c_{\hat{\mathbf{w}}_{*}}(\mathbf{x})=\operatorname*{arg\,max}_{k\in[K]}E_{\hat
{\mathbf{w}}_{*}}(\mathbf{x},\mathbf{y}_{k})=\operatorname*{arg\,max}_{k\in[K]
}p_{\hat{\mathbf{w}}_{*}}(\mathbf{y}_{k}\mid\mathbf{x})$. As long as the probability mass $\mathrm{Pr}(k\mid\mathbf{x})$ on class $k$ is proportional to the probability density $p(\mathbf{y}_{k}\mid\mathbf{x})$ on the example $\mathbf{y}_{k}$ of class $k$, the zero-one loss $\ell_{0/1}(\mathbf{x},c(\mathbf{x});\hat{\mathbf{w}}_{*})=\mathbb{I}[c_{\hat{
\mathbf{w}}_{*}}(\mathbf{x})\neq c(\mathbf{x})]$ on the data-label pair $(\mathbf{x},c(\mathbf{x}))$ of the downstream classification approaches zero when $\hat{\mathbf{w}}_{*}\stackrel{{\scriptstyle p}}{{\rightarrow}}\mathbf{w}_{*}$.

## Appendix D Proof of Theorem 1

The structure of our proof is as follows:

- •
Section D.1 presents necessary lemmas for our generalization analysis;
- •
Section D.2 decomposes the generalization error into two parts, which are handled by Section D.3 and Section D.4, respectively;
- •
Section D.5 provides bounds for Rademacher complexities of function classes parameterized by deep neural networks.

The main theorem can be proved by combining (20), (21), (22), (23), (24), (27), (30), (31).

### D.1 Lemmas

The following two lemmas provide contraction lemmas on Rademacher complexities. Lemma 1 considers the class of real-valued functions, and Lemma 2 considers the class of vector-valued functions [Maurer, , Lei et al., ].
Let $\epsilon_{i}$ and $\epsilon_{i,j}$ be independent Rademacher variables, i.e., they take values in $\{+1,-1\}$ with the same probability.

###### Lemma 1 (Contraction Lemma, Thm 11.6 in Boucheron et al. [ 2013 ] ) .

Let $\tau:\mathbb{R}_{+}\mapsto\mathbb{R}_{+}$ be convex and nondecreasing. Suppose $\psi:\mathbb{R}\mapsto\mathbb{R}$ is contractive ($|\psi(t)-\psi(\tilde{t})|\leq G|t-\tilde{t}|$) and $\psi(0)=0$. Then for any $\widetilde{\mathcal{F}}$ we have

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}}\tau\bigg{(}\sup_{f\in \widetilde{\mathcal{F}}}\sum_{i=1}^{n}\epsilon_{i}\psi\big{(}f(x_{i})\big{)} \bigg{)}\leq\mathbb{E}_{\boldsymbol{\epsilon}}\tau\bigg{(}G\sup_{f\in \widetilde{\mathcal{F}}}\sum_{i=1}^{n}\epsilon_{i}f(x_{i})\bigg{)}.$ $$

We say that a function $\psi:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is $G$-Lipschitz continuous w.r.t. $\left\|\cdot\right\|_{2}$ if $|\psi(x)-\psi(\mathbf{x})|\leq G\left\|\mathbf{x}-\mathbf{x}^{\prime}\right\|_
{2}$ for a $G>0$ and any $\mathbf{x},\mathbf{x}^{\prime}\in\mathbb{R}^{d}$.

###### Lemma 2 .

Let $\mathcal{F}$ be a class of bounded functions $f:\mathcal{Z}\mapsto\mathbb{R}^{d}$ which contains the zero function.
Let $\tau:\mathbb{R}_{+}\to\mathbb{R}_{+}$ be a continuous, non-decreasing, and convex function.
Assume $\tilde{g}_{1},\ldots,\tilde{g}_{n}:\mathbb{R}^{d}\to\mathbb{R}$ are $G$-Lipschitz continuous w.r.t. $\|\cdot\|_{2}$ and satisfy $\tilde{g}_{i}(\boldsymbol{0})=0$.
Then

$$ $\mathbb{E}_{\boldsymbol{\epsilon}\sim\{\pm 1\}^{n}}\tau\Big{(}\sup_{f\in \mathcal{F}}\sum_{i=1}^{n}\epsilon_{i}\tilde{g}_{i}(f(\mathbf{x}_{i}))\Big{)} \leq\mathbb{E}_{\boldsymbol{\epsilon}\sim\{\pm 1\}^{nd}}\tau\Big{(}G\sqrt{2} \sup_{f\in\mathcal{F}}\sum_{i=1}^{n}\sum_{j=1}^{d}\epsilon_{i,j}f_{j}(\mathbf{ x}_{i})\Big{)}.$ (19) $$

The following lemma estimates the moment generating function of a Rademacher chaos variable of order $2$ [De la Pena and Giné, ].

###### Lemma 3 ( De la Pena and Giné 2012 ) .

Let $\epsilon_{i},i\in[n]$ be independent Rademacher variables. Let $a_{i,j}\in\mathbb{R},i,j\in[n]$.
Then for $Z=\sum_{1\leq i<j\leq n}\epsilon_{i}\epsilon_{j}a_{ij}$ we have

$$ $\mathbb{E}_{\boldsymbol{\epsilon}}\exp\Big{(}|Z|/(4es)\Big{)}\leq 2,\quad\text {where }s^{2}:=\sum_{1\leq i<j\leq n}a_{i,j}^{2}.$ $$

The following lemma is a version of Talagrand’s contraction lemma.

###### Lemma 4 (Lemma 8 in Mohri and Medina [ 2014 ] ) .

Let $\mathcal{H}$ be a hypothesis set of functions mapping $\mathcal{X}$ to $\mathbb{R}$ and $\psi$ is $G$-Lipschitz functions for some
$G>0$. Then, for any sample $S$ of $n$ points $x_{1},\dotsc,x_{n}\in\mathcal{X}$, the following inequality holds:

$$ $\displaystyle\frac{1}{n}\mathbb{E}_{\epsilon_{1:n}}\left[\sup_{h\in\mathcal{H} }\sum_{i=1}^{n}\epsilon_{i}\psi(h(x_{i}))\right]\leq\frac{G}{n}\mathbb{E}_{ \epsilon_{1:n}}\left[\sup_{h\in\mathcal{H}}\sum_{i=1}^{n}\epsilon_{i}h(x_{i}) \right].$ $$

###### Lemma 1 (Contraction Lemma, Thm 11.6 in Boucheron et al. [ 2013 ] ) .

###### Lemma 2 .

###### Lemma 3 ( De la Pena and Giné 2012 ) .

###### Lemma 4 (Lemma 8 in Mohri and Medina [ 2014 ] ) .

### D.2 Error Decomposition

Considering $\log_{e}x\leq x-1$ for any $x>0$, we have

$$ $\displaystyle\hat{\mathcal{L}}(\mathbf{w};\tilde{\mathbf{q}},\hat{\mathbf{S}}) -\mathcal{L}(\mathbf{w})$ $\displaystyle=\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x}, \mathbf{y})]-\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y} _{i})+\frac{1}{n}\sum_{i=1}^{n}\tau\log(\tilde{g}(\mathbf{w};\mathbf{x}_{i}, \hat{\mathbf{Y}}))-\mathbb{E}_{\mathbf{x}}\left[\tau\log g(\mathbf{w};\mathbf{ x},\mathcal{Y})\right]$ $\displaystyle=\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x}, \mathbf{y})]-\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y} _{i})+\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\mathbf{x}}\left[\tau\log\frac{\sum _{j=1}^{n}\frac{1}{\tilde{q}^{(j)}}\exp((E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf {y}_{j})-1)/\tau)}{\int_{\mathcal{Y}}\exp((E_{\mathbf{w}}(\mathbf{x},\mathbf{y })-1)/\tau)d\mu(\mathbf{y})}\right]$ $\displaystyle\leq\underbrace{\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}} (\mathbf{x},\mathbf{y})]-\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i },\mathbf{y}_{i})}_{\text{I}}$ $\displaystyle\quad\quad\quad\quad+\underbrace{\frac{\underline{C}}{n}\sum_{i=1 }^{n}\sum_{j=1}^{n}\frac{1}{\tilde{q}^{(j)}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{ x}_{i},\mathbf{y}_{j}))-\underline{C}\mathbb{E}_{\mathbf{x}}\left[\int_{ \mathcal{Y}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y}))d\mu(\mathbf{y}) \right]}_{\text{II}},$ (20) $$

where we define $\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y})\coloneqq\frac{E_{\mathbf{w}}(
\mathbf{x},\mathbf{y})-1}{\tau}\in[-2/\tau,0]$ such that $\exp(\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y}))\in[\exp(-2/\tau),1]$. Besides, and $\overline{C}\coloneqq\sup_{\mathbf{x}\in\mathcal{X}}\frac{\tau}{\int_{\mathcal
{Y}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y}))d\mu(\mathbf{y})}$. Due to Assumption 2, $\overline{C}\leq\frac{\tau\exp(2/\tau)}{\mu(\mathcal{Y})}<\infty$. In practice, $\overline{C}$ could be much smaller than the worst-case value $\frac{\tau\exp(2/\tau)}{\mu(\mathcal{Y})}$. Similarly, we have

$$ $\displaystyle\mathcal{L}(\mathbf{w})-\hat{\mathcal{L}}(\mathbf{w};\tilde{ \mathbf{q}},\hat{\mathbf{S}})$ $\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y}_{i})-\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x}, \mathbf{y})]$ (21) $\displaystyle\quad+\overline{C}^{\prime}\mathbb{E}_{\mathbf{x}}\left[\int_{ \mathcal{Y}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y}))d\mu(\mathbf{y}) \right]-\frac{\overline{C}^{\prime}}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{ \tilde{q}^{(j)}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j})),$ $$

where $\overline{C}^{\prime}=\frac{\tau\left\|\tilde{q}\right\|_{\infty}}{n}\exp(2/\tau)$.

### D.3 Bounding Term I

Define the function class $\mathcal{E}\coloneqq\{(\mathbf{x},\mathbf{y})\mapsto E_{\mathbf{w}}(\mathbf{x}
,\mathbf{y})\mid\mathbf{w}\in\mathcal{W}\}$. Since $(\mathbf{x}_{1},\mathbf{y}_{1}),\dotsc,(\mathbf{x}_{n},\mathbf{y}_{n})$ are i.i.d. and Assumption 1 ($E_{\mathbf{w}}(\mathbf{x},\mathbf{y})\in[-1,1]$ for any $\mathbf{w}\in\mathcal{W}$), we can apply the McDiarmid’s inequality to $\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x},\mathbf{y})]-
\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i})$ and utilize the symmeterization argument following Theorem 3.3 in Mohri et al. []. With probability at least $1-\frac{\delta}{4}$,

$$ $\displaystyle\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x}, \mathbf{y})]\leq\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf {y}_{i})+2\mathfrak{R}_{n}(\mathcal{E})+6\sqrt{\frac{\log(8/\delta)}{2n}},$ $$

where $\mathfrak{R}_{n}(\mathcal{E})\coloneqq\mathbb{E}_{\hat{\mathbf{X}},\hat{
\mathbf{Y}}}[\hat{\mathfrak{R}}^{+}_{n}(\mathcal{E})]$, $\hat{\mathfrak{R}}^{+}_{n}(\mathcal{E})\coloneqq\mathbb{E}_{\epsilon_{1:n}}
\left[\sup_{e\in\mathcal{E}}\frac{1}{n}\sum_{i=1}^{n}\epsilon_{i}E_{\mathbf{w}
}(\mathbf{x}_{i},\mathbf{y}_{i})\right]$ is the empirical Rademacher complexity of $\mathcal{E}$ on the sample $\hat{\mathbf{X}}\times\hat{\mathbf{Y}}$, and $\epsilon_{1},\dotsc,\epsilon_{n}$ are Rademacher random variables. Similarly, we can also apply McDiarmid’s inequality to $\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i})-\mathbb
{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x},\mathbf{y})]$ and then use the symmetrization argument. With probability at least $1-\frac{\delta}{4}$,

$$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y} _{i})\leq\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x},\mathbf{ y})]+2\mathfrak{R}_{n}(\mathcal{E})+6\sqrt{\frac{\log(8/\delta)}{2n}},$ $$

Thus, with probability at least $1-\frac{\delta}{2}$, we have

$$ $\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}E_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y}_{i})-\mathbb{E}_{\mathbf{x},\mathbf{y}}[E_{\mathbf{w}}(\mathbf{x}, \mathbf{y})]\right|\leq 2\mathfrak{R}_{n}(\mathcal{E})+6\sqrt{\frac{\log(8/ \delta)}{2n}}.$ (22) $$

### D.4 Bounding Term II

We decompose the term II in (20) as follows.

$$ $\displaystyle\mathrm{II}$ $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{\tilde{q}^{(j)}} \exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j}))-\mathbb{E}_{\mathbf{ x}}\left[\int_{\mathcal{Y}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y}))d \mu(\mathbf{y})\right]$ $\displaystyle=\underbrace{\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\left(\frac{1 }{\tilde{q}^{(j)}}-\frac{1}{q^{(j)}}\right)\exp(\bar{E}_{\mathbf{w}}(\mathbf{x }_{i},\mathbf{y}_{j}))}_{\mathrm{II.a}}$ $\displaystyle\quad\quad+\underbrace{\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n} \frac{1}{q^{(j)}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j}))- \mathbb{E}_{\mathbf{x}}\left[\int_{\mathcal{Y}}\exp(\bar{E}_{\mathbf{w}}( \mathbf{x},\mathbf{y}))d\mu(\mathbf{y})\right]}_{\mathrm{II.b}}.$ (23) $$

Thus, we have $|\mathrm{II}|\leq|\mathrm{II.a}|+|\mathrm{II.b}|$.
Since $\exp(\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y}))=\exp((E_{\mathbf{w}}(\mathbf
{x},\mathbf{y})-1)/\tau)\leq 1$ for any $\mathbf{x}\in\mathcal{X},\mathbf{y}\in\mathcal{Y}$, we have

$$ $\displaystyle|\mathrm{II.a}|$ $\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\left|\frac{1}{\tilde{ q}^{(j)}}-\frac{1}{q^{(j)}}\right|\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y}_{j}))\leq\sum_{j=1}^{n}\left|\frac{1}{\tilde{q}^{(j)}}-\frac{1}{q^{ (j)}}\right|.$ (24) $$

We define $\Gamma(\hat{\mathbf{X}},\hat{\mathbf{Y}})\coloneqq\sup_{\mathbf{w}}\left\{
\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{q^{(j)}}\exp(\bar{E}_{\mathbf{
w}}(\mathbf{x}_{i},\mathbf{y}_{j}))-\mathbb{E}_{\mathbf{x}}\left[\int_{
\mathcal{Y}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y}))d\mu(\mathbf{y})
\right]\right\}$.
We denote that $\hat{\mathbf{X}}_{\ell}=(\hat{\mathbf{X}}\backslash\{\mathbf{x}_{\ell}\})\cup
\{\mathbf{x}_{\ell}^{\prime}\}$, $\hat{\mathbf{Y}}_{\ell}=(\hat{\mathbf{Y}}\backslash\{\mathbf{y}_{\ell}\})\cup
\{\mathbf{y}_{\ell}^{\prime}\}$, where $(\mathbf{x}_{1}^{\prime},\mathbf{y}_{1}^{\prime}),\dotsc,(\mathbf{x}_{n}^{
\prime},\mathbf{y}_{n}^{\prime})$ are i.i.d. to $(\mathbf{x}_{1},\mathbf{y}_{1}),\dotsc,(\mathbf{x}_{n},\mathbf{y}_{n})$. We denote that $q(\mathbf{y};\hat{\mathbf{X}})\coloneqq\sum_{\mathbf{x}\in\hat{\mathbf{X}}}p(
\mathbf{y}\mid\mathbf{x})$ such that $q^{(j)}=q(\mathbf{y}_{j};\hat{\mathbf{X}})$. If $q^{(j)}=\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{j}\mid\mathbf{x}_{j^{\prime}})
\geq\Omega(n)$ almost surely, we have

$$ $\displaystyle|\Gamma(\hat{\mathbf{X}},\hat{\mathbf{Y}})-\Gamma(\hat{\mathbf{X} }_{\ell},\hat{\mathbf{Y}})|$ $\displaystyle=\left|\sup_{\mathbf{w}}\frac{1}{n}\sum_{j=1}^{n}\frac{1}{q^{(j)} }\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{\ell},\mathbf{y}_{j}))-\sup_{\mathbf{w} }\frac{1}{n}\sum_{j=1}^{n}\frac{1}{q(\mathbf{y}_{j};\hat{\mathbf{X}}_{\ell})} \exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{\ell}^{\prime},\mathbf{y}_{j}))\right| \leq O(1/n),$ $\displaystyle|\Gamma(\hat{\mathbf{X}},\hat{\mathbf{Y}})-\Gamma(\hat{\mathbf{X} },\hat{\mathbf{Y}}_{\ell})|$ $\displaystyle=\left|\sup_{\mathbf{w}}\frac{1}{n}\sum_{i=1}^{n}\frac{1}{q( \mathbf{y}_{\ell};\hat{\mathbf{X}})}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y}_{\ell}))-\sup_{\mathbf{w}}\frac{1}{n}\sum_{i=1}^{n}\frac{1}{q( \mathbf{y}_{\ell}^{\prime};\hat{\mathbf{X}})}\exp(\bar{E}_{\mathbf{w}}(\mathbf {x}_{i},\mathbf{y}_{\ell}^{\prime}))\right|\leq O(1/n).$ $$

Since $\mathbf{x}_{i}$ and $\mathbf{y}_{j}$ are mutually dependent only when $i=j$, we then apply the McDiarmid-Type inequalities for graph-dependent variables (Theorem 3.6 in Zhang et al. []) to the term II.b and $-$II.b. With probability at least $1-\frac{\delta}{4}$, $\delta\in(0,1)$, we have

$$ $\displaystyle\mathrm{II.b}\leq\mathbb{E}\left[\sup_{\mathbf{w}}\mathrm{II.b} \right]+O\left(\sqrt{\frac{10\log(4/\delta)}{n}}\right).$ (25) $$

Similarly, with probability at least $1-\frac{\delta}{4}$, $\delta\in(0,1)$, we have

$$ $\displaystyle-\mathrm{II.b}\leq\mathbb{E}\left[\sup_{\mathbf{w}}\left\{- \mathrm{II.b}\right\}\right]+O\left(\sqrt{\frac{10\log(4/\delta)}{n}}\right).$ (26) $$

Let $(\mathbf{x}_{1}^{\prime},\mathbf{y}_{1}^{\prime}),\dotsc,(\mathbf{x}_{n}^{
\prime},\mathbf{y}_{n}^{\prime})$ be a virtual sample i.i.d. to $(\mathbf{x}_{1},\mathbf{y}_{1}),\dotsc,(\mathbf{x}_{n},\mathbf{y}_{n})$. Denote that $\hat{\mathbf{X}}^{\prime}\coloneqq\{\mathbf{x}_{1}^{\prime},\dotsc,\mathbf{x}_
{n}^{\prime}\}$, $\hat{\mathbf{Y}}^{\prime}\coloneqq\{\mathbf{y}_{1}^{\prime},\dotsc,\mathbf{y}_
{n}^{\prime}\}$. Due to (17), we have

$$ $\displaystyle\mathbb{E}_{\mathbf{x}}\left[\int_{\mathcal{Y}}\exp(\bar{E}_{ \mathbf{w}}(\mathbf{x},\mathbf{y}))d\mu(\mathbf{y})\right]=\mathbb{E}_{\hat{ \mathbf{X}}^{\prime},\hat{\mathbf{Y}}^{\prime}}\left[\frac{1}{n}\sum_{i=1}^{n} \sum_{j=1}^{n}\frac{1}{q(\mathbf{y}_{j}^{\prime};\hat{\mathbf{X}}^{\prime})} \exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i}^{\prime},\mathbf{y}_{j}^{\prime})) \right].$ $$

We can rewrite and decompose the $\mathbb{E}\left[\sup_{\mathbf{w}}\mathrm{II.b}\right]$ term as

$$ $\displaystyle\mathbb{E}\left[\sup_{\mathbf{w}}\mathrm{II.b}\right]=\mathbb{E} \left[\sup_{\mathbf{w}}\left\{\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{ q^{(j)}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j}))-\mathbb{E}_{ \mathbf{x}}\left[\int_{\mathcal{Y}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}, \mathbf{y}))d\mu(\mathbf{y})\right]\right\}\right]$ $\displaystyle=\mathbb{E}\left[\sup_{\mathbf{w}}\left\{\frac{1}{n}\sum_{i=1}^{n }\sum_{j=1}^{n}\frac{1}{q^{(j)}}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i}, \mathbf{y}_{j}))-\mathbb{E}_{\hat{\mathbf{X}}^{\prime},\hat{\mathbf{Y}}^{ \prime}}\left[\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{q(\mathbf{y}_{j} ^{\prime};\hat{\mathbf{X}}^{\prime})}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{i}^ {\prime},\mathbf{y}_{j}^{\prime}))\right]\right\}\right]$ $\displaystyle\leq\mathbb{E}_{\hat{\mathbf{X}},\hat{\mathbf{Y}},\hat{\mathbf{X} }^{\prime},\hat{\mathbf{Y}}^{\prime}}\left[\sup_{\mathbf{w}}\left\{\frac{1}{n} \sum_{i=1}^{n}\frac{1}{q(\mathbf{y}_{i};\hat{\mathbf{X}})}\exp(\bar{E}_{ \mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{i}))-\frac{1}{n}\sum_{i=1}^{n}\frac{1}{ q(\mathbf{y}_{i}^{\prime};\hat{\mathbf{X}}^{\prime})}\exp(\bar{E}_{\mathbf{w}} (\mathbf{x}_{i}^{\prime},\mathbf{y}_{i}^{\prime}))\right\}\right]$ $\displaystyle\quad+\mathbb{E}_{\hat{\mathbf{X}},\hat{\mathbf{Y}},\hat{\mathbf{ X}}^{\prime},\hat{\mathbf{Y}}^{\prime}}\left[\sup_{\mathbf{w}}\left\{\frac{1}{ n}\sum_{i=1}^{n}\sum_{j\neq i}\frac{1}{q^{(j)}}\exp(\bar{E}_{\mathbf{w}}( \mathbf{x}_{i},\mathbf{y}_{j}))-\frac{1}{n}\sum_{i=1}^{n}\sum_{j\neq i}\frac{1 }{q(\mathbf{y}_{j}^{\prime};\hat{\mathbf{X}}^{\prime})}\exp(\bar{E}_{\mathbf{w }}(\mathbf{x}_{i}^{\prime},\mathbf{y}_{j}^{\prime}))\right\}\right]$ $\displaystyle\leq O(1/n)+\mathbb{E}\left[\sup_{\mathbf{w}}\left\{\frac{1}{n} \sum_{i=1}^{n}\sum_{j\neq i}\frac{1}{q^{(j)}}\exp(\bar{E}_{\mathbf{w}}(\mathbf {x}_{i},\mathbf{y}_{j}))-\frac{1}{n}\sum_{i=1}^{n}\sum_{j\neq i}\frac{1}{q( \mathbf{y}_{j}^{\prime};\hat{\mathbf{X}}^{\prime})}\exp(\bar{E}_{\mathbf{w}}( \mathbf{x}_{i}^{\prime},\mathbf{y}_{j}^{\prime}))\right\}\right],$ $$

the last step is due to the assumption $q(\mathbf{y}_{i};\hat{\mathbf{X}})=\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{i}\mid
\mathbf{x}_{j^{\prime}})\geq\Omega(n)$. Next, we adapt the proof technique in Theorem 6 of Waida et al. []. W.l.o.g., we assume that $n$ is even (If $n$ is odd, we can apply the following analysis to the first $n-1$ terms in the summation, where $n-1$ is even. The last term in the summation is a $O(1/n)$ term, which does not change the result). Suppose that $S_{n}$ is the set of all permutations (the symmetric group of degree $n$). Then, for each $s\in S$, pairs $(\mathbf{x}_{s(2i-1)}),\mathbf{y}_{s(2i)})$ ($i=1,\dotsc,n/2$) are mutually independent. Consider the alternative expression of a U-statistics of order 2  (See Appendix 1 in Clémencon et al. []):

$$ $\displaystyle\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j\neq i}\frac{1}{q^{(j)}}\exp (\bar{E}_{\mathbf{w}}(\mathbf{x}_{i},\mathbf{y}_{j}))=\frac{1}{n!(n/2)}\sum_{s \in S_{n}}\sum_{i=1}^{n/2}\frac{1}{q(\mathbf{y}_{s(2i)};\hat{\mathbf{X}})}\exp (\bar{E}_{\mathbf{w}}(\mathbf{x}_{s(2i-1)},\mathbf{y}_{s(2i)})).$ $$

It then follows that

$$ $\displaystyle\mathbb{E}\left[\sup_{\mathbf{w}}\mathrm{II.b}\right]\leq O(1/n)+ \frac{n\!-\!1}{n/2}\mathbb{E}\left[\sup_{\mathbf{w}}\frac{1}{n!}\sum_{s\in S_{ n}}\sum_{i=1}^{n/2}\left(\frac{\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{s(2i-1)}, \mathbf{y}_{s(2i)}))}{q(\mathbf{y}_{s(2i)};\hat{\mathbf{X}})}-\frac{\exp(\bar{ E}_{\mathbf{w}}(\mathbf{x}^{\prime}_{s(2i-1)},\mathbf{y}^{\prime}_{s(2i)}))}{q (\mathbf{y}_{s(2i)}^{\prime};\hat{\mathbf{X}}^{\prime})}\right)\right]$ $\displaystyle\leq O(1/n)+\frac{n\!-\!1}{n/2}\frac{1}{n!}\sum_{s\in S_{n}} \mathbb{E}\left[\sup_{\mathbf{w}}\sum_{i=1}^{n/2}\left(\frac{\exp(\bar{E}_{ \mathbf{w}}(\mathbf{x}_{s(2i-1)},\mathbf{y}_{s(2i)}))}{q(\mathbf{y}_{s(2i)}; \hat{\mathbf{X}})}-\frac{\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}^{\prime}_{s(2i-1 )},\mathbf{y}^{\prime}_{s(2i)}))}{q(\mathbf{y}_{s(2i)}^{\prime};\hat{\mathbf{X }}^{\prime})}\right)\right]$ $\displaystyle=O(1/n)+\frac{n\!-\!1}{n/2}\mathbb{E}\left[\sup_{\mathbf{w}}\sum_ {i=1}^{n/2}\left(\frac{\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{2i-1},\mathbf{y}_ {2i}))}{q(\mathbf{y}_{2i};\hat{\mathbf{X}})}-\frac{\exp(\bar{E}_{\mathbf{w}}( \mathbf{x}^{\prime}_{2i-1},\mathbf{y}^{\prime}_{2i}))}{q(\mathbf{y}_{2i}^{ \prime};\hat{\mathbf{X}}^{\prime})}\right)\right]$ $\displaystyle=O(1/n)+\frac{n\!-\!1}{n/2}\mathbb{E}\left[\sup_{\mathbf{w}}\sum_ {i=1}^{n/2}\epsilon_{i}\left(\frac{\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{2i-1} ,\mathbf{y}_{2i}))}{q(\mathbf{y}_{2i};\hat{\mathbf{X}})}-\frac{\exp(\bar{E}_{ \mathbf{w}}(\mathbf{x}^{\prime}_{2i-1},\mathbf{y}^{\prime}_{2i}))}{q(\mathbf{y }_{2i}^{\prime};\hat{\mathbf{X}}^{\prime})}\right)\right]$ $\displaystyle\leq O(1/n)+\frac{2(n-1)}{n/2}\mathbb{E}\left[\sup_{\mathbf{w}} \sum_{i=1}^{n/2}\frac{\epsilon_{i}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{2i-1}, \mathbf{y}_{2i}))}{q(\mathbf{y}_{2i};\hat{\mathbf{X}})}\right],$ $$

where we have used the symmetry between the permutations in $S_{n}$ and $(\mathbf{x}_{i},\mathbf{y}_{i}),(\mathbf{x}_{i}^{\prime},\mathbf{y}_{i}^{
\prime})$. By Lemma 4 and the assumption $q(\mathbf{y}_{2i};\hat{\mathbf{X}})=\sum_{j^{\prime}=1}^{n}p(\mathbf{y}_{2i}
\mid\mathbf{x}_{j^{\prime}})\geq\Omega(n)$, we further get

$$ $\mathbb{E}\left[\sup_{\mathbf{w}}\mathrm{II.b}\right]\leq O(1/n)+O(1/n)\mathbb {E}\left[\sup_{\mathbf{w}}\sum_{i=1}^{n/2}\epsilon_{i}\exp(\bar{E}_{\mathbf{w} }(\mathbf{x}_{2i-1},\mathbf{y}_{2i}))\right].$ $$

Define the function class $\bar{\mathcal{G}}=\{(\mathbf{x},\mathbf{y})\mapsto\exp(\bar{E}_{\mathbf{w}}(
\mathbf{x},\mathbf{y}))\mid\mathbf{w}\in\mathcal{W}\}$. Then, we define the following empirical Rademacher complexity

$$ $\hat{\mathfrak{R}}^{-}_{n/2}(\bar{\mathcal{G}};s)\coloneq\frac{2}{n}\mathbb{E} _{\epsilon_{1:n/2}}\left[\sup_{\mathbf{w}}\sum_{i=1}^{n/2}\epsilon_{i}\exp( \bar{E}_{\mathbf{w}}(\mathbf{x}_{s(2i-1)},\mathbf{y}_{s(2i)}))\right].$ $$

We further define the Rademacher complexity $\mathfrak{R}^{-}_{n/2}(\bar{\mathcal{G}})\coloneqq\max_{s\in S_{n}}\mathbb{E}_
{\hat{\mathbf{X}},\hat{\mathbf{Y}}}[\hat{\mathfrak{R}}^{-}_{n/2}(\bar{\mathcal
{G}};s)]$. We can also apply the symmetrization argument above to bound $\mathbb{E}[\sup_{\mathbf{w}}\{-\mathrm{II.b}\}]$. Due to Assumption 1, we can bound the II.b term as: With probability $1-\frac{\delta}{2}$, $\delta\in(0,1)$, we have

$$ $\displaystyle|\text{II.b}|$ $\displaystyle\leq O(1)\hat{\mathfrak{R}}^{-}_{n/2}(\bar{\mathcal{G}};s)+O\left (\frac{1}{n}+\sqrt{\frac{10\log(4/\delta)}{n}}\right).$ (27) $$

### D.5 Bounding Rademacher Complexities

We consider the specific similarity function:

$$ $\displaystyle E_{\mathbf{w}}(\mathbf{x},\mathbf{y})=E_{1}(\mathbf{w}_{1}; \mathbf{x})^{\top}E_{2}(\mathbf{w}_{2};\mathbf{y}).$ $$

We consider $L$-layer neural networks

$$ $\displaystyle E_{1}(\mathbf{w}_{1};\mathbf{x})\in\mathcal{F}_{1,L}=\{\mathbf{x }\rightarrow\sigma(W_{1,L}\sigma(W_{1,L-1}\dotsc\sigma(W_{1,1}\mathbf{x}))): \left\|W_{1,l}\right\|_{F}\leq B_{l}\},$ $\displaystyle E_{2}(\mathbf{w}_{2};\mathbf{y})\in\mathcal{F}_{2,L}=\{\mathbf{y }\rightarrow\sigma(W_{2,L}\sigma(W_{2,L-1}\dotsc\sigma(W_{2,1}\mathbf{y}))): \left\|W_{2,l}\right\|_{F}\leq B_{l}\}.$ $$

Suppose that $W_{1,l}\in\mathbb{R}^{d_{1,l}\times d_{1,l-1}}$, $W_{2,l}\in\mathbb{R}^{d_{2,l}\times d_{2,l-1}}$ and $d_{1,0}=d_{1}$, $d_{2,0}=d_{2}$, $d_{1,L}=d_{2,L}=d_{L}$. Define $W_{l}^{\top}=(W_{l}^{(1)},\dotsc,W_{1}^{(d_{l})})$, where $W_{l}^{(\iota)}$ is the $\iota$-th row of matrix $W_{l}$. The following results are adaptions of the results in Golowich et al. [].

#### D.5.1 Bounding ℜ n ⁢ ( ℰ ) subscript ℜ 𝑛 ℰ \mathfrak{R}_{n}(\mathcal{E}) fraktur_R start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( caligraphic_E )

Define $h:\mathbb{R}^{2d}\rightarrow\mathbb{R}$ as $h(\mathbf{v})=\mathbf{v}_{1}^{\top}\mathbf{v}_{2}$, where $\mathbf{v}=\begin{pmatrix}\mathbf{v}_{1}\\
\mathbf{v}_{2}\end{pmatrix}$ and $\mathbf{v}_{1},\mathbf{v}_{2}\in\mathbb{R}^{d}$. It is clear that $E_{\mathbf{w}}(\mathbf{x},\mathbf{y})=h(E_{1}(\mathbf{w}_{1};\mathbf{x}),E_{
\mathbf{w}}(\mathbf{w}_{2};\mathbf{y}))$. Due to Assumption 2, we have $\left\|E_{1}(\mathbf{w}_{1};\mathbf{x})\right\|_{2}\leq 1$ and $\left\|E_{2}(\mathbf{w}_{2};\mathbf{y})\right\|_{2}\leq 1$. For any $\mathbf{v}=\begin{pmatrix}\mathbf{v}_{1}\\
\mathbf{v}_{2}\end{pmatrix},\mathbf{v}^{\prime}=\begin{pmatrix}\mathbf{v}_{1}^
{\prime}\\
\mathbf{v}_{2}^{\prime}\end{pmatrix}$ and $\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{v}_{1}^{\prime},\mathbf{v}_{2}^{\prime}
\in[0,1]^{d}$, we have

$$ $\displaystyle\left(h(\mathbf{v})-h(\mathbf{v}^{\prime})\right)^{2}$ $\displaystyle\leq 2(\mathbf{v}_{1}^{\top}(\mathbf{v}_{2}-\mathbf{v}_{2}^{ \prime}))^{2}+2((\mathbf{v}_{1}-\mathbf{v}_{1}^{\prime})^{\top}\mathbf{v}_{2}^ {\prime})^{2}\leq 2\left\|\mathbf{v}-\mathbf{v}^{\prime}\right\|_{2}^{2},$ $$

where we have used $(a+b)^{2}\leq 2a^{2}+2b^{2}$ and the decomposition $\mathbf{v}_{1}^{\top}\mathbf{v}_{2}-(\mathbf{v}_{1}^{\prime})^{\top}\mathbf{v}
_{2}^{\prime}=\mathbf{v}_{1}^{\top}(\mathbf{v}_{2}-\mathbf{v}_{2}^{\prime})+(
\mathbf{v}_{1}-\mathbf{v}_{1}^{\prime})^{\top}\mathbf{v}_{2}^{\prime}$. Thus, we can conclude that $h$ is 1-Lipschitz continuous to $\mathbf{v}$ and apply Lemma 2 to the function $E_{\mathbf{w}}(\mathbf{x},\mathbf{y})=h(E_{1}(\mathbf{w}_{1};\mathbf{x}),E_{2}
(\mathbf{w}_{2};\mathbf{y}))$:

$$ $\displaystyle\hat{\mathfrak{R}}^{+}_{n}(\mathcal{E})$ $\displaystyle=\mathbb{E}_{\epsilon_{1:n}}\left[\sup_{e\in\mathcal{E}}\frac{1}{ n}\sum_{i=1}^{n}\epsilon_{i}E(\mathbf{x}_{i},\mathbf{y}_{i})\right]$ $\displaystyle\leq\frac{1}{n}\mathbb{E}_{\boldsymbol{\epsilon}_{1},\boldsymbol{ \epsilon}_{2}\in\{\pm 1\}^{nd_{L}}}\left[\sup_{\mathbf{w}}\sum_{i=1}^{n}\sum_{ \iota=1}^{d_{L}}\left(\epsilon_{1}^{(i,\iota)}E_{1}^{(\iota)}(\mathbf{w}_{1}, \mathbf{x}_{i})+\epsilon_{2}^{(i,\iota)}E_{2}^{(\iota)}(\mathbf{w}_{2},\mathbf {y}_{i})\right)\right]$ $\displaystyle\leq\frac{1}{n}\mathbb{E}_{\boldsymbol{\epsilon}_{1}\in\{\pm 1\}^ {nd_{L}}}\left[\sup_{\mathbf{w}}\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\epsilon_{ 1}^{(i,\iota)}E_{1}^{(\iota)}(\mathbf{w}_{1},\mathbf{x}_{i})\right]+\frac{1}{n }\mathbb{E}_{\boldsymbol{\epsilon}_{1},\boldsymbol{\epsilon}_{2}\in\{\pm 1\}^{ nd_{L}}}\left[\sup_{\mathbf{w}}\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\epsilon_{2 }^{(i,\iota)}E_{2}^{(\iota)}(\mathbf{w}_{2},\mathbf{y}_{i})\right]$ $\displaystyle=\frac{1}{n}\mathbb{E}_{\boldsymbol{\epsilon}_{1}\in\{\pm 1\}^{nd _{L}}}\left[\sup_{W_{1,L},f_{1,L-1}\in\mathcal{F}_{1,L-1}}\sum_{i=1}^{n}\sum_{ \iota=1}^{d_{L}}\epsilon_{1}^{(i,\iota)}\sigma(f_{1,L-1}(\mathbf{x}_{i})^{\top }W_{1,L}^{(\iota)})\right]$ $\displaystyle\quad\quad+\frac{1}{n}\mathbb{E}_{\boldsymbol{\epsilon}_{2}\in\{ \pm 1\}^{nd_{L}}}\left[\sup_{W_{2,L},f_{2,L-1}\in\mathcal{F}_{2,L-1}}\sum_{i=1 }^{n}\sum_{\iota=1}^{d_{L}}\epsilon_{2}^{(i,\iota)}\sigma(f_{2,L-1}(\mathbf{y} _{i})^{\top}W_{2,L}^{(\iota)})\right].$ $$

For simplicity, we can only consider one of the terms above and neglect the index of embedding networks (1 or 2). Let $\mathbf{x}_{i}$ be one of $\mathbf{x}_{i}$ and $\mathbf{y}_{i}$. Cauchy-Schwarz and $(\sup x)^{2}\leq\sup x^{2}$ imply

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}\in\{\pm 1\}^{nd_{L}}}\left[\sup _{W_{L},f\in\mathcal{F}_{L-1}}\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{ \epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota)})\right]$ $\displaystyle\leq\left(\mathbb{E}_{\boldsymbol{\epsilon}\in\{\pm 1\}^{nd_{L}}} \left[\left(\sup_{W_{L},f\in\mathcal{F}_{L-1}}\sum_{i=1}^{n}\sum_{\iota=1}^{d_ {L}}\boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{( \iota)})\right)^{2}\right]\right)^{\frac{1}{2}}$ $\displaystyle\leq\left(\mathbb{E}_{\boldsymbol{\epsilon}\in\{\pm 1\}^{nd_{L}}} \left[\sup_{W_{L},f\in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_ {L}}\boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{( \iota)})\right)^{2}\right]\right)^{\frac{1}{2}}.$ (28) $$

For a $\lambda>0$, Jensen’s inequality implies that

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}\in\{\pm 1\}^{nd_{L}}}\left[\sup _{W_{L},f\in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}} \boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota) })\right)^{2}\right]$ $\displaystyle=\frac{1}{\lambda}\log\exp\left(\lambda\mathbb{E}_{\boldsymbol{ \epsilon}}\left[\sup_{W_{L},f\in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{ \iota=1}^{d_{L}}\boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{ \top}W_{L}^{(\iota)})\right)^{2}\right]\right)$ $\displaystyle\leq\frac{1}{\lambda}\log\left(\mathbb{E}_{\boldsymbol{\epsilon}} \exp\left(\lambda\sup_{W_{L},f\in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{ \iota=1}^{d_{L}}\boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{ \top}W_{L}^{(\iota)})\right)^{2}\right)\right).$ (29) $$

We utilize the following facts: (i) $\sup_{x}x^{2}\leq\max\{(\sup_{x}x)^{2},(\sup_{x}(-x))^{2}\}$ and for a Rademacher random variable $\epsilon$, we have $\epsilon,-\epsilon$ are i.i.d.; (ii) Lemma 1 with $\tau(t)=\exp(\lambda t^{2})$ and $\sigma$ is 1-Lipschitz; (iii) $(\sup x)^{2}\leq\sup x^{2}$:

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_{L},f \in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{ \epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota)})\right)^{2 }\right)$ $\displaystyle\stackrel{{\scriptstyle\text{(i)}}}{{\leq}}2\mathbb{E}_{ \boldsymbol{\epsilon}}\exp\left(\lambda\left(\sup_{W_{L},f\in\mathcal{F}_{L-1} }\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{\epsilon}^{(i,\iota)}\sigma(f (\mathbf{x}_{i})^{\top}W_{L}^{(\iota)})\right)^{2}\right)$ $\displaystyle\stackrel{{\scriptstyle\text{(ii)}}}{{\leq}}2\mathbb{E}_{ \boldsymbol{\epsilon}}\exp\left(\lambda\left(\sup_{W_{L},f\in\mathcal{F}_{L-1} }\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{\epsilon}^{(i,\iota)}f( \mathbf{x}_{i})^{\top}W_{L}^{(\iota)}\right)^{2}\right)$ $\displaystyle\stackrel{{\scriptstyle\text{(iii)}}}{{\leq}}2\mathbb{E}_{ \boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_{L},f\in\mathcal{F}_{L-1}}\left (\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{\epsilon}^{(i,\iota)}f( \mathbf{x}_{i})^{\top}W_{L}^{(\iota)}\right)^{2}\right).$ $$

Due to (iv) $\left\|W_{l}\right\|_{F}\leq B_{l}$ for each $l\in[L]$, we further have

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_{L},f \in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{ \epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota)})\right)^{2 }\right)$ $\displaystyle\leq 2\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_ {L},f\in\mathcal{F}_{L-1}}\left(\sum_{\iota=1}^{d_{L}}\left\|\sum_{i=1}^{n} \boldsymbol{\epsilon}^{(i,\iota)}f(\mathbf{x}_{i})\right\|_{2}\left\|W_{L}^{( \iota)}\right\|_{2}\right)^{2}\right)$ $\displaystyle\leq 2\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_ {L},f\in\mathcal{F}_{L-1}}\left\|W_{L}\right\|_{F}^{2}\sum_{\iota=1}^{d_{L}} \left\|\sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i,\iota)}f(\mathbf{x}_{i})\right \|_{2}^{2}\right)$ $\displaystyle\stackrel{{\scriptstyle\text{(iv)}}}{{\leq}}2\mathbb{E}_{ \boldsymbol{\epsilon}}\exp\left(\lambda B_{L}^{2}\sup_{f\in\mathcal{F}_{L-1}} \sum_{\iota=1}^{d_{L}}\left\|\sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i,\iota)}f( \mathbf{x}_{i})\right\|_{2}^{2}\right)$ $\displaystyle=2\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda B_{L}^{2} \sup_{W_{L-1},f\in\mathcal{F}_{L-2}}\sum_{\iota=1}^{d_{L}}\left\|\sum_{i=1}^{n }\boldsymbol{\epsilon}^{(i,\iota)}\sigma(W_{L-1}f(\mathbf{x}_{i}))\right\|_{2} ^{2}\right).$ $$

Due to the positive-homogeneous property of the activation function $\sigma(\cdot)$, we have

$$ $\displaystyle\sum_{\iota=1}^{d_{L}}\left\|\sum_{i=1}^{n}\boldsymbol{\epsilon}^ {(i,\iota)}\sigma(W_{L-1}f(\mathbf{x}_{i}))\right\|_{2}^{2}=\sum_{\iota=1}^{d_ {L}}\left\|\begin{pmatrix}\sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i,\iota)} \sigma(f(\mathbf{x}_{i})^{\top}W_{L-1}^{(1)})\\ \vdots\\ \sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top} W_{L-1}^{(d_{L-1})})\end{pmatrix}\right\|_{2}^{2}$ $\displaystyle=\sum_{\iota=1}^{d_{L}}\sum_{r=1}^{d_{L-1}}\left(\sum_{i=1}^{n} \boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L-1}^{(r)}) \right)^{2}=\sum_{r=1}^{d_{L-1}}\left\|W_{L-1}^{(r)}\right\|_{2}^{2}\sum_{ \iota=1}^{d_{L}}\left(\sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i,\iota)}\sigma \left(f(\mathbf{x}_{i})^{\top}\frac{W_{L-1}^{(r)}}{\left\|W_{L-1}^{(r)}\right \|_{2}}\right)\right)^{2}$ $\displaystyle\leq\left\|W_{L-1}\right\|_{F}^{2}\max_{r\in[d_{L-1}]}\sum_{\iota =1}^{d_{L}}\left(\sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i,\iota)}\sigma\left(f( \mathbf{x}_{i})^{\top}\frac{W_{L-1}^{(r)}}{\left\|W_{L-1}^{(r)}\right\|_{2}} \right)\right)^{2}$ $\displaystyle\leq B_{L-1}^{2}\sup_{\mathbf{w}:\left\|\mathbf{w}\right\|_{2} \leq 1}\sum_{\iota=1}^{d_{L}}\left(\sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i, \iota)}\sigma\left(f(\mathbf{x}_{i})^{\top}\mathbf{w}\right)\right)^{2}.$ $$

Thus, we can obtain

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_{L},f \in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{ \epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota)})\right)^{2 }\right)$ $\displaystyle\leq 2\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda B_{L}^{ 2}B_{L-1}^{2}\sup_{\left\|\mathbf{w}\right\|_{2}\leq 1,f\in\mathcal{F}_{L-2}} \sum_{\iota=1}^{d_{L}}\left(\sum_{i=1}^{n}\boldsymbol{\epsilon}^{(i,\iota)} \sigma(f(\mathbf{x}_{i})^{\top}\mathbf{w})\right)^{2}\right)$ $\displaystyle\leq 2\mathbb{E}_{\epsilon_{1:n}}\exp\left(d_{L}\lambda B_{L}^{2} B_{L-1}^{2}\sup_{\left\|\mathbf{w}\right\|_{2}\leq 1,f\in\mathcal{F}_{L-2}} \left(\sum_{i=1}^{n}\epsilon_{i}\sigma(f(\mathbf{x}_{i})^{\top}\mathbf{w}) \right)^{2}\right).$ $$

Applying Lemma 1 with $\tau_{\lambda}(t)=\exp(d_{L}\lambda B_{L}^{2}B_{L-1}^{2}t^{2})$ gives

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_{L},f \in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{ \epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota)})\right)^{2 }\right)\leq 2\mathbb{E}_{\epsilon_{1:n}}\left[\tau_{\lambda}\left(\sup_{\left \|\mathbf{w}\right\|_{2}\leq 1,f\in\mathcal{F}_{L-2}}\left|\sum_{i=1}^{n} \epsilon_{i}\sigma(f(\mathbf{x}_{i})^{\top}\mathbf{w})\right|\right)\right]$ $\displaystyle\leq 2\mathbb{E}_{\epsilon_{1:n}}\left[\tau_{\lambda}\left(\sup_{ \left\|\mathbf{w}\right\|_{2}\leq 1,f\in\mathcal{F}_{L-2}}\sum_{i=1}^{n} \epsilon_{i}\sigma(f(\mathbf{x}_{i})^{\top}\mathbf{w})\right)\right]+2\mathbb{ E}_{\epsilon_{1:n}}\left[\tau_{\lambda}\left(\sup_{\left\|\mathbf{w}\right\|_{ 2}\leq 1,f\in\mathcal{F}_{L-2}}-\sum_{i=1}^{n}\epsilon_{i}\sigma(f(\mathbf{x}_ {i})^{\top}\mathbf{w})\right)\right]$ $\displaystyle=4\mathbb{E}_{\epsilon_{1:n}}\left[\tau_{\lambda}\left(\sup_{ \left\|\mathbf{w}\right\|_{2}\leq 1,f\in\mathcal{F}_{L-2}}\sum_{i=1}^{n} \epsilon_{i}\sigma(f(\mathbf{x}_{i})^{\top}\mathbf{w})\right)\right]\leq 4 \mathbb{E}_{\epsilon_{1:n}}\left[\tau_{\lambda}\left(\sup_{\left\|\mathbf{w} \right\|_{2}\leq 1,f\in\mathcal{F}_{L-2}}\sum_{i=1}^{n}\epsilon_{i}f(\mathbf{x }_{i})^{\top}\mathbf{w}\right)\right]$ $\displaystyle\leq 4\mathbb{E}_{\epsilon_{1:n}}\left[\tau_{\lambda}\left(\sup_{ W_{L-2},f\in\mathcal{F}_{L-3}}\left\|\sum_{i=1}^{n}\epsilon_{i}\sigma(W_{L-2}f (\mathbf{x}_{i}))\right\|_{2}\right)\right]\leq 4\mathbb{E}_{\epsilon_{1:n}} \left[\tau_{\lambda}\left(B_{L-2}\sup_{\left\|\mathbf{w}\right\|_{2}\leq 1,f \in\mathcal{F}_{L-3}}\left|\sum_{i=1}^{n}\epsilon_{i}f(\mathbf{x}_{i})^{\top} \mathbf{w}\right|\right)\right],$ $$

where in the last step we have used the positive-homogeneous property of $\sigma(\cdot)$ (e.g., analysis similar to handling the supremum over $W_{L},f\in\mathcal{F}_{L-1}$).
Applying the inequality above recursively over the layers leads to

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}}\exp\left(\lambda\sup_{W_{L},f \in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}}\boldsymbol{ \epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota)})\right)^{2 }\right)\leq 2^{L}\mathbb{E}_{\epsilon_{1:n}}\left[\tau_{\lambda}\left(\prod_{ l=1}^{L-2}B_{l}\left\|\sum_{i=1}^{n}\epsilon_{i}\mathbf{x}_{i}\right\|_{2} \right)\right].$ $$

Plug the inequality above into (29):

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}\in\{\pm 1\}^{nd_{L}}}\left[\sup _{W_{L},f\in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}} \boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota) })\right)^{2}\right]\leq\frac{1}{\lambda}\log\left(2^{L}\mathbb{E}_{\epsilon_{ 1:n}}\exp\left(d_{L}\lambda\left(\prod_{l=1}^{L}B_{l}^{2}\right)\left\|\sum_{i =1}^{n}\epsilon_{i}\mathbf{x}_{i}\right\|_{2}^{2}\right)\right).$ $$

Let $\tilde{\lambda}=d_{L}\lambda\left(\prod_{l=1}^{L}B_{l}^{2}\right)$ and choose $\lambda=\frac{1}{8esd_{L}\left(\prod_{l=1}^{L}B_{l}^{2}\right)}$, $s=\left(\sum_{1\leq i\leq\tilde{i}\leq n}(\mathbf{x}_{i}^{\top}X_{\tilde{i}})^
{2}\right)^{\frac{1}{2}}$. Then, $\tilde{\lambda}=1/(8es)$ and we can apply Lemma 3 to show $\mathbb{E}_{\epsilon_{1:n}}\left[\exp\left(2\tilde{\lambda}\sum_{1\leq i\leq
\tilde{i}\leq n}\epsilon_{i}\epsilon_{\tilde{i}}\mathbf{x}_{i}^{\top}X_{\tilde
{i}}\right)\right]\leq 2$ such that

$$ $\displaystyle\mathbb{E}_{\epsilon_{1:n}}\exp\left(\tilde{\lambda}\left\|\sum_{ i=1}^{n}\epsilon_{i}\mathbf{x}_{i}\right\|_{2}^{2}\right)=\mathbb{E}_{\epsilon _{1:n}}\left[\exp\left(\tilde{\lambda}\sum_{i=1}^{n}\left\|\mathbf{x}_{i} \right\|_{2}^{2}+2\tilde{\lambda}\sum_{1\leq i\leq\tilde{i}\leq n}\epsilon_{i} \epsilon_{\tilde{i}}\mathbf{x}_{i}^{\top}X_{\tilde{i}}\right)\right]$ $\displaystyle=\exp\left(\tilde{\lambda}\sum_{i=1}^{n}\left\|\mathbf{x}_{i} \right\|_{2}^{2}\right)\mathbb{E}_{\epsilon_{1:n}}\left[\exp\left(2\tilde{ \lambda}\sum_{1\leq i\leq\tilde{i}\leq n}\epsilon_{i}\epsilon_{\tilde{i}} \mathbf{x}_{i}^{\top}X_{\tilde{i}}\right)\right]\leq 2\exp\left(\tilde{\lambda }\sum_{i=1}^{n}\left\|\mathbf{x}_{i}\right\|_{2}^{2}\right).$ $$

Since $\lambda=\frac{1}{8esd_{L}\left(\prod_{l=1}^{L}B_{l}^{2}\right)}$ and $s^{2}\leq\sum_{1\leq i\leq\tilde{i}\leq n}\left\|\mathbf{x}_{i}\right\|_{2}^{2
}\left\|\mathbf{x}_{\tilde{i}}\right\|_{2}^{2}\leq\left(\sum_{i=1}^{n}\left\|
\mathbf{x}_{i}\right\|_{2}^{2}\right)^{2}$, we can obtain

$$ $\displaystyle\mathbb{E}_{\boldsymbol{\epsilon}\in\{\pm 1\}^{nd_{L}}}\left[\sup _{W_{L},f\in\mathcal{F}_{L-1}}\left(\sum_{i=1}^{n}\sum_{\iota=1}^{d_{L}} \boldsymbol{\epsilon}^{(i,\iota)}\sigma(f(\mathbf{x}_{i})^{\top}W_{L}^{(\iota) })\right)^{2}\right]\leq\frac{1}{\lambda}\log\left(2^{L+1}\exp\left(\tilde{ \lambda}\sum_{i=1}^{n}\left\|\mathbf{x}_{i}\right\|_{2}^{2}\right)\right)$ $\displaystyle=\frac{(L+1)\log 2}{\lambda}+d_{L}\left(\prod_{l=1}^{L}B_{l}^{2} \right)\sum_{i=1}^{n}\left\|\mathbf{x}_{i}\right\|_{2}^{2}\leq d_{L}\left( \prod_{l=1}^{L}B_{l}^{2}\right)\left(8(L+1)e\log 2+1\right)\sum_{i=1}^{n}\left \|\mathbf{x}_{i}\right\|_{2}^{2}.$ $$

Due to (28), we can obtain

$$ $\displaystyle\hat{\mathfrak{R}}^{+}_{n}(\mathcal{E})$ $\displaystyle=\mathbb{E}_{\epsilon_{1:n}}\left[\sup_{e\in\mathcal{E}}\frac{1}{ n}\sum_{i=1}^{n}\epsilon_{i}E(\mathbf{x}_{i},\mathbf{y}_{i})\right]\leq\frac{1 }{\sqrt{n}}\sqrt{d_{L}\left(\prod_{l=1}^{L}B_{l}^{2}\right)\left(8(L+1)e\log 2 +1\right)}(c_{1}+c_{2}).$ (30) $$

#### D.5.2 Bounding ℜ n / 2 − ⁢ ( 𝒢 ¯ ) subscript superscript ℜ 𝑛 2 ¯ 𝒢 \mathfrak{R}^{-}_{n/2}(\bar{\mathcal{G}}) fraktur_R start_POSTSUPERSCRIPT - end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_n / 2 end_POSTSUBSCRIPT ( over¯ start_ARG caligraphic_G end_ARG )

We define the dataset $\hat{\mathbf{D}}_{s}\coloneqq\{(\mathbf{x}_{s(1)},\mathbf{y}_{s(2)}),\dotsc,(
\mathbf{x}_{s(n-1)},\mathbf{y}_{s(n)})\}$. Consider $\mathcal{E}\coloneqq\{(\mathbf{x},\mathbf{y})\mapsto E_{\mathbf{w}}(\mathbf{x}
,\mathbf{y})\mid\mathbf{w}\in\mathcal{W}\}$ and the following two function classes

$$ $\displaystyle\bar{\mathcal{E}}\coloneqq\{(\mathbf{x},\mathbf{y})\mapsto\bar{E} _{\mathbf{w}}(\mathbf{x},\mathbf{y})\mid\mathbf{w}\in\mathcal{W}\},\quad\bar{ \mathcal{G}}=\{(\mathbf{x},\mathbf{y})\mapsto\exp(\bar{E}_{\mathbf{w}}(\mathbf {x},\mathbf{y}))\mid\mathbf{w}\in\mathcal{W}\}.$ $$

The empirical Rademacher complexities of $\bar{\mathcal{E}},\bar{\mathcal{G}}$ on $\hat{\mathbf{D}}_{s}$ can be defined as

$$ $\displaystyle\hat{\mathfrak{R}}^{-}_{n/2}(\bar{\mathcal{E}};s)$ $\displaystyle=\mathbb{E}_{\epsilon_{1:n/2}}\left[\frac{2}{n}\sup_{\mathbf{w}} \sum_{i=1}^{n/2}\epsilon_{i}\bar{E}_{\mathbf{w}}(\mathbf{x}_{s(2i-1)},\mathbf{ y}_{s(2i)})\right],$ $\displaystyle\hat{\mathfrak{R}}^{-}_{n/2}(\bar{\mathcal{G}};s)$ $\displaystyle=\mathbb{E}_{\epsilon_{1:n/2}}\left[\frac{2}{n}\sup_{\mathbf{w}} \sum_{i=1}^{n/2}\epsilon_{i}\exp(\bar{E}_{\mathbf{w}}(\mathbf{x}_{s(2i-1)}, \mathbf{y}_{s(2i)}))\right].$ $$

Note that $t\mapsto\exp(t)$ is $1$-Lipschitz when $t\leq 0$. Due to Lemma 4 and $\bar{E}_{\mathbf{w}}(\mathbf{x},\mathbf{y})=(E_{\mathbf{w}}(\mathbf{x},\mathbf
{y})-1)/\tau$,

$$ $\displaystyle\hat{\mathfrak{R}}^{-}_{n/2}(\bar{\mathcal{G}};s)\leq\hat{ \mathfrak{R}}^{-}_{n/2}(\bar{\mathcal{E}};s)=\frac{1}{\tau}\hat{\mathfrak{R}}^ {-}_{n/2}(\mathcal{E};s).$ (31) $$

Then, we can bound $\hat{\mathfrak{R}}^{-}_{n/2}(\mathcal{E};s)$ in the way similar to bounding $\hat{\mathfrak{R}}^{+}_{n}(\mathcal{E})$ in Section D.5.1.

## Appendix E Proof of Theorem 2

###### Proof.

The problem in (8) is equivalent to

$$ $\displaystyle\min_{\boldsymbol{\zeta}\in\mathbb{R}^{n}}\left\{\frac{1}{n}\sum_ {i=1}^{n}\tau\log\left(\sum_{j=1}^{n}\exp((E(\mathbf{x}_{i},\mathbf{y}_{j})-E( \mathbf{x}_{i},\mathbf{y}_{i})-\zeta^{(j)})/\tau)\right)+\frac{1}{n}\sum_{j=1} ^{n}\zeta^{(j)}\right\}.$ (32) $$

We define that $\Phi(\boldsymbol{\zeta})\coloneqq\frac{1}{n}\sum_{i=1}^{n}\tau\log\left(\sum_{
j=1}^{n}\exp((E(\mathbf{x}_{i},\mathbf{y}_{j})-E(\mathbf{x}_{i},\mathbf{y}_{i}
)-\zeta^{(j)})/\tau)\right)+\frac{1}{n}\sum_{j=1}^{n}\zeta^{(j)}$. Due to the first-order optimality
condition, any optimal solution $\boldsymbol{\zeta}_{*}$ to (8) satisfies

$$ $\displaystyle\exp(\zeta_{*}^{(j)}/\tau)=\sum_{j^{\prime}=1}^{n}\frac{\exp(E( \mathbf{x}_{j^{\prime}},\mathbf{y}_{j})/\tau)}{\sum_{i^{\prime}=1}^{n}\exp((E( \mathbf{x}_{j^{\prime}},\mathbf{y}_{i^{\prime}})-\zeta_{*}^{i^{\prime}})/\tau) },\quad\forall j\in[n].$ $$

Then, we can obtain (9) by changing the variable $\bar{q}^{(j)}=\exp(\zeta_{*}^{(j)}/\tau)$ for any $j\in[n]$.

Due to the property of the log-sum-exp function and $E(\mathbf{x}_{i},\mathbf{y}_{j})\in[-1,1]$, we have

$$ $\displaystyle\Phi(\boldsymbol{\zeta})\geq\frac{1}{n}\sum_{i=1}^{n}\max_{j\in[n ]}\left\{E(\mathbf{x}_{i},\mathbf{y}_{j})-E(\mathbf{x}_{i},\mathbf{y}_{i})- \zeta^{(j)}\right\}+\frac{1}{n}\sum_{j=1}^{n}\zeta^{(j)}\geq-2-\min_{j\in[n]} \zeta^{(j)}+\frac{1}{n}\sum_{j=1}^{n}\zeta^{(j)}\geq-2.$ $$

Thus, $\Phi(\boldsymbol{\zeta})$ is proper convex. Besides, each line parallel to the diagonal line $\mathfrak{d}_{n}=\{\boldsymbol{\zeta}\mid\boldsymbol{\zeta}=z\mathbf{1}_{n},z
\in\mathbb{R}\}$ can be expressed as $\bar{\mathfrak{d}}_{n}(\mathbf{b})=\left\{\boldsymbol{\zeta}\Big{\mid}
\boldsymbol{\zeta}=z\mathbf{1}_{n}+\begin{bmatrix}0\\
\mathbf{b}\end{bmatrix},z\in\mathbb{R}\right\}$ with some unique $\mathbf{b}\in\mathbb{R}^{n-1}$. For example, in the case $n=2$, a line $\bar{\mathfrak{d}}_{2}(b)$ with some $b\in\mathbb{R}$ can be rewritten as $\bar{\mathfrak{d}}_{2}(b)=\{\boldsymbol{\zeta}\mid\zeta^{(2)}=\zeta^{(1)}+b\}$, which is parallel to the diagonal line $\mathfrak{d}_{2}=\{\boldsymbol{\zeta}\mid\zeta^{(2)}=\zeta^{(1)}\}$. For any point $\boldsymbol{\zeta}$ on a line $\bar{\mathfrak{d}}_{n}(\mathbf{b})$,

$$ $\displaystyle\Phi(\boldsymbol{\zeta})$ $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\tau\log\left(\sum_{j=1}^{n}\exp((E( \mathbf{x}_{i},\mathbf{y}_{j})-E(\mathbf{x}_{i},\mathbf{y}_{i})-z+b^{(j)})/ \tau)\right)+z+\frac{1}{n}\sum_{j=1}^{n}b^{(j)}$ $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\tau\log\left(\exp(-z/\tau)\sum_{j=1}^{ n}\exp((E(\mathbf{x}_{i},\mathbf{y}_{j})-E(\mathbf{x}_{i},\mathbf{y}_{i})-b^{( j)})/\tau)\right)+z+\frac{1}{n}\sum_{j=1}^{n}b^{(j)}$ $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\tau\log\left(\sum_{j=1}^{n}\exp((E( \mathbf{x}_{i},\mathbf{y}_{j})-E(\mathbf{x}_{i},\mathbf{y}_{i})-b^{(j)})/\tau) \right)+\frac{1}{n}\sum_{j=1}^{n}b^{(j)},$ $$

where the expression on the R.H.S is fixed when $z$ varies. Then, the value of $\Phi(\boldsymbol{\zeta})$ does not change along the line $\bar{\mathfrak{d}}_{n}(\mathbf{b})$. Note that every point $\boldsymbol{\zeta}\in\mathbb{R}^{n}$ is uniquely located on one line $\bar{\mathfrak{d}}_{n}(\mathbf{b})$ parallel to the diagonal line $\mathfrak{d}_{n}$. Thus, if $\boldsymbol{\zeta}_{*}=z_{*}\mathbf{1}_{n}+\begin{bmatrix}0\\
\mathbf{b}_{*}\end{bmatrix}$ with specific $z_{*}\in\mathbb{R}$ and $\mathbf{b}_{*}\in\mathbb{R}^{n-1}$ is a minimum of $\Phi(\boldsymbol{\zeta})$, then any point on the line $\bar{\mathfrak{d}}_{n}(\mathbf{b}_{*})$ is a minimum of $\Phi(\boldsymbol{\zeta})$.

There may exist uncountably infinite many $\mathbf{b}_{*}\in\mathbb{R}^{n-1}$ and every point on $\bar{\mathfrak{d}}_{n}(\mathbf{b}_{*})$ minimizes $\Phi(\boldsymbol{\zeta})$. However, we can rule out such a case since the set of minima of a proper convex function is convex and $\Phi(\boldsymbol{\zeta})$ is strictly convex along any direction other than the diagonal and parallel lines(^6^66Each log-sum-exp function is strictly convex along any direction other than diagonal and parallel lines.). Thus, there is only a unique $\mathbf{b}_{*}\in\mathbb{R}^{n-1}$ and any point on $\bar{\mathfrak{d}}_{n}(\mathbf{b}_{*})=\left\{\boldsymbol{\zeta}\Big{\mid}
\zeta=z\mathbf{1}_{n}+\begin{bmatrix}0\\
\mathbf{b}_{*}\end{bmatrix},z\in\mathbb{R}\right\}$ minimizes $\Phi(\boldsymbol{\zeta})$, i.e., the minimum of $\Phi(\boldsymbol{\zeta})$ is unique up to an arbitrary scalar additive term $z\in\mathbb{R}$.

∎

###### Proof.

## Appendix F More Discussions on NUCLR in Algorithm 1

### F.1 Computational and Memory Overheads of NUCLR

Compared to the $O(Bd)$ per-iteration computational cost of the CLIP, the proposed NUCLR leads to a computational overhead $O(B)$ similar to SogCLR for updating the scalars $\{u^{(i)}\}_{i\in\mathcal{B}_{t}}$ and $\{\zeta^{(i)}\}_{i\in\mathcal{B}_{t}}$. This extra $O(B)$ cost can be ignored since $d$ is extremely large in modern deep neural networks. Owing to the moving average estimator $\mathbf{u}$, our NUCLR does not require a huge batch size $B$ for good performance, unlike SimCLR/CLIP. Thus, NUCLR is also more memory-efficient, making it suitable for environments with limited GPU resources, similar to SogCLR.

Compared to SogCLR, NUCLR needs to store one extra $n$-dimensional vector $\boldsymbol{\zeta}$. Maintaining $\boldsymbol{\zeta}$ in GPU only requires less than 100MB for 12 million data points, which is negligible compared to the GPU memory required for backpropagation. Moreover, we may instead maintain the vector $\boldsymbol{\zeta}$ in CPU and only transfer those needed $\{\zeta^{(j)}\}_{j\in\mathcal{B}_{t}}$ to GPU in each iteration. The overhead can be further reduced by overlapping the communication and computation.

### F.2 Margin Interpretation of NUCLR

Cross-entropy and contrastive losses with a positive additive margin have been widely studied in the literature [Li et al., , Liu et al., , Wang et al., , Cao et al., , Li et al., , Zhu et al., ], which can be viewed as a smooth version of the hinge loss to separate the matching (positive) pair $(\mathbf{x}_{i},\mathbf{y}_{i})$ from negative pairs $\{(\mathbf{x}_{i},\mathbf{y}_{j})\mid\mathbf{y}_{j}\neq\mathbf{y}_{i},\mathbf{
y}_{j}\in\mathcal{Y}\}$. In supervised learning tasks such as face verification and multi-class classification, using a relatively large positive margin has been shown to be beneficial [Wang et al., , Cao et al., ]. However, the “false negative” issue is more pronounced in self-supervised learning. Determining the appropriate (positive or negative) margin becomes more difficult, as aggressively and uniformly pushing negative pairs away from positive pairs may hurt the performance [Xie et al., ]. As shown in the objective in (11), our NUCLR algorithm adopts an *individualized* negative margin $\zeta^{(j)}$ for each negative data $\mathbf{y}_{j}$ when updating the model parameter $\mathbf{w}$. Rather than relying on an expensive grid search for individualized margins, our method learns them in a principled way. Recall $\zeta^{(j)}$ can serve as a measure of the popularity since $\tilde{q}^{(j)}\propto\exp(\zeta^{(j)}/\tau)$ when $\zeta^{(j)}$ is optimized. As a result, NUCLR may help *tolerate* potential false negatives because the negative margin $\zeta^{(j)}$ between pairs $(\mathbf{x}_{i},\mathbf{y}_{i})$ and $(\mathbf{x}_{i},\mathbf{y}_{j})$ is larger when $\mathbf{y}_{j}$ is popular, as it is more likely to be a false negative.

## Appendix G More Experimental Results

### G.1 Additional Plot of the Toy Experiment in Section 3.3

Figure: Figure 6: Comparing the error term $\mathcal{E}(\tilde{\mathbf{q}},\mathbf{q},\hat{\mathbf{S}})$ of our method and GCL across various $n$.
Refer to caption: x14.png

In addition to the figure of generalization error shown in Figure 2, we also plot the error term $\mathcal{E}(\tilde{\mathbf{q}},\mathbf{q},\hat{\mathbf{S}})$ as defined in Theorem 1 of our method and GCL (uniform approximation).

### G.2 More Bimodal Experimental Results

#### G.2.1 Comparison with iSogCLR

In Section 4 of the main paper, both our algorithm NUCLR and baselines (CLIP, SigLIP, CyCLIP, DCL, SogCLR) use a shared temperature parameter $\tau$ for all data pairs. In contrast, the iSogCLR algorithm [Qiu et al., ] learns an individual temperature parameter for each pair of data. Here we compare the testing performance of our NUCLR algorithm to that of iSogCLR. The results show that our NUCLR outperforms iSogCLR on most metrics except for the zero-shot classification on ImageNet1k when pretrained on the CC12M dataset. Lastly, we note that the individual temperature learned by the approach in iSogCLR could also be incorporated into our algorithm.

**Table 2: Test performance of NUCLR and iSogCLR.**
| Dataset | Algorithm | MSCOCO | Flickr30k | CIFAR100 | ImageNet1k | Mean |
| --- | --- | --- | --- | --- | --- | --- |
| CC3M | iSogCLR | 28.50 $\pm$ 0.21 | 51.56 $\pm$ 0.38 | 35.57 $\pm$ 0.99 | 40.18 $\pm$ 0.28 | 38.95 $\pm$ 0.30 |
| NUCLR | 29.55 $\pm$ 0.26 | 53.55 $\pm$ 0.22 | 37.45 $\pm$ 0.45 | 40.49 $\pm$ 0.30 | 40.26 $\pm$ 0.19 |  |
| CC12M | iSogCLR | 34.09 $\pm$ 0.25 | 59.42 $\pm$ 0.41 | 27.78 $\pm$ 1.75 | 50.23 $\pm$ 0.18 | 42.88 $\pm$ 0.38 |
| NUCLR | 34.36 $\pm$ 0.13 | 60.45 $\pm$ 0.03 | 28.16 $\pm$ 1.35 | 49.82 $\pm$ 0.23 | 43.20 $\pm$ 0.39 |  |

#### G.2.2 Effect of the Initial Value ζ 0 subscript 𝜁 0 \zeta_{0} italic_ζ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT

**Table 3: Test performance of NUCLR with different initial values of $\zeta_{0}$.**
| Dataset | Algorithm | MSCOCO | Flickr30k | CIFAR100 | ImageNet1k | Mean |
| --- | --- | --- | --- | --- | --- | --- |
| CC3M | SogCLR | 28.54 $\pm$ 0.25 | 52.20 $\pm$ 0.64 | 35.50 $\pm$ 1.71 | 40.40 $\pm$ 0.12 | 39.16 $\pm$ 0.33 |
| NUCLR ($\zeta_{0}=0.0$) | 29.51 $\pm$ 0.12 | 53.10 $\pm$ 0.35 | 37.17 $\pm$ 1.27 | 40.21 $\pm$ 0.24 | 40.00 $\pm$ 0.26 |  |
| NUCLR ($\zeta_{0}=-0.05$) | 29.55 $\pm$ 0.26 | 53.55 $\pm$ 0.22 | 37.45 $\pm$ 0.45 | 40.49 $\pm$ 0.30 | 40.26 $\pm$ 0.19 |  |
| CC12M | SogCLR | 33.91 $\pm$ 0.26 | 59.28 $\pm$ 0.07 | 26.10 $\pm$ 0.88 | 49.82 $\pm$ 0.14 | 42.28 $\pm$ 0.27 |
| NUCLR ($\zeta_{0}=0.0$) | 34.36 $\pm$ 0.13 | 60.45 $\pm$ 0.03 | 28.16 $\pm$ 1.35 | 49.82 $\pm$ 0.23 | 43.20 $\pm$ 0.39 |  |
| NUCLR ($\zeta_{0}=-0.05$) | 34.35 $\pm$ 0.14 | 59.69 $\pm$ 0.16 | 26.30 $\pm$ 2.35 | 49.90 $\pm$ 0.28 | 42.56 $\pm$ 0.59 |  |

We also investigate the performance of NUCLR with two different initial values $\zeta_{0}$ of the auxiliary variable $\boldsymbol{\zeta}$: 1) The natural choice $\zeta_{0}=0$; and 2) The value $\zeta_{0}=-0.05$ according to EqCo [Zhu et al., ]. As seen in Table 3, $\zeta_{0}=-0.05$ performs better on the CC3M dataset while $\zeta_{0}=0$ yields better results on the CC12M dataset. It is worth noting that our NUCLR algorithm, with either initialization $\zeta_{0}$, leads to overall better performance than SogCLR.

#### G.2.3 Histograms of 𝜻 𝜻 \boldsymbol{\zeta} bold_italic_ζ after the Final Epoch

Figure: Figure 7: Histograms of the learned $\boldsymbol{\zeta}$ after the final epoch on the CC3M and CC12M datasets.
Refer to caption: x15.png

In Figure 7, we visualize the final distributions of the learned $\boldsymbol{\zeta}$ for images and texts in the CC3M and CC12M datasets. These distributions exhibit long tails, suggesting that many images or texts have low popularities in web-based datasets like CC3M and CC12M.

#### G.2.4 Images from CC3M dataset with large and small learned popularity

Figure 8 and 9 provide more images from CC3M with small and large learned popularity $\tilde{q}^{\prime}$.

Figure: Figure 8: Images (and their captions) from the CC3M dataset with large $\tilde{q}^{\prime}$ (high learned popularity).
Refer to caption: x19.png

Figure: Figure 9: Images (and their captions) from the CC3M dataset with small $\tilde{q}^{\prime}$ (low learned popularity).
Refer to caption: x20.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: please refer to the weakness part.
soundness: 2
strengths: 1.	It seems a novel and interesting perspective to improve self-supervised contrastive learning with discriminative probabilistic modeling.
2.	It provides theoretical insights into the generalization error limitation in self-supervised contrastive learning and presents a reasonable solution to address this problem.
3.	Experiments on CC3M and CC12M are conducted to verify the effectiveness of the proposed method.
4.	The paper is generally well-structured and easy to follow.
summary: This paper explores discriminative probabilistic modeling in unimodal and multimodal self-supervised learning. By applying multiple importance sampling to Monte Carlo integration, it recovers the contrastive loss as a special case. The paper also identifies a generalization error limitation inherent in contrastive learning frameworks. To address this, it proposes a non-parametric method aimed at mitigating this error. The effectiveness of the approach is validated through experiments on multimodal datasets.
weaknesses: -	This paper claims that the vanilla contrastive learning objective (e.g., GCL) results in non-diminishing error terms due to a misalignment between the uniform distribution and the true data distribution on y. However, it is unclear if these error terms are substantial and whether they might worsen under long-tailed data distributions. Would incorporating long-tailed contrastive learning methods [1,2], which use automatic or prior population distribution estimation, help address this issue?
-	The paper lacks quantitative or qualitative results demonstrating that the proposed method effectively reduces the non-diminishing error terms (as seen in Eq. 6) on real-world datasets like CC3M and CC12M.
-	The proposed method may introduce additional computational overhead. It would be beneficial to discuss and compare the computational cost relative to baseline methods. 
-	According to Figure 4 in the ablation study, the proposed method appears to require specific design tricks for optimal performance, raising questions about its generalizability and applicability across varied scenarios. Clarifying whether substantial hyperparameter tuning is needed for different applications would be helpful.
-	The theoretical framework suggests applicability to unimodal contrastive learning. The paper would be further strengthened by demonstrating its effectiveness on common SSL benchmarks for unimodal data. 
-	The empirical evaluations focus on classification and retrieval tasks across four datasets (MSCOCO, Flickr, CIFAR100, and ImageNet). Including more challenging datasets and a wider variety of tasks would enhance the robustness of the results.

[1] The Hidden Uniform Cluster Prior in Self-Supervised Learning.
[2] Combating Representation Learning Disparity with Geometric Harmonization.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: Given the method is to handle discriminative probabilistic modeling problem on a continuous domain for (multimodal) self-supervised representation learning, the scope of the method that can be applied seems to be narrow. I am wondering if the method can be extended to other cases outside of self-supervised representation learning such as unsupervised learning.

It will be also useful to provide examples in the experimental results to show the method can be applied to image and text, audio and image with the connection to the proposed theory to show the advantages. So that the readers can deeply understand how to apply the theory to the real examples.
soundness: 3
strengths: The paper addresses an interesting and important task of discriminative probabilistic modeling problem on a continuous domain for (multimodal) self-supervised representation learning. The idea appears to be novel in this field.

The theoretical analysis is thorough and with a good depth. The generalization error analysis is conducted to show the limitation of the existing algorithm. 

The experimental evaluation and comparison are convincing and sufficient to show the advantages of the proposed algorithm. Overall, it is a solid paper.
summary: The manuscript studies the discriminative probabilistic modeling problem on a continuous domain for (multimodal) self-supervised representation learning. To address the challenge of computing the integral in the partition function for each anchor data, the multiple importance sampling (MIS) technique is considered for robust Monte Carlo integration, which can recover InfoNCE-based contrastive loss as a special case. Typically, the generalization error analysis is conducted to reveal the limitation of current InfoNCE-based contrastive loss for self-supervised representation learning and derive insights for developing better approaches by reducing the error of Monte Carlo integration.

In order to resolve the problem, a novel non-parametric method is proposed for approximating the sum of conditional densities required by MIS through convex optimization, yielding a new contrastive objective for self-supervised representation learning. Moreover,  an efficient algorithm is designed for solving the proposed objective. The algorithm is compared with representative baselines on the contrastive image-language pretraining task. Experimental results on the CC3M and CC12M datasets demonstrate the superior overall performance of the proposed algorithm.
weaknesses: The weakness of the paper includes:

(1) As it is a theory heavy paper, it will be better to provide more motivation of the proposed method so that the readers can deeply understand why MIS is utilized and how is the convergence of MIS. In particulr, multiple importance sampling is a routine and well-known method. I believe that the uniqueness is probably the application of MIS in the particular framework. However, the author needs to specify and highlight what exact the uniqueness it is. (Is it a different math formulation etc)

(2) It is necessary to highlight the complexity and running time of the proposed method. I suggest to add a subsection for discussion the complexity of the proposed algorithm and comparing the complexity and running time with competing methods.

(3) The scope of the paper seems to be narrow (Can only be applied in  a continuous domain for (multimodal) self-supervised representation learning). For instance, can the proposed method be extend to unsupervised learning or semi-supervised learning ? Moreover, since it is multimodal based solution, can the proposed method be extended to text and image, image and audio etc ?

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: Besides the above weaknesses, other questions are
1. How to determine the involved sizes m and n and is there an optimal relation between them? 
2. Should the m and n be related to dimensions of x and y spaces? What relation is to the bootstrap?
3. Formally, (1) can replace as p_v(x|y)!  When the dimensions of x and y are imbalanced, e.g., y's dimension is far greater than that of x, whether does paired sampling influence the quality of analysis?
4. Assumption 2 requires the same d_L, I still want to know if the large difference of x's and y's dimensions impact their corresponding E1 and E2.
5. (7) has some flavor of the attention mechanism to great extent, it is unclear whether this can be interpreted.
6. Is the objective or loss measure (10) robust? e,g,, x_i or y_i is noisy!
7. In Proposition ,can its (iii) condition satisfied? can the authors give some concrete cases?
soundness: 3
strengths: 1. For the DMP problem on a continuous domain for (multimodal) self-supervised representation learning, the authors first point out a challenge in computing the integral in the partition function for each anchor data, then adopt the MIS method to realize robust Monte
Carlo integration. Particular;y, the authors make generalization error analysis.
2. For current limitation of current InfoNCE-based contrastive loss for self-supervised representation learning, the authors propose a new non-parametric method for approximating the sum of conditional probability densities required by MIS through convex optimization, as a result, yielding a new contrastive objective which is optimized by so-designed an efficient algorithm called NUCLR.
3. So-conducted coparative results with baselines get improved performance.
summary: The authors study the DPM problem on a continuous domain for (multimodal) self-supervised contrastive representation learning by leveraging the MIS method and proposing a novel nonparametric method for approximating the sum of conditional probability densities, establish a new contrastive loss for self-supervised representation learning and optimize it by an efficient algorithm called NUCLR. Experimental results on bimodal pretraining confirm the improvement in performance of their method compared to baseline approaches
on downstream tasks.
weaknesses: In my opinion, the manuscript is basically complete, just some weaknesses are listed as follows:
1. Lack an analysis of robustness of their results;
2. Assumption 2 requires the same d_L, I want to know if the large difference of x's and y's dimensions impact their corresponding E1 and E2.
3. In Proposition ,its (iii) condition seems difficult to be satisfied.

Meta-review: The paper introduces a non-parametric method leveraging multiple importance sampling (MIS) to enhance self-supervised representation learning, tackling the limitations associated with the InfoNCE loss. By utilizing multiple importance sampling for Monte Carlo integration, the approach retrieves the contrastive loss as a specific instance. Additionally, the paper highlights a limitation related to generalization error that is inherent in contrastive learning frameworks. This approach has demonstrated superior performance over existing baselines in image-language pretraining tasks.

The reviewers unanimously agree that the paper addresses an interesting and important task regarding discriminative probabilistic modeling problems, and the ideas presented in this field appear to be novel. Additionally, the paper provides theoretical analyses that further expand the interpretability of the proposed method. The overall review of the paper is positive, leading to a recommendation for acceptance at this time.

Decision: Accept (Poster)