\section{Introduction}

The natural language processing community has made tremendous progress
in using pre-trained language models to improve predictive accuracy \citep{Devlin2019BERTPO,Raffel2019ExploringTL}. Models have now surpassed human performance on language understanding benchmarks such as SuperGLUE \citep{Wang2019SuperGLUEAS}. However, studies have shown that these results are partially driven by these models detecting superficial cues that correlate well with labels but which may not be useful for the intended underlying task \citep{Jia2017AdversarialEF,Schwartz2017TheEO}. This brittleness leads to overestimating model performance on the artificially constructed tasks and poor performance in out-of-distribution or adversarial examples.

A well-studied example of this phenomenon is the natural language inference dataset MNLI \citep{N18-1101}. The generation of this dataset led to spurious surface patterns that correlate noticeably with the labels. \cite{Poliak2018HypothesisOB} highlight that negation words (``not'', ``no'', etc.) are often associated with the contradiction label. \cite{Gururangan2018AnnotationAI,Poliak2018HypothesisOB,Tsuchiya2018PerformanceIC} show that a model trained solely on the hypothesis, completely ignoring the intended signal, reaches strong performance. We refer to these surface patterns as \textit{dataset biases} since the conditional distribution of the labels given such biased features is likely to change in examples outside the training data distribution (as formalized by \cite{He2019UnlearnDB}).

A major challenge in representation learning for NLP is to produce models that are robust to these dataset biases. Previous work \citep{He2019UnlearnDB,Clark2019DontTT,Mahabadi2020EndtoEndBM} has targeted removing dataset biases by explicitly factoring them out of models. These works explicitly construct a biased model, for instance, a hypothesis-only model for NLI experiments, and use it to improve the robustness of the main model. The core idea is to encourage the main model to find a different explanation where the biased model is wrong. During training, products-of-experts ensembling \citep{Hinton2002TrainingPO} is used to factor out the biased model.

While these works show promising results, the assumption of knowledge of the underlying dataset bias is quite restrictive. Finding dataset biases in established datasets is a costly and time-consuming process, and may require access to private details about the annotation procedure, while actively reducing surface correlations in the collection process of new datasets is challenging given the number of potential biases \citep{Zellers2019HellaSwagCA,Sakaguchi2020WINOGRANDEAA}.

In this work, we explore methods for learning from biased datasets which do not require such an explicit formulation of the dataset biases. We first show how a model with limited capacity, which we call a \textit{weak learner}, trained with a standard cross-entropy loss learns to exploit biases in the dataset. We then investigate the biases on which this weak learner relies and show that they match several previously manually identified biases. Based on this observation, we leverage such limited capacity models in a product of experts ensemble to train a more robust model and evaluate our approach in various settings ranging from toy datasets up to large crowd-sourced benchmarks: controlled synthetic bias setup \citep{He2019UnlearnDB,Clark2019DontTT}, natural language inference \citep{McCoy2019RightFT} and extractive question answering \citep{Jia2017AdversarialEF}.

Our contributions are the following: (a) we show that weak learners are prone to relying on shallow heuristics and highlight how they rediscover previously human-identified dataset biases; (b) we demonstrate that we do not need to explicitly know or model dataset biases to train more robust models that generalize better to out-of-distribution examples; (c) we discuss the design choices for weak learners and show trade-offs between higher out-of-distribution performance at the expense of the in-distribution performance.

\section{Related work}

Many studies have reported dataset biases in various settings. Examples include visual question answering \citep{Jabri2016RevisitingVQ,Zhang2016YinAY}, story completion \citep{Schwartz2017TheEO}, and reading comprehension \citep{Kaushik2018HowMR,Chen2016ATE}. Towards better evaluation methods, researchers have proposed to collect ``challenge'' datasets that account for surface correlations a model might adopt \citep{Jia2017AdversarialEF,McCoy2019RightFT}. Standard models without specific robust training methods often drop in performance when evaluated on these challenge sets.

While these works have focused on data collection, another approach is to develop methods allowing models to ignore dataset biases during training. Several active areas of research tackle this challenge by adversarial training \citep{Belinkov2019DontTT,Stacey2020ThereIS}, example forgetting \citep{Yaghoobzadeh2019RobustNL} and dynamic loss adjustment \citep{Cadne2019RUBiRU}. Previous work \citep{He2019UnlearnDB,Clark2019DontTT,Mahabadi2020EndtoEndBM} has shown the effectiveness of product of experts to train un-biased models.  In our work, we show that we do not need to explicitly model biases to apply these de-biasing methods and can use a more general setup than previously presented.

Orthogonal to these evaluation and optimization efforts, data augmentation has attracted interest as a way to reduce model biases by explicitly modifying the dataset distribution \citep{Min2020SyntacticDA,Belinkov2018SyntheticAN}, either by leveraging human knowledge about dataset biases such as swapping male and female entities \citep{Zhao2018GenderBI} or by developing dynamic data collection and benchmarking \citep{nie2020adversarial}. Our work is mostly orthogonal to these efforts and alleviates the need for a human-in-the-loop setup which is common to such data-augmentation approaches.

Large pre-trained language models have contributed to improved out-of-distribution generalization \citep{Hendrycks2020PretrainedTI}. However, in practice, that remains a challenge in natural language processing \citep{Linzen2020HowCW,Yogatama2019LearningAE} and our work aims at out-of-distribution robustness without significantly compromising in-distribution performance.

Finally, as we were preparing this manuscript for submission, we became aware of a parallel work \citep{utama2020debiasing} which presents a related de-biasing method leveraging \textit{shallow models}'s mistakes without the need to explicitly model dataset biases. Our approach is different in several ways, in particular we advocate for using limited capacity weak learner while \cite{utama2020debiasing} uses the same architecture as the robust model trained on a few thousands examples. We investigated the trade-off between learner's capacity and resulting performances as well as the resulting few-shot learning regime in the limit of a high capacity weak model.

\section{Method}

\subsection{Overview}

Our approach utilizes \textit{product of experts} \citep{Hinton2002TrainingPO} to factor dataset biases out of a learned model. We have access to a training set $(x_i, y_i)_{1 \leq i \leq N}$ where each example $x_i$ has a label $y_i$ among $K$ classes. We use two models $f_W$ (weak) and $f_M$ (main) which produce respective logits vectors $\mathbf{w}$ and $\mathbf{m} \in \mathbb{R}^K$. The product of experts ensemble of $f_W$ and $f_M$ produces logits vector $\mathbf{e}$
\begin{equation}
    \label{eq:poe}
    \forall 1 \leq j \leq K, e^j = w^j + m^j
\end{equation}
Equivalently, we have $\mathrm{softmax}(\mathbf{e}) \propto \mathrm{softmax}(\mathbf{w}) \odot \mathrm{softmax}(\mathbf{m})$ where $\odot$ is the element-wise multiplication.

Our training approach can be decomposed in two successive stages: (a) training the weak learner  $f_W$ with a standard cross-entropy loss (CE) and (b) training a main (robust) model $f_M$ via product of experts (PoE) to learn from the errors of the weak learner. The core intuition of this method is to encourage the robust model to learn to make predictions that take into account the weak learner's mistakes.

We do not make any assumption on the biases present (or not) in the dataset and rely on letting the weak learner discover them during training. Moreover, in contrast to prior work \citep{Mahabadi2020EndtoEndBM,He2019UnlearnDB,Clark2019DontTT} in which the weak learner had a hand-engineered bias-specific structure, our approach does not make any specific assumption on the weak learner such as its architecture, capacity, pre-training, etc. The weak learner $f_W$ is trained with standard cross-entropy.

The final goal is producing main model $f_M$. After training, the weak model $f_W$ is frozen and used only as part of the product of experts. Since the weak model is frozen, only the main model $f_M$ receives gradient updates during training. This is similar to \cite{He2019UnlearnDB,Clark2019DontTT} but differs from \cite{Mahabadi2020EndtoEndBM} who train both weak and main models jointly. For convenience, we refer to the cross-entropy of the prediction $\mathbf{e}$ of Equation~\ref{eq:poe} as the \textit{PoE cross-entropy}.

\subsection{Analysis: The robust model learns from the errors of the weak learner}
\label{sec:gradients_interpretation}

