\section{Introduction}
Boosting is a powerful technique for producing highly accurate voting classifiers by combining less accurate base learners.
Boosting algorithms are typically easy to fine tune and obtain state-of-the-art performance on many learning tasks.
Boosting dates back to the seminal work introducing the AdaBoost algorithm~\cite{freund_jcsc} and much work has gone into understanding and developing better boosting algorithms.
The best performing boosting algorithms are typically variants of gradient boosters~\cite{Friedman00greedyfunction}, such as LightGBM~\cite{lightgbm} and XGBoost~\cite{xgboost}, using Regression Trees as base learners.

Classic experiments~\cite{SFBL98} showed that boosting algorithms tend to improve their test accuracy even when training past the point of perfectly classifying the training data.
This may seem counter-intuitive, as adding more base learners, results in a more complex model, that hence might be more prone to overfitting.
This phenomenon is often explained by observed improvements in margins.
 For binary classification with a sample space $\mathcal{X}$, labels in $\{-1,1\}$ and a class of base learners $\H \subseteq \mathcal{X} \to [-1,1]$, a voting classifier $f : \mathcal{X} \to \{-1,1\}$ has the form $f(x) = \sign(\sum_{h \in \H} \alpha_h h(x))$ with all $\alpha_h \geq 0$. A voting classifier thus takes a weighted ``vote'' among the base learners to obtain its prediction.
When speaking of margins, we assume $\sum_h \alpha_h = 1$, which can always be achieved by rescaling the $\alpha$'s by their sum without changing $f$.
The margin of a training point $(x,y)$ with $x \in \mathcal{X}$ and $y \in \{-1,1\}$ is then defined as $y \sum_h \alpha_h h(x)$. The margin is thus a value in $[-1,1]$ which is positive when $f(x) = y$ and negative otherwise.
Intuitively, large (positive) margins mean that $f$ is not only correct but very certain in its predictions.
Margin theory, starting with the work of Schapire et al.~\cite{freund_jcsc}, formalized this by proving generalization bounds demonstrating that large margins imply better generalization. It was also shown that the theoretical generalization bounds fit very well with the observed behavior of AdaBoost that tends to keep improving margins even when training past the point of perfectly classifying the training data ~\cite{SFBL98}.

However, shortly after ~\cite{freund_jcsc} and~\cite{SFBL98} was published, Breiman \cite{arc_gv} proved a generalization bound based on the minimal margin (the smallest margin achieved by a training point) that was sharper than the generalization bound in Schapire et al.~\cite{freund_jcsc}. He then designed a new boosting algorithm, named \textit{Arc-GV}, that provably optimizes the minimal margin, which AdaBoost does not (see~\cite{gronlund_icml_19} for the full story of maximizing the minimal margin).
In the same paper, Breiman experimentally showed that Arc-GV produced not just a better minimal margin, but better margins overall, than AdaBoost. However, AdaBoost still obtained a better generalization and test error.
This seemed to contradict margin theory, as according to margin theory, all other things being equal, then larger margins should imply better generalization.
Later it was shown by Reyzin and Schapire~\cite{reyzin_2006} that Breiman's experiments did not accurately take into account the complexity of the base learner trees created by AdaBoost and Arc-GV, as repeating the experiments showed that Arc-GV produced trees of larger depth than AdaBoost, and deeper trees may be more prone to overfitting.
Reyzin and Schapire then considered the same experiments using  stumps as base learners, forcing identical depth trees between the algorithms, and in this case, AdaBoost produced better margin distributions than Arc-GV and also generalized better.
These findings support the view that better margins provide better generalization as presented in \cite{freund_jcsc, SFBL98}.

Later, \cite{emargin, koltchinskii2002, GZ13} showed improved generalization bounds that subsumed both the generalization bounds by Schapire et al., and Breiman, providing further theoretical support for margin theory.
The current strongest generalization bounds are as follows. Let $\mathcal{D}$ be any distribution over $\mathcal{X} \times \{-1,1\}$ and define $\mathcal{L}_\mathcal{D}(f) = \Pr_{(x,y) \sim \mathcal{D}}[f(x) \neq y]$ as the out-of-sample error of a voting classifier $f$. Also, for a set $S = \{(x_i,y_i)\}_{i=1}^m$ of $m$ labeled samples drawn i.i.d. from $\mathcal{D}$, define $\mathcal{L}_S^\theta(f) = \Pr_{(x,y) \sim S}[yf(x) < \theta]$ as the fraction of points in $S$ with margin less than $\theta$ (the notation $(x,y) \sim S$ denotes a uniform random point $(x,y)$ in $S$). With this notation, there are two strongest current generalization bounds. The first~\cite{koltchinskii2002} uses Rademacher complexity to show that with high probability over the sample set $S$, it holds for every margin $\theta \in (0,1]$ and every voting classifier $f$ that:
\begin{equation}
  \label{eq:rademacher}
\mathcal{L}_\mathcal{D}(f) \leq \mathcal{L}_S^\theta(f)+ O\left(\sqrt{\frac{\lg |\H|}{\theta^2 m}}\right).
  \end{equation}
The $k$'th margin bound by Gao and Zhou~\cite{GZ13} improves  this for $\mathcal{L}_S^\theta(f) = o(1/\lg m)$ and is as follows:
\begin{equation}
  \label{eq:kthmargin}
\mathcal{L}_\mathcal{D}(f) \leq \mathcal{L}_S^\theta(f)+ O\left(\frac{\lg |\H| \lg m}{\theta^2 m} + \sqrt{\mathcal{L}_S^\theta(f) \cdot \frac{\lg |\H| \lg m}{\theta^2 m}}\right).
\end{equation}
The $k$'th margin bound subsumes both Breiman's min margin generalization bound and the original generalization bound by Schapire et al.
For infinite $\H$, one may replace $\lg |\H|$ in the above bounds with the VC-dimension of $\H$ times a $\lg m$ factor (as is standard). For simplicity, we focus on the case of finite $\H$ throughout the paper.
Moreover, recent work by Gr\o nlund et al.~\cite{gronlund_nips19} shows that the margin bounds above are near-tight. Formally, they show that for (almost) all margins $\theta$, there exists a data distribution $\mathcal{D}$ and a set of base learners $\H$, such that with constant probability over the sample set $S$, there is a voting classifier $f$ such that
\begin{equation}
  \label{eq:lower}
\mathcal{L}_\mathcal{D}(f) \geq \mathcal{L}_S^\theta(f) + \Omega\left( \frac{\lg |\H| \lg m}{\theta^2 m} + \sqrt{\mathcal{L}_S^\theta(f) \cdot \frac{\lg |\H|}{\theta^2 m}}\right).
\end{equation}
Moreover, the lower bound holds for any  value of $\mathcal{L}_S^\theta(f) \leq 49/100$ and any value of $\lg|\H|$~\cite{gronlund_nips19}.

\textbf{Remark.} Many boosting algorithms produce classifiers $f = \sum_h \alpha_h h$ where $\sum_h \alpha_h \neq 1$ or where base learners output values in  $\mathbb{R}$ rather than $[-1,1]$.
To apply margin theory, following \cite{singer_confidence}, such classifiers are rescaled as follows: For each $h$ with output range $[a_h,b_h]$ and coefficient $\alpha_h$,  divide all outputs of $h$ by $\Delta_h = \max\{|a_h|,|b_h|\}$, multiply $\alpha_h$ by $\Delta_h$, ans then divide all $\alpha_h$ by $\sum_h \alpha_h$.

\subsection{Our contribution.}
\paragraph{A new margin lower bound:}
Comparing the current best upper and lower bounds, we see that~\eqref{eq:kthmargin} and~\eqref{eq:lower} match when $\mathcal{L}_S^\theta(f)$ approaches $0$. Similarly, we see that~\eqref{eq:kthmargin} and~\eqref{eq:rademacher} match as $\mathcal{L}_S^{\theta}(f)$ approaches a constant. But what is the true behavior in-between? The $k$'th margin bound~\eqref{eq:kthmargin} gained the factor $\mathcal{L}_S^\theta(f)$ inside the $\sqrt{\cdot}$ but lost a factor $\lg m$ compared to~\eqref{eq:rademacher}. Can the $\lg m$ factor be removed? What is the correct behavior as $\mathcal{L}_S^\theta(f)$ goes from $0$ towards $1$? In this work, we show an improved generalization lower bound of:
\begin{equation}
  \label{eq:our}
\mathcal{L}_\mathcal{D}(f) \geq \mathcal{L}_S^\theta(f) + \Omega\left( \frac{\lg |\H| \lg m}{\theta^2 m} + \sqrt{\mathcal{L}_S^\theta(f) \cdot \frac{\lg |\H| \lg(\mathcal{L}_S^\theta(f)^{-1})}{\theta^2 m}}\right).
\end{equation}
Our lower bound shows that the $\lg m$ factor inside the $\sqrt{\cdot }$ has to show up as $\mathcal{L}_S^\theta(f)$ drops to $m^{-\varepsilon}$ for any constant $\varepsilon > 0$. Moreover, our new lower bound completely settles the generalization performance of boosting in terms of margins whenever $\mathcal{L}_S^\theta(f)$ is outside the range $m^{-o(1)}$ to $o(1)$. It also nicely interpolates between the $\mathcal{L}_s^\theta(f) = 0$ and $\mathcal{L}_S^\theta(f) = 1$ case. We conjecture that the lower bound gives the correct margin-based tradeoff, i.e. that it is possible to improve the upper bounds~\eqref{eq:rademacher} and~\eqref{eq:kthmargin} to match~\eqref{eq:our}. Our proof is based on the work in \cite{gronlund_nips19}, and the recent near-tight generalization lower bound proof for Support Vector Machines shown in \cite{gronlund_icml_20}.
\paragraph{A new refined margin generalization bound:}
The main part of our paper considers a new refined margin based generalization bound for voting classifiers (boosting algorithms).
First, we present experiments showing that the classic margin bounds alone fail to explain the performance of state-of-the art gradient boosting algorithms. More concretely, we show that gradient boosters actually may produce smaller and smaller margins when run for many iterations, despite the test accuracy staying the same or even improving. We additionally demonstrate that the classic version of AdaBoost may produce significantly better margins than gradient boosters, despite gradient boosters obtaining similar or even better test accuracy and generalization error than AdaBoost. To explain this inconsistency, we observe experimentally that the trees produced by gradient boosters return very small values on all but a few training points, thus making minimal changes to most predictions when added to the voting classifier. We then use this insight to prove a new margin-based generalization bound for boosting algorithms which also take into account the magnitude of predictions by base learners. Finally, we run experiments demonstrating that our refined generalization bounds in fact succeed in explaining and predicting the performance of boosting algorithms.
In addition to achieving a better theoretical understanding of boosting algorithms, in particular gradient boosters, these new insights may potentially lead to new algorithms with better accuracy by using regularization inspired by our new generalization bound or more directly optimizing it.

\section{Insufficiency of current margin bounds}
\label{sec:insufficient}
From the margin-based upper and lower bounds, it may seem that we have all the theory necessary for understanding the generalization performance of boosters.
To confirm the theory, we ran experiments with AdaBoost and the state-of-the-art gradient booster LightGBM on standard data sets with the same size trees as base learners.
For all experiments we only change the tree size and learning rate of the LightGBM hyperparameters.
For AdaBoost we allow the same tree size, unlimited depth, as well as forcing a minimum number of elements in each tree learner to be 20 as is default in LightGBM.

