# Robust Partial Least Squares Using Low Rank and Sparse Decomposition

Farwa Abbas Hussain Ahmad Department of Electrical and Electronic Engineering
School of Computer and Mathematical Sciences Imperial College London
University of Adelaide, Australia<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper proposes a framework for simultaneous dimensionality reduction and
regression in the presence of outliers in data by applying low-rank and sparse
matrix decomposition. For multivariate data corrupted with outliers, it is
generally hard to estimate the true low dimensional manifold from corrupted
data. The objective of the proposed framework is to find a robust estimate of
the low dimensional space of data to reliably perform regression. The
effectiveness of the proposed algorithm is demonstrated experimentally for
simultaneous regression and dimensionality reduction in the presence of
outliers in data.

###### Index Terms:

Robust, partial least squares, outliers, low-rank, sparse, regression

##  I Introduction

In classical Multivariate Linear Regression (MLR), the objective is to model
relationship between response variables ${\boldsymbol{Y}}$ and predictors
${\boldsymbol{X}}$ by fitting a linear map $\mathbf{\Theta}$ to observed
samples such that
${\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}}={\boldsymbol{X}}\mathbf{\Theta}$
where ${\boldsymbol{\Delta}_{y}}$ is the measurement error. The solution to
this problem i.e.
$\hat{\mathbf{\Theta}}=({\boldsymbol{X}}^{\mathrm{T}}{\boldsymbol{X}})^{-1}{\boldsymbol{X}}^{\mathrm{T}}{\boldsymbol{Y}}$
can be found by solving the following program.

$\displaystyle\min_{\boldsymbol{{\boldsymbol{\Delta}_{y}}},\boldsymbol{\Theta}}\|{\boldsymbol{\Delta}_{y}}\|_{F}^{2}$
(1) $\displaystyle\text{ s.t.
}{\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}}={\boldsymbol{X}}\mathbf{\Theta}.$

However, if the model is also prone to error then accounting for both the
measurement and modelling errors in the equation
${\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}}=({\boldsymbol{X}}-{\boldsymbol{\Delta}_{x}})\mathbf{\Theta}$
using Total Least Squares (TLS) corresponds to solving the objective mentioned
below.

$\displaystyle\min_{{\boldsymbol{\Delta}_{y}},{\boldsymbol{\Delta}_{x}},\boldsymbol{\Theta}}\|{\boldsymbol{\Delta}_{y}}\|_{F}^{2}+\|{\boldsymbol{\Delta}_{x}}\|_{F}^{2}$
(2) $\displaystyle\text{ s.t.
}{\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}}=({\boldsymbol{X}}-{\boldsymbol{\Delta}_{x}})\mathbf{\Theta}.$

The above program tries to find a linear map from ${\boldsymbol{X}}$ to
${\boldsymbol{Y}}$ to model the relationship between predictors and responses
while penalizing the errors in both ${\boldsymbol{X}}$ and ${\boldsymbol{Y}}$.
In most practical cases, there are strong correlations between subsets of
predictors. Therefore, they tend to span a low dimensional subspace and hence
the covariance matrix turns out to be singular, making the linear regression
problem ill-conditioned [1]. To resolve this problem, one possible solution is
to project each sample onto the low-dimensional subspace spanned by
uncorrelated predictors and form a smaller set of new orthogonal predictors,
each being the linear combination of original set of predictors. This can be
done by finding singular value decomposition of ${\boldsymbol{X}}$ such that
${\boldsymbol{X}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}$
where ${\boldsymbol{Q}}$ and ${\boldsymbol{\Lambda}_{x}}$ are the scores and
loadings matrices respectively. The new set of orthogonal predictors in
${\boldsymbol{Q}}$ can be regressed to ${\boldsymbol{Y}}$ by regression
equation
${\boldsymbol{Y}}={\boldsymbol{Q}}\mathbf{\Theta}+{\boldsymbol{\Delta}_{y}}$
and the solution to this new problem is
$\tilde{\mathbf{\Theta}}=({\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Q}})^{-1}{\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Y}}$.
It must be noted that our new covariance matrix is now invertible even in the
presence of strong correlations between predictors. This is the main idea
behind Principal Component Regression (PCR) [2]. However PCA, being an
unsupervised dimensionality reduction approach, computes latent components
while completely disregarding the information obtained from responses and
hence, the predictive information can probably be discarded and noise may be
preserved inducing a poor regression model [3].

In order to estimate the low dimensional subspace by taking into account the
information provided by both predictors and responses, we aim to find a low
rank approximation of both ${\boldsymbol{X}}$ and ${\boldsymbol{Y}}$ by
solving the following optimization problem.

$\displaystyle\min_{{\boldsymbol{\Delta}_{y}},{\boldsymbol{\Delta}_{x}},\boldsymbol{\Theta}}\|{\boldsymbol{\Delta}_{y}}\|_{F}^{2}+\|{\boldsymbol{\Delta}_{x}}\|_{F}^{2}+\|{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\|_{*}+\|{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}\|_{*}$
(3) $\displaystyle\text{ s.t.
}{\boldsymbol{X}}-{\boldsymbol{\Delta}_{x}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}},{\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}.\hskip
28.45274pt$

We can argue that Partial Least Squares (PLS) [4] is an iterative algorithm
that fundamentally tries to solve the above mentioned objective by maximizing
the covariance between predictors and responses. It aims to find the first
latent component ${\vec{t}}={\boldsymbol{X}}{\vec{w}}$ by finding weights
${\vec{w}}$ while maximizing the covariance
$\operatorname{cov}({\boldsymbol{X}}{\vec{w}},{\boldsymbol{Y}})$ subject to
${\vec{w}}^{\mathrm{T}}{\vec{w}}=1$. The subsequent latent components are
found by singular value decomposition of residual matrix that is updated
iteratively. Unlike conventional PLS method that computes the scores matrix as
${\boldsymbol{Q}}={\boldsymbol{X}}{\boldsymbol{W}}({\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}{\boldsymbol{W}})^{-1}$,
we observe that it can rather be estimated in a flexible way by optimizing the
objective in (3) without imposing any explicit form. Unlike PCA that induces
an orthogonal projection
${\boldsymbol{Q}}({\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Q}})^{-1}{\boldsymbol{Q}}^{\mathrm{T}}$
matrix, PLS algorithm instead induces an oblique projection
${\boldsymbol{W}}({\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}{\boldsymbol{W}})^{-1}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}$
matrix in a supervised way [5]. We argue that the optimization program in (3)
tries to find the low dimensional manifold similar to the one found by PLS
iterative algorithm without finding any oblique projections. Moreover, we
claim that we can achieve the same objective without requiring a weight matrix
${\boldsymbol{W}}$ for regression or dimensionality reduction.

##  II Related Work

PLS algorithm has widely been used in chemometrics [6, 7], clinical chemistry
[8, 9], industrial process control [10, 11, 12] and many other applications
[13, 14, 15]. However, its effectiveness is restricted due its extreme
sensitivity to outliers in data [16]. Several methods have been proposed in
past decades to develop a robust PLS model that include [16, 17, 18, 19, 20].
Primarily, these methods have proposed to use a robust estimate of covariance
matrix to mitigate the effect of outliers in data. Two most commonly used
methods to estimate a robust covariance matrix are Orthogonalized
Gnanadesikan-Kettenring (OGK) [21] and Minimum Covariance Determinant (MCD)
[22]. However, the utility of these methods is restricted to the cases where
the number of samples are very large as compared to the number of predictors
which is unrealistic in many situations [23].

It can be observed that, instead of minimizing
$\|{\boldsymbol{\Delta}_{x}}\|_{F}^{2}$ and
$\|{\boldsymbol{\Delta}_{y}}\|_{F}^{2}$ as in conventional PLS method, we can
minimize instead $\|{\boldsymbol{\Delta}_{x}}\|_{1}$ and
$\|{\boldsymbol{\Delta}_{y}}\|_{1}$ to get rid of outliers in data. Motivated
by this observation, we propose to employ a new low-rank and sparse
decomposition framework for simultaneously decomposing both predictor and
response matrices. We refer to our proposed algorithm as Robust PLS (RPLS). We
have shown experimentally that by allowing sparse errors in both predictors
and responses, the low-dimensional space can be discovered reliably in the
presence of outliers.