To better explore the impact of PoE training with a weak learner, we consider the special case of binary classification with logistic regression. Here $w$ and $m$ are scalar logits and the softmax becomes a sigmoid.  The loss of the product of experts for a single positive example is:
\begin{equation}
    \label{eq:toy_loss}
    \mathcal{L}_{PoE,binary} = - m - w + \log \big(1 + \exp(m+w) \big)
\end{equation}
Logit $w$ is a fixed value since the weak learner is frozen. We also define the entropy of the weak learner as  $\mathcal{H}_w = - p \log(p) - (1-p) \log(1-p)$ where $p=\sigma(w)$ as our measure of certainty.

Different values of $w$ from the weak learner induce different gradient updates in the main model. Figure~\ref{fig:gradients} shows the gradient update of the main model logit $m$. Each of the three curves corresponds to a different value of $w$ the weak model.
\begin{itemize}
    \item \textbf{Weak Model is Certain / Incorrect}: the first case (in blue) corresponds to low values of $w$. The entropy is low and the loss of the weak model is high. The main model receives gradients even when it is classifying the point correctly ($\approx m=5$) which encourages $m$ to compensate for the weak model's mistake.
    \item \textbf{Weak Model is Uncertain}: the second case (in red) corresponds to $w=0$ which means the weak model's entropy is high (uniform probability over all classes). In this case, product of experts is equal to the main model, and the gradient is equal to the one obtained with cross-entropy.
    \item \textbf{Weak Model is Certain / Correct}: the third case (in green) corresponds to high values of $w$. The entropy is low and the loss of the weak model is low. In this case, $m$'s gradients are ``cropped'' early on and the main model receives less gradients on average. When $w$ is extremely high, $m$ receives no gradient (and the current example is simply ignored).
\end{itemize}

Put another way, the logit values for which $m$ receives gradients are shifted according the correctness and certainty of the weak model. Figure~\ref{fig:2d_plots} shows the concentration of training examples of MNLI \citep{N18-1101} projected on the 2D coordinates (correctness, certainty) from a trained weak learner (described in Section~\ref{sec:qualitative_examples}). We observe that there are many examples for the 3  cases. More crucially, we verify that the group \textbf{certain / incorrect} is not empty since the examples in this group encourage the model to not rely on the dataset biases.

\begin{figure}[h]
  \centering
  \subfloat[Gradient update of $m$ for different values of $w$ on binary classification.]{
    \includegraphics[width=0.48\columnwidth]{gradients.png}
    \label{fig:gradients}
  }\hfill
  \subfloat[2D projection of MNLI examples from a trained weak learner. Colors indicate the concentration and are in log scale.]{
    \includegraphics[width=0.46\columnwidth]{ent_loss_tiny_mnli_concentration.png}
    \label{fig:2d_plots}
  }
  \caption{The analysis of the gradients reveals 3 regimes where the gradient is shifted by the certainty and correctness of the weak learner. These 3 regions are present in real dataset such as MNLI.}
  \label{fig:understanding_the_maths}
\end{figure}

\subsection{Connection to distillation}
\label{sec:distillation}

Our product of experts setup bears similarities with \textit{knowledge distillation} \citep{Hinton2015DistillingTK} where a student network is trained to mimic the behavior of a frozen teacher model.

In our PoE training, we encourage the main model $f_M$ (analog to the student network) to learn from the errors of the weak model $f_W$ (analog to the teacher network): instead of mimicking, it learns an orthogonal behavior when the teacher is incorrect. To recognize errors from the weak learner, we use the gold labels which alleviates the need to use pseudo-labelling or data augmentation as it is commonly done in distillation setups \citep{Furlanello2018BornAN,Xie2020SelfTrainingWN}.

Similarly to \citet{Hinton2015DistillingTK}, our final loss is a linear combination of the original cross-entropy loss (CE) and the PoE cross-entropy. We refer to this multi-loss objective as \textit{PoE + CE}.

\section{Experiments}

We consider several different experimental settings that explore the use of a weak learner to isolate and train against dataset biases. All the experiments are conducted on English datasets, and follow the standard setup for BERT training. Our main model is BERT-base \citep{Devlin2019BERTPO} with 110M parameters. Except when indicated otherwise, our weak learner is a significantly smaller pre-trained masked language model known as  TinyBERT \citep{Turc2019WellReadSL} with 4M parameters (2 layers, hidden size of 128). The weak learner is fine-tuned on exactly the same data as our main model. For instance, when trained on MNLI, it gets a 67\% accuracy on the matched development set (compared to 84\% for BERT-base).

Part of our discussion relies on natural language inference, which has been widely studied in this area.  The classification task is to determine whether a \textit{hypothesis} statement is true (entailment), false (contradiction) or undetermined (neural) given a \textit{premise} statement. MNLI \citep{N18-1101} is the canonical large-scale English dataset to study this problem with 433K labeled examples. For evaluation, it features matched sets (examples from domains encountered in training) and mismatched sets (domains not-seen during training).

Experiments first examine qualitatively the spurious correlations picked up by the method. We then verify the validity of the method on a synthetic experimental setup. Finally, we verify the impact of our method by evaluating robust models on several out-of-distribution sets and discuss the choice of the weak learner.

\subsection{Weak learners rediscover previously reported dataset biases}
\label{sec:qualitative_examples}

\begin{wraptable}{R}{0.45\textwidth}
    \caption{Breakdown of the 1,000 top \textbf{certain / incorrect} training examples.}
    \label{tab:qualitative_breakdown}
    \centering
    \vspace{-5pt}
    \begin{tabular}{p{4.6cm}c}
        \toprule
        Category & (\%) \\
        \midrule
        Predicted \textit{Contradiction} & 46 \\
        ~~~~ Neg. in the hyp. & 43 \\
        Predicted \textit{Entailment} & 51 \\
        ~~~~ High word overlap prem./hyp. & 43 \\
        Predicted \textit{Neutral} & 3 \\
        \bottomrule
    \end{tabular}
    \vspace{-10pt}
\end{wraptable}

Most approaches for circumventing dataset bias require modeling the bias explicitly, for example using a model limited to only the hypothesis in NLI \citep{Gururangan2018AnnotationAI}. These approaches are effective, but require isolating specific biases present in a dataset. Since this process is costly, time consuming and error-prone, it is unrealistic to expect such analysis for all new datasets.  On the contrary, we hypothesize that weak learners might operate like \textit{rapid surface learners} \citep{Zellers2019HellaSwagCA}, picking up on dataset biases without specific signal or input curation and being rather certain of their biased errors (high certainty on the biased prediction errors).

We first investigate whether our weak learner re-discovers two well-known dataset biases reported on NLI benchmarks: (a) the presence of negative word in the hypothesis is highly correlated with the \textit{contradiction} label \citep{Poliak2018HypothesisOB,Gururangan2018AnnotationAI}, (b) high word overlap between the premise and the hypothesis is highly correlated with the \textit{entailment} label \citep{McCoy2019RightFT}.

To this aim, we fine-tune a weak learner on MNLI \citep{N18-1101}. Hyper-parameters can be found in Appendix~\ref{sec:hyperparams}. We extract and manually categorize 1,000 training examples wrongly predicted by the weak learner (with a \textbf{high loss} and a \textbf{high certainty}). Table~\ref{tab:qualitative_breakdown} breaks them down per category. Half of these incorrect examples are wrongly predicted as  \textit{Contradiction} and almost all of these contain a negation\footnote{We use the following list of negation words: \textit{no, not, none, nothing, never, aren't, isn't, weren't, neither, don't, didn't, doesn't, cannot, hasn't, won't}.} in the hypothesis. Another half of the examples are incorrectly predicted as \textit{Entailment}, a majority of these presenting a high lexical overlap between the premise and the hypothesis (5 or more words in common). The weak learner thus appears to predict with high-certainty a \textit{Contradiction} label whenever the hypothesis contains a negative word and with high-certainty an \textit{Entailment} label whenever there is a strong lexical overlap between premise/hypothesis. Table~\ref{tab:qualitative_biases} in Appendix~\ref{sec:qualitative_biases} presents qualitative examples of dataset biases identified by the fine-tuned weak learner.