Figure~\ref{fig:simple_margin} shows a plot of the margin distributions for the two boosters trained on the Forest Cover dataset.
From this plot, it is obvious that AdaBoost achieves significantly better margins than LightGBM. Indeed, the $k$'th smallest margin of AdaBoost, is much larger than the $k$'th smallest margin of LightGBM for all $k$ where at least one of the two margins are non-negative. Thus, from the generalization bounds \eqref{eq:rademacher} and \eqref{eq:kthmargin}, AdaBoost should have a much smaller out-of-sample error than LightGBM. However, the corresponding test errors in Figure~\ref{fig:simple_test} show a very different story, with LightGBM slightly outperforming AdaBoost. Furthermore, as shown in Section \ref{sec:refined}, the trees produced by LightGBM are in fact deeper than the trees produced by AdaBoost. This gives rise to some concerns regarding the explanatory power of margins.
\begin{figure}[h]
  \begin{subfigure}[t]{0.5\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{perf_plot_ds--forest_nl-256-ns-200.pdf}
    \caption{Mean training and test error over five runs. The standard deviation of the final test error is 0.00037 for AdaBoost and smaller for LightGBM.}
    \label{fig:simple_test}
  \end{subfigure}%
  $\quad$
  \begin{subfigure}[t]{0.5\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{quantile_plot_ds-forest_nl-256-ns-200.pdf}
    \caption{Sorted margin values.}
      \label{fig:simple_margin}
  \end{subfigure}
  \caption{Accuracy and margin plots for AdaBoost and LightGBM on the Forest Cover data set.}
\end{figure}
To further underline the theoretical inconsistency, we examine the two generalization bounds \eqref{eq:rademacher} and \eqref{eq:kthmargin}. When applying the generalization bounds to AdaBoost and LightGBM, then for any choice of $p = \mathcal{L}_S^\theta(f) \in [0,1]$, the only parameter that vary between AdaBoost and LightGBM is $\theta^{-2}$. That is, if we choose $\theta$ as the $(pm)$'th smallest margin, i.e. fix $\mathcal{L}_S^\theta(f) = p$, then only the value of $\theta$ differ between the two boosters and the generalization error grows as $\theta^{-2}$. Figure~\ref{fig:thetasquare} shows a plot of $\theta^{-2}$ as a function of $\mathcal{L}_S^\theta(f)$ for the two boosters. Clearly the penalty in the generalization error is much smaller for AdaBoost, suggesting that AdaBoost should perform much better than LightGBM, despite the test errors in Figure~\ref{fig:simple_test} showing that LightGBM outperforms AdaBoost.
\begin{figure}[h]
  \begin{subfigure}[t]{0.5\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{gen_error_plot_ds-forest_nl-256-ns-200.pdf}
    \caption{Plot of $\theta^{-2}$ when choosing $\theta$ as the $(pm)$'th smallest margin for $p \in [0,1]$. The margins are those also shown in Figure~\ref{fig:simple_margin}.}
         \label{fig:thetasquare}
  \end{subfigure}
    $\quad$
    \begin{subfigure}[t]{0.5\textwidth}    \centering
    \includegraphics[width=.9\linewidth]{margin_dev_plot-forest_nl-256-ns-200.pdf}
    \caption{Development in margin distributions for AdaBoost and LightGBM.}
    \label{fig:margin_dev}
  \end{subfigure}%
  \caption{Generalization penalties and margin distributions on the Forest Cover data set.}
\end{figure}
To investigate this phenomenon further, we have plotted the margin distribution of the two boosters after $t=10,20$ and $50$ iterations of training, see Figure~\ref{fig:margin_dev}. It is clear from this plot that the margins of the gradient booster, learned by LightGBM, deteriorate quickly with the number of training iterations.  To explain why the margins quickly drops towards $0$ for the gradient booster, we take a closer look at the  trees produced by LightGBM compared to AdaBoost. Figure~\ref{fig:histogram} shows a histogram of the predictions made by the trees produced by LightGBM. It is very striking from this histogram that the trees making up the LightGBM gradient booster makes very small (in absolute value) predictions on most data points, whereas AdaBoost always makes predictions in $\{-1,1\}$. Note that each tree always has its largest prediction among $\{-1,1\}$. Thus, LightGBM produces  trees that only significantly change the predictions of very few data points, while leaving almost all others unchanged. When training more and more trees, this causes the margins to diminish. To see this, consider as an example a training point $(x,1)$ and assume the first trained  tree $h$ makes a (correct) prediction of $h(x)=1$ and is assigned a weight of $\alpha_h=1$. After the first training iteration, the margin of $(x,1)$ is $1$. However, as training progresses, many more trees may be produced that all predict $0$ on $x$ while being assigned a weight of $1$. Since margins are normalized, $\sum_{h \in \H} \alpha_h = 1$, this means that the margin of $x$ drops to $1/t$ after $t$ rounds of training. The drop in predicted accuracy by the generalization bounds~\eqref{eq:rademacher} and~\eqref{eq:kthmargin} seem unreasonable if we think about the data point $x$ (the error is expected to grow as $t^2$ or $t$).
\begin{figure}[h]
    \centering
    \includegraphics[width=.7\linewidth]{pred_hist_plot_ds-forest_nl-256-ns-200.pdf}
    \caption{Histogram of base learner predictions for LightGBM on the Forest Cover data set. Only about 1 in 5000 predictions are larger than 0.95 in absolute value.}
    \label{fig:histogram}
\end{figure}
A possible explanation of the shortcomings of current generalization bounds is thus that they simply treat base learners as arbitrary functions in $\mathcal{X} \to [-1,1]$. That is, they pay no attention to the fact that base learners trained by gradient boosters make very small predictions on almost all data points.
To further support this claim, we note that the proof of the previous generalization lower bound~\eqref{eq:lower} as well as our improved bound~\eqref{eq:our} construct a set of base learners $\H$ where all $h \in \H$ make predictions among $\{-1,1\}$, i.e. they make no predictions of small magnitude. This further supports the belief that an explanation based on the magnitude of predictions may be found, which is the focus of the next section.
We have used a tree size of 256 as large tree sizes are used in practice and provide better test errors. Furthermore, the phenomena we are studying is clearer for large tree sizes.
In Section \ref{sec:refined} we show results for both large trees and stumps.
We note that base learners with real valued predictions were first considered by Schapire and Singer \cite{singer_confidence} that generalized the generalization bound of Schapire et al. \cite{SFBL98} to work with real values but without otherwise changing the bound.

\section{Refined margin bounds}
\label{sec:refined}
Motivated by the empirical observations in the previous section, we prove a more refined margin based generalization bound for voting classifiers. Define from a voting classifier $f$ the notation $\Delta(x,h) := |f(x) - h(x)|$. Intuitively, if a voting classifier $f$ has a small margin on a training point $x$, but this is the result of using mostly base learners $h$ that make small predictions (in absolute value), then $\Delta(x,h)$ will be small for most $h$ in $f$. Also define from a voting classifier $f = \sum_h \alpha_h h$ the distribution $\mathcal{Q}(f)$ over base learners, which simply returns $h$ with probability $\alpha_h$. With this notation, our new generalization bound states that for any distribution $\mathcal{D}$ over $\mathcal{X} \times \{-1,1\}$ and for any margin $\theta$, it holds with high probability over a set $S \sim \mathcal{D}^m$ that all voting classifiers $f$ satisfy:
\begin{equation}
  \label{eq:newupper}
\mathcal{L}_{\mathcal{D}}(f) \le \mathcal{L}_S^{\theta}(f) + O\left(\frac{N\lg|\H| \lg m}{m} + \sqrt{\mathcal{L}_S^{\theta}(f) \cdot \frac{N\lg|\H| \lg m}{m}}\right)\;,
\end{equation}
where $N = \max\{\theta^{-2}\cdot \left(\E_{(x,y)\sim S}\left[\E_{h \sim \mathcal{Q}(f)}\left[\Delta(x,h)^2\right]^{(\lg(16m))/2}\right]\right)^{2/(\lg(16m)} , \theta^{-1}\}$.

\paragraph{Never worse.}
Comparing our bound to the $k$'th margin bound~\eqref{eq:kthmargin}, we see that~\eqref{eq:newupper} equals the $k$'th margin bound when $N = \Theta(\theta^{-2})$. First, we argue that we always have $N = O(\theta^{-2})$, i.e.~\eqref{eq:newupper} is never worse than the $k$'th margin bound. To see this, observe that $\Delta(x,h) \leq 2$ since all $h \in \H$ produce values in $[-1,1]$. Thus, $\Delta(x,h)^2 \leq 4$ and $\E_{h \sim \mathcal{Q}(f)}\left[\Delta(x,h)^2\right] \leq 4$. This implies $\left(\E_{(x,y)\sim S}\left[\E_{h \sim \mathcal{Q}(f)}\left[\Delta(x,h)^2\right]^{(\lg(16m))/2}\right]\right)^{2/\lg(16m)} \leq 4$, hence we always have $N = O(\theta^{-2})$.

\paragraph{Potentially much better.}
Next, we demonstrate that our new bound may be significantly better than previous generalization bounds for very natural voting classifiers. For any desired margin $\theta \in (0,1]$, consider an example of a voting classifier $f(x) = \sum_{i=1}^{1/\theta}\theta h_i(x)$ such that for each training point $(x,y)$, there is exactly one hypothesis $h_i$ with $h_i(x) = y$ and all others have $h_j(x)=0$. This example is quite similar to the empirical performance of LightGBM seen in Section~\ref{sec:insufficient}, where most hypotheses make small predictions on most training points. The voting classifier $f$ has a margin of $\theta$ on all training points and thus the $k$'th margin bound predicts a generalization error of $O(\lg |\H| \lg m /(m \theta^2))$ (since $\mathcal{L}_S^\theta(f) = 0$ when all points have margin $\theta$). Let us now estimate $N$ in~\eqref{eq:newupper}. First, fix an $(x,y) \in S$ and consider the expression $\E_{h \sim \mathcal{Q}(f)}\left[\Delta(x,h)^2\right]^{(\lg(16m))/2} = \left(\sum_{i=1}^{1/\theta}\theta \cdot \Delta(x,h_i)^2 \right)^{(\lg(16m))/2} =\left(\theta \cdot (1-\theta)^2 + (1-\theta)\theta^2 \right)^{(\lg(16m))/2} < \theta^{(\lg(16m))/2}$. Since this holds for every $(x,y)$, we have $\left(\E_{(x,y)\sim S}\left[\E_{h \sim \mathcal{Q}(f)}\left[\Delta(x,h)^2\right]^{(\lg(16m))/2}\right]\right)^{2/\lg(16m)} < \theta$. Plugging that into the definition of $N$, we see that $N \leq \max\{\theta^{-2} \cdot \theta, \theta^{-1}\} = \theta^{-1}$. That is, the dependency on the margin has improved by a factor $\theta$ and our new generalization bound predicts $\mathcal{L}_\mathcal{D}(f) = O(\lg |\H| \lg m/(m \theta))$.

\begin{figure}[htb]
    \centering
    \includegraphics[width=.7\linewidth]{new_gen_error_plot_ds-forest_nl-256-ns-200.pdf}
    \caption{Generalization penalty $N$ on the Forest Cover data set when choosing $\theta$ as the $(pm)$'th smallest margin for $p \in [0,1]$.}
    \label{fig:new_gen_error}
\end{figure}
\paragraph{Comparison to earlier work.}
In recent work, Cortes et al. \cite{cortes_regboost}, also proved refined generalization bounds for gradient boosters.
Their works shows, that if the $q$-norm of the vector of leaf predictions for each tree trained by a gradient booster is small, then the trees have smaller VC-dimension and hence the voting classifier has better generalization performance (by using previous generalization bounds).
Note that their bound only depends on the leaf predictions and does not take into account the number of training points in each leaf. Our experiment in Figure \ref{fig:histogram} shows that for each base learner, only a tiny fraction (about 1 in 5000) of training points end in a leaf with large prediction, which our bound takes into account.

\begin{table}[ht]
  \caption{Comparing AdaBoost with LightGBM.
    In this experiment the trees used as bare learners are of increasing size relative to the data size.
    Each value shown is the average over several runs and each run use 200 rounds of boosting.
    Moment is  $\left(\E_{(x,y)\sim S}\left[\E_{h \sim \mathcal{Q}(f)}\left[\Delta(x,h)^2\right]^{(\lg(m))/2}\right]\right)^{2/\lg(m)}$.
  }
  \begin{center}
    \centering
    \begin{tabular}{lllllllll}
  \toprule
  Data Set & Alg. &  Train Err & Test Err &  Mean Margin & Max Depth & Mean Depth & Moment\\
  \midrule
  \multirow{2}{*}{Forest}
    & ada  &  0.0001 & 0.0331 & 0.1696 & 22.0 & 12.4 & 0.969 \\
  & lgb  &  0.0002 & 0.0291 & 0.0280 & 23.7 & 13.9 & 0.025 \\

  \midrule
  \multirow{2}{*}{Boone}   & ada  & 0.00009 & 0.0589 & 0.311 & 17.5 & 10.2 & 0.917 \\
    & lgb  & 0.00009 & 0.0552 & 0.0818 & 17.6 & 10.4 & 0.0564 \\
  \midrule
   \multirow{2}{*}{Higgs}  & ada  &  0.178 & 0.277 & 0.0747 & 24.9 & 13.5 & 0.99 \\
    & lgb  &  0.185 & 0.251 & 0.018 & 26 & 14.7 & 0.0289 \\
  \midrule
  \multirow{2}{*}{Diabetes}  & ada  &  0 & 0.268 & 0.148 & 3.5 & 2.63 & 0.973 \\
    & lgb  &  0.0264 & 0.26 & 0.142 & 3.5 & 2.63 & 0.214 \\
  \bottomrule
\end{tabular}
  \end{center}
\label{tab:large_trees}
\end{table}
\textbf{Empirical evaluation.}
Our new generalization bound carefully takes the magnitude of predictions made by the base learners into account, thus there is hope that~\eqref{eq:newupper} may better explain the experiments in the previous section. To test this, we have run the experiments again, this time plotting the value of $N$ as a function of $p = \mathcal{L}_S^\theta(f)$. That is, we notice that for two voting classifiers produced by AdaBoost and LightGBM, respectively, the only thing that varies in~\eqref{eq:newupper} when choosing the $(pm)$'th smallest margin, i.e. $p = \mathcal{L}_S^\theta(f)$, is the value of $N$. Thus smaller values of $N$ imply better generalization according to the theory. Figure~\ref{fig:new_gen_error} shows the result of the experiment. Quite remarkably, the relative ordering of AdaBoost and LightGBM match the observed test errors from Figure~\ref{fig:simple_test} much better, i.e. LightGBM slightly outperforms AdaBoost. We have repeated the same experiment on more data sets and summarized the results in Table~\ref{tab:large_trees}. The parameters for the experiments are shown in Table~\ref{tab:sets}

\begin{table}[htb]
    \caption{Data sets, all freely available, and parameters considered in the experiments. LR means learning rate as used in LightGBM.
    For each experiment we randomly split the data set in half to get a training set and a test set of equal size.
    For the Higgs dataset of size 11 million, we sample a subset of 2 million data points that we randomly split evenly into train and test set.
    For Forest Cover only the first two classes are used to make it into a binary classification problem.
  }
  \centering
  \begin{tabular}{llllll}
    \toprule
    Data Set &  Data Size & Tree Size & LR  & Stumps LR & Runs \\
    \midrule
    Diabetes & 768 & 5 & 0.1 & 0.1 &  100 \\
    \midrule
    Boone & 65032 & 96 & 0.2 & 0.6 & 5 \\
    \midrule
    Forest Cover & 495141 & 256 & 0.3 & 0.3 & 5\\
    \midrule
    Higgs & 2000000 & 512 & 0.3 & 0.3 & 5  \\
    \bottomrule
  \end{tabular}
  \label{tab:sets}
\end{table}

In all experiments, the margin distribution, here represented by the mean margin, is much worse for the LightGBM classifier, while the height of the trees used, both the max height and the mean height, is larger. Still the  LightGBM classifier generalizes at least as well (in fact,  slightly better) than the AdaBoost classifier.
Table \ref{tab:large_trees} also shows that the moment value from our generalization bound is significantly better for the LightGBM classifier.
When we consider our new generalization bound, the theory nicely matches the observed test errors in the same way as was shown in Figure \ref{fig:new_gen_error} for all data sets.
While not final proof that this is the real or only explanation, it suggests that the success of gradient boosters, despite having poor margins,
may be explained by the many small predictions made by the base learner trees.
The standard deviations of the test statistics are left out since they are extremely small for the three large data sets (and we have run 100 iterations of the small Diabetes data set).
For completeness we have included the same experiment replacing the large trees with stumps and shown the results in Table~\ref{tab:stumps}.
The results for stumps match those from the larger trees, just with a smaller difference in margins and moment values.
\begin{table}[htb]
  \caption{Experiments with stumps as base learners. Same setup as in Table~\ref{tab:large_trees}. }

  \centering
  \begin{tabular}{lllllllll}
    \toprule
    Data Set & Alg. &  Train Err & Test Err &  Mean Margin & Moment\\
    \midrule
    \multirow{2}{*}{Forest}  & ada  &  0.223 & 0.224 & 0.0754 & 0.987 \\
    & lgb  &  0.217 & 0.218 & 0.0225 & 0.0986 \\
    \midrule
  \multirow{2}{*}{Boone} & ada  &  0.0781 & 0.0817 & 0.138 & 0.975 \\
  & lgb  &  0.0669 & 0.0744 & 0.0422 &  0.239 \\
  \midrule
  \multirow{2}{*}{Higgs}   & ada  &  0.309 & 0.31 & 0.059  & 0.986 \\
  & lgb  &  0.301 & 0.302 & 0.0309 &  0.329 \\
  \midrule
  \multirow{2}{*}{Diabetes}  & ada  &  0.161 & 0.246 & 0.108 & 0.976 \\
    & lgb  &  0.176 & 0.238 & 0.138  & 0.299 \\
  \bottomrule
  \end{tabular}
\label{tab:stumps}
\end{table}

\section{Generalization Bound Proof} \label{sec:upperProof}
This section is devoted to the proof of our refined margin based generalization bound for voting classifiers, presented hereafter as Theorem~\ref{th:upperMain}. First we recollect some notation. Let $\mathcal{X}$ be some ground set, ${\cal D}$ a distribution over $\mathcal{X} \times [-1,1]$, ${\cal H} \subseteq \mathcal{X} \to [-1,1]$, and $\mathcal{C}=\mathcal{C}(\H)$ be the convex hull of $\H$.
Fix a voting classifier $f$, then there exists a sequence $\left\langle \alpha_h \right\rangle_{h \in {\H}} \in \mathbb{R}_+^{\H}$ such that $\sum_{h \in {\cal H}}{\alpha_h}=1$ and $f = \sum_{h \in {\cal H}}{\alpha_h \cdot h}$. Thus $f$ implicitly defines a distribution $\mathcal{Q} = \mathcal{Q}(f)$ over ${\cal H}$, where $\Pr_{h \sim \mathcal{Q}}[h=h'] = \alpha_{h'}$ for all $h' \in {\cal H}$. Finally, let $\Delta : \mathcal{X} \times \H \to \mathbb{R}$ be defined by $\Delta(x,h) := |f(x) - h(x)|$ for every $x \in \mathcal{X}$, $h \in \H$. We show the following.

\begin{theorem}\label{th:upperMain}
Let $\mathcal{D}$ be a distribution over $\mathcal{X}\times\{-1,1\}$ where $\mathcal{X}$ is some ground set, and let ${\cal H} \subseteq \mathcal{X} \to [-1,1]$. For every $\delta > 0$, it holds with probability at least $1-\delta$ over a set of $m$ samples $S \sim {\cal D}^m$, that for every voting classifier $f \in \mathcal{C}(\H)$ and every margin $\theta > 0$, we have
\begin{equation}
\mathcal{L}_{\mathcal{D}}(f) \le \mathcal{L}_S^{\theta}(f) + O\left(\frac{N\lg|\H| + \lg(1/\delta)}{m} + \sqrt{\frac{N\lg|\H|+\lg(1/\delta)}{m}\mathcal{L}_S^{\theta}(f)}\right)\;,
\label{eq:upperMain}
\end{equation}
where $N = O\left(\max\{\theta^{-2}\cdot \left(\E_{(x,y)\sim S}\left[\E_{h \sim \mathcal{Q}(f)}\left[\Delta(x,h)^2\right]^{(\lg(16m))/2}\right]\right)^{2/(\lg(16m)} , \theta^{-1}\}\lg m\right)$.
\end{theorem}

Denote by ${\cal E}= {\cal E}(\delta)$ the event that for every voting classifier $f$ and every margin $\theta > 0$, the bound in \eqref{eq:upperMain} holds with $N$ as defined in Theorem \ref{th:upperMain}.
In these notations we prove that $\Pr_{S \sim \mathcal{D}^m}[{\cal E}] \ge 1-\delta$.

\textbf{Proof overview.} Inspired by techniques presented by Schapire {\em et al.\ }\xspace \cite{SFBL98} and employed by Gao and Zhou \cite{GZ13}, our proof incorporates a discretization of the set of all voting classifiers over $\H$ to a discrete {\em net} of classifiers, such that, loosely speaking, every voting classifier over $\H$ can be approximated by a classifier that belongs to the net, and in addition, the size of the net is not too big, and thus union bounding over the net yields the desired probability bounds. Thus, intuitively speaking, by randomly rounding every voting classifier $f$ to the net we get an upper bound on the out of sample error for $f$.
More specifically, $N \in \mathbb{N}^+$ be some positive integer. We define a net ${\cal C}_N$ of voting classifier by
${\cal C}_N := \left\{\frac{1}{N}\sum_{j \in [N]}{h_j} : \left\langle h_j \right\rangle_{j \in [N]} \in {\cal H}^N\right\}$.
For every voting classifier $f$ over $\H$, we then give a randomized rounding scheme that essentially associates a random net element $g \in \mathcal{C}_N$ with $f$, and show that with high probability the out of sample error with respect to $g$ well-approximates that of $f$. By choosing $N$ carefully and union bounding over $\mathcal{C}_N$ we get an upper bound on the out of sample error for all voting classifiers $f$.
The crux of the proof lies in carefully choosing the size of the net, namely $N$. Loosely speaking, the net size $N$ has to be large enough, so that the net is rich enough to approximate every voting classifier well, but on the other hand small enough, so that union bounding over the net does not incur too large a cost for the probability bound. By subtly choosing $N$ and proving refined bounds on the rounding scheme we get the bound in Theorem~\ref{th:upperMain}.

Formally we define for every $N \in \mathbb{N}^+$, the event ${\cal E}_N$ to be the set of all samples $S \in (\mathcal{X}\times\{-1,1\})^m$ satisfying that for all voting classifiers $g \in C_N$ and integer $\ell \in [0,N]$ it holds that

\begin{equation*}
\mathcal{L}_\mathcal{D}^{\ell/N}(g) \le \mathcal{L}_S^{\ell/N}(g) + \frac{8\ln(2\delta^{-1}N(N+1)^2|\H|^N)}{m} + 4\sqrt{\frac{\ln(N(N+1)^2|\H|^N/\delta)}{m}\mathcal{L}_{S}^{\ell/N}(g)} \; ;
\end{equation*}
and
\begin{equation*}
\Pr_{(x,y)\sim\mathcal{D}}[\;|f(x)-g(x)| > \ell/N ] \le 2\Pr_{(x,y)\sim S}[\;|f(x)-g(x)| > \ell/N ] + \frac{8\ln(4\delta^{-1}N(N+1)^2|\H|^N)}{m} \;.
\end{equation*}
Intuitively speaking, for $S \in {\cal E}_N$, the first bound ensures a good generalization bound for every voting classifier $g$ in the net, whereas the second bound shows that $g$ approximates $f$ over $\mathcal{D}$ almost as well as it approximates $f$ over $S$. In turn these two bounds imply that the behavior of $f,g$ over $S$ predicts their behavior over $\mathcal{D}$.
As $\sum_{N=1}^\infty{\frac{1}{N(N+1)}}=1$, the following lemma implies Theorem~\ref{th:upperMain} by applying a union bound.
\begin{lemma}\label{l:mainTechnical}
For every $N \in \mathbb{N}^+$ we have $\Pr\limits_{S \sim\mathcal{D}^m}[{\cal E}_N] \ge 1 - \frac{\delta}{N(N+1)}$, and moreover, $\bigcap\limits_{N \in \mathbb{N}^+}{\cal E}_N \subseteq {\cal E}$.
\end{lemma}
The proof of the lemma is quite involved technically, and most of the proof is thus deferred to the appendix. Our main novelty lies in showing that for our choice of $N=N(f,\theta)$, for every sample set $S \in \supp(\mathcal{D}^m)$, with very high probability over the choice of a point $x \in \mathcal{X}$ and a net-classifier $g \in \mathcal{C}_N$, $g$ approximates $f$. In turn, this implies that if $S \in \bigcap_{N \in \mathbb{N}^+}{\cal E}_N$, then for every voting classifier $f$ and $\theta > 0$, $f$ is well-approximated by a randomized rounding to the net $\mathcal{C}_N$. Formally we show the following for every $f$ and $\theta$.
\begin{lemma}\label{l:gridApproxWHP}
$\Pr_{\substack{(x,y)\sim S \\ g \sim \mathcal{Q}^N}}[\;\Delta(x,g) > 49\theta/100] \le \frac{1}{m^2}$, where
\begin{equation*}
\begin{split}
N = N(f,\theta) := \lg(16m) \cdot \max\{&256 \theta^{-1}\|\Delta(x,h)\|_{\lg(16m)}, 100/\theta \;,\\
&128e \theta^{-2}\cdot \left(\E_{(x,y)\sim S}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^2\right]^{(\lg(16m))/2}\right]\right)^{2/(\lg(16m)} \}\;.
\end{split}
\end{equation*}
\end{lemma}
\begin{proof}
Let $Z = \Delta(x,g)$, then for every integer $r \ge 1$ we conclude from Markov's inequality that
\begin{equation}
\Pr_{\substack{(x,y)\sim S \\ g \sim \mathcal{Q}^N}}[\; Z > 49\theta/100] =  \Pr_{\substack{(x,y)\sim\mathcal{D} \\ g \sim \mathcal{Q}^N}}[Z^r > (49\theta/100)^r] \le \left(\frac{100}{49\theta}\right)^r\|Z\|_r^r\;.
\label{eq:momentBoundZ}
\end{equation}
It is therefore enough to show $\|Z\|_r^r \le \left(\frac{49\theta}{100}\right)^r m^{-2}$ for some positive integer $r \ge 1$. Let $r = 2 \cdot \left\lceil \lg(4m) / 2\right\rceil$, then $r$ is an even integer, satisfying $\lg(4m) = 2\lg(4m) / 2 \le r \le \lg(4m) + 2 \le N$. Since $r$ is even, then for $g = \frac{1}{N}\sum_{j \in [N]}{h_j}$ we get that
\begin{equation*}
Z^r=Z(x,g)^r = \left(\frac{1}{N}\sum_{j \in [N]}{(f(x)-h_j(x))}\right)^r = \frac{1}{N^r}\sum_{T = (j_i)_{i\in[r]}\in [N]^r}\prod_{i \in [r]}(f(x)-h_{j_i}(x))\;.
\end{equation*}
For every $T= (j_i)_{i\in[r]} \in [N]^r$ let $D(T) := \{j \in [N] : \exists i \in [r]. j_i=j\}$ be the set of distinct indices occurring in $T$, and for every $j \in [N]$, let $c_T(j) := |\{i \in [r] : j_i=j\}|$ be the number of times $j$ occurs in $T$. Then in these notations we have
$$Z^r = \frac{1}{N^r}\sum_{T \in [N]^r}\prod_{j \in D(T)}(f(x)-h_j(x))^{c_T(j)}\;.$$
As $h_1,\ldots,h_N$ are chosen independently, we get that
$$\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}[Z^r] = \frac{1}{N^r}\sum_{T \in [N]^r}\prod_{j \in D(T)}\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}\left[(f(x)-h_j(x))^{c_T(j)}\right]\;.$$
Let $T \in [N]^r$, and assume that for some $j \in D(T)$ we have $c_T(j)=1$, then
$$\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}\left[(f(x)-h_j(x))^{c_T(j)}\right] = \E_{h \sim \mathcal{Q}}\left[f(x)-h(x)\right] = f(x) - \E_{h \sim \mathcal{Q}}\left[h(x)\right] = f(x) - \sum_{h \in \H}{\alpha_hh(x)}=0 \;,$$
Denote $\mathcal{T} := \{T \in [N]^r : \forall j \in D(T). \; c_T(j) > 1\}$, then
\begin{equation}
\begin{split}
\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}[Z^r] &= \frac{1}{N^r}\sum_{T \in \mathcal{T}}\prod_{j \in D(T)}\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}\left[(f(x)-h_j(x))^{c_T(j)}\right] \\
&= \frac{1}{N^r}\sum_{T \in \mathcal{T}}\prod_{j \in D(T)}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{c_T(j)}\right]\;.
\label{eq:momentBound}
\end{split}
\end{equation}
By Lyapunov's Theorem (see, e.g. \cite{MOA11}), $\E_{h\sim\mathcal{Q}}[\Delta(x,h)^\xi]$ is logarithmic convex for $\xi \in [1,+\infty)$, and as $c_T(j) \ge 2$ for all $j \in D(T)$ we get that
$$\prod_{j \in D(T)}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{c_T(j)}\right] \le \E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{|D(T)|-1}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{r-2|D(T)|+2}\right] \;.$$ Plugging into \eqref{eq:momentBound} we get that
\begin{equation}
\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}[Z^r] \le \frac{1}{N^r}\sum_{T \in \mathcal{T}}{\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{|D(T)|-1}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{r-2|D(T)|+2}\right]}\;.
\label{eq:momentBoundSimplified}
\end{equation}
For every $d \in \mathbb{N}$ denote $\mathcal{T}_d := \{T \in \mathcal{T} : |D(T)|=d\}$. Since for every $T \in \mathcal{T}$ and every $j \in D(T)$, we know that $c_T(j) \ge 2$, then for every $d > r/2$ we get that $\mathcal{T}_d = \emptyset$. Therefore $\mathcal{T} = \charfusion[\mathop]{\bigcup}{\cdot}_{d \in [r/2]}{|\mathcal{T}_d|}$. Moreover, for every $d \in [r/2]$ and every $T \in \mathcal{T}_d$, we have
$$\E_{h \sim \mathcal{Q}}\left[|h(x)|^{2}\right]^{|D(T)|-1}\E_{h \sim \mathcal{Q}}\left[|h(x)|^{r-2|D(T)|+2}\right] = \E_{h \sim \mathcal{Q}}\left[|h(x)|^{2}\right]^{d-1}\E_{h \sim \mathcal{Q}}\left[|h(x)|^{r-2d+2}\right] \;.$$ We therefore refine \eqref{eq:momentBoundSimplified} to get
\begin{equation}
\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}[Z^r] \le \frac{1}{N^r}\sum_{d \in [r/2]}{|\mathcal{T}_d|\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{d-1}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{r-2d+2}\right]}\;.
\label{eq:momentBoundAlmostFinal}
\end{equation}
\begin{claim}
For every $d \in [r/2]$, $|\mathcal{T}_d| \le r^r\sqrt{2e\pi r}\left(\frac{Ne}{r}\right)^d $.
\end{claim}
\begin{proof}
Fix some $d \in [r/2]$. There are at most $\binom{N}{d}$ ways to choose a subset $Y \subseteq [N]$ such that $|Y|=d$. Once such a set $Y$ is fixed, there are at most $\binom{d+(r-2d)-1}{r-2d}$ solution to the equation $\sum_{j \in Y}{y_j} = r$ under the constraint that $y_j \in \mathbb{N}\setminus\{0,1\}$ for all $j \in Y$. Moreover, once $\{y_j\}_{j \in Y}$ is fixed, there are $r! \cdot \prod_{j \in Y}(y_j!)^{-1}$ ways to form a sequence $T$ satisfying that $D(T)=Y$, $c_T(j) = y_j$ for all $j \in Y$ and $c_T(j)=0$ otherwise. Note that $\prod_{j \in Y}(y_j!) \ge ((r/d)!)^d$ for every choice of $\{y_j\}_{j \in Y}$, and therefore
\begin{equation*}
\begin{split}
|\mathcal{T}_d| &\le \binom{N}{d}\cdot\binom{r-d-1}{r-2d} \cdot \frac{r!}{((r/d)!)^d} \le \left(\frac{Ne}{d}\right)^d\cdot 2^{r-d} \cdot \frac{\sqrt{2e\pi r}(r/e)^{r}}{(\sqrt{2 \pi (r/d)}(r/(ed))^{r/d})^d} \\
&\le \sqrt{2e\pi r}\left(Ne\right)^d \cdot r^{r-d} \le r^r\sqrt{2e\pi r}\left(\frac{Ne}{r}\right)^d
\end{split}
\end{equation*}
\end{proof}
Plugging into \eqref{eq:momentBoundAlmostFinal} we conclude that
\begin{equation*}
\begin{split}
\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}[Z^r] &\le \frac{1}{N^r}\sum_{d \in [r/2]}{r^r\sqrt{2e\pi r}\left(\frac{Ne}{r}\right)^d\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{d-1}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{r-2d+2}\right]}\\
&=\sqrt{2e\pi r}\left(\frac{r}{N}\right)^r\sum_{d \in [r/2]}{\left(\frac{Ne}{r}\right)^d\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{d-1}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{r-2d+2}\right]}
\label{eq:momentBoundNoT}
\end{split}
\end{equation*}
As $\left(\frac{Ne}{r}\right)^\xi, \E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{\xi-1}, \E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{r-2\xi+2}\right]$ are all logarithmic convex for $\xi \in [1,r/2]$, their product is also logarithmic convex over that range, and thus gets its maximum on either $1$ or $r/2$. Concluding we get that
\begin{equation*}
\begin{split}
\E_{(h_k)_{k \in [N]} \sim \mathcal{Q}^N}[Z^r]  \le \frac{r}{2} \cdot \sqrt{2e\pi r}\left(\frac{r}{N}\right)^r\left(\left(\frac{Ne}{r}\right)\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{r}\right] + \left(\frac{Ne}{r}\right)^{r/2}\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{r/2}\right) \;.
\label{eq:momentBoundFinal}
\end{split}
\end{equation*}
Taking the expectation over $(x,y)\sim\mathcal{D}$ gives
\begin{equation}
\begin{split}
\|Z\|_r^r \le \frac{r}{2}\sqrt{2e\pi r}\left(\frac{r}{N}\right)^r\left(\left(\frac{Ne}{r}\right)\|\Delta(x,h)\|_r^r + \left(\frac{Ne}{r}\right)^{r/2}\E_{(x,y)\sim\mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^2\right]^{r/2}\right]\right)
\end{split}
\label{eq:ZMomentBound}
\end{equation}
To finish the proof of Lemma~\ref{l:gridApproxWHP}, we show that our bound on $N$ implies that $\|Z\|_r^r \le \left(\frac{49\theta}{100}\right)^r m^{-2}$.
Denote
\begin{equation*}
\begin{split}
\Psi_1 &= \frac{r}{2} \cdot \sqrt{2e\pi r}\left(\frac{r}{N}\right)^r \cdot \left(\frac{Ne}{r}\right)\|\Delta(x,h)\|_r^{r} = \frac{r}{2} \cdot \sqrt{2e\pi r}\left(\frac{r\|\Delta(x,h)\|_r}{N }\right)^r \cdot \left(\frac{Ne}{r}\right)\\
\Psi_2 &= \frac{r}{2} \cdot \sqrt{2e\pi r}\left(\frac{r}{N}\right)^r\left(\left(\frac{Ne}{r}\right)^{r/2}\E_{(x,y) \sim \mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{r/2}\right]\right)\\
\end{split}
\end{equation*}
Plugging into \eqref{eq:ZMomentBound} we get that $\|Z\|_r^r \le  \Psi_1 + \Psi_2$.