In contrast to iterative PLS scheme, we propose a new way to find the
decompositions
${\boldsymbol{X}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}+{\boldsymbol{\Delta}_{x}}$
and
${\boldsymbol{Y}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}+{\boldsymbol{\Delta}_{y}}$
that does not involve any power method or oblique projections at all. We
introduce a novel framework to perform regression by finding projections onto
the latent space estimated by a robust method. More precisely, we make the
assumption that after removing errors in both ${\boldsymbol{X}}$ and
${\boldsymbol{Y}}$, the low rank matrices
${\boldsymbol{X}}-{\boldsymbol{\Delta}_{x}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}$
and
${\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}$
can be related through regression equation as
${\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\mathbf{\Theta}$.
Moreover, if ${\boldsymbol{Q}}$ is constrained to have orthonormal columns
i.e. ${\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Q}}={\boldsymbol{I}}$ then our
regression equation reduces to
${\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}={\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\mathbf{\Theta}$.
Our method is straightforward and also handles outliers in data. We show the
usefulness of proposed framework for simultaneous regression and
dimensionality reduction. We, further, show that by performing regression
after projecting onto the low dimensional space, the response variable can be
predicted in a reliable and robust way.

###  II-A Our Contributions

In this paper, we propose a novel method to formulate PLS regression problem
by introducing a low-rank and sparse decomposition framework. To the best of
our knowledge, this is the first work that solves the PLS regression problem
in an optimization framework by effectively penalizing the outliers in data.
In particular, we make the following contributions:

  1. 1.

We propose a novel way to formulate PLS regression problem as a low rank and
sparse decomposition problem by minimizing the proposed objective using
Alternating Direction Method of Multipliers (ADMM).

  2. 2.

We propose a regression approach based on projection onto a latent space.

  3. 3.

Finally, we demonstrate experimentally the superior performance of our
proposed method on both synthetic and real world datasets corrupted with
outliers.

##  III Mathematical Formulation

Motivated by the approach discussed in [24], we develop a robust PLS model in
low-rank and sparse decomposition framework. More formally, we aim to find
${\boldsymbol{Q}}\in\mathbb{R}^{n\times k}$ (such that
${\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Q}}={\boldsymbol{I}}$),
${\boldsymbol{\Lambda}_{x}}\in\mathbb{R}^{p\times k}$ and
${\boldsymbol{\Lambda}_{y}}\in\mathbb{R}^{r\times k}$ such that
${\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\in\mathbb{R}^{n\times
p}$ and
${\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}\in\mathbb{R}^{n\times
r}$ are low-rank matrices where $k$ is an upper bound on the rank of
${\boldsymbol{X}}$ and ${\boldsymbol{Y}}$ i.e.,
$\text{rank}({\boldsymbol{X}})\leq k$ and $\text{rank}({\boldsymbol{Y}})\leq
k$. We aim to find a low rank approximation for ${\boldsymbol{X}}$ and
${\boldsymbol{Y}}$ while ensuring the error matrices
${\boldsymbol{\Delta}_{x}}$ and ${\boldsymbol{\Delta}_{y}}$ to be sparse. The
problem can be mathematically formulated as follows.

$\displaystyle\min_{{\boldsymbol{Q}},{\boldsymbol{\Lambda}_{x}},{\boldsymbol{\Lambda}_{y}},{\boldsymbol{\Delta}_{x}},{\boldsymbol{\Delta}_{y}}}\|{\boldsymbol{\Delta}_{x}}\|_{1}+\|{\boldsymbol{\Delta}_{y}}\|_{1}+\lambda_{1}\|{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\|_{*}+\lambda_{2}\|{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}\|_{*}$
$\displaystyle\text{ subject to
}{\boldsymbol{X}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}+{\boldsymbol{\Delta}_{x}}\text{
,
}{\boldsymbol{Y}}={\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}+{\boldsymbol{\Delta}_{y}}\text{
and }{\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Q}}={\boldsymbol{I}}.$ (4)

In this paper, we restrict our attention to solve the optimization program in
(III). The partial augmented Lagrangian function of the program in (III) can
be written as follows.

$\displaystyle\mathcal{L}({\boldsymbol{Q}},{\boldsymbol{\Lambda}_{x}},{\boldsymbol{\Lambda}_{y}},{\boldsymbol{\Delta}_{x}},{\boldsymbol{\Delta}_{y}},{\boldsymbol{L}},{\boldsymbol{M}})=\|{\boldsymbol{\Delta}_{x}}\|_{1}+\lambda_{1}\|{\boldsymbol{\Lambda}_{x}}\|_{*}+\|{\boldsymbol{\Delta}_{y}}\|_{1},$
$\displaystyle+\lambda_{2}\|{\boldsymbol{\Lambda}_{y}}\|_{*}+\langle{\boldsymbol{L}},{\boldsymbol{X}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{x}}\rangle+\frac{\alpha_{1}}{2}\|{\boldsymbol{X}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{x}}\|_{F}^{2},$
$\displaystyle+\langle{\boldsymbol{M}},{\boldsymbol{Y}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{y}}\rangle+\frac{\alpha_{2}}{2}\|{\boldsymbol{Y}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{y}}\|_{F}^{2},$
(5)

where ${\boldsymbol{L}}$ and ${\boldsymbol{M}}$ are the matrices of Lagrange
multipliers and $\alpha_{1}$ and $\alpha_{2}$ are penalty parameters for first
two constraints. The last constraint regarding the orthogonality of
${\boldsymbol{Q}}$ is enforced separately while computing the update for
${\boldsymbol{Q}}$ . $\lambda_{1}$ and $\lambda_{2}$ are weights of the second
and fourth terms respectively.

The objective function defined above is non-convex therefore we resort to
alternative minimization through non-convex ADMM. By keeping all other
matrices to be fixed, the update for each individual matrix in the objective
function can be obtained. First, the update for ${\boldsymbol{Q}}$ can be
computed while keeping all other matrices fixed.

$\displaystyle{\boldsymbol{Q}}$
${}^{(k+1)}=\arg\min_{\boldsymbol{Q}}\mathcal{L}({\boldsymbol{Q}},{\boldsymbol{\Lambda}_{x}}^{(k)},{\boldsymbol{\Lambda}_{y}}^{(k)},{\boldsymbol{\Delta}_{x}}^{(k)},{\boldsymbol{\Delta}_{y}}^{(k)},{\boldsymbol{L}}^{(k)},{\boldsymbol{M}}^{(k)})$
$\displaystyle=\arg\min_{\boldsymbol{Q}}\langle{\boldsymbol{L}},{\boldsymbol{X}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{x}}\rangle+\frac{\alpha_{1}}{2}\|{\boldsymbol{X}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{x}}\|_{F}^{2}$
$\displaystyle+\langle{\boldsymbol{M}},{\boldsymbol{Y}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{y}}\rangle+\frac{\alpha_{2}}{2}\|{\boldsymbol{Y}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{y}}\|_{F}^{2}$
$\displaystyle=\arg\min_{\boldsymbol{Q}}\frac{\alpha_{1}}{2}\Big{\|}\frac{{\boldsymbol{L}}}{\alpha_{1}}+({\boldsymbol{X}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{x}})\Big{\|}_{F}^{2}$
$\displaystyle+\frac{\alpha_{2}}{2}\Big{\|}\frac{{\boldsymbol{M}}}{\alpha_{2}}+({\boldsymbol{Y}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}-{\boldsymbol{\Delta}_{y}})\Big{\|}_{F}^{2}$
$\displaystyle=\arg\min_{\boldsymbol{Q}}\frac{\alpha_{1}}{2}\Big{\|}(\frac{{\boldsymbol{L}}}{\alpha_{1}}+{\boldsymbol{X}}-{\boldsymbol{\Delta}_{x}})-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\Big{\|}_{F}^{2}$
$\displaystyle+\frac{\alpha_{2}}{2}\Big{\|}(\frac{{\boldsymbol{M}}}{\alpha_{2}}+{\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}})-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}\Big{\|}_{F}^{2}$
$\displaystyle=\arg\min_{\boldsymbol{Q}}\frac{\alpha_{1}}{2}\|{\boldsymbol{B}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\|_{F}^{2}+\frac{\alpha_{2}}{2}\|{\boldsymbol{A}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}\|_{F}^{2}.$

By substituting
${\boldsymbol{B}}=\frac{{\boldsymbol{L}}}{\alpha_{1}}+{\boldsymbol{X}}-{\boldsymbol{\Delta}_{x}}$
and
${\boldsymbol{A}}=\frac{{\boldsymbol{M}}}{\alpha_{2}}+{\boldsymbol{Y}}-{\boldsymbol{\Delta}_{y}}$
above we get the last equation. By expanding the squared terms we get the
following.

$\displaystyle{\boldsymbol{Q}}^{(k+1)}$
$\displaystyle=\arg\min_{\boldsymbol{Q}}-\alpha_{1}\langle{\boldsymbol{B}},{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\rangle-\alpha_{2}\langle{\boldsymbol{A}},{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}\rangle$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\alpha_{1}\langle{\boldsymbol{B}},{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}\rangle+\alpha_{2}\langle{\boldsymbol{A}},{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}\rangle$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\langle\alpha_{1}{\boldsymbol{B}}{\boldsymbol{\Lambda}_{x}},{\boldsymbol{Q}}\rangle+\langle\alpha_{2}{\boldsymbol{A}}{\boldsymbol{\Lambda}_{y}},{\boldsymbol{Q}}\rangle$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\langle\alpha_{1}{\boldsymbol{B}}{\boldsymbol{\Lambda}_{x}}+\alpha_{2}{\boldsymbol{A}}{\boldsymbol{\Lambda}_{y}},{\boldsymbol{Q}}\rangle$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\langle{\boldsymbol{D}},{\boldsymbol{Q}}\rangle,$

where
${\boldsymbol{D}}=\alpha_{1}{\boldsymbol{B}}{\boldsymbol{\Lambda}_{x}}+\alpha_{2}{\boldsymbol{A}}{\boldsymbol{\Lambda}_{y}}$.
Now we want to minimize the following objective while simultaneously enforcing
${\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Q}}={\boldsymbol{I}}$.

$\displaystyle{\boldsymbol{Q}}^{(k+1)}$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\langle{\boldsymbol{D}},{\boldsymbol{Q}}\rangle$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\langle{\boldsymbol{U}}{\boldsymbol{S}}{\boldsymbol{V}}^{\mathrm{T}},{\boldsymbol{Q}}\rangle$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\langle{\boldsymbol{S}},{\boldsymbol{U}}^{\mathrm{T}}{\boldsymbol{Q}}{\boldsymbol{V}}\rangle$
$\displaystyle=\arg\max_{\boldsymbol{Q}}\langle{\boldsymbol{S}},{\boldsymbol{J}}\rangle\hskip
36.98866pt\text{ s.t.
}{\boldsymbol{Q}}^{\mathrm{T}}{\boldsymbol{Q}}={\boldsymbol{I}}.$

In the set of equations above
${\boldsymbol{D}}={\boldsymbol{U}}{\boldsymbol{S}}{\boldsymbol{V}}^{\mathrm{T}}$
is the singular value decomposition of matrix ${\boldsymbol{D}}$. The matrix
${\boldsymbol{J}}={\boldsymbol{U}}^{\mathrm{T}}{\boldsymbol{Q}}{\boldsymbol{V}}$
is an orthonormal matrix, as it is a product of orthonormal matrices, and thus
the expression is maximised when ${\boldsymbol{J}}$ equals identity matrix
${\boldsymbol{I}}$. Thus

$\displaystyle{\boldsymbol{J}}$
$\displaystyle={\boldsymbol{U}}^{\mathrm{T}}{\boldsymbol{Q}}{\boldsymbol{V}}={\boldsymbol{I}}$
$\displaystyle{\boldsymbol{Q}}^{(k+1)}$
$\displaystyle={\boldsymbol{U}}{\boldsymbol{V}}^{\mathrm{T}}.$ (6)

Similarly, the updates for ${\boldsymbol{\Lambda}_{x}}$ and
${\boldsymbol{\Lambda}_{y}}$ can be computed as follows.

$\displaystyle{\boldsymbol{\Lambda}_{x}}^{(k+1)},{\boldsymbol{\Lambda}_{y}}^{(k+1)}$
$\displaystyle=\arg\min_{{\boldsymbol{\Lambda}_{x}},{\boldsymbol{\Lambda}_{y}}}\frac{\lambda_{1}}{\alpha_{1}}\|{\boldsymbol{\Lambda}_{x}}\|_{*}+\frac{1}{2}\|{\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{Q}}-{\boldsymbol{\Lambda}_{x}}\|_{F}^{2}$
$\displaystyle+\frac{\lambda_{2}}{\alpha_{2}}\|{\boldsymbol{\Lambda}_{y}}\|_{*}+\frac{1}{2}\|{\boldsymbol{A}}^{\mathrm{T}}{\boldsymbol{Q}}-{\boldsymbol{\Lambda}_{y}}\|_{F}^{2}$
$\displaystyle{\boldsymbol{\Lambda}_{x}}^{(k+1)}$
$\displaystyle=\tilde{{\boldsymbol{U}}}\mathcal{S}_{\lambda_{1}/\alpha_{1}}(\tilde{\mathbf{\Sigma}})\tilde{{\boldsymbol{V}}}^{\mathrm{T}}$
(7) $\displaystyle{\boldsymbol{\Lambda}_{y}}^{(k+1)}$
$\displaystyle=\widehat{{\boldsymbol{U}}}\mathcal{S}_{\lambda_{2}/\alpha_{2}}(\widehat{\mathbf{\Sigma}})\widehat{{\boldsymbol{V}}}^{\mathrm{T}},$
(8)

where
${\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{Q}}=\tilde{{\boldsymbol{U}}}\tilde{\mathbf{\Sigma}}\tilde{{\boldsymbol{V}}^{\mathrm{T}}}$
and
${\boldsymbol{A}}^{\mathrm{T}}{\boldsymbol{Q}}=\widehat{{\boldsymbol{U}}}\widehat{\mathbf{\Sigma}}\widehat{{\boldsymbol{V}}}^{\mathrm{T}}$
are the singular value decompositions of
${\boldsymbol{B}}^{\mathrm{T}}{\boldsymbol{Q}}$ and
${\boldsymbol{A}}^{\mathrm{T}}{\boldsymbol{Q}}$ respectively. The above result
is obtained from the proximal operator of nuclear norm minimization problem
where $\mathcal{S}_{\epsilon}$ is the soft thresholding operator as defined in
(9).

$\displaystyle\mathcal{S}_{\epsilon}(\boldsymbol{K})$
$\displaystyle=\begin{cases}K_{ij}-\epsilon&K_{ij}>\epsilon\\\
K_{ij}+\epsilon&K_{ij}<\epsilon\\\ 0&\text{otherwise}\end{cases}$ (9)

By keeping all the other matrices to be fixed we can find the update equations
for ${\boldsymbol{\Delta}_{x}}$ and ${\boldsymbol{\Delta}_{y}}$ by using the
proximal operator of $\ell_{1}$ norm minimization problem.

$\displaystyle{\boldsymbol{\Delta}_{x}}^{(k+1)},{\boldsymbol{\Delta}_{y}}^{(k+1)}=\arg\min_{{\boldsymbol{\Delta}_{x}},{\boldsymbol{\Delta}_{y}}}+\frac{1}{2}\Big{\|}({\boldsymbol{X}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{x}}^{\mathrm{T}}+\frac{{\boldsymbol{L}}}{\alpha_{1}})-{\boldsymbol{\Delta}_{x}}\Big{\|}_{F}^{2}$
$\displaystyle+\frac{1}{2}\Big{\|}({\boldsymbol{Y}}-{\boldsymbol{Q}}{\boldsymbol{\Lambda}_{y}}^{\mathrm{T}}+\frac{{\boldsymbol{M}}}{\alpha_{2}})-{\boldsymbol{\Delta}_{y}}\Big{\|}_{F}^{2}+\frac{1}{\alpha_{1}}\|{\boldsymbol{\Delta}_{x}}\|_{1}+\frac{1}{\alpha_{2}}\|{\boldsymbol{\Delta}_{y}}\|_{1}$
$\displaystyle{\boldsymbol{\Delta}_{x}}^{(k+1)}$
$\displaystyle=\mathcal{S}_{1/\alpha_{1}}\left({\boldsymbol{X}}-{\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{x}}^{(k+1)\mathrm{T}}+\frac{{\boldsymbol{L}}^{(k)}}{\alpha_{1}}\right)$
(10) $\displaystyle{\boldsymbol{\Delta}_{y}}^{(k+1)}$
$\displaystyle=\mathcal{S}_{1/\alpha_{2}}\left({\boldsymbol{Y}}-{\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{y}}^{(k+1)\mathrm{T}}+\frac{{\boldsymbol{M}}^{(k)}}{\alpha_{2}}\right)$
(11)

In the end the matrices of Lagrange multipliers are updated as follows.

$\displaystyle{\boldsymbol{L}}^{(k+1)}$
$\displaystyle={\boldsymbol{L}}^{(k)}+\alpha_{1}^{(k)}({\boldsymbol{X}}-{\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{x}}^{(k+1)\mathrm{T}}-{\boldsymbol{\Delta}_{x}}^{(k+1)})$
(12) $\displaystyle{\boldsymbol{M}}^{(k+1)}$
$\displaystyle={\boldsymbol{M}}^{(k)}+\alpha_{2}^{(k)}({\boldsymbol{Y}}-{\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{y}}^{(k+1)\mathrm{T}}-{\boldsymbol{\Delta}_{y}}^{(k+1)}).$
(13)

The complete algorithm is shown in Algorithm 1. Here
$\hat{{\boldsymbol{X}}}^{(k+1)}$ and $\hat{{\boldsymbol{Y}}}^{(k+1)}$ are the
approximations of ${\boldsymbol{X}}$ and ${\boldsymbol{Y}}$ at iteration $k+1$
i.e.
$\hat{{\boldsymbol{X}}}^{(k+1)}={\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{x}}^{(k+1)\mathrm{T}}+{\boldsymbol{\Delta}_{x}}^{(k+1)}$
and
$\hat{{\boldsymbol{Y}}}^{(k+1)}={\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{y}}^{(k+1)\mathrm{T}}+{\boldsymbol{\Delta}_{y}}^{(k+1)}$.

Algorithm 1 RPLS via ADMM

Input: $\lambda_{1},\lambda_{2},tol$
Initialization:
$\alpha_{1}^{(0)},\alpha_{2}^{(0)},\alpha_{max},\rho,{\boldsymbol{Q}}^{(0)}=\operatorname{eye}(n,k),$
${\boldsymbol{\Lambda}_{x}}^{(0)}=\boldsymbol{0},{\boldsymbol{\Lambda}_{y}}^{(0)}=\boldsymbol{0},{\boldsymbol{\Delta}_{x}}^{(0)}=\boldsymbol{0},{\boldsymbol{\Delta}_{y}}^{(0)}=\boldsymbol{0}$
Output:
${\boldsymbol{Q}}$,${\boldsymbol{\Lambda}_{x}}$,${\boldsymbol{\Lambda}_{y}}$,${\boldsymbol{\Delta}_{x}}$,${\boldsymbol{\Delta}_{y}}$

while not converged do

Update ${\boldsymbol{Q}}^{(k+1)}$ by (6)

Update ${\boldsymbol{\Lambda}_{x}}^{(k+1)}$ by (7)

Update ${\boldsymbol{\Lambda}_{y}}^{(k+1)}$ by (8)

Update ${\boldsymbol{\Delta}_{x}}^{(k+1)}$ by (10)

Update ${\boldsymbol{\Delta}_{y}}^{(k+1)}$ by (11)

Update ${\boldsymbol{L}}^{(k+1)}$ by (12)

Update ${\boldsymbol{M}}^{(k+1)}$ by (13)

Update $\alpha_{1}^{(k+1)}$ by

$\alpha_{1}^{(k+1)}=\min(\rho\alpha_{1}^{(k)},\alpha_{max})$

Update $\alpha_{2}^{(k+1)}$ by

$\alpha_{2}^{(k+1)}=\min(\rho\alpha_{2}^{(k)},\alpha_{max})$

Check convergence by

$\|{\boldsymbol{X}}-\hat{{\boldsymbol{X}}}^{(k+1)}\|_{F}+\|{\boldsymbol{Y}}-\hat{{\boldsymbol{Y}}}^{(k+1)}\|_{F}<$
$tol$

end while

##  IV COMPUTATION COMPLEXITY ANALYSIS

In the above algorithm the computational cost in all the updates is primarily
due to svd and matrix multiplications. In equation (6), for matrix
${\boldsymbol{D}}\in\mathbb{R}^{n\times k}$ the computational complexity of
full svd is $\mathcal{O}(nk^{2})$ $(n>k)$. Similarly the next two updates (7)
, (8) involve svd having costs $\mathcal{O}(pk^{2})$, and
$\mathcal{O}(rk^{2})$ respectively. The cost of matrix multiplications
${\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{x}}^{(k+1)\mathrm{T}}$ in (10)
and (12) is $\mathcal{O}(nkp)$. Similarly, the cost associated with
${\boldsymbol{Q}}^{(k+1)}{\boldsymbol{\Lambda}_{y}}^{(k+1)\mathrm{T}}$ in (11)
and (13) is $\mathcal{O}(nkr)$. Hence the total cost is
$\mathcal{O}(T*((k^{2}(n+p+r)+2nk(p+r))))$ where $T$ is total number of
iterations.

##  V EXPERIMENTAL DETAILS

In this section we compare the performance of our proposed algorithm on both
synthetic and real world datasets. The experiments on synthetic data show the
performance on dimensionality reduction and the real dataset is evaluated for
both dimensionality reduction and regression. All the experiments have been
performed on MATLAB 2018a with Core i5 CPU (2.50 GHz), 8 GB RAM and Windows 10
operating system.

###  V-A Results on Synthetic dataset

For randomly generated synthetic data, the predictor matrix is generated
randomly having dimensions $150\times 40$ i.e. we have $150$ samples and $40$
predictors. The number of orthogonal components is chosen to be $5$. To
simulate multicollinearity in data, out of $40$ predictors, some of the
predictors are linearly combined to make other predictors. Sample size is
deliberately made larger as compared to the number of predictors so that
existing methods for robust PLS can be applied. We compare our results with
traditional PLS and PLS models robustified by using frequently used robust
estimates of covariance matrix i.e. PLS-OGK and PLS-MCD. Four responses are
being generated from linear functions of randomly selected predictors and are
then stacked in a $150\times 4$ dimensional response matrix.

To simulate the effect of presence of outliers, we have introduced two
different types of outliers in data. First ones are the sparse outliers, that
are added randomly to some of the predictors and responses. Second type of
outliers perturb lowest $10\%$ values in the response to $10$ times their
value. To increase reliability, this experiment has been repeated four times
and results for all four experiments are shown in Figure 1. By introducing
outliers, the traditional PLS algorithm fails to correctly discover the low
dimensional structure in data as depicted in Figure 1. The plotted ellipses
show the region of $95\%$ confidence determined by different algorithms.

Figure 1: Results of PLS, PLS-OGK, PLS-MCD and RPLS on four times randomly
generated synthetic data that is corrupted by outliers. RPLS is able to
correctly determine the low-dimensional space in all cases.

###  V-B Results on Real dataset

Partial least squares regression has extensively been applied in spectroscopy
to reveal the spectral content of various chemical compositions. In past
decades, PLS has achieved most promising results to predict the octane number
of gasoline using near infrared (NIR) spectroscopic techniques [25, 26]. We
test the performance of our algorithm on a dataset containing the NIR spectra
of 60 gasoline samples with known octane numbers and covers the range from 900
to 1700 nm in 2 nm intervals. The size of complete predictor matrix is
$60\times 401$ and number of orthogonal components are chosen to be $10$.

Since the number of samples are very less as compared to the dimensionality of
data, therefore, PLS-OGK and PLS-MCD could not be applied. The data has been
corrupted with second type of outliers i.e. outliers introduced by adding a
perturbation on $10\%$ lowest values and results produced by both PLS and RPLS
are compared in Figure 2.

Figure 2: RPLS has correctly identified the underlying manifold representing
the spectral measurements of gasoline samples.

After learning the PLS model with $80\%$ of the data, the remaining $20\%$
data is tested to predict the response. Results for regression are compared
for MLR, PCR, and PLSR by using their original regression equations. The
proposed regression results i.e. regression by projection are also compared
for both PLS and RPLS. All these results are summarized in Table I. It can be
seen that MLR is unable to produce reliable estimate of the response. PCR
produces better results but it completely disregards information in response
while finding projections. PLS projects the data in a supervised way but its
estimates are degraded due to outliers. RPLS performs well as compared to all
three since it finds robust projections in a supervised manner. NMSE defined
by
$\|{\boldsymbol{Y}}_{true}-{\boldsymbol{Y}}_{est}\|_{2}/\|{\boldsymbol{Y}}_{true}\|_{2}$
can be compared for predictions obtained below.

Octane Num | MLR | PCR | PLS | PLSR | RPLS
---|---|---|---|---|---
88.45 | 74.83 | 87.97 | 81.61 | 94.73 | 90.12
88.70 | 111.46 | 89.23 | 88.48 | 95.41 | 89.09
88.10 | 32.23 | 91.80 | 110.69 | 89.62 | 88.22
87.60 | 80.89 | 86.17 | 85.28 | 83.17 | 87.20
88.35 | 65.94 | 90.63 | 87.10 | 89.13 | 88.43
85.10 | 112.51 | 87.16 | 91.89 | 86.98 | 85.05
85.10 | 108.17 | 87.78 | 76.13 | 87.75 | 85.27
84.70 | 69.04 | 84.02 | 72.18 | 82.19 | 83.92
87.20 | 103.74 | 94.09 | 81.34 | 91.61 | 86.67
86.60 | 101.65 | 87.13 | 85.44 | 81.50 | 85.57
89.60 | 106.25 | 88.93 | 91.66 | 81.61 | 88.84
87.10 | 79.30 | 86.46 | 94.23 | 82.25 | 87.66
NMSE | 0.271670 | 0.029946 | 0.101360 | 0.053032 | 0.0081

TABLE I: Predictions obtained by MLR, PCR, PLS, PLSR, and RPLS.

##  VI Conclusion

This paper proposed a novel framework for dimensionality reduction and
regression by finding a low-rank and sparse decompositions for both predictors
and responses. Experiments have shown that the proposed algorithm performs
better as compared to the other techniques in presence of outliers and
corruptions. We further show that our method approximates the regression
function better than the existing algorithms and hence achieves smallest error
while making predictions on unseen data in the presence of outliers.

## References

  * [1] Herman Wold.  Nonlinear estimation by iterative least squares procedures.  F. N. David and J. Neyman, Eds., Research Papers in Statistics, Festschrift for J. Neyman, 1966.

  * [2] Sung H Park.  Collinearity and optimal restrictions on regression parameters for estimating responses.  Technometrics, 23(3):289–295, 1981.
  * [3] Paul Geladi and Bruce R Kowalski.  Partial least-squares regression: a tutorial.  Analytica chimica acta, 185:1–17, 1986.
  * [4] Michel Tenenhaus.  La régression PLS: théorie et pratique.  Editions technip, 1998.
  * [5] Aloke Phatak and Sijmen De Jong.  The geometry of partial least squares.  Journal of Chemometrics: A Journal of the Chemometrics Society, 11(4):311–338, 1997.
  * [6] Ahmad Mani-Varnosfaderani, Samaneh Ehsani, and Yadollah Yamini.  Investigating the effects of chemical composition of motor oils on their viscosity indices using gas chromatography and chemometrics techniques.  Petroleum Science and Technology, pages 1–9, 2019.
  * [7] Pedro Sousa Sampaio, Andreia Soares, Ana Castanho, Ana Sofia Almeida, Jorge Oliveira, and Carla Brites.  Optimization of rice amylose determination by nir-spectroscopy using pls chemometrics algorithms.  Food chemistry, 242:196–204, 2018.
  * [8] Jeffrey W Hall and Alan Pollard.  Near-infrared spectrophotometry: a new dimension in clinical chemistry.  Clinical chemistry, 38(9):1623–1631, 1992.
  * [9] M Grafen, K Nalpantidis, D Ihrig, HM Heise, and A Ostendorf.  Performance testing of a mid-infrared spectroscopic system for clinical chemistry applications utilising an ultra-broadband tunable ec-qcl radiation source.  In Optical Diagnostics and Sensing XVI: Toward Point-of-Care Diagnostics, volume 9715, page 971517. International Society for Optics and Photonics, 2016.
  * [10] Ao Chen, Hongpeng Zhou, Yujian An, and Wei Sun.  Pca and pls monitoring approaches for fault detection of wastewater treatment process.  In 2016 IEEE 25th International Symposium on Industrial Electronics (ISIE), pages 1022–1027. IEEE, 2016.
  * [11] Jie Dong, Kai Zhang, Ya Huang, Gang Li, and Kaixiang Peng.  Adaptive total pls based quality-relevant process monitoring with application to the tennessee eastman process.  Neurocomputing, 154:77–85, 2015.
  * [12] Tao Li, Yihai He, and Chunling Zhu.  Big data oriented macro-quality index based on customer satisfaction index and pls-sem for manufacturing industry.  In 2016 International Conference on Industrial Informatics-Computing Technology, Intelligent Technology, Industrial Information Integration (ICIICII), pages 181–186. IEEE, 2016.
  * [13] V Casal, PJ Martin-Alvarez, and T Herraiz.  Comparative prediction of the retention behaviour of small peptides in several reversed-phase high-performance liquid chromatography columns by using partial least squares and multiple linear regression.  Analytica Chimica Acta, 326(1-3):77–84, 1996.
  * [14] Manuel J Sanchez-Franco.  Exploring the influence of gender on the web usage via partial least squares.  Behaviour & Information Technology, 25(1):19–36, 2006.
  * [15] JW Einax, O Kampe, and D Truckenbrodt.  Assessing the deposition and remobilisation behaviour of metals between river water and river sediment using partial least squares regression.  Fresenius’ journal of analytical chemistry, 361(2):149–154, 1998.
  * [16] Juan A Gil and Rosario Romera.  On robust partial least squares (pls) methods.  Journal of Chemometrics: A Journal of the Chemometrics Society, 12(6):365–378, 1998.
  * [17] Mia Hubert and K Vanden Branden.  Robust methods for partial least squares regression.  Journal of Chemometrics: A Journal of the Chemometrics Society, 17(10):537–549, 2003.
  * [18] Javier González, Daniel Peña, and Rosario Romera.  A robust partial least squares regression method with applications.  Journal of Chemometrics: A Journal of the Chemometrics Society, 23(2):78–90, 2009.
  * [19] IN Wakelinc and HJH Macfie.  A robust pls procedure.  Journal of Chemometrics, 6(4):189–198, 1992.
  * [20] Irene Hoffmann, Peter Filzmoser, Sven Serneels, and Kurt Varmuza.  Sparse and robust pls for binary classification.  Journal of Chemometrics, 30(4):153–162, 2016.
  * [21] Ricardo A Maronna and Ruben H Zamar.  Robust estimates of location and dispersion for high-dimensional datasets.  Technometrics, 44(4):307–317, 2002.
  * [22] Peter J Rousseeuw and Katrien Van Driessen.  A fast algorithm for the minimum covariance determinant estimator.  Technometrics, 41(3):212–223, 1999.
  * [23] Kris Boudt, Peter J Rousseeuw, Steven Vanduffel, and Tim Verdonck.  The minimum regularized covariance determinant estimator.  Statistics and Computing, pages 1–16, 2017.
  * [24] Fanhua Shang, Yuanyuan Liu, James Cheng, and Hong Cheng.  Recovering low-rank and sparse matrices via robust bilateral factorization.  In 2014 IEEE International Conference on Data Mining, pages 965–970. IEEE, 2014.
  * [25] Gisele Mendes, Helga G Aleme, and Paulo JS Barbeira.  Determination of octane numbers in gasoline by distillation curves and partial least squares regression.  Fuel, 97:131–136, 2012.
  * [26] Durmuş Özdemir.  Determination of octane number of gasoline using near infrared spectroscopy and genetic multivariate calibration methods.  Petroleum science and technology, 23(9-10):1139–1152, 2005.


# Centrality bin size dependence of multiplicity correlation
in central Au+Au collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV

Yu-Liang Yan1, Dai-Mei Zhou2, Bao-Guo Dong1,3, Xiao-Mei Li1, Hai-Liang Ma1,
Ben-Hao Sa1,2,4111Corresponding author<EMAIL_ADDRESS>1 China Institute of
Atomic Energy, P.O. Box 275(18), Beijing 102413, China
2 Institute of Particle Physics, Huazhong Normal University, Wuhan 430079,
China
3 Center of Theoretical Nuclear Physics, National Laboratory of Heavy Ion
Collisions, Lanzhou 730000,China
4 CCAST (World Laboratory), P. O. Box 8730 Beijing 100080, China

###### Abstract

We have studied the centrality bin size dependence of charged particle
forward-backward multiplicity correlation strength in 5%, 0-5%, and 0-10% most
central Au+Au collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV with a parton and
hadron cascade model, PACIAE based on PYTHIA. The real (total), statistical,
and NBD (Negative Binomial Distribution) correlation strengths are calculated
by the real events, the mixed events, and fitting the charged particle
multiplicity distribution to the NBD, respectively. It is turned out that the
correlation strength increases with increasing centrality bin size
monotonously. If the discrepancy between real (total) and statistical
correlation strengths is identified as dynamical one, the dynamical
correlation may just be a few percent of the total (real) correlation.

###### pacs:

24.10.Lx, 24.60.Ky, 25.75.Gz

##  I INTRODUCTION

The study of fluctuations and correlations has been suggested as a useful
means for revealing the mechanism of particle production and Quark-Gluon-
Plasma (QGP) formation in Relativistic Heavy Ion Collisions hwa2 ; naya .
Correlations and fluctuations of the thermodynamic quantities and/or the
produced particle distributions may be significantly altered when the system
undergoes phase transition from hadronic matter to quark-gluon matter because
the degrees of freedom in two matters is very different.

The experimental study of fluctuations and correlations becomes a hot topic in
relativistic heavy ion collisions with the availability of high multiplicity
event-by-event measurements at the CERN-SPS and BNL-RHIC experiments. An
abundant experimental data have been reported star2 ; phen2 ; phob  where a
lot of new physics arise and are urgent to be studied. A lot of theoretical
investigations have been reported as well paja ; hwa1 ; konc ; brog ; fu ; yan
; bzda .

Recently STAR collaboration have measured the charged particle forward-
backward multiplicity correlation strength $b$ in Au+Au collisions at
$\sqrt{s_{\rm{NN}}}$=200 GeV star3 ; star4 . The outstanding features of STAR
data are:

  * •

In most central collisions, the correlation strength $b$ is approximately flat
across a wide range in $\Delta\eta$ which is the distance between the centers
of forward and backward (pseudo)rapidity bins.

  * •

This trend disappears slowly with decreasing centrality and approaches a
exponential function of $\Delta\eta$ at the peripheral collisions.

That has stimulated a lot of theoretical interests konc ; brog ; fu ; bzda .

In Ref. fu , a statistical model was proposed to calculate the charged
particle forward-backward multiplicity correlation strength $b$ in 0-10% most
central Au+Au collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV. One outstanding
feature of STAR data, the $b$ as a function of $\Delta\eta$ is approximately
flat, was well reproduced. The calculated value of $b\approx 0.44$ was
compared with STAR data of $b\approx 0.60$ star3 ; star4 .

However, in this statistical model fu  the Negative Binomial Distribution
(NBD) is assumed for the charged multiplicity distribution and the NBD
parameters of $\mu$ and $k$ (see later) are extracted from fit in with PHENIX
charged particle multiplicity distribution phen3 . It is turned out in Ref.
yan  that the experimental $\eta$ and $p_{T}$ acceptances have large
influences on the correlation strength $b$. The STAR experimental acceptances
are quite different from PHENIX, thus the inconsistency, using PHENIX
multiplicity data to explain STAR correlation data, involved in fu  have to be
studies further. Meanwhile, what is the discrepancy between $b\approx 0.60$
(STAR datum) and $b\approx 0.44$ (NBD) also needs to be answered.

In this paper we use a parton and hadron cascade model PACIAE sa , to
investigate the centrality bin size dependence of charged particle
multiplicity correlation in 5, 0-5, and 0-10% most central Au+Au collisions at
$\sqrt{s_{\rm{NN}}}$=200 GeV. Following Ref. yan  we generate the real events
(6000) by the PACIAE model, construct the mixed events according to real
events one by one, and extract the NBD parameters ($\mu$ and $k$) from fitting
the real events charged particle multiplicity distribution to the NBD. Then
the charged particle forward-backward multiplicity correlation strength $b$ is
calculated for the real events (real correlation strength), the mixed events
(statistical correlation strength), and the NBD (NBD correlation strength),
respectively. They are all nearly flat across a widw range in $\Delta\eta$.
Their magnitude in 0-10% most central Au+Au collisions are about 0.63, 0.59,
and 0.52, respectively. So the corresponding STAR data is well reproduced. It
is turned out that the real (total), statistical, and NBD correlation
strengths increase with increasing centrality bin size monotonously. If the
discrepancy between real (total) and statistical correlation strengths is
identified as dynamical one, then the dynamical correlation strength may just
be a few percent of the total (real) correlation strength.

##  II PACIAE MODEL

The parton and hadron cascade model, PACIAE sa , is based on PYTHIA soj2
which is a model for hadron-hadron ($hh$) collisions. The PACIAE model is
composed of four stages: parton initialization, parton evolution
(rescattering), hadronization, and hadron evolution (rescattering).

###  II.1 Parton initialization

In the PACIAE model a nucleon-nucleon ($NN$) collision is described with
PYTHIA model, where a $NN$ ($hh$) collision is decomposed into the parton-
parton collisions. The hard parton-parton collision is described by the
lowest-leading-order (LO) pQCD parton-parton cross section comb  with
modification of parton distribution function in the nucleon. And the soft
parton-parton interaction is considered empirically. The semihard, between
hard and soft, QCD $2\rightarrow 2$ processes are also involved in PYTHIA
(PACIAE) model. Because of the initial- and final-state QCD radiation added to
the above processes, the PYTHIA (PACIAE) model generates a multijet event for
a $NN$ ($hh$) collision. That is followed, in the PYTHIA model, by the string-
based fragmentation scheme (Lund model and/or Independent Fragmentation
model), thus a hadronic state is reached for a $NN$ ($hh$) collision. However,
in the PACIAE model above string fragmentation is switched off temporarily, so
the result is a multijet event (composed of quark pairs, diquark pairs and
gluons) instead of a hadronic state. If the diquarks (anti-diquarks) are split
forcibly into quarks (anti-quarks) randomly, the consequence of a $NN$ ($hh$)
collision is its initial partonic state composed of quarks, anti-quarks, and
gluons.

A nucleus-nucleus collision, in the PACIAE model, is decomposed into the
nucleon-nucleon collisions based on the collision geometry. A nucleon in the
colliding nucleus is randomly distributed in the spatial coordinate space
according to the Woods-Saxon distribution ($r$) and the 4$\pi$ uniform
distribution ($\theta$ and $\phi$). The beam momentum is given to $p_{z}$ and
$p_{x}=p_{y}=0$ is assumed for each nucleon in the colliding nucleus. A
closest approaching distance of two assumed straight line trajectories is
calculated for each $NN$ pair. If this distance is less than or equal to
$\displaystyle{\sqrt{\sigma_{\rm{tot}}/\pi}}$, then it is considered as a
collision pair. Here $\sigma_{\rm{tot}}$ refers to the total cross section of
$NN$ collision assumed to be 40 mb. The corresponding collision time of this
collision pair is then calculated. So the particle list and the $NN$ collision
(time) list can be constructed. A $NN$ collision pair with smallest collision
time is selected from the $NN$ collision (time ) list and performed by the
method in former paragraph. After upgrading the particle list and collision
(time) list we select and perform a new $NN$ collision pair again. Repeat
these processes until the collision (time) list is empty we obtain a initial
partonic state for a nucleus-nucleus collision.

###  II.2 Parton evolution (rescattering)

The next step, in the PACIAE model, is parton evolution (partonic
rescattering). Here the $2\rightarrow 2$ LO-pQCD differential cross sections
comb  are employed. The differential cross section for a subprocess
$ij\rightarrow kl$ reads

$\frac{d\sigma_{ij\rightarrow
kl}}{d\hat{t}}=K\frac{\pi\alpha_{s}^{2}}{\hat{s}}\sum_{ij\rightarrow kl},$ (1)

where the $K$ factor is introduced for higher order corrections and the non-
perturbative QCD correction as usual. Take the process $q_{1}q_{2}\rightarrow
q_{1}q_{2}$ as an example, one has

$\sum_{q_{1}q_{2}\rightarrow
q_{1}q_{2}}=\frac{4}{9}\frac{\hat{s}^{2}+\hat{u}^{2}}{\hat{t}^{2}},$ (2)

which can be regularized as

$\sum_{q_{1}q_{2}\rightarrow
q_{1}q_{2}}=\frac{4}{9}\frac{\hat{s}^{2}+\hat{u}^{2}}{(\hat{t}-m^{2})^{2}}$
(3)

by introducing the parton colour screen mass, $m$=0.63 GeV. In above equation
$\hat{s}$, $\hat{t}$, and $\hat{u}$ are the Mandelstam variables and
$\alpha_{s}$= 0.47 stands for the running coupling constant. The total cross
section of the parton collision $i+j$ is

$\sigma_{ij}(\hat{s})=\sum_{k,l}\int_{-\hat{s}}^{0}d\hat{t}\frac{d\sigma_{ij\to
kl}}{d\hat{t}}.$ (4)

With these total and differential cross sections the parton evolution
(rescattering) can be simulated by the Monte Carlo method until the parton-
parton collision is ceased (partonic freeze-out).

###  II.3 Hadronization

In the PACIAE model the partons can be hadronized with the string-based
fragmentation scheme or by the coalescence (recombination) models biro ; csiz
; grec ; hwo ; frie . The Lund string fragmentation regime, involved in the
PYTHIA model, is adopted for hadronization in this paper, see soj2  for the
details.

Meanwhile, we have proposed a simulant coalescence (recombination) model which
can be briefly explained as follows:

  1. 1.

The Field-Feynman parton generation mechanism ff1  is first applied to
deexcite the energetic parton and thus to increase the parton multiplicity.
This deexcitation of an energetic parton plays a similar role as string
multiple fragmentation in the Lund model and1 .

  2. 2.

The gluons are forcibly split into $q\bar{q}$ pair randomly.

  3. 3.

In the program there is a hadron table composed of mesons and baryons. The
pseudoscalar and vector mesons made of u, d, s, and c quarks, as well as
$B^{+}$, $B^{0}$, $B^{*0}$, and $\Upsilon$ are considered. The SU(4)
multiplets of baryons made of u, d, s, and c quarks (except those with double
c quarks) as well as $\Lambda^{0}_{b}$ are considered.

  4. 4.

Two partons can coalesce into a meson and three partons into a baryon
(antibaryon) according to the flavor, momentum, and spatial coordinates of
partons and the valence quark structure of hadron.

  5. 5.

When the coalescing partons can form either a pseudoscalar meson or a vector
meson (e. g. $u\bar{d}$ can form either a $\pi^{+}$ or a $\rho^{+}$) a
judgment of less discrepancy between the invariant mass of coalescing partons
and the mass of coalesced hadron is invoked to select one from two mesons
above. In the case of baryon, e. g. both $p$ and $\Delta^{+}$ are composed of
$uud$, the same judgment is invoked to select one baryon from both of
$\frac{1}{2}^{+}$ and $\frac{3}{2}^{+}$ baryons.

  6. 6.

Four momentum conservation is required.

  7. 7.

There is a phase space condition

$\frac{16\pi^{2}}{9}\Delta r^{3}\Delta p^{3}=\frac{h^{3}}{d},$ (5)

where $h^{3}/d$ is the volume occupied by a single hadron in the phase space,
$d$=4 refers to the spin and parity degeneracies, $\Delta r$ and $\Delta p$
stand for the spatial and momentum distances between coalescing partons,
respectively.

###  II.4 Hadron evolution (rescattering)

We obtain a configuration of hadrons in spatial and momentum coordinate spaces
for a nucleus-nucleus collision after the hadronization. If one only considers
the rescattering among
$\pi,k,p,n,\rho(\omega),\Delta,\Lambda,\Sigma,\Xi,\Omega,J/\Psi$ and their
antiparticles, the particle list is then constructed by the above hadrons. A
closest approaching distance of two assumed straight line trajectories is
calculated for each $hh$ pair. If this distance is less than or equal to
$\displaystyle{\sqrt{\sigma_{\rm{tot}}^{hh}/\pi}}$ sa1 , then it is considered
as a collision pair. Here $\sigma_{\rm{tot}}^{hh}$ refers to the total cross
section of $hh$ collision. The corresponding collision time of this collision
pair is then calculated. So the $hh$ collision (time) list can be constructed.
A $hh$ collision pair with smallest collision time is selected from the
collision (time) list and performed by the usual two-body collision method sa1
. After upgrading the particle list and collision (time) list we select and
perform a new $hh$ collision pair again. Repeat these processes until the
collision (time) list is empty (hadronic freeze-out).

A isospin averaged parametrization formula is used for the $hh$ cross section
koch ; bald . However, we also provide a option of constant total, elastic,
and inelastic cross sections sa1 : $\sigma_{\rm{tot}}^{NN}=40$ mb,
$\sigma_{\rm{tot}}^{\pi N}=25$ mb, $\sigma_{\rm{tot}}^{kN}=35$ mb,
$\sigma_{\rm{tot}}^{\pi\pi}=10$ mb, and the assumed ratio of inelastic to
total cross section equals 0.85. We also assume

$\sigma_{pp}=\sigma_{pn}=\sigma_{nn}=\sigma{\Delta N}=\sigma{\Delta\Delta}.$
(6)

The cross section of $\pi\bar{N}$ and $k\bar{N}$, for instance, is assumed to
be equal to the cross section of $\pi N$ and $kN$, respectively.

The momentum of scattered particles in a $hh$ elastic collision is simulated
according to that the $hh$ differential cross section,
$d\sigma_{\rm{tot}}^{hh}/dt$, is assumed to be an exponential function of $t$
which is squared momentum transfer sa1 . As it is impossible to include all
inelastic channels, we consider only a part of them ($\approx$ 600) which have
noticeable effects on the hadronic final state, and the rest is attributed to
the elastic scattering. Take incident channel $\pi N$ as an example, if there
are possible final channels of $\pi N\rightarrow\pi\Delta$, $\pi
N\rightarrow\rho N$, and $\pi N\rightarrow k\Lambda$, their relative
probabilities are then used to select one among above three channels. The
momentum of scattered particles in a $hh$ inelastic collision is simulated
according to the usual two-body kinematics sa1 ; pdg .

##  III CALCULATION AND RESULT

Following cape  the charged particle forward-backward multiplicity correlation
strength $b$ is defined as

$\ b=\frac{\langle n_{f}n_{b}\rangle-\langle n_{f}\rangle\langle
n_{b}\rangle}{\langle n_{f}^{2}\rangle-\langle
n_{f}\rangle^{2}}=\frac{{\rm{cov}}(n_{f},n_{b})}{{\rm{var}}(n_{f})},$ (7)

where $n_{f}$ and $n_{b}$ are, respectively, the number of charged particles
in forward and backward pseudorapidity bins defined relatively and
symmetrically to a given pseudorapidity $\eta$. $\langle n_{f}\rangle$ refers
to the mean value of $n_{f}$ for instance. cov($n_{f}$,$n_{b}$) and
var($n_{f}$) are the forward-backward multiplicity covariance and forward
multiplicity variance, respectively.

Table 1: Total charged particle multiplicity in three $\eta$ fiducial ranges in 0-6% most central Au+Au collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV. | $N_{\rm{ch}}(|\eta|<4.7)$ | $N_{\rm{ch}}(|\eta|<5.4)$ |  $N_{\rm{ch}}$(total)
---|---|---|---
PHOBOSa | 4810 $\pm$ 240 | 4960 $\pm$ 250 | 5060 $\pm$ 250
PACIAE | 4819 | 4983 | 5100
a The experimental data are taken from phob2 .

In the calculations the default values given in the PYTHIA model are adopted
for all model parameters except the parameters $K$ and $b_{s}$ (in the Lund
string fragmentation function). The $K$=3 is assumed and the $b_{s}$=6 is
fixed by fitting the charged particle multiplicity to the corresponding PHOBOS
data in 0-6% most central Au+Au collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV
phob2  as shown in Tab. 1. Therefore in the generation of real events there is
no free parameters. The mapping relation sa2  between the centrality
definition in theory and experiment

$b_{i}=\sqrt{g}b_{i}^{\rm{max}},\qquad b_{i}^{\rm{max}}=R_{A}+R_{B},$ (8)

is employed. In the above equation $b_{i}$ (in fm) refers to the theoretical
impact parameter and $g$ stands for the percentage of geometrical (total)
cross section used in experiment to define the centrality.
$R_{A}=1.12A^{1/3}+0.45$ fm, for instance, is the radius of nucleus $A$. Thus
the 0-10, 0-6, 0-5, and 5% most central collisions, for instance, are mapped
to $0<b_{i}<4.46$, $0<b_{i}<3.53$, $0<b_{i}<3.20$, and $b_{i}=3.20$ fm,
respectively.

In this paper we propose a mixed event method where the mixed events are
generated according to real events one by one. We first assume the charged
particle multiplicity $n$ in a mixed event is the same as one corresponding
real event. However, $n$ particles of this mixed event are sampled randomly
from the particle reservoir composed of all particles in all real events.
Therefore, there is no dynamical relevance among the particles in a mixed
event. So the correlation calculated by mixed events is reasonably to be
identified as the statistical correlation yan .

It is known that the statistical correlation can also be studied by the NBD
method, because the charged particle multiplicity distribution in high energy
heavy-ion collisions is close to NBD phen3 . For an integer $n$ the NBD reads

$\ P(n;\mu,k)=\begin{pmatrix}n+k-1\\\
k-1\end{pmatrix}\frac{(\mu/k)^{n}}{(1+\mu/k)^{n+k}},$ (9)

where $\mu\equiv\langle n\rangle$ is a parameter, $P(n;\mu,k)$ is normalized
in $0\leq n\leq\infty$, and $k$ is another parameter responsible for the shape
of the distribution. As proved in yan  the correlation strength can be
expressed as

$\ b=\frac{\langle n_{f}\rangle}{\langle n_{f}\rangle+k},$ (10)

where the parameter $k$ is fixed by fitting the charged particle multiplicity
to the NBD usually.

Figure 1: Charged particle pseudorapidity distribution in Au+Au collision at
$\sqrt{s_{\rm{NN}}}$=200 GeV: (a) 0-6% most central collision and (b) 0-10,
0-5, and 5% most central collision. The experimental data are taken from phob2
.

We compare the theoretical charged particle pseudorapidity distribution (open
circles) in 0-6% most central Au+Au collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV
with the corresponding PHOBOS data (solid squares) phob2  in Fig. 1 (a). One
sees here that the PHOBOS data are well reproduced. In Fig. 1 (b), we compare
the charged particle pseudorapidity distributions in 0-5% (open circles) and
5% (open triangles) most central Au+Au collisions with the 0-10% one (open
squares). We see in Fig. 1 (b) that the pseudorapidity distribution in 5% most
central collision is quite close to the 0-10% one, because the 5% centrality
is nearly equal to the average centrality of 0-10% centrality bin.

Figure 2: Charged particle forward-backward multiplicity correlation strength
$b$ in 0-10, 0-5, and 5% most central Au+Au collisions at
$\sqrt{s_{\rm{NN}}}$=200 GeV. The experimental data are taken from star3 .

In Fig. 2 we compare the calculated real (total) correlation strength $b$
(open squares) as a function of $\Delta\eta$ in 0-10% most central Au+Au
collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV with the corresponding STAR data
(solid squares) star3 . The STAR data feature of correlation strength $b$ is
approximately flat across a wide range in $\Delta\eta$ are well reproduced.
For comparison we also give the real (total) correlation strength in 0-5 and
5% most central collisions by open circles and triangles, respectively. One
sees here that the real (total) correlation strength decreases with decreasing
centrality bin size monotonously, because the charged particle multiplicity
fluctuation decreases from 0-10 to 0-5 and to 5% monotonously, as one will see
in Fig. 3. This first result of the correlation strength increases with
increasing centrality bin size monotonously given in the transport model
remains to be proved experimentally.

Figure 3: Charged particle multiplicity distributions in 0-10 (open squares),
0-5 (open circles), and 5% (open triangles) most central Au+Au collisions at
$\sqrt{s_{\rm{NN}}}$=200 GeV. The dotted, dashed, and solid lines are the
corresponding NBD fits, respectively.

The calculated charged particle multiplicity distributions in 0-10, 0-5, and
5% most central Au+Au collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV are given in
Fig. 3, respectively, by the open squares, circles and triangles. The
corresponding NBD fits are shown by dotted, dashed, and solid lines,
respectively. One sees in Fig. 3 that the charged particle multiplicity
fluctuation is increased and the NBD fit is worsened with increasing
centrality bin size monotonously.

Figure 4: The calculated charged particle total (real), statistical, and NBD
correlation strengths in 0-10, 0-5 and, 5% most central Au+Au collisions at
$\sqrt{s_{\rm{NN}}}$=200 GeV.

In Fig. 4 we compare the calculated charged particle real (solid symbols),
statistical (open symbols), and NBD (lines) correlation strengths as a
function of $\Delta\eta$ in 0-10, 0-5, and 5% most central Au+Au collisions at
$\sqrt{s_{\rm{NN}}}$=200 GeV. The solid squares, open squares, and dotted line
are for 0-10% most central collisions, solid circles, open circles, and dashed
line for 0-5%, and solid triangles, open triangles, and solid line for 5%,
respectively. We see in Fig. 4 that the behavior of correlation strength
increases with increasing centrality bin size monotonously is not only existed
in the real correlation strength but also in the statistical and NBD ones.

If the discrepancy between real (total) and statistical correlation strengths
is identified as the dynamical correlation strength, one then sees in Fig. 4
that the dynamical correlation strength may just be a few percent of the total
(real) correlation strength. The dynamical correlation strength in 0-10% most
central collision is close to the one in 5% most central collision globally
speaking. That is because the later centrality is nearly the average of the
former one. The dynamical correlation strength in 0-10% most central
collisions is globally less than 0-5% most central collision. That is because
the interactions (represented by the collision number for instance) in the
former collisions is weaker than the later one. We also see in Fig. 4 that the
statistical correlation strength is nearly the same as the NBD one in the 5%
most central collision, that is consistent with the results in $p+p$
collisions at the same energy yan . However the discrepancy between
statistical and NBD correlation strengths seems to be increased with
increasing centrality bin size monotonously. That is mainly because the NBD
fitting to the charged particle multiplicity distribution becomes worse with
increasing centrality bin size monotonously.

##  IV CONCLUSION

In summary, we have used a parton and hadron cascade model, PACIAE, to study
the centrality bin size dependence of charged particle forward-backward
multiplicity correlation strength in 5, 0-5, and 0-10% most central Au+Au
collisions at $\sqrt{s_{\rm{NN}}}$=200 GeV. The real (total), statistical, and
NBD correlation strengths are calculated by real events, mixed events, and NBD
method, respectively. The corresponding STAR data feature of the correlation
strength $b$ is approximately flat across a wide range in $\Delta\eta$ in most
central Au+Au collisions is well reproduced. It is turned out that the
correlation strength increases with increasing centrality bin size
monotonously. This first result, given in the transport model, remains to be
proved experimentally. If the discrepancy between real (total) and statistical
correlation strengths is identified as dynamical one yan , then the dynamical
correlation may be just a few percent of the total (real) correlation. As a
next step, we will investigate the relation between correlation strength $b$
and the centrality bin size in the mid-central and peripheral collisions, and
the STAR data feature of $b$ approaches an exponential function of
$\Delta\eta$ at the peripheral collisions.

ACKNOWLEDGMENT

Finally, the financial support from NSFC (10635020, 10605040, and 10705012) in
China is acknowledged

## References

  * (1) R. C. Hwa, Int. J. Mod. Phys. E 16, 3395 (2008).
  * (2) T. K. Nayak, J. of Phys. G 32, S187 (2006).
  * (3) J. Adams, et al., STAR Collaboration, Phys. Rev. C 75, 034901 (2007).
  * (4) A. Adare, et al., PHENIX Collaboration, Phys. Rev. Lett. 98, 232302 (2007).
  * (5) Zheng-Wei Chai , et al., PHOBOS Collaboration, J. of Phys.: Conference Series 27, 128 (2005).
  * (6) N. S. Amelin, N. Armesto, M. A. Braun, E. G. Ferreiro, and C. Pajares, Phys. Rev. Lett. 73, 2813 (1994); N. Armesto, M. A. Braun, and C. Pajares, Phys. Rev. C 75, 054902 (2007).
  * (7) R. C. Hwa and C. B. Yang , nucl-th/0705.3073.
  * (8) V. P. Konchakovski, M. I. Gorenstein, and E. L. Bratkovskaya, Phys. Rev. C 76, 031901(R) (2007); V. P. Konchakovski, M. Hauer, G. Torrieri, M. I. Gorenstein, and E. L. Bratkovskaya, nucl-th/0812.3967 .
  * (9) P.Brogueira, J. Dias de Deus, and J. G. Milhano, Phys. Rev. C 76, 064901 (2007).
  * (10) Jinghua Fu, Phys. Rev. C 77, 027902 (2008).
  * (11) Yu-Liang Yan, Bao-Guo Dong, Dai-Mei Zhou, Xiao-Mei Li, and Ben-Hao Sa, Phys. Lett. B 660, 478 (2008).
  * (12) A. Bzdak, hep-ph/0902.2639.
  * (13) T. Tarnowsky, STAR Collaboration, arXiv:0711.1175v1 (PoS CP0D07 (2007) 019); Int. J. Mod. Phys. E 16, 3363 (2008).
  * (14) B. K. Srivastavs, STAR Collaboration, Int. J. Mod. Phys. E 16, 3371 (2008).
  * (15) S. S. Adler et al., PHENIX Collaboration, Phys. Rev. C, 76, 034903 (2007).
  * (16) Dai-Mei Zhou, Xiao-Mei Li, Bao-Guo Dong, and Ben-Hao Sa, Phys. Lett. B 638, 461 (2006); Ben-Hao Sa, Xiao-Mei Li, Shou-Yang Hu, Shou-Ping Li, Jing Feng, and Dai-Mei Zhou, Phys. Rev. C 75, 054912 (2007).
  * (17) T. Sj$\ddot{o}$strand, Comput. Phys. Commun. 82, 74 (1994).
  * (18) B. L. Combridge, J. Kripfgang, and J. Ranft, Phys. Lett. B 70, 234 (1977).
  * (19) T. S. Biró, P. Lévai, and J. Zimányi, Phys. Rev. C 59, 1547 (1999).
  * (20) P. Csizmadia and P. Lévai, Phys. Rev. C 61, 031903(R) (2000).
  * (21) V. Greco, C. M. Ko, and P. Lévai, Phys. Rev. Lett. 90, 202302 (2003).
  * (22) R. C. Hwa and C. B. Yang, Phys. Rev. C 67, 034902 (2003).
  * (23) R. J. Fries, B. Müller, and C.Nonaka, Phys. Rev. Lett. 90, 202303 (2003).
  * (24) R. D. Field and R. P. Feynman, Phys. Rev. D 15, 2590 (1977); Nucl. Phys. B 138, 1(1978); R. P. Feynman, R. D. Field, and G. C. Fox, Phys. Rev. D 18, 3320 (1978).
  * (25) B. Andersson, G. Gustafson, G. Ingelman, T. Sj$\ddot{o}$strand, Phys. Rep. 97, 33 (1983); B. Andersson, G. Gustafson, B. S$\ddot{o}$derberg, Nucl. Phys. B 264, 29 (1986).
  * (26) Ben-Hao Sa and Tai An, Comput. Phys. Commun. 90, 121 (1995); Tai An and Ben-Hao Sa, Comput. Phys. Commun. 116, 353 (1999).
  * (27) P. Koch, B. Müller, and J. Rafelski Phys. Rep. 142, 167 (1986).
  * (28) A. Baldini, et al., “Total cross sections for reactions of high energy particles”, Springer-Verlag, Berlin, 1988.
  * (29) PDG, “Particle Physics Booklet”, Extracted from C. Amsler, et al., Phys. Lett. B 667, 1 (2008).
  * (30) A. Capella, U. Sukhatme, C.-I. Tan, and J. Tran Thanh Van, Phys. Rep. 236, 225 (1994).
  * (31) B. B. Back, et al., PHOBOS Collaboration, Phys. Rev. Lett. 91, 052303 (2003).
  * (32) Ben-Hao Sa, A. Bonasera, An Tai and Dai-Mei Zhou, Phys. Lett. B 537, 268 (2002).