This analysis is based on a set of biases referenced in the literature and does not exclude the possibility of other biases being detected by the weak learner. For instance, during this investigation we notice that the presence of ``negative sentiment'' words (for instance: \textit{dull, boring}) in the hypothesis appears to be often indicative of a \textit{Contradiction} prediction. We leave further investigation on such behaviors to future work.

\subsection{Synthetic experiment: cheating feature}

\begin{wrapfigure}{R}{0.5\textwidth}
  \vspace{-20pt}
  \begin{center}
    \includegraphics[width=0.48\textwidth]{synthetic_cheating.png}
  \end{center}
  \vspace{-10pt}
  \caption{Accuracy on MNLI matched development set for models with a cheating feature. The model trained with PoE (\textit{Main Model}) is less sensitive to this synthetic bias.}
  \label{fig:synthetic}
  \vspace{-15pt}
\end{wrapfigure}

We consider a controlled synthetic experiment described in \cite{He2019UnlearnDB,Clark2019DontTT} that simulates bias. We modify 20,000 MNLI training examples by injecting a \textit{cheating feature} which encodes an example's label with probability $p_{cheat}$ and a random label selected among the two incorrect labels otherwise. For simplicity, we consider the first 20,000 examples. On the evaluation sets, the cheating feature is random and does not convey any useful information. In the present experiment, the cheating feature takes the form of a prefix added to the hypothesis (``0'' for \textit{Contradiction}, ``1'' for \textit{Entailment}, ``2'' for \textit{Neutral}). We train the weak and main models on these 20,000 examples and evaluate their accuracy on the matched development set.\footnote{We observe similar trends on the mismatched development set.} We expect a biased model to rely mostly on the cheating feature thereby leading to poor evaluation performance.

Figure~\ref{fig:synthetic} shows the results. As the proportion of examples containing the bias increases, the evaluation accuracy of the weak learner quickly decreases to reach 0\% when $p_{cheat}=0.9$. The weak learner detects the cheating feature during training and is mainly relying on the synthetic bias which is not directly indicative of the gold label.

Both \citet{He2019UnlearnDB} and \citet{Clark2019DontTT} protect against the reliance on this cheating feature by ensembling the main model with a biased model that only uses the hypothesis (or its first token). We instead train the main model in the product of experts setting, relying on the weak learner to identify the bias. Figure~\ref{fig:synthetic} shows that when a majority of the training examples contain the bias ($p_{cheat} \geq 0.6$), the performance of the model trained with cross-entropy drops faster than the one trained in PoE. PoE training leads to a more robust model by encouraging it to learn from the mistakes of the weak learner. As $p_{cheat}$ comes close to 1, the model's training enters a ``few-shot regime'' where there are very few incorrectly predicted biased examples to learn from (examples where following the biased heuristic lead to a wrong answer) and the performance of the model trained with PoE drops as well.

\subsection{Adversarial datasets: NLI and QA}
\label{sec:adv_experiments}

\paragraph{NLI}

The HANS adversarial dataset~\citep{McCoy2019RightFT} was constructed by writing templates to generate examples with a high premise/hypothesis word overlap to attack models that rely on this bias. In one template the word overlap generates entailed premise/hypothesis pairs (\textit{heuristic-entailed examples}), whereas in another the examples contradict the heuristic (\textit{non-heuristic-entailed}). The dataset contains 30K evaluation examples equally split between both.

Table~\ref{tab:mnli_hans} shows that the weak learner exhibits medium performance on the in-distribution sets (MNLI) and that on out-of-distribution evaluation (HANS), it relies heavily on the word overlap heuristic. Product of experts training is effective at reducing the reliance on biases and leads to significant gains on the heuristic-non-entailed examples when compared to a model trained with standard cross-entropy leading to an improvement of +24\%.

The small degradation on in-distribution data is likely because product of experts training does not specialize for in-distribution performance but focuses on the weak model errors \citep{He2019UnlearnDB}. The linear combination of the original cross-entropy loss and the product of experts loss (PoE + CE) aims at counteracting this effect. This multi-loss objective trades off out-of-distribution generalization for in-distribution accuracy. A similar trade-off between accuracy and robustness has been reported in adversarial training \citep{Zhang2019TheoreticallyPT,Tsipras2019RobustnessMB}. In Appendix~\ref{sec:alpha_multi_loss}, we detail the influence of this multi-loss objective.

\begin{table*}
    \caption{MNLI matched dev accuracies, HANS accuracies and MNLI matched hard test set. All numbers are averaged on 6 runs (with standard deviations). Detailed results on HANS are given Appendix~\ref{sec:hans_details}. Reported results are indicated with *. $\clubsuit$ \cite{utama2020debiasing} is a concurrent work where they use a BERT-base fine-tuned on 2000 random examples from MNLI as a weak learner and ``PoE + An.'' refers to the annealing mechanism proposed by the authors.
    }
    \label{tab:mnli_hans}
    \centering
    \begin{tabular}{lc|ccc|c}
        \toprule
         & Loss & MNLI & \multicolumn{2}{c}{HANS} & Hard \\
         & & & Ent & Non-Ent & \\
        \midrule
        \cite{Clark2019DontTT}* & PoE & 82.97 & 64.67 & 71.16 & - \\
        \cite{Mahabadi2020EndtoEndBM}* & PoE & 84.19 & 95.99 & 33.30 & 76.81  \\
        \cite{utama2020debiasing}*$\clubsuit$  & PoE & 80.70 & 86.13 & 55.20 & - \\
        \cite{utama2020debiasing}*$\clubsuit$  & PoE + An. & 81.90 & 88.40 & 47.13 & - \\
        \midrule
        BERT-base & CE & \textbf{84.52}$\pm$0.27 & 98.12$\pm$0.62 & 26.74$\pm$6.15 & 76.96$\pm$0.38 \\
        TinyBERT - Weak & CE & 66.93$\pm$0.12 & \textbf{99.80}$\pm$0.09 &  0.44$\pm$0.26 & 46.65$\pm$0.48 \\
        \midrule
        BERT-base - Main & PoE & 81.35$\pm$0.40 & 81.13$\pm$8.1 &  \textbf{56.41}$\pm$5.91 & 76.54$\pm$0.56 \\
        BERT-base - Main & PoE + CE & 83.32$\pm$0.24 & 94.51$\pm$0.82 &  41.35$\pm$8.25 & \textbf{77.63}$\pm$0.49 \\
        \bottomrule
    \end{tabular}
\end{table*}

We also evaluate our method on MNLI's hard test set \citep{Gururangan2018AnnotationAI} which is expected to be less biased than MNLI's standard split. These examples are selected such that a hypothesis-only model cannot predict the label accurately. Table~\ref{tab:mnli_hans} shows the results of this experiment. Our method surpasses the performance of a PoE model trained with a hypothesis-only biased learner. Results on the mismatched set are given in Appendix~\ref{sec:hans_details}.

\paragraph{QA}

Question answering models often rely on heuristics such as type and keyword-matching \citep{Weissenborn2017MakingNQ} that can do well on benchmarks like SQuAD \citep{Rajpurkar2016SQuAD10}. We evaluate on the Adversarial SQuAD dataset \citep{Jia2017AdversarialEF} built by appending distractor sentences to the passages in the original SQuAD. Distractors are constructed such that they look like a plausible answer to the question while not changing the correct answer or misleading humans.

Results on SQuAD v1.1 and Adversarial SQuAD are listed in Table~\ref{tab:adv_qa}. The weak learner alone has low performance both on in-distribution and adversarial sets. PoE training improves the adversarial performance (+1\% on AddSent) while sacrificing some in-distribution performance. A multi-loss optimization closes the gap and even boosts adversarial robustness (+3\% on AddSent and +2\% on AddOneSent). In contrast to our experiments on MNLI/HANS, multi-loss training thus leads here to better performance on out-of-distribution as well. We hypothesize that in this dataset, the weak learner picks up more useful information and removing it entirely might be non-optimal. Multi-loss in this case allows us to strike a balance between learning from, or removing, the weak learner.