We will show that $\max\{\Psi_1,\Psi_2\} \le \left(\frac{49\theta}{100}\right)^r \cdot \frac{1}{2m^2}$, which proves the claim.
To bound $\Psi_1$, note first that $\Psi_1$ decreases as a function of $N$ (since $r \ge 2$). Since $N \ge  256\theta^{-1}\lg(16m) \cdot \|\Delta(x,h)\|_{\lg(16m)}$ we get that
\begin{equation*}
\Psi_1 \le \frac{r}{2} \cdot \sqrt{2e\pi r}\left(\frac{r \cdot \|\Delta(x,h)\|_r}{256\theta^{-1}\lg(16m) \cdot \|\Delta(x,h)\|_{\lg(16m)}}\right)^r \cdot \left(\frac{256\theta^{-1}\lg(16m) \cdot \|\Delta(x,h)\|_{\lg(16m)} \cdot e}{r}\right)
\label{eq:BoundPsi1}
\end{equation*}
Since $r < \lg(16m)$, and by monotonicity of norms,
$\|\Delta(x,h)\|_r \le \|\Delta(x,h)\|_{\lg(16m)} \le 2$,
where the last inequality is due to the fact that $|f(x) - h(x)| \le 2$ for all $h \in \H$, $x \in \mathcal{X}$. Moreover, $\lg(4m) \le r \le \lg(16m) \le 2(\lg(4m))$, therefore
\begin{equation*}
\begin{split}
\Psi_1 &\le \frac{r}{2} \cdot \sqrt{2e\pi r}\left(\frac{\theta}{256}\right)^r \cdot 1024e\theta^{-1}\\
&\le \left(\frac{49\theta}{100}\right)^r \cdot 3r^{3/2}125^{-r} \cdot \left(1024e \theta^{-1}\right) \le \left(\frac{49\theta}{100}\right)^r \cdot  3r^{3/2}64^{-\lg m}125^{-2} \cdot \left(1024e \theta^{-1}\right) \\
&\le \left(\frac{49\theta}{100}\right)^r \cdot \frac{1}{5}\lg^{3/2}(4m) \cdot m^{-6}\theta^{-1}\le \left(\frac{49\theta}{100}\right)^r \cdot \frac{1}{2m^2} \cdot \frac{1}{2}(\lg(4m)/m)^{3/2}(m^{5/2}\theta)^{-1}
\end{split}
\end{equation*}
For large enough $m$, we have that $\lg(4m)/m \le 5/8$, and therefore $(\lg(4m)/m)^{3/2} \le 1/2$. Since $\theta \ge 1/m$ we get that $\Psi_1 \le \left(\frac{49\theta}{100}\right)^r \cdot \frac{1}{2m^2}$.
We now turn to bound $\Psi_2$. Recall that $N \ge  128e \theta^{-2}\lg(16m)\cdot \left(\E_{(x,y)\sim\mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^2\right]^{\lg(16m)/2}\right]\right)^{2/\lg(16m)}$, and therefore
\begin{equation*}
\begin{split}
\Psi_2 &\le 3r^{3/2}\left(\frac{er\E_{(x,y) \sim \mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{r/2}\right]^{2/r}}{128 e \theta^{-2}\lg(16m)\left(\E_{(x,y)\sim\mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^2\right]^{\lg(16m)/2}\right]\right)^{2/\lg(16m)}}\right)^{r/2}\\
&\le \left(\frac{49\theta}{100}\right)^r \cdot 3r^{3/2} \left(\frac{r\E_{(x,y) \sim \mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^{2}\right]^{r/2}\right]^{2/r}}{30 \lg(16m)\left(\E_{(x,y)\sim\mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^2\right]^{\lg(16m)/2}\right]\right)^{2/\lg(16m)}}\right)^{r/2}\\
\end{split}
\label{eq:BoundPsi2}
\end{equation*}
Since $r< \log(16m)$, and by monotonicity of norms of random variables, we get that
\begin{equation*}
\E_{(x,y)\sim\mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^2\right]^{r/2}\right]^{2/r} \le \E_{(x,y)\sim\mathcal{D}}\left[\E_{h \sim \mathcal{Q}}\left[\Delta(x,h)^2\right]^{\log(16m)/2}\right]^{2/\log(16m)}\;.
\end{equation*}
Therefore
\begin{equation*}
\begin{split}
\Psi_2&\le \left(\frac{49\theta}{100}\right)^r \cdot 3r^{3/2}\left(30\right)^{-r/2} \le \left(\frac{49\theta}{100}\right)^r \cdot 3r^{3/2}\left(30\right)^{-(\lg m)/2 - 1} \le \left(\frac{49\theta}{100}\right)^r \cdot \frac{1}{2m^2} \cdot \frac{1}{5} r^{3/2}m^{-2/5}
\end{split}
\end{equation*}
Similarly to before, for large enough $m$, $\lg^{3/2}(4m)\cdot m^{-2/5} \le 5$, and therefore we conclude that $\Psi_2 \le \left(\frac{49\theta}{100}\right)^r \cdot \frac{1}{2m^2}$, which completes the proof of the lemma.
\end{proof}