\begin{table*}
    \caption{F1 Scores on SQuAD and Adversarial QA. The AddOneSent set is model agnostic while we use the AddSent set obtained using an ensemble of BiDAF models \citep{Seo2017BidirectionalAF}. * are reported results.}
    \label{tab:adv_qa}
    \centering
    \begin{tabular}{lc|ccc}
        \toprule
         & Loss & SQuAD & \multicolumn{2}{c}{Adversarial QA} \\
         & & & AddSent & AddOneSent \\
        \midrule
        \multirow{2}{*}{\shortstack{\cite{Clark2019DontTT}*\\BiDAF}} & CE & 80.61 & 42.54 & 53.91 \\
        & PoE & 78.63 & 57.64 & 57.17 \\
        \midrule
        BERT-base & CE & \textbf{88.68} & 53.98 & 58.84 \\
        TinyBERT - Weak & CE & 41.08 & 16.02 & 18.63 \\
        \midrule
        BERT-base - Main & PoE & 83.11 & 54.92  & 58.44 \\
        BERT-base - Main & PoE + CE & 86.49 & \textbf{56.80} & \textbf{61.04} \\
        \bottomrule
    \end{tabular}
\end{table*}

\section{Analysis}

\subsection{Reducing bias: correlation analysis}
\label{sec:correlations}

\begin{wrapfigure}{R}{0.52\textwidth}
  \vspace{-20pt}
  \begin{center}
    \includegraphics[width=0.45\textwidth]{correlations.png}
  \end{center}
  \vspace{-10pt}
  \caption{Pearson correlation between the losses (on evaluation sets) of the biased model and different training methods. The PoE training is effective at reducing the correlation with the biased model.}
  \label{fig:correlations}
  \vspace{-5pt}
\end{wrapfigure}