\section{Generalization lower bound} \label{sec:lowerProof}
In this section we state and prove our new generalization lower bound, presented as Theorem~\ref{th:lowerEx}.
\begin{theorem} \label{th:lowerEx}
For every large enough integer $N$, every $\theta \in \left(1/N, 1/40\right)$, $\tau \in [0,1]$ and every $\left(\theta^{-2}\ln N\right)^{1 + \Omega(1)} \le m \le 2^{N^{O(1)}}$, if $\frac{\ln N \ln m}{m\theta^2} \le \tau \le 1$, then there exist a set $\mathcal{X}$, a hypothesis set $\H$ over $\mathcal{X}$ and a distribution $\mathcal{D}$ over $\mathcal{X} \times \{-1,1\}$ such that $\ln |\H|= \Theta(\ln N)$ and with probability at least $1/100$ over the choice of samples $S \sim \mathcal{D}^m$ there exists a voting classifier $f_S \in C(\H)$  such that
\begin{enumerate}
	\item $\mathcal{L}_S^{\theta}(f_S) \le \tau$; and
	\item $\mathcal{L}_{\mathcal{D}}(f_S) \ge \mathcal{L}_S^{\theta}(f_S) + \Omega\left(\frac{\ln|\H|\ln m}{m \theta^2} + \sqrt{\tau\ln(\tau^{-1}) \cdot \frac{\ln|\H|}{m \theta^2}}\right)$.
\end{enumerate}
\end{theorem}

Our proof is inspired by the constructions in \cite{gronlund_icml_20, gronlund_nips19} and  makes use of the following lemma, whose proof can be found in \cite{gronlund_nips19}.
\begin{lemma} \label{l:hypoDist}
For every $\theta \in (0,1/40)$, $\delta \in(0,1)$ and integers $d \le u$, there exists a distribution $\mu=\mu(u,d,\theta,\delta)$ over hypothesis sets $\H \subset \mathcal{X} \to \{-1,1\}$, where $\mathcal{X}$ is a set of size $u$, such that the following holds for $N = \Theta\left(\theta^{-2}\ln d \ln (\theta^{-2}d\delta^{-1})e^{\Theta(\theta^2 d)}\right)$.
\begin{enumerate}
	\item For all $\H \in \supp(\mu)$, we have $|\H|=N$; and
	\item For every labeling $\ell \in \{-1,+1\}^u$, if no more
  than $d$ points $x \in \mathcal{X}$ satisfy $\ell(x) = -1$, then
\[
\Pr_{\H \sim \mu}[\exists f \in \mathcal{C}(\H) : \forall x \in \mathcal{X}. \;
\ell(x)f(x) \geq \theta] \geq 1-\delta \;,
\]
\end{enumerate}
\end{lemma}

We start by describing the outlines of the proofs.
To this end fix some integer $N$, and fix $\theta \in \left(1/N, 1/40\right)$. Let $u$ be an integer, and let $\mathcal{X} = \{\xi_1,\ldots,\xi_u\}$ be some set with $u$ elements.
The distribution $\mathcal{D}$ over $\mathcal{X} \times \{-1,1\}$, is simply the uniform distribution over $\mathcal{X}\times\{1\}$. That is for every $i \in [u]$ and $y \in \{-1,1\}$, $\Pr_{\mathcal{D}}[(\xi_i,y)] = \frac{1+y}{2u}$. The following claim is straightforward.

\begin{claim}\label{c:boundGenPsi}
For every $f:\mathcal{X}\to \mathbb{R}$ we have $\Pr\limits_{(x,y)\sim\mathcal{D}}[yf(x) < 0] = \frac{1}{u}\sum_{i \in [u]}{\mathbbm{1}_{f(\xi_i)<0}}$.
\end{claim}

We will show that with some constant probability over a random choice $S \sim \mathcal{D}^m$, an adversarial voting classifier has a high generalization probability.
We additionally show existence of a hypothesis set $\hat{\mathcal{H}}$ such that with very high (constant) probability over a random choice of $\ell \in \{-1,1\}^u$, $C(\hat{\mathcal{H}})$ contains a voting classifier that attains high margins with $\ell$ over the entire set $\mathcal{X}$.
Finally, we conclude that with positive probability over a random choice of $S \sim \mathcal{D}^m$ both properties are satisfied.

To prove existence of a "rich" yet small enough hypothesis set $\hat{\mathcal{H}}$ we apply Lemma~\ref{l:hypoDist} together with Yao's minimax principle. In order to ensure that the hypothesis sets constructed using Lemma~\ref{l:hypoDist} is small enough, and specifically has size $N^{O(1)}$, we need to focus our attention on sparse labelings $\ell \in \{-1,1\}^u$ only. That is, the labelings cannot contain more than $\frac{\ln N}{\theta^2}$ entries equal to $-1$. To this end we will focus on $d$-sparse vectors. More formally, we define a set of labelings of interest $\L(u,d)$ as follows.
\begin{equation}
\L(u,d) := \{\ell \in \{-1,1\}^u : |\{i \in [u]: \ell_i = -1\}|\le d\}\;.
\label{eq:labelingSet}
\end{equation}
We next show that there exists a small enough (with respect to $N$) hypothesis set $\hat{\mathcal{H}}$ that is rich enough. That is, with high probability over $\ell \in \L(u,d)$, there exists a voting classifier $f \in C(\hat{\mathcal{H}})$ that attains high minimum margin with $\ell$ over the entire set $\mathcal{X}$. Note that the following result, similarly to Lemma~\ref{l:hypoDist} does not depend on the size of $\mathcal{X}$, but only on the sparsity of the labelings in question.
\begin{claim} \label{c:HspExistsEx}
If $u \le 2^{N^{O(1)}}$ and $d \le \frac{\ln N}{\theta^2}$ then there exists a hypothesis set $\hat{\mathcal{H}}$ such that $\ln|\hat{\mathcal{H}}| = \Theta\left(\ln N\right)$ and
\[\Pr_{\ell \in_R \L(u,d)}[\exists f \in \mathcal{C}(\hat{\mathcal{H}}) : \forall i \in [u]. \; \ell_if(\xi_i) \geq \theta] \geq 1-1/N \;.\]
\end{claim}
\begin{proof}
Let $\mu = \mu(u,d,\theta, 1/N)$, be the distribution whose existence is guaranteed in Lemma~\ref{l:hypoDist}. Then for every labeling $\ell \in \L(u,d)$, with probability at least $99/100$ over $\H \sim \mu$, there exists a voting classifier $f \in C(\H)$ that has minimal margin of $\theta$. That is, for every $i \in [u]$, $\ell_i f(\xi_i) \ge \theta$. By Yao's minimax principle, there exists a hypothesis set $\hat{\mathcal{H}} \in \supp(\mu)$ such that
\begin{equation*}
\Pr_{\ell \in_R \L(u,d)}[\exists f \in \mathcal{C}(\hat{\mathcal{H}}) : \forall i \in [u]. \; \ell_if(x_i) \geq \theta] \geq 1-1/N \;.
\end{equation*}
Moreover, since $\hat{\mathcal{H}} \in \supp(\mu)$, then $|\hat{\mathcal{H}}| = \Theta\left(\theta^{-2}\ln u \cdot \ln(N\theta^{-2}\ln u)\cdot e^{\Theta(\theta^2 d)}\right)$. Since $\theta \ge 1/N$, $\ln u \le N^{O(1)}$, and $d \le \frac{\ln N}{\theta^2}$, and thus $e^{\theta^2 d} = N$ we get that there exists some universal constant $C>0$ such that $|\hat{\mathcal{H}}| = \Theta(N^C)$, and thus $\ln |\hat{\mathcal{H}}| = \Theta(\ln N)$.
\end{proof}

Let $u = \frac{\ln N}{16 \tau \theta^2}$, and let $d = \frac{\ln N}{16e^{28}\theta^2}$. We next introduce some notation. With every set $T \subseteq [u]$ we associate the classifier $h_T : \mathcal{X} \to \{-1,1\}$ satisfying that for every $x \in \mathcal{X}$, $h_T(x)= -1$ if and only if $x \in T$. For every $m$-point sample $S \in (\mathcal{X}\times\{1\})^m$ and every $i \in [u]$, let $b_i^S$ be the number of times $\xi_i$ is sampled into $S$. If the set $S$ is clear from context, we simply denote $b_i$. In these notations, $\mathcal{L}_S(h_T) = \frac{1}{m}\sum_{i \in T}b_i^S$ for every $T \subseteq [u]$. Given a sample set $S$ Let $T^*=T^*(S) \subseteq [u]$ be a random set of size $d$ that minimizes $\mathcal{L}_S(h_{T^*(S)}) = \sum_{i \in T^*(S)}{b_i^S}$. We will show the following.

\begin{lemma} \label{l:existsLabelEx}
With probability at least $1/100$ over the choice of sample $S \sim \mathcal{D}^m$, the following holds.
\begin{enumerate}
	\item There exists a voting classifier $f_S \in C(\hat{\mathcal{H}})$ such that $f_S(\xi_i)h_{T^*(S)}(\xi_i) \ge \theta$ for all $i \in [u]$; and
	\item $ \mathcal{L}_S(h_{T^*(S)}) \le \frac{d}{u}\left(1 - \sqrt{\frac{\ln(u/2d)}{9m/u}}\right) \;.$
\end{enumerate}
\end{lemma}
Note that as $\tau \ge \frac{\ln N \ln m}{m \theta^2}$ we know that $u = \frac{\ln N}{16 \tau \theta^2} \le \frac{m}{16 \ln m}$ and therefore $\frac{\ln(u/2d)}{9m/u} \le \frac{u\ln(e^{28}/\tau)}{9m} \le \frac{\ln(e^{28}/\tau)}{144 \ln m} \le \frac{1}{2}$ for large enough $N$, and therefore the bound in the second part of Lemma~\ref{l:existsLabelEx} is meaningful. We first show that the lemma implies Theorem~\ref{th:lowerEx}.
\begin{proof}[Proof of Theorem~\ref{th:lowerEx}]
Fix some $\frac{\ln N \ln m}{m \theta^2} \le \tau \le 1$. From Lemma~\ref{l:existsLabelEx} with probability $1/100$ over the choice of a sample $S \sim \mathcal{D}^m$ there exists a voting classifier $f_S \in C(\hat{\mathcal{H}})$ such that $f_S(\xi_i)h_{T^*(S)}(\xi_i) \ge \theta$ for all $i \in [u]$ and moreover $\mathcal{L}_S(h_{T^*(S)}) \le \tau$. Consider $f_S$, and note first that
\[
\mathcal{L}_{\mathcal{D}}(f_S) = \frac{1}{u}\sum_{i \in [u]}{\mathbbm{1}_{f_S(\xi_i)<0}} = \frac{1}{u}\sum_{i \in [u]}{\mathbbm{1}_{h_{T^*(S)}(\xi_i)<0}} = \frac{|T^*(S)|}{u} = \frac{d}{u} \;.
\]
Additionally, since for every $i \in [u]$, $f_S(\xi_i) \le 0$ if and only if $f_S(\xi_i) \le \theta$, then
\[
\mathcal{L}_S^\theta(f_S) = \mathcal{L}_S(f_S) = \mathcal{L}_S(h_{T^*(S)}) \le \frac{d}{u}\left(1 - \sqrt{\frac{\ln(u/2d)}{9m/u}}\right) \le \frac{d}{2u} \le \tau\;.
\]
Summing up we get also that
\[
\mathcal{L}_{\mathcal{D}}(f_S) - \mathcal{L}_S^\theta(f_S) \ge \frac{d}{u}\sqrt{\frac{\ln(u/2d)}{9m/u}} = \Omega\left (\tau\sqrt{\frac{u\ln(\tau^{-1})}{m}} \right) = \Omega\left (\sqrt{\frac{\ln N \tau \ln(\tau^{-1})}{m\theta^2}} \right) \;.
\]
\end{proof}