To investigate the behavior of the ensemble of the weak and main learner, we compute the Pearson correlation between the element-wise loss of the weak (biased) learner and the loss of the trained models following \citet{Mahabadi2020EndtoEndBM}. A correlation of 1 indicates a linear relation between the per-example losses (the two learners make the same mistakes), and 0 indicates the absence of linear correlation (models' mistakes are uncorrelated). Figure~\ref{fig:correlations} shows that models trained with a linear combination of the PoE cross-entropy and the standard cross-entropy have a higher correlation than when trained solely with PoE. This confirms that PoE training is effective at reducing biases uncovered by the weak learner and re-emphasizes that adding standard cross-entropy leads to a trade-off between the two.

\subsection{How weak do the weak learners need to be?}

We consider parameter size as a measure of the capacity or "weakness" of the weak learner. We fine-tune different sizes of BERT~\citep{Turc2019WellReadSL} ranging from 4.4 to 41.4 million parameters and use these as weak models in a PoE setting. Figure~\ref{fig:poe_weak_analysis} shows the accuracies on MNLI and HANS of the weak learners and main models trained with various weak learners.

Varying the capacity of the weak models affects both in-distribution and out-of-distribution performance. Out-of-distribution performance of the main model increases as the weak model becomes stronger (more parameters) up to a certain point while in-distribution performances drop slightly at first and then more strongly. When trained jointly with the larger MediumBERT weak learner (41.4 million parameters), the main model gets 97\% accuracy on HANS's heuristic-non-entailed set but a very low accuracy on the in-distribution examples (28\% on MNLI and 3\% on the heuristic-entailed examples).

\begin{figure}[t]
  \centering
  \subfloat[Performance of weak learners (CE)]{
    \includegraphics[width=0.5\columnwidth]{weak_analysis.png}
    \label{fig:weak_analysis}
  }
  \subfloat[Performance of main models (PoE)]{
    \includegraphics[width=0.5\columnwidth]{poe_weak_analysis.png}
    \label{fig:poe_weak_analysis}
  }
  \caption{Weaker learners assure a good balance between out-of-distribution and in-distribution while stronger learners encourage out-of-distribution generalization at the expense of in-distribution performance. We indicate the number of parameters in parenthesis (in millions).}
  \label{fig:weak_to_strong}
\end{figure}

As a weak model grows in capacity, it becomes a better learner. The average loss decreases and the model becomes more confident in its predictions. As a result, the group \textbf{certain / correct} becomes more populated and the main model receives on average a smaller gradient magnitude per input. On the contrary, the \textbf{certain / incorrect} group (which generally align with out-of-distribution samples and induce higher magnitude gradient updates, encouraging generalization at the expense of in-distribution performance) becomes less populated. These results corroborate and complement insights from \cite{Yaghoobzadeh2019RobustNL}. This is also reminiscent of findings from \citet{Vodrahalli2018AreAT} and \citet{Shrivastava2016TrainingRO}: not all training samples contribute equally towards learning and in some cases, a carefully selected subset of the training set is sufficient to match (or surpass) the performance on the whole set.

We present complementary analyses in Appendix. To further show the effectiveness of our method, we included in Appendix~\ref{sec:fever} an additional experiment on facts verification \citep{Thorne2018FEVERAL,Schuster2019TowardsDF}. In Appendix~\ref{sec:transfer_nli}, we evaluate the ability of our method to generalize to other domains that do not share the same annotation artifacts. We highlight the trade-off between in-distribution performance and out-of-distribution robustness by quantifying the influence of the multi-loss objective in Appendix~\ref{sec:alpha_multi_loss} and draw a connection between our 3 groups of examples and recently introduced \textit{Data Maps} \citep{Swayamdipta2020DatasetCM}.

\section{Conclusion}
We have presented an effective method for training models robust to dataset biases. Leveraging a weak learner with limited capacity and a modified product of experts training setup, we show that dataset biases do not need to be explicitly known or modeled to be able to train models that can generalize significantly better to out-of-distribution examples. We discuss the design choices for such weak learner and investigate how using higher-capacity learners leads to higher out-of-distribution performance and a trade-off with in-distribution performance. We believe that such approaches capable of automatically identifying and mitigating datasets bias will be essential tools for future bias-discovery and mitigation techniques.


\section{Introduction}

This note is devoted to a quantitative analysis of the small correlations in chaotic states of a rarefied gas.
Particles have a little ``core'' of size $\varepsilon$, meaning that their interaction is repulsive at the origin, and negligible on scales larger than $\varepsilon$.
The number of particles is increased with $\varepsilon$ in such a way that, letting the system evolve, we would see
in average a finite number of interactions per unit of time, and a variation of the macroscopic density on
distances of order one (order of the mean free path).

This is known as Boltzmann-Grad regime. Its interest
arises in connection with validation issues in kinetic theory, as it is the regime in which the Boltzmann
equation has been proved to hold rigorously, at least for short range potentials and small enough times;
see \cite{La75} for a foundational paper and \cite{S2,CIP94,GG20, PS20} for surveys.
In such works, the famous ``Stosszahlansatz'' of Boltzmann is ensured by the randomness of the initial conditions: one starts with a {\em chaotic} measure such that the joint $j$-particle distributions factorize as $\varepsilon\to 0$, and shows that the factorization is propagated at positive times.

It is difficult to establish the minimal set of probabilistic hypotheses that allow the result to hold. Actually, one gives sufficient conditions, and verifies at the end that these conditions are valid for simple microscopic distributions. Therefore in practice, one starts with measures proportional to $f^{\otimes n} e^{-\beta U^\varepsilon_n}$, $\beta>0$ where $f$ is an assigned (nonequilibrium) $1-$particle density, and $U^{\varepsilon}_n$ is the potential energy of the $n$-particle configuration. As $U^\varepsilon_n$ is concentrated on configurations of particles having mutual distance $\varepsilon$, this measure can be regarded as an almost Poissonian, or ``maximally factorized'' state. The only source of correlation comes from the repulsion at short scale.

Our aim here is to quantify such correlations in detail,
and especially how they decay in space, passing from microscopic to macroscopic scales.
In order to do that,
the most powerful tool at hand is provided by cumulants, or {\em truncated} correlation functions. Roughly, the $j$th cumulant measures the collective correlation in a cluster of $j$ particles.
This is a standard notion in statistical mechanics;
see for instance \cite{DIS73,DIS74} where
strong bounds on truncated functions are used to prove analyticity of thermodynamic functionals; or \cite{D-CGR20} for a recent investigation of the Ising model.

In the kinetic limit of nonequilibrium processes,
truncated functions have been carefully studied over recent times.
It has been argued that, as $\varepsilon\to 0$,
the full collection of cumulants retains {\em all} the information about the correlations \cite{BGSS1} (see also \cite{GG18,LMN18}); in particular, $L^1$ estimates on cumulants have been used to provide a rigorous theory of fluctuations and large deviations of a hard sphere gas \cite{BGSS2}\footnote{It seems however hard to obtain estimates in more accurate norms.}.

The latter work, and most of the rigorous literature on the Boltzmann-Grad limit, focus on the simplest case of hard-sphere interactions (see \cite{GSRT,PSS13,Ay17} for exceptions). Nevertheless the main results are expected to hold as well for smooth short range interactions (for which they are proved in some cases), and even for power laws with fast enough decay (for which no result is available). Notice that, in the particular case of hard spheres, the maximally chaotic state (MCS from now on) is a pure product up to the excluded volume; while for long range interactions the product structure is perturbed by the tails.

In this paper we study grand canonical, maximally chaotic distributions for a general class of pair interactions. We shall not be concerned with the dynamical problem, but only with the MCS itself. For such states, we prove decay of the $j$th truncated function, for arbitrary $j$, in terms of a suitable quantity measuring the internal spatial separation in a cluster of $j$ particles. This translates into an exponential decay in $\varepsilon$ at macroscopic distances in the case of compactly supported interactions.

Besides the role of ``typical'' initial state in the rigorous derivation of the Boltzmann equation, we conclude by mentioning yet another application of MCS to the dynamical problem. Of course the dynamics prevents the MCS structure to be propagated in time. Nevertheless, these states have been proved useful as an argument to obtain finite density corrections to the Boltzmann equation \cite{Re78}. In the previous reference it is shown (restricting to hard spheres) that, if the ansatz of maximal chaoticity is made {\em at any time} (thus disregarding dynamical correlations), then the density evolves according to an Enskog equation, and an $H-$Theorem can be deduced for the associated entropy; see also \cite{LPS67,vBE73}. These remarks have not been object of rigorous investigation so far; see however \cite{GG12,PS17, PST18} for related discussions.

The rest of this work is organized in two parts. In Section \ref{sec:states} we set the assumptions on the interaction, introduce the maximally chaotic states, and state our main result (Theorem \ref{thm: bound_rtcf} below). Section \ref{sec:proof} is devoted to the proof. This combines classical cluster expansion (\cite{Ru69}) with estimates on long chains of interacting particles. Finally in the appendix, we report the proof of a crucial combinatorial argument from \cite{PU09}.

\section{Estimates on truncated correlations}
\label{sec:states}

\subsection{Potentials and kinetic limit}

We consider a system of identical point particles placed in an open connected set $\Lambda \subset \mathbb{R}^{d}$, $d \geq 1$ and we denote by
\begin{equation*}
z_i = (x_i,v_i)\in\Lambda \times \mathbb{R}^{d} \;,
\qquad
i \in\mathbb{N} \;,
\end{equation*}
the configuration of the $i$-th particle having position $x_i$ and velocity $v_{i}$.
With
\begin{equation*}
{\bm z}_j = (z_1,\ldots,z_j) \;,
\qquad
j\in\mathbb{N}
\end{equation*}
we indicate the configuration of the first $j$ particles.

The grand canonical phase space is
\begin{equation}
\Omega = \bigcup_{n\geq0} \Omega_{n} \;, \nonumber
\end{equation}
where the canonical $n$-particle phase space is
\begin{equation}
\Omega_{n} = \left\{
{\bm z}_{n}\in\left(\Lambda \times \mathbb{R}^{d}\right)^n\mid x_i \neq x_k\;,\  i,k=1,\ldots,n,\  i\neq k
\right\} \;.\nonumber
\end{equation}

The particles interact through a translation invariant pair potential given by a piecewise continuous function
$\phi: \mathbb{R}^{d}\rightarrow\mathbb{R} \cup \{\infty\}$.
As customary in statistical mechanics \cite{Ru69}, we assume:
\begin{itemize}
\item \textit{stability}: there exists a constant $B\geq0$ such that
\begin{equation}
\sum_{1\leq i<k \leq n}\phi(x_i - x_k) \geq -Bn \;,
\label{eq: phi_stab}
\end{equation}
for all $n\geq0$ and $x_1,\ldots,x_n\in\mathbb{R}^{d}$;
\item \textit{decay}:
\begin{equation}
C_{\beta}:= \int_{\mathbb{R}^{d}}\text{d}x \left|\zeta(x,0)\right| < + \infty \;,\qquad \beta >0
\label{eq: phi_reg}
\end{equation}
where $\zeta(x_{i},x_{k}):=\text{e}^{-\beta\phi(x_{i}-x_{k})} - 1 $.
\end{itemize}

Stability is usually required to ensure the existence of a partition function.
In our setting, the stability property is motivated by \cite{PSS13}, where the validity theorem for the Boltzmann equation has been extended to forces including an attractive part provided that condition \eqref{eq: phi_stab} holds.

Since the potential is bounded below, \eqref{eq: phi_reg} corresponds to a fast enough decrease at infinity.
Although we will be actually interested in finite range interactions (for which \eqref{eq: phi_reg} is trivially true), our result could be extended to any potential satisfying the above decay assumption.

We are interested here in a low density regime where particles interact strongly at a short distance $\varepsilon>0$. We will use then the potential in the rescaled form $ \phi\left( x / \varepsilon \right)$. In the Boltzmann-Grad scaling of a constant mean free path, the average number of particles should diverge
as
$$
\mu_\varepsilon := \varepsilon^{-(d-1)}\;.
$$
We introduce next a class of measures describing this asymptotic.

\subsection{Grand canonical chaotic states}

By definition, a {\em state} of the system is
an absolutely continuous measure on $\Omega$
admitting a set of densities $\{W^{\varepsilon}_{n}\}_{n\geq0}$, where
$W^{\varepsilon}_{n}: \Omega_{n}\rightarrow\mathbb{R}^{+}$
are positive Borel functions. For each $n$, $W^{\varepsilon}_{n}$ is assumed to be invariant under permutations of the particle labels.
The quantity
$(1/n!)W^{\varepsilon}_{n}({\bm z}_{n})$
represents the probability density of finding $n$ particles in the configuration
${\bm z}_{n}=(z_{1},\ldots,z_{n})$. In particular the distribution of the number of particles
is $(1/n!)\int  W^{\varepsilon}_{n}$ and the
normalization condition reads
\begin{equation*}
\sum_{n=0}^{\infty}\frac{1}{n!}
\int_{\Omega_n}W^{\varepsilon}_{n}({\bm z}_n)\,\text{d}{\bm z}_{n} = 1\;.
\end{equation*}

The collection of {\em correlation functions}
$\{\rho^{\varepsilon}_{j}\}_{j\geq0}$ with $\rho^{\varepsilon}_{j}:\Omega_{j}\rightarrow\mathbb{R}^{+}$ is given by
\begin{equation}
\label{eq: cf_def}
\rho_{j}^{\varepsilon}({\bm z}_{j}) :=
\sum_{n=0}^{\infty}\frac{1}{n!}
\int_{\Omega_n}\
W^{\varepsilon}_{j+n}({\bm z}_{j+n})\, \text{d}z_{j+1}\ldots \text{d}z_{j+n} \;.
\end{equation}
Note that $\rho_{0}^{\varepsilon} = 1$.
We say that the state admits correlation functions when the series in the right hand side is convergent, together with the series in the inverse formula
$$
W^{\varepsilon}_{j}({\bm z}_{j}) =
\sum_{n=0}^{\infty}\frac{(-1)^n}{n!}
\int_{\Omega_n}\
\rho^{\varepsilon}_{j+n}({\bm z}_{j+n})\, \text{d}z_{j+1}\ldots \text{d}z_{j+n} \;.
$$
In this case the correlation functions are an alternative description of all the statistical properties of the system.
After suitable normalization, the function $\rho^{\varepsilon}_{j}$ provides the joint distribution of $j$ distinct particles.

Since the average number of particles is equal to the integral of $\rho_{1}^{\varepsilon}$,
the {\em Boltzmann-Grad} condition can be formulated by requiring that
\begin{equation*}
\lim_{\varepsilon\rightarrow0}\,\mu_\varepsilon^{-1}
\int_{\Omega_1} \rho^{\varepsilon}_{1}(z)\, \text{d}z = 1\;.
\label{eq: BGlim}
\end{equation*}
Therefore the correlation functions $\rho_{j}^{\varepsilon}$ will diverge as $\mu_\varepsilon^{j}$.

In order to state the aforementioned validity theorem for the Boltzmann equation (\cite{La75,GSRT,PSS13}), it is essential to assume that the rescaled correlation functions $\mu_\varepsilon^{-j}\rho_{j}^{\varepsilon}$ are bounded, uniformly in $\varepsilon$, by the correlation functions of an equilibrium state.
In fact, the invariance of the energy is used to rule out the emergence of too strong correlations.
Moreover, a strong chaos property is assumed to the extent that
\begin{equation}\label{eq:PC}
\lim_{\varepsilon\rightarrow0} \mu_\varepsilon^{-j}\rho_{j}^{\varepsilon}  = f^{\otimes j}\;,
\end{equation}
uniformly on compact sets of $\Omega_j$. Here $f^{\otimes j}({\bm z}_{j}) =  f(z_1)f(z_2)\ldots f(z_j)$ and $f$ is a continuous probability density on $\Omega_1$, playing the role of the Boltzmann density at a given time.

To fulfill such hypotheses, one resorts to the following grand canonical nonequilibrium prescription.
\begin{defi} \label{def:MCS}
Let $\varepsilon>0$ and $f$ be a continuous probability density on $\Omega_1$ with spatial density $\rho(x):= \int_{\mathbb{R}^{d}} \ f(x,v)\text{d}v$ such that
$
\bar\rho := \|\rho\|_{L^\infty} < \infty \;.
$
A ``maximally chaotic'' state (MCS) is a grand canonical measure
with densities
\begin{equation}
\frac{1}{n!}W^{\varepsilon}_{n}({\bm z}_n):=\frac{1}{{\cal Z}_\varepsilon} \frac{\mu_\varepsilon^n}{n!}
f^{\otimes n}({\bm z}_n) \,\psi^\varepsilon_n({\bm x}_n) \;,  \qquad n \geq 0
\label{eq: Wn}
\end{equation}
where ${\cal Z}_\varepsilon$ is the partition function
\begin{equation} \label{eq:FP}
{\cal Z}_\varepsilon:=\sum_{n\geq 0} \frac{\mu_{\varepsilon}^n}{n!}\int_{\Omega_{n}} \
f^{\otimes n}({\bm z}_n)\,\psi^\varepsilon_n({\bm x}_n) \text{d}{\bm z}_n
\end{equation}
and $\{\psi^\varepsilon_n\}_{n \geq 0}$ are the Boltzmann factors
\begin{equation}
\psi^\varepsilon_n({\bm x}_n) :=
\prod_{1\leq \alpha < \alpha' \leq n}
\text{e}^{ -\beta\phi(\frac{x_\alpha-x_{\alpha'}}{\varepsilon})}\;, \qquad \beta > 0
\label{eq: Boltzmann}
\end{equation}
$(\psi^\varepsilon_0=\psi^\varepsilon_1=1)$.
\end{defi}

The equilibrium is recovered by the choice $f(x,v) = (z/|\Lambda|)M_{\beta}(v)$, where $z>0$ and $M_\beta$ is the normalized Maxwellian with inverse temperature $\beta$.
Hence Definition \ref{def:MCS} resembles an equilibrium distribution for particles interacting by means of $\phi$ in an external field $-\beta^{-1}\log \rho( x)$. This allows to apply known cluster expansion methods to analyze its properties in the limit $\varepsilon \to 0$.

\subsection{Main result}

We first introduce several notations.

\begin{defi}
(i) $\mathcal{G}_{k}$ is the set  of  all graphs $G$ (unoriented, with no loops and no multiple edges) with vertices $V(G)=\{1,\ldots, k\}$ and edges $E(G) \subset\{ \{\alpha,\alpha'\}\ |\ \alpha,\alpha'\in V(G)\,, \alpha \neq \alpha'\}$.\\
(ii) ${\cal C}_{k}\subset \mathcal{G}_{k}$ is the set of all connected graphs on $\{1,\ldots, k\}$, i.e.\,graphs $G$ such that
for any $\alpha,\alpha'\in V(G)$ with $\alpha \neq \alpha'$ there exist $\alpha_1,\cdots,\alpha_\ell$ with $\{\alpha,\alpha_1\}, \{\alpha_1,\alpha_2\},\cdots\{\alpha_\ell,\alpha'\}$ in $E(G)$.\\
(iii) ${\cal T}_{k}\subset{\cal C}_{k}$ is the set of trees on $\{1,\ldots, k\}$, i.e.\,connected graphs with no cycles.
\end{defi}

We shall adopt the following abbreviations. For $S=\{i_{1},i_2,\ldots, i_{s}\}$ a set of indices of cardinality $s = |S|$, and $g_s$ a generic function of $s$ variables,
we write
\begin{eqnarray}
&& {\bm y}_S = (y_{i_1},\cdots,y_{i_s})\;,\nonumber\\
&&g_S= g_{s} ({\bm y}_S)=g_{s}(y_{i_1}, \cdots, y_{i_s})\;.\nonumber
\\
&& g^{ \otimes S}= g(y_{i_1})g(y_{i_2}) \cdots  g(y_{i_s})\;.
\label{eq:notintro}
\end{eqnarray}
Moreover, we indicate by $(S_{1},\ldots,S_k)_{S}$
a generic non trivial partition of $S$ into subsets $S_{1},\ldots,S_{k}$.
We shall use later on the superscript zero, $(S_{1},\ldots,S_k)_{S}^{0}$,
for partitions of $S$ into possibly trivial subsets.

\begin{defi}
The truncated correlation (or cluster, cumulant) functions $\{\rho^{\varepsilon,T}_{j}\}_{j\geq0}$
are defined recursively on $\Omega_j$ by $\rho^{\varepsilon,T}_0 = 0$,
\begin{equation}
\begin{cases}
\displaystyle
\rho_{1}^{\varepsilon,T}:= \rho_{1}^{\varepsilon} \; , \\
\displaystyle
\rho_{J}^{\varepsilon,T} := \rho_{J}^{\varepsilon} -
\sum_{k=2}^{j}\sum_{(J_{1},\ldots,J_{k})_{J}} \prod_{i=1}^{k} \rho_{J_i}^{\varepsilon,T} \;,
&\qquad j\geq2 \;.
\end{cases}
\label{eq: rtcf_rec}
\end{equation}
\end{defi}

Our main result is that truncated correlations behave like ``stretched trees''. More precisely, their magnitude is controlled by the length of a tree graph constructed on the space configuration $(x_1,\cdots,x_j)$ and measuring the collective separation between particles; which we define next.

Given a tree $T\in {\cal T}_j$ and a configuration  of particles ${\bm x}_j$, the {\em length} of $T$ on ${\bm x}_j$ is
\begin{equation}\label{eq:LLTj}
{\cal L}_T({\bm x}_j) := \sum_{\{\alpha,\alpha'\}\in E(T)} |x_\alpha - x_{\alpha'}|\;.
\end{equation}
The length of the cluster ${\cal L}({\bm x}_{j})$ is the minimal length of all tree graphs connecting the vertices
$x_1,\ldots x_j$, and possibly other vertices:
\begin{equation} \label{eq:LLj}
{\cal L}({\bm x}_j) := \min_{n} \min_{T \in {\cal T}_{j+n}} \min_{x_{j+1},\cdots,x_{j+n}} {\cal L}_T({\bm x}_{j+n})\;.
\end{equation}

Due to the scaling \eqref{eq:PC}, we are interested in estimating the truncated functions multiplied by the factor $\mu_\varepsilon^{-j}$.
We then obtain the following result.

\begin{thm}
\label{thm: bound_rtcf}
Suppose that the pair potential satisfies \eqref{eq: phi_stab} and has support $\{|x|<1/2\}$.
There exist positive constants $\varepsilon_0=\varepsilon_0(\phi,\bar\rho,\beta,B)$ and $A$ such that the truncated
correlation functions of the maximally chaotic state \eqref{eq: Wn} satisfy the bounds:
\begin{equation}
\frac{|\rho^{\varepsilon,T}_j({\bm z}_j)|}{\mu_\varepsilon^j} \leq \left(A f\right)^{\otimes j}\
 j^{j-2} \frac{\left(\varepsilon /\varepsilon_0\right)^{\left( {\cal L} \left({\bm x}_{j}\right) / \varepsilon - j+1 \right)_+}}{1-\left(\varepsilon /\varepsilon_0\right)}
 \;,
\label{eq: rtcf_est}
\end{equation}
for all $j$, ${\bm z}_j \in \Omega_j$ and $\varepsilon < \varepsilon_0$.
\end{thm}

In particular if the particles are at macroscopic distance, the minimal distance between them is $a > 0$ and the estimate provides the exponential decay:
$$
\frac{|\rho^{\varepsilon,T}_j({\bm z}_j)|}{\mu_\varepsilon^j} \leq \left(A  f\right)^{\otimes j} j^{j-2} e^{-(j-1)a' \varepsilon^{-1} \log\varepsilon^{-1}}
$$
for some $a'>0$.

Similar results have been provided in \cite{DIS73,DIS74}, under the name of `strong cluster property', for the thermodynamic limit of equilibrium states.

\section{Proof of Theorem \ref{thm: bound_rtcf}} \label{sec:proof}

The correlation functions \eqref{eq: cf_def} of an MCS (see Def.\,\ref{def:MCS}) are given
by the following explicit expression:
\begin{equation} \label{eq:exprcf}
\frac{\rho_{j}^{\varepsilon}({\bm z}_{j})}{\mu_{\varepsilon}^{j}} =
\frac{f^{\otimes j}({\bm z}_j)}{{\cal Z}_\varepsilon}
\sum_{n=0}^{\infty}\frac{\mu_\varepsilon^n}{n!}
\int_{\Lambda^{n}}
\psi^\varepsilon_{j+n}({\bm x}_{j+n})\,
\text{d}\rho(x_{j+1})\cdots \text{d}\rho(x_{j+n})
\end{equation}
where
$$\text{d}\rho(x) = \rho(x)\text{d}x\;.$$
We introduce a rescaled decay function as in Eq.\,\eqref{eq: phi_reg} by
\begin{equation} \label{eq:rdf}
\zeta^\varepsilon(x_{i},x_{k}) := \text{e}^{-\beta\phi(\frac{x_{i}-x_{k}}{\varepsilon})} - 1 \;.
\end{equation}
Expanding the Boltzmann factor
\begin{equation*}
\psi^\varepsilon_k({\bm x}_k)  =
\prod_{1\leq \alpha < \alpha' \leq k} \left(1 + \zeta^\varepsilon(x_{\alpha},x_{\alpha'})\right)
\label{eq: Boltzmann}
\end{equation*}
one finds a sum over all possible graphs on $k$ vertices:
\begin{equation*}
\psi^\varepsilon_k({\bm x}_k) =
\sum_{G\in\mathcal{G}_{k}}\,\prod_{\{\alpha,\alpha'\}\in E(G)}\zeta^\varepsilon(x_\alpha,x_{\alpha'}) \;.
\end{equation*}

Truncated functions have a similar expression, the main difference being that generic graphs
are replaced by connected graphs. We introduce Ursell functions $\{u^\varepsilon_k\}_{k \geq 0}$ by $u^\varepsilon_0=0$, $u^\varepsilon_1=1$ and
\begin{equation*}
u^\varepsilon_k({\bm x}_k) :=
\sum_{G \in {\cal C}_{k}} \prod_{\{\alpha,\alpha'\}\in E(G)}\zeta^\varepsilon(x_\alpha,x_{\alpha'}) \;,\qquad k > 1\;.
\label{eq: Ursell}
\end{equation*}
Then we have the following classical result.
\begin{lem}\label{eq: Zce*}
Suppose that conditions \eqref{eq: phi_stab} and \eqref{eq: phi_reg} are satisfied.
Then the MCS has partition function
\begin{equation}
\mathcal{Z}_{\varepsilon} =
\exp{\left\{\mu_\varepsilon\sum_{m=1}^{\infty}
\frac{\mu_\varepsilon^{m-1}}{m!}
\int\,u^\varepsilon_m({\bm x}_m)\,\text{d}\rho^{\otimes m}({\bm x}_{m}) \right\}} \;,
\label{eq: Zce}
\end{equation}
and the truncated correlation functions are given by
\begin{equation}
\frac{\rho_{j}^{\varepsilon,T}({\bm z}_{j})}{ \mu_{\varepsilon}^{j}} =
f^{\otimes j}({\bm z}_j)
\sum_{n=0}^{\infty}\frac{\mu_{\varepsilon}^{n}}{n!}
\int_{\Lambda^n}
u^\varepsilon_{j+n}({\bm x}_{j+n})
\,\emph{d}\rho(x_{j+1})\ldots \emph{d}\rho(x_{j+n}) \ \;,
\label{eq: rtcf}
\end{equation}
where the series and integrals are absolutely convergent, uniformly for
\begin{equation} \label{eq:radconv}
\varepsilon<  \frac{1}{\bar\rho\, C_\beta e^{2\beta B+1}}\;.
\end{equation}
\label{prop: prop_ce}
\end{lem}

The result follows directly from the following subtle estimate, known as {\em tree-graph inequality}.
\begin{lem} \label{eq:treein}
Suppose that condition \eqref{eq: phi_stab} is satisfied.
Then one has
\begin{equation}
\left|\sum_{G \in {\cal C}_{k}} \prod_{\{\alpha,\alpha'\}\in E(G)}\zeta^\varepsilon(x_\alpha,x_{\alpha'})\right| \leq
\emph{e}^{2k\beta B}\sum_{T \in {\cal T}_{k}} \prod_{\{\alpha,\alpha'\}\in E(T)}\left|\zeta^\varepsilon(x_\alpha,x_{\alpha'})\right| \;.
\label{eq: trees}
\end{equation}
\end{lem}

This bound goes back to work of Penrose for positive interactions \cite{Pe67}, later extended to stable interactions \cite{PU09} (see also \cite{Ue17} for an improvement). Due to its importance, we report the proof in appendix.

\bigskip
{\em Proof of Lemma \ref{eq: Zce*}.} Given Lemma \ref{eq:treein}, the proof reduces to simple combinatorics.
Since
$$
\psi^\varepsilon_n({\bm x}_n) = \sum_{k=1}^{n} \sum_{(J_{1},\ldots, J_{k})_{J}}
\prod_{i=1}^{k}
u^\varepsilon_{J_{i}}
$$
where $J = \{1,\cdots,j\}$ (and using the notations introduced in \eqref{eq:notintro})
one gets from \eqref{eq:FP} that
$$
\mathcal{Z}_{\varepsilon} = 1 + \sum_{n\geq1}
\sum_{k=1}^{n}\frac{1}{k!}
\sum_{\substack{m_1,\ldots,m_k >0 \\ m_1 + \ldots + m_k = n}}
\prod_{i=1}^{k}
\left(\frac{\mu_\varepsilon^{m_i}}{m_i!}
\int_{\Lambda^{m_i}}u^\varepsilon_{m_i}({\bm x}_{m_i})\,\text{d}\rho^{\otimes m_{i}}({\bm x}_{m_i})  \right)  \;,
$$
from which interchanging the sums we find
\begin{align}
\mathcal{Z}_{\varepsilon} & = 1 +
\sum_{k=1}^{n}\frac{1}{k!}
\left(
\sum_{m=1}^{\infty}
\frac{\mu_\varepsilon^m}{m!}
\int_{\Lambda^{m}}u^\varepsilon_{m}({\bm x}_{m})\,\text{d}\rho^{\otimes m}({\bm x}_{m}) \right)^k \;.
\nonumber
\end{align}
Thus \eqref{eq: Zce} holds.

Similarly, given a partition $(J_{1},\ldots,J_{k})_J$ of $J=\{1,\ldots,j\}$,
a subset $I_0\subset J^c=\{j+1,\ldots,j+n\}$ and a partition of $J^{c}\setminus I_0$
into (possibly trivial) subsets $I_{1},\ldots,I_{k}$,  we expand $\psi^\varepsilon_{j+n}$ in \eqref{eq:exprcf}
into corresponding connected components
$$
\psi^\varepsilon_{j+n}({\bm x}_{j+n}) =
\sum_{I_0\subset J^c} \psi^{\varepsilon}_{I_{0}}\,
\sum_{k=1}^{j}
\sum_{(J_{1},\ldots,J_{k})_{J}}
\sum_{(I_{1},\ldots,I_{k})^{0}_{J^c\setminus I_0}}
\prod_{i=1}^{k}u^\varepsilon_{J_i \cup I_i}
$$
and obtain
\begin{align*}
\frac{\rho_{j}^{\varepsilon}({\bm z}_{j})}{\mu_{\varepsilon}^{j}} & =
\frac{f^{\otimes j}({\bm z}_j)}{{\cal Z}_\varepsilon}
\sum_{k=1}^{j}
\sum_{(J_{1},\ldots,J_{k})_{J}}
\sum_{n=0}^{\infty}\frac{\mu_\varepsilon^n}{n!}
\sum_{(I_0,I_{1},\ldots,I_{k})^{0}_{J^c}}
\left(\int_{\Lambda^{i_0}} \psi^{\varepsilon}_{I_{0}}\, \text{d}\rho^{\otimes I_0}\right)
\prod_{\ell=1}^{k}\left(\int_{\Lambda^{i_\ell}}u^\varepsilon_{J_\ell\cup I_\ell} \, \text{d}\rho^{\otimes I_\ell} \right)
\nonumber \\ & =
\frac{f^{\otimes j}({\bm z}_j)}{{\cal Z}_\varepsilon}
\sum_{k=1}^{j}
\sum_{(J_{1},\ldots,J_{k})_{J}}
\sum_{n=0}^{\infty}
\sum_{\substack{m_0+\ldots+m_{k} = n \\ m_0,\ldots,m_{k}\geq0}}
\left(\frac{\mu_\varepsilon^{m_0}}{m_0!}
\int_{\Lambda^{m_0}} \psi^{\varepsilon}_{m_{0}}({\bm y}_{m_0})\, \text{d}\rho^{\otimes m_0}({\bm y}_{m_0})\right)
\nonumber \\ &
\qquad\qquad\qquad\qquad\qquad
\prod_{\ell=1}^{k}\left(\frac{\mu_\varepsilon^{m_\ell}}{m_\ell !}
\int_{\Lambda^{m_\ell}}u^\varepsilon_{j_\ell + m_\ell}({\bm x}_{J_\ell},{\bm y}_{m_\ell}) \, \text{d}\rho^{\otimes m_\ell}({\bm y}_{m_\ell})\right) \;.
\end{align*}
By comparing with \eqref{eq: rtcf_rec}, we get \eqref{eq: rtcf}.

Absolute convergence in the region \eqref{eq:radconv} follows from \eqref{eq: trees}, recalling that $|{\cal T}_k| = k^{k-2}$ (Cayley's formula), that $\mu_\varepsilon \varepsilon^{d} = \varepsilon$, and using \eqref{eq: phi_reg} to control the decay functions. A more explicit bound for the truncated functions is provided in the rest of this section.
\qed

\bigskip
We turn to the proof of \eqref{eq: rtcf_est}.

>From \eqref{eq: rtcf} and \eqref{eq: trees} we get
\begin{equation}
\left|\frac{\rho_{j}^{\varepsilon,T}({\bm z}_{j})}{ \mu_{\varepsilon}^{j}}\right|  \leq
f^{\otimes j}({\bm z}_j) \emph{e}^{2j\beta B}
\sum_{n=0}^{\infty}\frac{\mu_{\varepsilon}^{n}}{n!}\emph{e}^{2n\beta B}
\int_{\Lambda^n}
\bar u^\varepsilon_{j+n}({\bm x}_{j+n})
\,\text{d}\rho(x_{j+1})\ldots \text{d}\rho(x_{j+n}) \
\label{eq: rtcf_tree}
\end{equation}
where
\begin{equation}
\bar u^\varepsilon_{j+n}({\bm x}_{j+n}) := \sum_{T \in {\cal T}_{j+n}} \prod_{\{\alpha,\alpha'\}\in E(T)}\left|\zeta^\varepsilon(x_\alpha,x_{\alpha'})\right| \;.
\end{equation}

For given $n$ and $T \in {\cal T}_{j+n}$, we study the integral
$$
I^\varepsilon_{T}({\bm x}_j) :=\int_{\Lambda^n}
 \prod_{\{\alpha,\alpha'\}\in E(T)}\left|\zeta^\varepsilon(x_\alpha,x_{\alpha'})\right|
\,\text{d}\rho(x_{j+1})\ldots \text{d}\rho(x_{j+n}) \;.
$$
Since the potential is compactly supported
in the ball of radius 1, we can restrict the integration region  to the set
$$
S^\varepsilon_{T}({\bm x}_j) := \Big\{|x_\alpha - x_{\alpha'}| < \varepsilon\ \ \forall \{\alpha,\alpha'\}\in E(T)\Big\}\;.
$$
In particular $I^\varepsilon_{T}({\bm x}_j) =0$ certainly holds unless $n \geq n_0 = n_0({\bm x}_j)$, where $n_0$
is defined as the minimum number of open balls of radius $\varepsilon$, $B_\varepsilon(x_{j+1}),\cdots,B_\varepsilon(x_{j+n})$, necessary to realize a connection between the points $(x_1,\cdots,x_j)$; meaning that $\cup_{1 \leq i \leq j+n}B_\varepsilon(x_i)$  is a connected set.
By \eqref{eq:LLTj}-\eqref{eq:LLj}, we observe that
$
n_0
$
cannot be smaller than ${\cal L} \left({\bm x}_{j}\right) / \varepsilon - j$.

Given ${\bm x}_j$, $n \geq n_0$ and $T \in {\cal T}_{j+n}$, we apply Fubini and perform the integrations in the following ordered way:

(i) first, we integrate with respect to all the available leaves (vertices of degree 1), thus obtaining a smaller tree;

(ii) we iterate the procedure and prune the tree graph leaf by leaf, until we are left with a union of $j-1$ linear chains joining $x_1,\cdots,x_j$;

(iii) we integrate with respect to the remaining variables forming the chains.

\noindent
Each one of the elementary integrals in steps (i) and (ii) is bounded by (cf.\,\eqref{eq: phi_reg})
$$
\int  \left|\zeta^\varepsilon(0,x)\right| \text{d}\rho(x) \leq \varepsilon^d \bar\rho\, C_\beta\;.
$$
Similarly using \eqref{eq: phi_stab} (with $n=2$), a chain made of $m+1$ edges is bounded by $(1+e^{2\beta B}) \left(\varepsilon^d \bar\rho\, C_\beta\right)^m$.
We conclude that
$$
I^\varepsilon_{T}({\bm x}_j) \leq (1+e^{2\beta B})^{j-1} \left(\varepsilon^d \bar\rho\, C_\beta\right)^n\;.
$$

Thus
$$
\left|\frac{\rho_{j}^{\varepsilon,T}({\bm z}_{j})}{ \mu_{\varepsilon}^{j}}\right|  \leq
f^{\otimes j}({\bm z}_j) \emph{e}^{2j\beta B}(1+e^{2\beta B})^{j-1}
\sum_{n=n_0}^{\infty}\frac{\varepsilon^{n}}{n!}\emph{e}^{2n\beta B}
(j+n)^{j+n-2}
\left(\bar\rho\, C_\beta\right)^n
$$
where we used again $\mu_\varepsilon \varepsilon^{d} = \varepsilon$.
As
$$
\frac{(j+n)^{j+n-2}}{n!} \leq j^{j-2} (A')^j (2e)^{n}
$$
for some pure constant $A'>0$, we obtain the  final result by choosing $A = A' e^{2\beta B}(1+e^{2\beta B})$
and $$\varepsilon_0 = \frac{1}{2\bar\rho\, C_\beta e^{2\beta B+1}}\;,$$.

\qed