For the rest of the section we therefore prove Lemma~\ref{l:existsLabelEx}.
First note that since $\mathcal{D}$ is uniform over $\mathcal{X} \times \{1\}$, and since given $S \sim \mathcal{D}^m$, $T^*$ is sampled uniformly over all subsets $T \in \binom{[u]}{d}$ such that the sum $\sum_{i \in T}{b_i^S}$ is minimized, we get that for every $T \in \binom{[u]}{d}$, $\Pr_{S \sim \mathcal{D}^m}[T^*(S) = T] = \binom{u}{d}^{-1}$. In other words, for every $h \in \L(u,d)$, $\Pr_{S \sim \mathcal{D}^m}[h_{T^*(S)} = h] = \binom{u}{d}^{-1}$. Therefore $h_{T^*(S)}$ is uniformly distributed over $\L(u,d)$. From claim~\ref{c:HspExistsEx} it follows that for large enough $N$,
the probability over the choice of $S \sim \mathcal{D}^m$ that there exists $f_S \in C(\hat{\mathcal{H}})$ such that $f_S(\xi)h_{T^*(S)}(\xi_i) \ge \theta$ for all $i \in [u]$ is at least $99/100$. In order to prove Lemma~\ref{l:existsLabelEx}, it is therefore enough to show that with probability at least $1/50$ over the choice of $S \sim \mathcal{D}^m$, $\mathcal{L}_S(h_{T^*(S)}) \le \frac{d}{u}\left(1 - \sqrt{\frac{\ln(u/2d)}{9m/u}}\right)$. We will show that with probability at least $1/50$ over the choice of $S$ there exist $i_1,\ldots,i_d \in [u]$ such that for every $j \in [d]$, $b_{i_j}^S \le \frac{m}{u}\left(1 - \sqrt{\frac{\ln(u/2d)}{9m/u}}\right)$. Since $T^*(S)$ minimizes $\sum_{i \in T^*(S)}b_i^S$, it follows that
\[
\mathcal{L}_S(h_{T^*(S)}) = \frac{1}{m}\sum_{i \in T^*(S)}b_i^S \le \frac{1}{m}\sum_{j \in [d]}b_{i_j}^S \le \frac{d}{u}\left(1 - \sqrt{\frac{\ln(u/2d)}{9m/u}}\right) \;.
\]
To this end, fix some $i \in [u]$. For every $j \in [m]$, let $I_j^S$ be an indicator for the event that the $j$th element selected into $S$ is $(\xi_i,1)$. Then $b_i^S = \sum_{j \in [m]}{I_j^S}$, and as $\mathcal{D}$ is uniform, we get that $\mathbb{E}[b_i^S] = \sum_{j \in [m]}\mathbb{E}[I_j^S] = m/u$. We will use the following reverse Chernoff bound and show that with good enough probability, $b_i^S$ is far from its expectation.

\begin{lemma}\label{l:revChernoff}
Let $m \in \mathbb{N}^+$ and let $I_1,\ldots,I_m$ be independent indicator random variables with success probability $1/u$. Then for every $\sqrt{3/(m/u)} \le \delta \le 1/2$ we have
\[
\Pr\left[\sum_{j \in [m]}{I_j} \le (1-\delta)mp\right] \ge e^{-9m\delta^2/u} \;.
\]
\end{lemma}
Denote $\delta := \sqrt{\frac{\ln(u/2d)}{9m/u}}$. As we have shown earlier, $\delta \le 1/2$. Moreover, since $\frac{u}{2d} \ge e^{27}\tau^{-1} \ge e^{27}$, we get that $\delta \ge \sqrt{\frac{27}{9m/u}} = \sqrt{\frac{3}{m/u}}$. We can therefore conclude from Lemma~\ref{l:revChernoff} that
\[
\Pr[b_i^S \ge (1-\delta)m/u] \ge e^{-9m\delta^2/u} = e^{-\ln(u/2d)} = \frac{2d}{u} \;.
\]
Let $B_i^S$ be the indicator for the event $b_i^S \ge (1-\delta)m/u$, then $\mathbb{E}[B_i^S] \ge \frac{2d}{u}$. Finally, let $B^S = \sum_{i \in [u]}B_i^S$, then $\mathbb{E}[B^S] \ge 2d$. We will show that with probability at least $1/8 \ge 1/50$ we have $B^S \ge d$. This implies that there exist $i_1,\ldots,i_d$ such that for every $j \in [d]$, $b_{i_j}^S \le \frac{m}{u}\left(1 - \delta\right) = \frac{m}{u}\left(1 - \sqrt{\frac{\ln(u/2d)}{9m/u}}\right)$. To show $B^S \ge d$ with reasonable probability, we use the Paley-Zigmund inequality.
\[
\Pr[B^S \ge d] = \Pr\left[B^S \ge \frac{1}{2}\mathbb{E}[B^S]\right] \ge \frac{\mathbb{E}[B^S]^2}{4\mathbb{E}[(B^S)^2]} \;.
\]
Since $B_1^S,\ldots,B_u^S$ are negatively correlated, we have that $\mathbb{E}[B_i^SB_j^S] \le \mathbb{E}[B_i^S]\mathbb[B_j^S] = \mathbb{E}[B_1^S]^2$ for every $i,j \in [u]$. Moreover, as $B_1^S,\ldots,B_u^S$ are indicators,  $\mathbb{E}[(B_i^S)^2] = \mathbb{E}[B_i^S]$ for all $i \in [u]$. Therefore
\begin{equation*}
\begin{split}
\mathbb{E}[(B^S)^2] &= \sum_{i,j \in [u]}{\mathbb{E}[B_i^SB_j^S]} \le (u^2-u)\mathbb{E}[B_1^S]^2 + u\mathbb{E}[B_i^S] \\
&\le u^2\mathbb{E}[B_1^S]^2 + \mathbb{E}[B^S] = \mathbb{E}[B^S]^2 + \mathbb{E}[B^S] \le 2\mathbb{E}[B^S]^2 \;,
\end{split}
\end{equation*}
where the last inequality is due to the fact that $\mathbb{E}[B^S] \ge 2d \ge 1$. We conclude that
\[
\Pr[B^S \ge d] \ge \frac{\mathbb{E}[B^S]^2}{4\mathbb{E}[(B^S)^2]} \ge \frac{1}{8}\;.
\]
The proof of the lemma, and therefore of Theorem~\ref{th:lowerEx} is now complete.

\newcommand{\etalchar}[1]{$^{#1}$}


\section{Introduction}\label{sec1}

High energy heavy-ion collisions performed in the laboratory have been used to study the systems with hadronic or partonic degrees of freedom under extremely high temperature and density. It is commonly agreed that a new dense and hot phase of nuclear matter denoted as quark-gluon plasma (QGP), predicted by Quantum-Chromodynamics (QCD), is created where the partonic degrees of freedom come into play during the dynamical evolution of the system \cite{Wong1994, Csernai1994}. But in experiment, only the final states of produced particles can be measured, including the charged particle pseudorapidity distributions ($\frac{dN_{ch}}{d\eta}$), particle transverse momentum spectra and so on. Given those information, we are able to study the properties of QGP, particle production mechanism, particle correlations etc.. The study of the charged particle pseudorapidity distribution over a wide range of pseudorapidity ($\eta$) and its dependence on the colliding system, energy as well as centrality, are benchmark tools to constrain the models \cite{Abbas2013, Adam20162, Adam20172, Acharya2019}. It is important to  improve our understanding of the particle production mechanisms and provide insight into the partonic structure of the interacting nuclei.  At the same time, the total charged particle multiplicities and the evolution of their distributions with centrality add significant information.

In experiments, abundant data of the pseudorapidity distributions of charged particles have been measured by different experimental collaborations \cite{Abbas2013,Adam20162,Adam20172, Acharya2019,Alner1986,Alver2011,Aamodt2010,Adam2017,Adam2016,Abe1990,Arnison1983,thome1977,Alver2009,Back2003,Back2005,Aad2014,Back20052}, ranging from the small collision systems such as p+p to large ones such as Pb+Pb at different collision energies and centralities. In order to describe $\frac{dN_{ch}}{d\eta}$, experimentalists have proposed several parameterizations and utilized them to extract the total charged particle multiplicities ($N_{ch}^{tot}$) produced in the reactions \cite{Abbas2013, Adam20162, Back2003, basu2020}. Some well established  models, such as HIJING \cite{Wang1991}, AMPT \cite{Lin2005} (with and without string melting), EPOS-LHC \cite{Pierog2015}, UrQMD \cite{Mitrovski2009}, CGC based model \cite{Albacete2011}, introduced in the field of high energy heavy-ion collisions, are also adopted to understand the pseudorapidity distributions of charged particles as well as other physics.

In Refs. \cite{Sun2013, Li2014}, a multi-source thermal model, the four sources being projectile and target cylinders and the projectile and target leading particles, has been applied to describe the centrality dependence of the charged particle pseudorapidity distribution in Pb+Pb collisions. In Ref. \cite{Gao2015}, a new revised Landau hydrodynamic model following the same philosophy as in the multi-source thermal model is proposed to systematically study $\frac{dN_{ch}}{d\eta}$ produced in heavy-ion collisions. The system is assumed to be consisted of three sources, i.e., the central, target and projectile sources. The central source is described by the Landau hydrodynamic model and the target and projectile sources are assumed to emit particles in their rest frame isotropically. In Refs. \cite{jiang22018,jiang2016}, a 1+1 dimensional hydrodynamics model, which is analytically solvable,  is adopted by taking into account the collective motion to study the pseudorapidity distributions in both AA and pp collisions at currently available energies. In this model, the final particles have been categorized as the particles governed by the hydrodynamics and the leading particles. In Ref. \cite{jiang20182}, a relativistic viscous hydrodynamics model considering the longitudinal acceleration effect has been used to study the pseudorapidity distribution of charged particles and estimate the energy density from the most central AA collisions at RHIC and LHC. These models have different scenarios, but their results are in good agreement with the available experimental data.

Recently, the fireball model based on Tsallis thermodynamics had been successfully applied to describe the charged particle (pseudo)rapidity distributions produced in heavy-ion collisions \cite{Marques2015,Gao2017}. In this paper, we will systematically analyze the pseudorapidity distributions of charged particles produced in p+p($\rm \overline{p}$) collisions at energies ranging from $\sqrt{s_{NN}}=23.6$ GeV to 13 TeV, d+Au collisions at $\sqrt{s_{NN}}=$200 GeV, p+Pb collisions at $\sqrt{s_{NN}}=$5.02 TeV and A+A collisions including Cu+Cu, Au+Au, Pb+Pb and Xe+Xe at RHIC and LHC energies and different centralities. We will also extract the total charged particle multiplicities from the model and compare with the experimental data and study their dependence on centrality as well as collision energy.

The paper is organized as follows: In section II, we briefly introduce the fireball model with Tsallis thermodynamics for the pseudorapidity distribution. Then we extend the model to asymmetric collision system. In section III, we show the model results along with the experimental data of the pseudorapidity distributions of charged particles produced in collision systems selected. The corresponding total charged particle multiplicities are also extracted and studied. Furthermore, the collision energy and centrality dependence of the central position $y_0$ and its width $\sigma$ of the fireball distribution are also discussed. In section IV, a brief conclusion is drawn.

\begin{figure}[ht]
	\centering
	\begin{tabular}{c}

		\includegraphics[scale=0.72]{figure1.pdf}
	\end{tabular}
	\caption{(Color online) The pseudorapidity distribution of charged particles produced in proton-(anti)proton collisions for energies ranging from $\sqrt{s_{NN}}=23.6$ GeV to 13 TeV. The symbols are experimental data taken from Refs.  \cite{Alner1986,Alver2011,Aamodt2010,Adam2017,Adam2016,Abe1990,Arnison1983,thome1977}. The curves are the results from Eqs. (\ref{8}) and (\ref{10}).}
	\label{p+(p_)p}
\end{figure}

\section{Theoretical framework and formulas}\label{sec2}

Recently, the Tsallis distribution has been extensively and successfully applied to describe the particle spectra produced in pp, pA and AA collisions. In Refs. \cite{Zheng2016, Zheng2015}, it has been shown that different versions of the Tsallis distribution can reproduce the particle spectra resulting in different fitting parameters. In the framework of the fireball model with Tsallis thermodynamics \cite{Marques2015, Cleymans2012}, the transverse momentum distribution of particles is taken as
\begin{eqnarray}
		\frac{d^{2} N}{2 \pi p_{T} d p_{T} d y}=g V \frac{m_{T} \cosh y}{(2 \pi)^{3}}\left[1+(q-1) \frac{m_{T} \cosh y-\mu}{T}\right]^{-\frac{q}{q-1}},\label{1}
\end{eqnarray}
where $V$ is the volume, $g$ is the particle state degeneracy,  $m_{T}=\sqrt{m^{2}_{0}+p_{T}^{2}}$ is the transverse mass and $m_{0}$ is the particle rest mass, $y$ is the rapidity, $\mu$ is the chemical potential, $T$ is the temperature and $q$ is the entropic factor which measures the nonadditivity of the entropy. When $q=1$, the Tsallis distribution becomes the Boltzmann distribution. There are four parameters $V$, $\mu$, $T$ and $q$ in Eq. (\ref{1}) to describe the particle spectrum. To reduce the number of free parameters and simplify the problem, $\mu$ is assumed to be 0 for all the collision systems in this study. In the mid-rapidity $y\approx0$ region, Eq. (\ref{1}) can be written as
\begin{eqnarray}
	\frac{d^{2} N}{2 \pi p_{T} d p_{T} d y}=g V \frac{m_{T}}{(2 \pi)^{3}}\left[1+(q-1) \frac{m_{T}}{T}\right]^{-\frac{q}{q-1}}, \label{2}
\end{eqnarray}
which will be used to extract the parameters $q$ and $T$ by fitting the particle transverse momentum spectrum.

In the framework of the fireball model with Tsallis thermodynamics \cite{Marques2015, Gao2017}, the rapidity axis is populated with fireballs which follow the distribution function $\nu(y_{f})$, where $y_f$ is the rapidity of the fireball. The particles emitted from the fireball follow the Tsallis distribution Eq. (\ref{1}). Therefore, the particle distribution measured can be written as
\begin{equation}
\frac{d^{2} N}{p_{T} d p_{T} d y}= \frac{N}{A} \int_{-\infty}^{\infty} \nu\left(y_{f}\right) \frac{m_{T} \cosh \left(y-y_{f}\right)}{(2 \pi)^{2}}\left[1+(q-1) \frac{m_{T} \cosh \left(y-y_{f}\right)}{T}\right]^{-\frac{q}{q-1}} d y_{f},
 \label{3}
\end{equation}
where $N$ is the total particle multiplicity and $A$ is the normalization constant that ensures
\begin{eqnarray}\label{4}
	\int_{-\infty}^{\infty} \int_{0}^{\infty} \frac{d^{2} N}{p_Td p_{T} d y} p_Td p_{T} d y=N.
\end{eqnarray}
To deduce the formula of the rapidity density distribution $\frac{dN}{dy}$ of charged particles, we integrate Eq. (\ref{3}) over the transverse momentum $p_{T}$ and obtain \cite{Marques2015, Gao2017}
\begin{eqnarray}
	\frac{dN}{dy}&=&\frac{N}{A} \int_{-\infty}^{\infty} dy_f \nu(y_f) T \left[1+m_0 (q-1)\frac{\cosh(y-y_f)}{T}\right]^{-\frac{1}{q-1}} \nonumber\\
	&&\times \frac{-(q-2)m_0^2+2m_0T \textrm{sech}(y-y_f)+2T^2 \textrm{sech}^2(y-y_f)}{4\pi^2(q-2)(2q-3)}.\label{5}
\end{eqnarray}

In order to compare with the experimental data, one should convert Eq. (\ref{5}) to pseudorapidity density distribution $\frac{dN}{d\eta}$ of charged particles because data are usually measured in pseudorapidity space. The relation between the pseudorapidity and rapidity is \cite{Wong1994}
\begin{eqnarray}
	\frac{d y}{d \eta}\left(\eta, p_{T}\right)=\sqrt{1-\frac{m_{0}^{2}}{m_{T}^{2} \cosh ^{2} y}}.\label{6}
\end{eqnarray}
Substituting Eq. (\ref{6}) into Eq. (\ref{3}), we obtain:
\begin{equation}
\frac{d^2N}{p_Tdp_Td\eta } =  \frac{N}{A}\sqrt{1-\frac{m_0^2}{m_T^2 \cosh^2 y}}\int_{-\infty}^{\infty} \nu(y_f)\frac{m_T\cosh(y-y_f)}{(2\pi)^2}\left[1+(q-1)\frac{m_T\cosh(y-y_f)}{T}\right]^{-\frac{q}{q-1}}dy_f.
	\label{7}
\end{equation}

Because of the extra term $\sqrt{1-\frac{m_0^2}{m_T^2 \cosh^2 y}}$, Eq. (\ref{7}) cannot be analytically integrated over $p_T$. Therefore, we perform the numerical integration for the  pseudorapidity distribution
\begin{eqnarray}
	\frac{dN}{d\eta}&=& \frac{N}{A}\int_{-\infty}^{\infty}dy_f \int_0^\infty dp_T p_T   \sqrt{1-\frac{m_0^2}{m_T^2 \cosh^2 y}}  \nonumber\\
	&&\times \nu(y_f)\frac{m_T\cosh(y-y_f)}{(2\pi)^2}\left[1+(q-1)\frac{m_T\cosh(y-y_f)}{T}\right]^{-\frac{q}{q-1}},
	\label{8}
\end{eqnarray}

where
\begin{eqnarray}
	y=\frac{1}{2} \ln \left[\frac{\sqrt{p_{T}^{2} \cosh ^{2} \eta+m_{0}^{2}}+p_{T} \sinh \eta}{\sqrt{p_{T}^{2} \cosh ^{2} \eta+m_{0}^{2}}-p_{T} \sinh \eta}\right].
	\label{9}
\end{eqnarray}

The distribution $\nu(y_f)$ of the  fireballs has been assumed to be the summation of two q-Gaussian functions for the symmetric collision systems in Refs.\cite{Marques2015, Gao2017} as

\begin{equation}\label{10}
	\nu\left(y_{f}\right)=\frac{1}{\sqrt{2 \pi} \sigma}\left[1+\left(q^{\prime}-1\right) \frac{\left(y_{f}-y_{0}\right)^{2}}{2 \sigma^{2}}\right]^{-\frac{1}{q^{\prime}-1}}+\frac{1}{\sqrt{2 \pi} \sigma}\left[1+\left(q^{\prime}-1\right) \frac{\left(y_{f}+y_{0}\right)^{2}}{2 \sigma^{2}}\right]^{-\frac{1}{q^{\prime}-1}},
\end{equation}
where $y_0$ and $\sigma$ are fitting parameters referring to the central position and its width of the fireball distribution respectively, which are determined by the experimental data.  $q^{\prime}$ is another fitting parameter which can be different from $q$ for the particle transverse momentum spectrum Eq. (\ref{1}) while $q'=q$ is assumed in Refs.  \cite{Marques2015,Gao2017}. Since the symmetric and asymmetric collision systems will be considered in this work, we extend the fireball model with Tsallis thermodynamics to the asymmetric collision system and modify $\nu(y_f)$ as
\begin{equation}\label{11}
\nu\left(y_{f}\right)=\frac{1}{\sqrt{2 \pi} \sigma_{01}}\left[1+\left(q^{\prime}-1\right) \frac{\left(y_{f}-y_{01}\right)^{2}}{2 \sigma_{01}^{2}}\right]^{-\frac{1}{q^{\prime}-1}}+\frac{x}{\sqrt{2 \pi} \sigma_{02}}\left[1+\left(q^{\prime}-1\right) \frac{\left(y_{f}+y_{02}\right)^{2}}{2 \sigma_{02}^{2}}\right]^{-\frac{1}{q^{\prime}-1}},
\end{equation}
where $y_{0i}$ and $\sigma_{0i}$ ($i=1, 2$) are the central positions and their widths of the fireball distributions. $x$ is the relative weight introduced for the asymmetric collision system. For the symmetric collision system, $y_{01}=y_{02}$, $\sigma_{01}=\sigma_{02}$ and $x=1$ and Eq. (\ref{10}) is recovered. In analogy to Refs. \cite{Marques2015,Gao2017}, $q'=q$ is assumed.
\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.5]{cucu62gev.pdf}
	\caption{(Color online) The pseudorapidity distribution of charged particles produced in Cu+Cu collisions at  $\sqrt{s_{NN}}=62.4$ GeV for different centralities. The symbols are experimental data taken from Ref. \cite{Alver2009}. The curves are the results from Eqs. (\ref{8}) and (\ref{10}).}\label{Cu+Cu}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.5]{auau200gev.pdf}
	\caption{(Color online) The pseudorapidity distribution of charged particles produced in Au+Au collisions at  $\sqrt{s_{NN}}=200$ GeV for different centralities. The symbols are experimental data taken from Ref. \cite{Back2003}. The curves are the results from Eqs. (\ref{8}) and (\ref{10}).}
	\label{Au+Au}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.5]{xexe5440gev.pdf}
	\caption{(Color online) The pseudorapidity distribution of charged particles produced in Xe+Xe collisions at  $\sqrt{s_{NN}}=5.44$ TeV for different centralities. The symbols are experimental data taken from Ref. \cite{Acharya2019}. The curves are the results from Eqs. (\ref{8}) and (\ref{10}).}
	\label{Xe+Xe}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.5]{ppb5020gev.pdf}
	\caption{(Color online) The pseudorapidity distribution of charged particle produced in p+Pb collisions at  $\sqrt{s_{NN}}=5.02$ TeV for different centralities. The symbols are experimental data taken from Ref. \cite{Aad2014}. The curves are the results from Eqs. (\ref{8}) and (\ref{11}).}
	\label{p+Pb}
\end{figure}

\section{Results}
As discussed in Sec. \ref{sec2}, we need to extract the parameters $q$ and $T$ by fitting the particle transverse momentum spectra at the mid-rapidity, given by the experimental collaborations \cite{Adam2016,ptAdam2003,ptKhachatryan2010,ptAcharya2018,ptAdam2015,ptAad2016,ptArnison1982,ptAbe1988,ptAdamczyk2017,ptBack2005,ptAdler2002,ptBack2004,ptAlver2006,ptBack2003,ptAdam2016,ptAcharya2019}, using Eq. (\ref{2}).  We had pointed out that the single Tsallis distribution only can reproduce the particle spectra at low and intermediate $p_T$ region, where most of the particles are produced, in Pb+Pb and Xe+Xe collisions at LHC. Therefore, a $p_T$ cut at 7.5 GeV/c is applied for both collision systems when we fit the charged particle spectra. Since there are plenty of works to study the particle spectra with Tsallis distribution \cite{Zheng2016, Zheng2015, Zhengprd2015, pt1,pt2,pt3,pt4,pt5}, we will not show the figures here. Then we can use these two parameters and the fireball model with Tsallis thermodynamics, Eqs. (\ref{8}) and (\ref{11}) (or (\ref{10}) for the symmetric collision systems) to study the pseudorapidity distribution of charged particles for the system investigated.

In Fig.\ref{p+(p_)p}, the results of the pseudorapidity distribution of charged particles produced in proton-(anti)proton collisions, using the fireball model, from different experiments with collision energies ranging from 23.6 GeV up to 13 TeV, as well as the corresponding experimental data are shown. As one can see, the experimental results are well reproduced, analogous to the results shown in Ref. \cite{Marques2015} but with more collision energies. One should notice that the scale of the charged particle pseudorapidity density distribution $\frac{dN_{ch}}{d\eta}$ is different for each row in Fig. \ref{p+(p_)p} in order to better visualize the results. One can observe that a valley at $\eta=0$ is developed and its depth increases with the collision energy for both p+p and p+$\rm\overline{p}$ collisions. This valley is commonly related to partial transparency and limited stopping power of the colliding (anti)proton \cite{Wolschin2015}. A difference of the pseudorapidity distributions between p+p collisions and p+$\rm\overline{p}$ collisions that the valley of the pseudorapidity distribution is developed at higher collision energy in p+$\rm\overline{p}$ collisions can be easily noticed. This may reflect the fact that p and $\rm\overline{p}$ can annihilate each other when they collide while p and p cannot, so this should take off more p$\rm\overline{p}$ from zero pseudorapidity. One also can notice that the pseudorapidity distributions only cover a very limited $\eta$ range at some collision energies, e.g., 0.9 TeV of p+p collisions, because of the limited acceptance of the detectors in the experiment. This results in the large uncertainties of the central position values $y_0$ and its widths $\sigma$ of the fireball distributions. Therefore, we will not adopt these collision energies to explore the energy dependence of the central position and its width of the fireball distribution discussed later in this paper.

In Figs. \ref{Cu+Cu}, \ref{Au+Au}, \ref{Xe+Xe}, the results of the pseudorapidity distributions of charged particles produced in AA collisions, i.e., Cu+Cu collisions at $\sqrt{s_{NN}}=62.4$ GeV, Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV, Xe+Xe collisions at $\sqrt{s_{NN}}=5.44$ TeV, at different centrality bins are shown.  As one can see, the experimental results at different centralities are well reproduced. We also studied the other nucleus-nucleus collision systems, i.e., Cu+Cu collisions at $\sqrt{s_{NN}}=200$ GeV, Au+Au collisions at $\sqrt{s_{NN}}=19.6, 62.4, 130$ GeV and Pb+Pb collisions at $\sqrt{s_{NN}}=2.76, 5.02$ TeV,  at RHIC and LHC. Similar results are found. Actually, the results for Au+Au collisions at $\sqrt{s_{NN}}=130, 200$ GeV, Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV and the prediction for the central collisions of Pb+Pb at $\sqrt{s_{NN}}=5.02$ TeV were shown in Ref. \cite{Gao2017}. It is worth noting that the peak and valley structure of the charged particle pseudorapidity distributions produced in AA collisions are typically shallower and wider than those produced in proton-(anti)proton collisions at the similar collision energies. This can be associated with the dense medium produced in AA collisions and its evolution dynamics.

In Fig.\ref{p+Pb}, we studied the asymmetric collision system p+Pb at $\sqrt{s_{NN}}=5.02$ TeV. Different from the symmetric collision case, the pseudorapidity distribution is no longer symmetric respect to $\eta=0$ but obviously has forward/backward asymmetry that more charged particles are produced in the p side than the Pb side. The extended fireball model, Eqs. (\ref{8}) and (\ref{11}), is adopted to take into account asymmetric fact in this case. The results of the pseudorapidity distributions of charged particles for eight centrality bins are shown. A good agreement between the model results and the experimental data is shown. One can see that the pseudorapidity distribution of the charged particles becomes more symmetric from central to peripheral collisions because the peripheral collisions are more like the p+p collisions for the asymmetric collision system. Similar behavior is also observed in d+Au collisions at $\sqrt{s_{NN}}=200$ GeV \cite{Back20052}.

\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.5]{totalN.pdf}
	\caption{(Color online) The total charged particle multiplicity produced in Au+Au collisions at $\sqrt{s_{NN}}=62.4, 200$ GeV, Pb+Pb at $\sqrt{s_{NN}}=2.76$ TeV and Xe+Xe at  $\sqrt{s_{NN}}=5.44$ TeV versus the number of participants. Experimental data (solid circles) are taken from Refs. \cite{Abbas2013,Adam20162,Acharya2019, Alver2011}. The open circles with a cross inside are the model results. The lines are the fitting results  and their functions are indicated in the legend.}
	\label{totalN}
\end{figure}
\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.85]{totalNNpart.pdf}
	\caption{(Color online) The total charged particle multiplicity per participant pair produced in AA collisions at different centralities versus $\sqrt{s_{NN}}$. The experimental data are from RHIC (Cu+Cu and Au+Au) \cite{Alver2011} and LHC (Pb+Pb and Xe+Xe) \cite{Abbas2013,Adam20162,Adam20172,Acharya2019}. The curves are the fitting results and the corresponding function is adopted from Ref. \cite{Abbas2013}.}
	\label{totalNNpart}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.5]{figure7.pdf}
	\caption{(Color online) The collision energy dependence of the central position $y_0$ of the fireball distribution and its width $\sigma$ in proton-(anti)proton collisions (a) and (b), and the most central nucleus-nucleus collisions (c) and (d). The open circles in (a) and (b) are the model results for the p+p collisions without two humps structure in their pseudorapidity distributions because of the detector acceptance. The AA collision systems  in (c) and (d) are Au+Au(0-6\% 19.6 GeV), Au+Au(0-6\% 62.4 GeV), Au+Au(0-6\% 130 GeV), Au+Au(0-6\% 200 GeV), Cu+Cu(0-6\% 62.4 GeV), Cu+Cu(0-6\% 200 GeV), Pb+Pb(0-5\% 2.76 TeV), Pb+Pb(0-5\% 5.02 TeV), Xe+Xe(0-5\% 5.44 TeV). The solid lines are the linear fitting results. }
	\label{parameter(sqrt(s))}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[scale=1.7]{figure8.pdf}
	\caption{(Color online) The centrality dependence of the $y_0$ and $\sigma$ in Cu+Cu collisions at $\sqrt{s_{NN}}=62.4, 200$ GeV, Au+Au collisions at $\sqrt{s_{NN}}=62.4, 130, 200$ GeV, Pb+Pb collisions at $\sqrt{s_{NN}}=2.76, 5.02$ TeV and Xe+Xe collisions at $\sqrt{s_{NN}}=5.44$ TeV. The solid lines are the linear fitting results.}
	\label{parameter(centrality)}
\end{figure}

Now we turn to analyze the relation between the centrality (collision energy) and the total charged particle multiplicity produced in different collision systems. By integrating Eq. (\ref{8}) over $\eta$, the total number of charged particles produced can be extracted from the fireball model. To compare with the experimental data, the $\eta$ cut has been adopted same as the experiments when available.  Here we only show the results for AA collisions. For the pp($\rm\overline{p}$) collisions, we refer the results to Ref. \cite{Acharya2019}. While for the p(d)A collisions, we do not have enough data to conduct the analysis. In Fig. \ref{totalN}, we plot the total charged particle multiplicity versus the number of participants for the selected collision systems. As one can see that the model results are almost indistinguishable from the experimental data. Similar results are also obtained for the other collision systems investigated in this paper. The analytical function (see the legend) adopted from Ref. \cite{Adam20172} has been used to fit the total charged particle multiplicity as a function of $N_{part}$, shown by the dashed lines in the figure. This function can well reproduce all the experimental data except Xe+Xe collisions where there is a deviation for the data at central collisions. We proposed an analytical function $a(e^{N_{part}/b}-1)$, where $a$ and $b$ are fitting parameters, and it works for all the collision systems. According to the analytical functions, we cannot distinguish whether the total charged particle multiplicity increases with the number of participants exponentially or with a power-law.

In Fig. \ref{totalNNpart}, we investigate the total charged particle multiplicity per participant pair produced in AA collisions at different centralities as a function of the collision energy. One can see that the model results are in good agreement with the experimental data which is not surprising as one has seen it from Fig. \ref{totalN}. It is found that the analytical function adopted from Ref. \cite{Abbas2013}, applied only to the central collision data, can describe the experimental data well for different centralities. The fitting power $p$ decreases from central to peripheral collisions, which qualitatively reflects the fact that the dense medium created in the more central collisions causes more collisions and leads to more charged particle production per participant pair.

\renewcommand\arraystretch{1.5}
\begin{table}
	\centering
	\caption{The linear fit of $y_0$ and $\sigma$ versus centrality represented by c.}
	\begin{tabular}{>{\centering}p{1.5cm}>{\centering}p{1.7cm}>{\centering}p{2.4cm}p{2.4cm}<{\centering}}
		\hline
		System&$\sqrt{s_{NN}}$&$y_0$&$\sigma$\\
		\hline
		\multirowcell{2}{Cu+Cu}&62.4 GeV&$0.11c+1.47$&$0.48c+1.23$\\
		&200 GeV&$0.22c+1.79$&$0.41c+1.63$\\
		\hline
		\multirowcell{3}{Au+Au}&62.4 GeV&$0.25c+1.47$&$0.68c+1.18$\\
		&130 GeV&$0.26c+1.64$&$0.57c+1.38$\\
		&200 GeV&$0.22c+1.76$&$0.63c+1.46$\\
		\hline
		\multirowcell{2}{Pb+Pb}&2.76 TeV&$0.47c+2.26$&-\\
		&5.02 TeV&$1.03c+2.10$&-\\
		\hline
		\multirowcell{1}{Xe+Xe}&5.44 TeV&$0.93c+2.34$&-\\
		\hline
	\end{tabular}\label{c1}
\end{table}

We also analyze the collision energy (centrality) dependence of the central position $y_0$ and its width $\sigma$ of the fireball distribution. For the same reason as stated before, the analysis of the collision energy dependence on $y_0$ and $\sigma$ for pA collision systems is not conducted. In Fig. \ref{parameter(sqrt(s))}(a) and (b), the results of $y_0$ and $\sigma$ versus collision energy ($\ln \sqrt{s_{NN}}$) of proton-(anti)proton collisions are shown respectively. For pp collisions, the pseudorapidity distributions of charged particles do not show the two humps structure because of the limited detector acceptance of the experiments for some energies, i.e., 0.9 TeV, 2.36 TeV, 13 TeV, see Fig. \ref{p+(p_)p}. Therefore, we show the model results, which are not well constrained by the experimental data for these energies, with open circles to see the trend.  The solid circles show that both $y_0$ and $\sigma$ increase with collision energy to a critical value and then saturate to a certain value. Although there are only three collision energies for p$\rm\overline{p}$ collisions, the same trend is observed (solid squares).  In Fig. \ref{parameter(sqrt(s))}(c) and (d), we show the model results for the most central AA collisions. Different from pp$(\rm\overline{p})$, both $y_0$ and $\sigma$ show a linear dependence of collision energy for the energies investigated which demonstrates our conjecture made in Ref. \cite{Gao2017}. A linear fit has been conducted and shown in the figure respectively. For $\sigma$,  we only do the linear fit up to 2.76 TeV. For the other centralities (We put the closest centralities together because the centrality classification is not exactly the same between RHIC and LHC), the similar behavior has been observed.

In Fig. \ref{parameter(centrality)}, the centrality dependence of the $y_0$ and $\sigma$ are shown for AA collisions. A nice linear relation between $y_0$ and centrality is found for the collision energy up to 2.76 TeV. While for the reactions Pb+Pb at $\sqrt{s_{NN}}=5.02$ TeV and Xe+Xe at $\sqrt{s_{NN}}=5.44$ TeV, the linear behavior is not so great. We performed linear fits shown with the lines in Fig. \ref{parameter(centrality)}(a). The fitting parameters are listed in Table. \ref{c1}. The positive slope indicates that the stopping power is decreasing from central to peripheral collisions because there are more collisions and dissipation processes at central collisions comparing with peripheral collisions. The linear behavior of the width $\sigma$ versus centrality is also observed at RHIC. The fitting results are shown by the lines Fig. \ref{parameter(centrality)}(b) and fitting parameters are listed in Table. \ref{c1}. For the collisions at LHC, $\sigma$ shows different dependence on centrality. For Pb+Pb at  $\sqrt{s_{NN}}=2.76$ TeV, $\sigma$ tends to be a constant for all the centralities and a constant line is plotted. While the collision energy goes up, i.e., Pb+Pb at $\sqrt{s_{NN}}=5.02$ TeV and Xe+Xe at $\sqrt{s_{NN}}=5.44$ TeV,  $\sigma$ shows a U shape. We also notice a different energy dependence behavior of $\sigma$ in Fig. \ref{parameter(sqrt(s))}(d). We do not have an explanation for these different behaviors. From the results shown in Fig. \ref{parameter(centrality)}, one also can conclude that the collision system size dependence of the $y_0$ and $\sigma$ at the RHIC is weak by comparing the results obtained from Cu+Cu and Au+Au collisions at the same collision energy. Combining the results shown in Figs. \ref{totalNNpart} and \ref{parameter(sqrt(s))}, it becomes possible that one can use the fireball model to predict the charged particle pseudorapidity distributions for the system size scan program proposed recently for the STAR experiment at RHIC. A machine learning technique can provide powerful tool to this application of the fireball model.

\section{Conclusions}
In this paper, the pseudorapidity distributions of charged particles produced in different collision systems at  RHIC and LHC energies have been studied in the fireball model based on Tsallis thermodynamics. We have extended the model to the asymmetric collision systems. The theoretical results are in good agreement with the experimental data for all the collision systems investigated. The total charged particle multiplicities are extracted from model and compared with experimental data. Whether the total charged particle multiplicity increases with the number of participants exponentially or with a power-law can not be distinguished. We also investigate the collision energy and centrality dependence of the central position $y_0$ and its width $\sigma$ of the fireball distribution. In AA collisions, both the $y_0$ and $\sigma$ have a linear relation with the collision energy and centrality, in particular at RHIC energy, respectively. With the results found in this paper, a possible application of the fireball model with Tsallis thermodynamics to predict the charged particle pseudorapidity distributions for the system size scan program proposed recently for the STAR experiment at RHIC is suggested.

\section*{Acknowledgments}

This work is supported by the National Natural Science Foundation of China (Grants No. 11905120), the Fundamental Research Funds for the Central Universities No. GK201903022 and No. GK202003019 and Natural Science Basic Research Plan in Shaanxi Province of China (program No. 2020JM-289).

\section*{References}