\section{Introduction}
Due to the recent advances in machine learning (ML), ML models are being more and more used in practice and applied to real-world scenarios for decision making~\cite{Goel2016PrecinctOP,creditscoresunfair,creditriskml}.  Essential demands for user acceptance as well as  legal regulations like the EU's "General Data Protection Right" (GDPR)~\cite{gdpr}, that contains a "right to an explanation", make it indispensable to explain the output and behavior of ML models in a comprehensible way.

As a consequence, many research approaches focused on the question how to realize  explainability and transparency in machine learning in recent years~\cite{explainingexplanations, explainingblackboxmodelssurvey, surveyxai, explainableartificialintelligence}. There exist diverse methods for explaining  ML models~\cite{explainingblackboxmodelssurvey,molnar2019}.
One family of methods are model-agnostic methods~\cite{explainingblackboxmodelssurvey,modelagnosticinterpretability}. Model-agnostic methods are flexible in the sense that they are not tailored to a particular model or representation. This makes model-agnostic methods (in theory) applicable to many different types of ML models. In particular, "truly" model-agnostic methods do not need access to the training data or model internals. It is sufficient to have an interface for passing data points to the model and obtaining the output/predictions of the model - the underlying model is viewed as a black-box.

Examples of model-agnostic methods are feature interaction methods~\cite{featureinteraction}, feature importance methods~\cite{featureimportance}, partial dependency plots~\cite{partialdependenceplots} and local methods that approximates the model locally by an explainable model (e.g. a decision tree)~\cite{lime2016,decisiontreecounterfactual}. These methods explain the models by using features as vocabulary.

A different class of model-agnostic explanations are example-based explanations where a prediction or behavior is explained by a (set of) data points~\cite{casebasedreasoning}. Instances of example-based explanations are prototypes \& criticisms~\cite{prototypescriticism} and influential instances~\cite{influentialinstances}. Another instance of example-based explanations are counterfactual explanations~\cite{counterfactualwachter}. A counterfactual explanation is a change of the original input that leads to a different (specific) prediction/behavior of the ML model - \textit{what has to be different in order to change the prediction of the model?} Such an explanation is considered to be fairly intuitive, human-friendly and useful because it tells people what to do in order to achieve a desired outcome~\cite{counterfactualwachter,molnar2019}.

Existing counterfactual explanations are mostly model agnostic methods, which are universally applicable but computationally expensive\cite{decisiontreecounterfactual,counterfactualsgeneticalgorithm,growingsphere}. In this work, we address the question how to efficiently compute counterfactuals for a popular class of models by referring to its specific structure.

Prototype based models such as learning vector quantization (LVQ) represent data by a set of representative samples~\cite{lvqreview}. LVQ models can be combined with metric learning and thereby increase the effectiveness of the model in case of few prototypes~\cite{gmlvq,mrslvq}.
Furthermore, LVQ models can be used in many settings like life-long learning~\cite{lifelonglvq}.

Here, we will consider the question how to efficiently compute counterfactual explanations of prototype-based classifiers, in particular LVQ models.
By exploiting the special structure of such models, we are able to
\begin{enumerate}
\item propose model- and regularization-dependent methods for efficiently computing counterfactual explanations of LVQ models,
\item and empirically demonstrate the efficiency of the modeling as regards speed-up as well as required amount of change in comparison to standard techniques.
\end{enumerate}
Further, the framework enables a straightforward incorporation of domain knowledge, which can be phrased as additional constraints. Such domain knowledge could be used for generating more plausible and feasible counterfactual explanations.\\\\
The remainder of this paper is structured as follows: First, we briefly review counterfactual explanations (section~\ref{sec:counterfactualexplanations}) and learning vector quantization models (section~\ref{sec:lvq}). Then, in section~\ref{sec:cf_lvq} we introduce our convex and non-convex programming framework for efficiently computing counterfactual explanations of different types of LVQ models - note that all derivations can be found in the appendix (section~\ref{sec:appendix}). Furthermore, we empirically evaluate the efficiency of our methods. Finally, section~\ref{sec:conclusion} summarizes the results of this paper.

\section{Counterfactual explanations}\label{sec:counterfactualexplanations}
Counterfactual explanations~\cite{counterfactualwachter} (often just called counterfactuals) are an instance of example-based explanations~\cite{casebasedreasoning}. Other instances of example-based explanations~\cite{molnar2019} are influential instances~\cite{influentialinstances} and prototypes \& criticisms~\cite{prototypescriticism}.

A counterfactual states a change to some features of a given input such that the resulting data point (called counterfactual) has a different (specified) prediction than the original input. The rational is considered to be intuitive, human-friendly and useful because it tells practitioners which minimum  changes can lead to a desired outcome~\cite{counterfactualwachter,molnar2019}.

A classical use case of counterfactual explanations is loan application~\cite{creditriskml,molnar2019}:
\textit{Imagine you applied for a credit at a bank. Unfortunately, the bank rejects your application. Now, you would like to know why. In particular, you would like to know what would have to be different so that your application would have been accepted. A possible explanation might be that you would have been accepted if you would earn 500\$ more per month and if you would not have a second credit card.}

Although counterfactuals constitute very intuitive explanation mechanisms, there do exist a couple of problems.

One problem is that there often exist more than one counterfactual - this is called \textit{Rashomon effect}~\cite{molnar2019}. If there are more than one possible explanation (counterfactual), it is not clear which one should be selected.

In our opinion, the main issue is feasibility and plausibility - in the sense that the counterfactual data point should be valid and plausible in the given domain. There exist first approaches like~\cite{counterfactualguidedbyprototypes, face} but in general it is still an open research question how to efficiently compute feasible and plausible counterfactuals.

An alternative - but very similar in the spirit - to counterfactuals~\cite{counterfactualwachter} is the Growing Spheres method~\cite{growingsphere}. However, this method suffers from the curse of dimensionality because it has to draw samples from the input space, which can become difficult if the input space is high-dimensional.

\begin{definition}[Counterfactual explanation]
Assume a prediction function $\ensuremath{h}$ is given. Computing a counterfactual $\ensuremath{\vec{x}'} \in \mathbb{R}^d$ for a given input $\ensuremath{\vec{x}} \in \mathbb{R}^d$ is phrased as an optimization problem~\cite{counterfactualwachter}:
\begin{equation}\label{eq:counterfactualoptproblem}
\underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d}{\arg\min}\; \loss\big(\ensuremath{h}(\ensuremath{\vec{x}'}), \ensuremath{y^{c}}\big) + C \cdot \regularization(\ensuremath{\vec{x}'}, \ensuremath{\vec{x}})
\end{equation}
where $\loss(\cdot)$ denotes the loss function,  $\ensuremath{y^{c}}$ the requested prediction, and  $\regularization(\cdot)$  a penalty term for deviations of $\ensuremath{\vec{x}'}$ from the original input $\ensuremath{\vec{x}}$. $C>0$ denotes the regularization strength.
\end{definition}

Two common regularizations are the weighted Manhattan distance and the generalized L2 distance.

The weighted Manhattan distance is defined as
\begin{equation}\label{eq:weighted_l1}
\regularization(\ensuremath{\vec{x}'}, \ensuremath{\vec{x}}) = \sum_j \alpha_j \cdot |(\ensuremath{\vec{x}})_j - (\ensuremath{\vec{x}'})_j|
\end{equation}
where $\alpha_j > 0$ denote the feature wise weights. A popular choice~\cite{counterfactualwachter} for $\alpha_j$ is the inverse median absolute deviation of the $j$-th feature median in the training data set $\set{D}$:
\begin{equation}
\begin{split}
& \alpha_j = \frac{1}{{\mad}_j} \\
& \text{where } {\mad}_j = \underset{\ensuremath{\vec{x}} \,\in\, \set{D}}{\median}\left(\left|(\ensuremath{\vec{x}})_j - \underset{\ensuremath{\vec{x}} \,\in\, \set{D}}{\median}\big((\ensuremath{\vec{x}})_j\big)\right|\right)
\end{split}
\end{equation}
The benefit of this choice is that it takes the (potentially) different variability of the features into account.  However, because we need access to the training data set $\set{D}$, this regularization is not a truly model-agnostic method - it is not usable if we only have access to a prediction interface of a black-box model.

The generalized L2 distance\footnote{also called Mahalanobis distance} is defined as
\begin{equation}\label{eq:general_l2}
\regularization(\ensuremath{\vec{x}'}, \ensuremath{\vec{x}}) = \pnorm{\ensuremath{\vec{x}} - \ensuremath{\vec{x}'}}_{\distmat}^2 = (\ensuremath{\vec{x}} - \ensuremath{\vec{x}'})^\top\distmat(\ensuremath{\vec{x}} - \ensuremath{\vec{x}'})
\end{equation}
where $\distmat$ denotes a symmetric positive semi-definite (s.psd) matrix. Note that the L2 distance can be recovered by setting $\distmat=\mat{\mathbb{I}}$. The generalized L2 distance can be interpreted as the Euclidean distance in a linearly distorted space.

Depending on the model and the choice of $\loss(\cdot)$ and $\regularization(\cdot)$, the final optimization problem might be differentiable or not. If it is differentiable, we can use a gradient-based optimization algorithm like conjugate gradients, gradient descent or (L-)BFGS. Otherwise, we have to use a black-box optimization algorithm like Downhill-Simplex (also called Nelder-Mead) or Powell. However, using a model and regularization specific optimization technique has the best chances to enable a simplification and computationally efficient solution. Unfortunately, there exit model specific optimizers for counterfactual explanations for few ML models only.

\section{Learning vector quantization}\label{sec:lvq}
In learning vector quantization (LVQ) models~\cite{lvqreview} we compute a set of labeled prototypes $\{(\vec{\prototype}_i, \ensuremath{o}_i)\}$ from a training data set of labeled real-valued vectors - we refer to the $i$-th prototype as $\vec{\prototype}_i$ and the corresponding label as $\ensuremath{o}_i$.
A new data point is classified according to the winner-takes-it-all scheme:
\begin{equation}\label{eq:lvq_predict}
\begin{split}
& \ensuremath{h}(\ensuremath{\vec{x}}) = \ensuremath{o}_i \quad \text{s.t. } \vec{\prototype}_i = \underset{\vec{\prototype}_j}{\arg\min}\;{\dist}(\ensuremath{\vec{x}}, \vec{\prototype}_j)
\end{split}
\end{equation}
where $\dist(\cdot)$ denotes a function for computing the distance between a data point and a prototype - in vanillas LVQ, this is independent chosen globally as the squared Euclidean distance:
\begin{equation}\label{eq:distfunction}
\dist(\ensuremath{\vec{x}}, \vec{\prototype}) = (\ensuremath{\vec{x}} - \vec{\prototype})^\top \mat{\mathbb{I}} (\ensuremath{\vec{x}} - \vec{\prototype})
\end{equation}
There exist matrix-LVQ models like GMLVQ, LGMLVQ~\cite{gmlvq}, MRSLVQ and LMRSLVQ~\cite{mrslvq} that learn a custom (class or prototype specific) distance matrix $\distmat_{\prototype}$ that is used instead of the identity $\mat{\mathbb{I}}$ when computing the distance between a data point and a prototype. This gives rise to the generalized L2 distance Eq.~\eqref{eq:general_l2}:
\begin{equation}\label{eq:lvq_generaldist}
\dist(\ensuremath{\vec{x}}, \vec{\prototype}) = (\ensuremath{\vec{x}} - \vec{\prototype})^\top{\distmat}_{\prototype}(\ensuremath{\vec{x}} - \vec{\prototype})
\end{equation}
Because ${\distmat}_{\prototype}$ must be a s.psd matrix, instead of learning ${\distmat}_{\prototype}$ directly, these LVQ variants learn a matrix $\distmatfactor_{\prototype}$ and compute the final distance matrix as:
\begin{equation}
{\distmat}_{\prototype}={\distmatfactor}_{\prototype}^\top{\distmatfactor}_{\prototype}
\end{equation}
By this, we can guarantee that the matrix $\distmat_{\prototype}$ is s.psd, whereas the model only has to learn an arbitrary matrix ${\distmatfactor}_{\prototype}$ - which is much easier than making some restriction on the definiteness. Training usually takes place by optimizing suitable cost functions as regards prototype locations and metric parameters. For counterfactual reasoning, the specific training method~\cite{gmlvq,mrslvq} is irrelevant and we refer to the final model only.

\section{Counterfactual explanations of LVQ models}\label{sec:cf_lvq}
We aim for an efficient explicit formulation how to find counterfactuals, for diverse  LVQ models.

\subsection{General approach}
Because a LVQ model assigns the label of the nearest prototype to a given input, we know that the nearest prototype of a counterfactual $\ensuremath{\vec{x}'}$ must be a prototype $\vec{\prototype}_i$ with $\ensuremath{o}_i=\ensuremath{y^{c}}$. In order to compute a counterfactual $\ensuremath{\vec{x}'}$ of a given input $\ensuremath{\vec{x}}$, it is sufficient to solve the following optimization problem for each prototype $\vec{\prototype}_i$ with $\ensuremath{o}_i=\ensuremath{y^{c}}$ and select the counterfactual $\ensuremath{\vec{x}'}$ yielding the smallest value of $\regularization(\ensuremath{\vec{x}'}, \ensuremath{\vec{x}})$:
\begin{subequations}
\label{eq:cflvq_general}
 \begin{align}
  &\underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d}{\arg\min}\;\regularization(\ensuremath{\vec{x}'}, \ensuremath{\vec{x}}) \label{eq:cflvq_general:regularization} \\
  & \quad \text{s.t. } \dist(\ensuremath{\vec{x}'}, \vec{\prototype}_i) + \epsilon \leq \dist(\ensuremath{\vec{x}'}, \vec{\prototype}_j) \quad \forall\, \vec{\prototype}_j\in\set{P}(\ensuremath{y^{c}}) \label{eq:cflvq_general:constraints}
 \end{align}
\end{subequations}
where $\set{P}(\ensuremath{y^{c}})$ denotes the set of all prototypes \emph{not} labeled as $\ensuremath{y^{c}}$ and $\epsilon>0$ is a small value preventing that the counterfactual lies exactly on the decision boundary. Note that the feasible region of Eq.~\eqref{eq:cflvq_general} is always non-empty - the prototype $\vec{\prototype}_i$ is always a feasible solution. Furthermore, in contrast to Eq.~\eqref{eq:counterfactualoptproblem}, the formulation does not include hyperparameters.

The pseudocode for computing a counterfactual of a LVQ model is described in Algorithm~\ref{algo:cflvq}.
\begin{algorithm}[!tb]
\caption{Computing a counterfactual of a LVQ model}\label{algo:cflvq}
\textbf{Input:} Original input $\ensuremath{\vec{x}}$, requested prediction $\ensuremath{y^{c}}$ of the counterfactual, the LVQ model\\
\textbf{Output:} Counterfactual $\ensuremath{\vec{x}'}$
\begin{algorithmic}[1]
 \State $\ensuremath{\vec{x}'} = \vec{0}$ \Comment{Initialize dummy solution}
 \State $z = \infty$
 \For{$\vec{\prototype}_i$ with $\ensuremath{o}_i=\ensuremath{y^{c}}$} \Comment{Try each prototype with the correct label}
 	\State Solving Eq.~\eqref{eq:cflvq_general} yields a counterfactual $\ensuremath{\vec{x}'}_{*}$
 	\If{$\regularization(\ensuremath{\vec{x}'}_{*}, \ensuremath{\vec{x}}) < z$} \Comment{Keep this counterfactual if it deviates less from the original input then the currently "best" counterfactual}
 		\State $z=\regularization(\ensuremath{\vec{x}'}_{*}, \ensuremath{\vec{x}})$
 		\State $\ensuremath{\vec{x}'} = \ensuremath{\vec{x}'}_{*}$
 	\EndIf
 \EndFor
\end{algorithmic}
\end{algorithm}
Note that the \texttt{for} loop in line 3 of Algorithm~\ref{algo:cflvq} can be easily parallelized.

For the weighted Manhattan distance as a regularization $\regularization(\cdot)$, the objective Eq.~\eqref{eq:cflvq_general:regularization} becomes linear in $\ensuremath{\vec{x}'}$, where
$\mat{\Upsilon}$ is the diagonal matrix with entries $\alpha_j$ and $\vec{\beta}$ is an auxiliary variable that can be discarded afterwards:
\begin{equation}\label{eq:objective:l1}
\begin{split}
& \underset{\ensuremath{\vec{x}'}, \vec{\beta} \,\in\, \mathbb{R}^d}{\min}\;\vec{1}^\top\vec{\beta} \\
& \text{s.t.} \;\mat{\Upsilon}\ensuremath{\vec{x}'} - \mat{\Upsilon}\ensuremath{\vec{x}} \leq \vec{\beta}, \quad -\mat{\Upsilon}\ensuremath{\vec{x}'} + \mat{\Upsilon}\ensuremath{\vec{x}} \leq \vec{\beta}, \quad\vec{\beta} \geq \vec{0}
\end{split}
\end{equation}
For the Euclidean distance Eq.~\eqref{eq:general_l2} as  regularization $\regularization(\cdot)$, the objective Eq.~\eqref{eq:cflvq_general:regularization} can be written in a convex quadratic form:
\begin{equation}\label{eq:objective:l2}
\underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d}{\min}\;\frac{1}{2}\ensuremath{\vec{x}'}^\top\ensuremath{\vec{x}'} - \ensuremath{\vec{x}'}^\top\ensuremath{\vec{x}}
\end{equation}
\mbox{}\\
In the subsequel  we explore Eq.~\eqref{eq:cflvq_general} for different regularizations $\regularization(\cdot)$ and LVQ models, and investigate how to solve it efficiently. Note that for the purpose of better readability and due to space constraints, we put all derivations in the appendix (section~\ref{sec:appendix}).

\subsection{(Generalized matrix) LVQ}
When using (generalized matrix) LVQ models - no class or prototype specific distance matrices - the optimization problem Eq.~\eqref{eq:cflvq_general} becomes a linear program (LP) when using the weighted Manhattan distance as a regularizer, and a convex quadratic program (QP) problem when using the Euclidean distance. Both programs can be solved efficiently and (up to equivalence) uniquely~\cite{Boyd2004}.
More precisely,
the constraints Eq.~\eqref{eq:cflvq_general:constraints} can be written as a set of linear inequality constraints:
\begin{equation}
\ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij} + \epsilon\leq 0 \quad \forall\, \vec{\prototype}_j\in\set{P}(\ensuremath{y^{c}})
\end{equation}
where
\begin{equation}
\vec{q}_{ij} = \distmat(\vec{\prototype}_{j} - \vec{\prototype}_i) \quad\quad r_{ij} =  \frac{1}{2}\left(\vec{\prototype}_{i}^\top\distmat\vec{\prototype}_{i} - \vec{\prototype}_{j}^\top\distmat\vec{\prototype}_{j}\right)
\end{equation}
We set $\distmat=\mat{\mathbb{I}}$ if the LVQ model uses the Euclidean distance.

Another benefit is that we can easily add additional convex quadratic constraints like box constraints, restricting the value range of linear interaction of features or freezing some features - these features are not allowed to be different in the counterfactual. We can use such additional constraints for computing a more plausible and feasible counterfactual.

\subsection{(Localized generalized matrix) LVQ}
In case of local matrix-LVQ models that learn a class or prototype specific distance matrix $\distmat_{\prototype}$, the optimization problem Eq.~\eqref{eq:cflvq_general} becomes a quadratically constrained quadratic program (QCQP) for Manhattan and Euclidean regularization, since
the constraints Eq.~\eqref{eq:cflvq_general:constraints} become a set of quadratic constraints:
\begin{equation}\label{eq:lgmlvq:constraints}
\frac{1}{2}\ensuremath{\vec{x}'}^\top\mat{Q}_{ij}\ensuremath{\vec{x}'} + \ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij} + \epsilon \leq 0 \quad \forall\, \vec{\prototype}_j\in\set{P}(\ensuremath{y^{c}})
\end{equation}
where
\begin{equation}
\mat{Q}_{ij} = {\distmat}_{i} - {\distmat}_{j}
\end{equation}
Unfortunately, we can not make any statement about the definiteness of $\mat{Q}_{ij}$. Because $\mat{Q}_{ij}$ is the difference of two s.psd matrices, all we know is that it is a symmetric matrix. Therefore, Eq.~\eqref{eq:lgmlvq:constraints} yields a non-convex QCQP and solving a non-convex QCQP is known to be NP-hard~\cite{park2017general}. However, there exist methods like the Suggest-Improve framework~\cite{park2017general} that can approximately solve a non-convex QCQP very efficiently - more details on how to apply this to Eq.~\eqref{eq:lgmlvq:constraints} can be found in the appendix (section~\ref{sec:appendix}).

\subsection{Experiments}
We empirically confirm the efficiency of our proposed methods in comparison to black-box mechanisms by means of the following experiments:
We use GLVQ, GMLVQ and LGMLVQ models with 3 prototypes per class for the "Breast Cancer Wisconsin (Diagnostic) Data Set"~\cite{breastcancer}, the "Optical Recognition of Handwritten Digits Data Set"~\cite{ocr} and the "Ames Housing dataset"~\cite{housepricedata}. Thereby, we use PCA-preprocessing
\footnote{This is for the purpose of better stability and better semantic meaning, since in the original domain already a small perturbation is sufficient for changing the class, since adversarial attacks exist even for linear functions in high dimensions if feature correlations are neglected. Since PCA can be approximately inverted, counterfactuals in PCA space can be lifted to the original data space.} to reduce the dimensionality of the digit data set to $10$ and of the breast cancer data set to $5$. We standardize the house data set and turned it into a binary classification problem\footnote{In addition, we select the following features: TotalBsmt, 1stFlr, 2ndFlr, GrLivA, WoodDeck, OpenP, 3SsnP, ScreenP and PoolA - When computing counterfactuals, we fix the last five features.} by setting the target to $1$ if the price is greater or equal to 160k\$ and $0$ otherwise.
We use the Suggest-Improve framework~\cite{park2017general} for  solving the non-convex QCQPs, where we pick the target prototype as initial  solution in the Suggest-step and we use the penalty convex-concave procedure (CCP) in the Improve-step~\cite{park2017general}.
For comparison, we use the optimizer for computing counterfactual explanations of LVQ models as implemented in ceml~\cite{ceml} - where the distance to the nearest prototype with the requested label $\ensuremath{y^{c}}$ is minimized by  Downhill-Simplex search or CMA-ES.

We report results for the Manhattan distance as regularizer - we used the Manhattan distance for enforcing a sparse solution.
For each possible combination of model, data set and method,  a 4-fold cross validation is conducted and the mean  distance is reported.
The results are listed in Table~\ref{table:experimentresults1}. In all cases, our method yields counterfactuals that are closer to the original data point than the one found by minimizing the original cost function Eq.~\eqref{eq:counterfactualoptproblem} with Downhill-Simplex search (DS) or CMA-ES.
In addition, our method is between $1.5$ and $158.0$ faster in comparison to  DS/CMA-ES method. Furthermore, Downhill-Simplex and CMA-ES did not always find a counterfactual when dealing with LGMLVQ models.
We would like to remark that our formulation can easily be extended by linear/quadratic constraints which can incorporate prior knowledge such as a maximum possible change of specific input features - see Table~\ref{table:pausiblecounterfactuals} for an example. Such extensions do not change the form of the optimization problem hence its complexity.

\begin{table}[t]
\centering
\footnotesize
\begin{tabular}{|c||c|c|c||c|c|c||c|c|c|}
 \hline
 \multicolumn{1}{|c||}{\textit{Data set}} & \multicolumn{3}{|c||}{Breast cancer} & \multicolumn{3}{|c||}{Handwritten digits} & \multicolumn{3}{|c|}{House prices}  \\
  \hline
 \textit{ Method} & DS & CMA & Ours & DS & CMA & Ours & DS & CMA & Ours \\
 \hline
 GLVQ   & 3.26 & 3.28 & \textbf{1.96} & 6.51 & 6.53 & \textbf{3.99} & 3.81 & 3.85 & \textbf{3.32} \\
 GMLVQ  & 2.71 & 6.49 & \textbf{2.46} & 21.34 & 11.63 & \textbf{4.40} & 5.06 & 8.63 & \textbf{3.78} \\
LGMLVQ & {\it 2.00}

 & {\it 1.61}

 & \textbf{1.57} & {\it 8.12}

 & {\it 7.88}

 & \textbf{7.53} & {\it 12.74}

 & {\it 12.59}

 & \textbf{8.20} \\
 \hline
\end{tabular}
\caption{Mean \emph{Manhattan distance} between the counterfactual and the original data point - best (smallest) values are \textbf{highlighted}.
For LGMLVQ with DS or CMA-ES (marked italic), in  5\% to 60\% of the cases no solution was found.}
\label{table:experimentresults1}
\end{table}
\begin{table}[htb]
\centering
\footnotesize
\begin{tabular}{|c||c|c|c|c|c|
 \hline
 \multicolumn{1}{|c||}{\textit{Data point}} & \multicolumn{1}{|c|}{TotalBsmt} & \multicolumn{1}{|c|}{1stFlr} & \multicolumn{1}{|c|}{2ndFlr} & \multicolumn{1}{|c|}{GrLivA}
 & \multicolumn{1}{|c|}{Label} \\
  \hline
 \textit{Original} & 0 & 1120 & 468 & 1588
 & 1 \\
 \hline
 \textit{Counterfactual} & 0 & 366 & 1824 & 2225
 & 0 \\
 \textit{Constrained Counterfactual}

 & 373 & 1454 & 1454 & 3125
 & 0 \\
 \hline
 \end{tabular}
 \caption{House prices, we obtain a "plausible" counterfactual by adding constrains, here the constraint "2ndFlr $\leq$ 1stFlr" is added.}
\label{table:pausiblecounterfactuals}
\end{table}
The implementation of our proposed method for computing counterfactual explanations is available online\footnote{\url{https://github.com/andreArtelt/efficient_computation_counterfactuals_lvq}}. All experiments were implemented in Python 3.6~\cite{van1995python} using the packages cvxpy~\cite{cvxpy} cvx-qcqp~\cite{qcqp}, sklearn-lvq~\cite{sklearn_lvq}, numpy~\cite{numpy}, scipy~\cite{scipy}, scikit-learn~\cite{scikit-learn} and ceml~\cite{ceml}.

\section{Conclusion}\label{sec:conclusion}
We proposed, and empirically evaluated, model- and regularization-dependent convex and non-convex programs for efficiently computing counterfactual explanations of LVQ models. We found that in many cases we get either a set of linear or convex quadratic programs which both can be solved efficiently. Only in the case of localized matrix-LVQ models we have to solve a set of non-convex quadratically constrained quadratic programs - we found that they can be efficiently approximately solved by using the Suggest-Improve framework.

\section{Appendix}\label{sec:appendix}
\subsection{Minimizing the Euclidean distance}
First, we expand the Euclidean distance Eq.~\eqref{eq:general_l2}:
\begin{equation}
\begin{split}
\pnorm{\ensuremath{\vec{x}'} - \ensuremath{\vec{x}}}_2^2 &= (\ensuremath{\vec{x}'} - \ensuremath{\vec{x}})^\top(\ensuremath{\vec{x}'} - \ensuremath{\vec{x}}) \\
&= \ensuremath{\vec{x}'}^\top\ensuremath{\vec{x}'} - \ensuremath{\vec{x}'}^\top\ensuremath{\vec{x}} - \ensuremath{\vec{x}}^\top\ensuremath{\vec{x}'} + \ensuremath{\vec{x}}^\top\ensuremath{\vec{x}} \\
&= \ensuremath{\vec{x}'}^\top\ensuremath{\vec{x}'} - 2\ensuremath{\vec{x}}^\top\ensuremath{\vec{x}'} + \ensuremath{\vec{x}}^\top\ensuremath{\vec{x}} \\
\end{split}
\end{equation}
Next, we note that that we can drop the constant $\ensuremath{\vec{x}}^\top\ensuremath{\vec{x}}$ when optimizing with respect to $\ensuremath{\vec{x}'}$:
\begin{equation}
\begin{split}
& \underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d}{\min}\; \pnorm{\ensuremath{\vec{x}'} - \ensuremath{\vec{x}}}_2^2 \\
& \Leftrightarrow \\
& \underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d}{\min}\; \frac{1}{2}\ensuremath{\vec{x}'}^\top\ensuremath{\vec{x}'} - \ensuremath{\vec{x}}^\top\ensuremath{\vec{x}'}
\end{split}
\end{equation}

\subsection{Minimizing the weighted Manhattan distance}
First, we transform the problem of minimizing the weighted Manhattan distance Eq.~\eqref{eq:weighted_l1} into epigraph form:
\begin{equation}
\begin{split}
& \underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d}{\min}\; \sum_j \alpha_j \cdot |(\ensuremath{\vec{x}'})_j - (\ensuremath{\vec{x}})_j| \\
& \Leftrightarrow \\
& \underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d, \beta\,\in\, \mathbb{R}}{\min}\; \beta \\
& \quad \text{s.t.} \quad \sum_j \alpha_j \cdot |(\ensuremath{\vec{x}'})_j - (\ensuremath{\vec{x}})_j| \leq \beta \\
& \quad \quad \quad \beta \geq 0
\end{split}
\end{equation}
Next, we separate the dimensions:
\begin{equation}
\begin{split}
& \underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d, \beta\,\in\, \mathbb{R}}{\min}\; \beta \\
& \quad \text{s.t.} \quad \sum_j \alpha_j \cdot |(\ensuremath{\vec{x}'})_j - (\ensuremath{\vec{x}})_j| \leq \beta \\
& \quad \quad \quad \beta \geq 0 \\
& \Leftrightarrow \\
& \underset{\ensuremath{\vec{x}'}, \vec{\beta} \,\in\, \mathbb{R}^d}{\min}\; \sum_j (\vec{\beta})_j \\
& \quad \text{s.t.} \quad \alpha_j \cdot |(\ensuremath{\vec{x}'})_j - (\ensuremath{\vec{x}})_j| \leq (\vec{\beta})_j \quad \forall\,j \\
& \quad \quad \quad (\vec{\beta})_j \geq 0 \quad \forall\,j
\end{split}
\end{equation}
After that, we remove the absolute value function:
\begin{equation}
\begin{split}
& \underset{\ensuremath{\vec{x}'}, \vec{\beta} \,\in\, \mathbb{R}^d}{\min}\; \sum_j (\vec{\beta})_j \\
& \quad \text{s.t.} \quad \alpha_j \cdot |(\ensuremath{\vec{x}'})_j - (\ensuremath{\vec{x}})_j| \leq (\vec{\beta})_j \quad \forall\,j \\
& \quad \quad \quad (\vec{\beta})_j \geq 0 \quad \forall\,j \\
& \Leftrightarrow \\
& \underset{\ensuremath{\vec{x}'}, \vec{\beta} \,\in\, \mathbb{R}^d}{\min}\; \sum_j (\vec{\beta})_j \\
& \quad \text{s.t.} \quad \alpha_j (\ensuremath{\vec{x}'})_j - \alpha_j (\ensuremath{\vec{x}})_j \leq (\vec{\beta})_j \quad \forall\,j \\
& \quad \quad \quad -\alpha_j (\ensuremath{\vec{x}'})_j + \alpha_j (\ensuremath{\vec{x}})_j \leq (\vec{\beta})_j \quad \forall\,j \\
& \quad \quad \quad (\vec{\beta})_j \geq 0 \quad \forall\,j
\end{split}
\end{equation}
Finally, we rewrite everything in matrix-vector notation:
\begin{equation}
\begin{split}
& \underset{\ensuremath{\vec{x}'}, \vec{\beta} \,\in\, \mathbb{R}^d}{\min}\;\vec{1}^\top\vec{\beta} \\
& \text{s.t.} \\
& \mat{\Upsilon}\ensuremath{\vec{x}'} - \mat{\Upsilon}\ensuremath{\vec{x}} \leq \vec{\beta} \\
& -\mat{\Upsilon}\ensuremath{\vec{x}'} + \mat{\Upsilon}\ensuremath{\vec{x}} \leq \vec{\beta} \\
& \vec{\beta} \geq \vec{0}
\end{split}
\end{equation}
where
\begin{equation}
\mat{\Upsilon} = \diag(\alpha_j)
\end{equation}

\subsection{Enforcing a specific prototype as the nearest neighbor}
By using the following set of strict inequalities, we can force the prototype $\vec{\prototype}_i$ to be the nearest neighbor of the counterfactual $\ensuremath{\vec{x}'}$ - which would cause $\ensuremath{\vec{x}'}$ to be classified as $\ensuremath{o}_i$:
\begin{equation}
\dist(\ensuremath{\vec{x}'}, \vec{\prototype}_i) < \dist(\ensuremath{\vec{x}'}, \vec{\prototype}_j) \quad \forall\, \vec{\prototype}_j\in\set{P}(\ensuremath{y^{c}})
\end{equation}
We consider a fixed pair of $i$ and $j$:
\begin{equation}\label{eq:nearestprototypeconstraint}
\begin{split}
& \dist(\ensuremath{\vec{x}'}, \vec{\prototype}_i) < \dist(\ensuremath{\vec{x}'}, \vec{\prototype}_j) \\
& \Leftrightarrow \pnorm{\ensuremath{\vec{x}'} - \vec{\prototype}_i}_{\distmat_i}^2 < \pnorm{\ensuremath{\vec{x}'} - \vec{\prototype}_j}_{\distmat_j}^2 \\
& \Leftrightarrow (\ensuremath{\vec{x}'} - \vec{\prototype}_i)^\top{\distmat}_i(\ensuremath{\vec{x}'} - \vec{\prototype}_i) < (\ensuremath{\vec{x}'} - \vec{\prototype}_j)^\top{\distmat}_j(\ensuremath{\vec{x}'} - \vec{\prototype}_j) \\
& \Leftrightarrow \ensuremath{\vec{x}'}^\top{\distmat}_i\ensuremath{\vec{x}'} - 2\ensuremath{\vec{x}'}^\top{\distmat}_i\vec{\prototype}_i + \vec{\prototype}_i^\top{\distmat}_i\vec{\prototype}_i < \ensuremath{\vec{x}'}^\top{\distmat}_j\ensuremath{\vec{x}'} - 2\ensuremath{\vec{x}'}^\top{\distmat}_j\vec{\prototype}_j + \vec{\prototype}_j^\top{\distmat}_i\vec{\prototype}_j \\
& \Leftrightarrow \ensuremath{\vec{x}'}^\top{\distmat}_i\ensuremath{\vec{x}'} - \ensuremath{\vec{x}'}^\top{\distmat}_j\ensuremath{\vec{x}'} - 2\ensuremath{\vec{x}'}^\top{\distmat}_i\vec{\prototype}_i + 2\ensuremath{\vec{x}'}^\top{\distmat}_j\vec{\prototype}_j + \vec{\prototype}_i^\top{\distmat}_i\vec{\prototype}_i - \vec{\prototype}_j^\top{\distmat}_i\vec{\prototype}_j < 0 \\
& \Leftrightarrow \ensuremath{\vec{x}'}^\top({\distmat}_i - {\distmat}_j)\ensuremath{\vec{x}'} + \ensuremath{\vec{x}'}^\top(-2{\distmat}_i\vec{\prototype}_i + 2{\distmat}_j\vec{\prototype}_j) + (\vec{\prototype}_i^\top{\distmat}_i\vec{\prototype}_i - \vec{\prototype}_j^\top{\distmat}_i\vec{\prototype}_j) < 0 \\
& \Leftrightarrow \frac{1}{2}\ensuremath{\vec{x}'}^\top({\distmat}_i - {\distmat}_j)\ensuremath{\vec{x}'} + \ensuremath{\vec{x}'}^\top({\distmat}_j\vec{\prototype}_j - {\distmat}_i\vec{\prototype}_i) + \frac{1}{2}(\vec{\prototype}_i^\top{\distmat}_i\vec{\prototype}_i - \vec{\prototype}_j^\top{\distmat}_i\vec{\prototype}_j) < 0 \\
& \Leftrightarrow \frac{1}{2}\ensuremath{\vec{x}'}^\top\mat{Q}_{ij}\ensuremath{\vec{x}'} + \ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij} < 0
\end{split}
\end{equation}
where
\begin{equation}
\mat{Q}_{ij} = {\distmat}_{i} - {\distmat}_{j}
\end{equation}
\begin{equation}
\vec{q}_{ij} = {\distmat}_j\vec{\prototype}_{j} - {\distmat}_i\vec{\prototype}_{i}
\end{equation}
\begin{equation}
r_{ij} = \frac{1}{2}\left(\vec{\prototype}_{i}\top{\distmat}_{i}\vec{\prototype}_{i} - \vec{\prototype}_{j}\top{\distmat}_{j}\vec{\prototype}_{j}\right)
\end{equation}
\mbox{}\\
If we only have a global distance matrix $\distmat$, we find that $\mat{Q}_{ij}=\mat{0}$ and the inequality Eq.~\eqref{eq:nearestprototypeconstraint} simplifies:
\begin{equation}
\begin{split}
& \dist(\ensuremath{\vec{x}}, \vec{\prototype}_i) < \dist(\ensuremath{\vec{x}}, \vec{\prototype}_j) \\
& \Leftrightarrow \ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij} < 0
\end{split}
\end{equation}
where
\begin{equation}
\vec{q}_{ij} = {\distmat}\left(\vec{\prototype}_{j} - \vec{\prototype}_{i}\right)
\end{equation}
\begin{equation}
r_{ij} = \frac{1}{2}\left(\vec{\prototype}_{i}^\top\distmat\vec{\prototype}_{i} - \vec{\prototype}_{j}^\top\distmat\vec{\prototype}_{j}\right)
\end{equation}
\mbox{}\\
If we do not use a custom distance matrix, we have $\distmat=\mat{\mathbb{I}}$ and Eq.~\eqref{eq:nearestprototypeconstraint} becomes:
\begin{equation}
\begin{split}
& \dist(\ensuremath{\vec{x}}, \vec{\prototype}_i) < \dist(\ensuremath{\vec{x}}, \vec{\prototype}_j) \\
& \Leftrightarrow \ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij} < 0
\end{split}
\end{equation}
where
\begin{equation}
\vec{q}_{ij} = \vec{\prototype}_{j} - \vec{\prototype}_{i}
\end{equation}
\begin{equation}
r_{ij} = \frac{1}{2}\left(\vec{\prototype}_{i}^\top\vec{\prototype}_{i} - \vec{\prototype}_{j}^\top\vec{\prototype}_{j}\right)
\end{equation}

\subsection{Approximately solving the non-convex QCQP}
Recall the non-convex quadratic constraint from Eq.~\eqref{eq:nearestprototypeconstraint}:
\begin{equation}\label{eq:nonconvexquadratic}
\frac{1}{2}\ensuremath{\vec{x}'}^\top\mat{Q}_{ij}\ensuremath{\vec{x}'} + \ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij} + \epsilon \leq 0
\end{equation}
where the matrix $\mat{Q}_{ij}$ is defined as the difference of two s.psd matrices:
\begin{equation}\label{eq:omegadiff}
\mat{Q}_{ij} = {\distmat}_{i} - {\distmat}_{j}
\end{equation}
By making use of Eq.~\eqref{eq:omegadiff}, we can rewrite Eq.~\eqref{eq:nonconvexquadratic} as:
\begin{equation}
\begin{split}
& \frac{1}{2}\ensuremath{\vec{x}'}^\top{\distmat}_{i}\ensuremath{\vec{x}'} + \ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij} + \epsilon - \frac{1}{2}\ensuremath{\vec{x}'}^\top{\distmat}_{j}\ensuremath{\vec{x}'} \leq 0 \\
& \Leftrightarrow f(\ensuremath{\vec{x}'}) - g(\ensuremath{\vec{x}'}) \leq 0
\end{split}
\end{equation}
where
\begin{equation}
f(\ensuremath{\vec{x}'}) = \frac{1}{2}\ensuremath{\vec{x}'}^\top{\distmat}_{i}\ensuremath{\vec{x}'} + \ensuremath{\vec{x}'}^\top\vec{q}_{ij} + r_{ij}  + \epsilon
\end{equation}
\begin{equation}
g(\ensuremath{\vec{x}'}) = \frac{1}{2}\ensuremath{\vec{x}'}^\top{\distmat}_{j}\ensuremath{\vec{x}'}
\end{equation}
Under the assumption that the regularization function $\regularization(\cdot)$ is a convex function\footnote{The weighted Manhattan distance and the Euclidean distance are convex functions!}, we can rewrite a generic version of the non-convex QCQP (induced by Eq.~\eqref{eq:lgmlvq:constraints}) as follows:
\begin{equation}\label{eq:dcp}
\begin{split}
& \underset{\ensuremath{\vec{x}'} \,\in\, \mathbb{R}^d}{\min}\;\regularization(\ensuremath{\vec{x}'}, \ensuremath{\vec{x}}) \\
& \text{s.t.} \\
& f(\ensuremath{\vec{x}'}) - g(\ensuremath{\vec{x}'}) \leq 0
\end{split}
\end{equation}
Because ${\distmat}_{i}$ and ${\distmat}_{j}$ are s.psd matrices, we know that $f(\ensuremath{\vec{x}'})$ and $g(\ensuremath{\vec{x}'})$ are convex functions. Therefore Eq.~\eqref{eq:dcp} is a difference-of-convex program (DCP).

This allows us to use the penalty convex-concave procedure (CCP)~\cite{park2017general} for computing an approximate solution of Eq.~\eqref{eq:dcp} that is equivalent to the original non-convex QCQP induced by Eq.~\eqref{eq:lgmlvq:constraints}. For using the penalty CCP, we need the first order Taylor approximation of $g(\ensuremath{\vec{x}'})$ around a current point $\ensuremath{\vec{x}}_k$:
\begin{equation}
\begin{split}
\hat{g}(\ensuremath{\vec{x}'})_{\ensuremath{\vec{x}}_k} &= g(\ensuremath{\vec{x}}_k) + (\nabla_{\ensuremath{\vec{x}'}}g)(\ensuremath{\vec{x}}_k)^\top(\ensuremath{\vec{x}'} - \ensuremath{\vec{x}}_k) \\
&= \frac{1}{2}\ensuremath{\vec{x}}_k^\top{\distmat}_{j}\ensuremath{\vec{x}}_k + ({\distmat}_{j}\ensuremath{\vec{x}}_k)^\top(\ensuremath{\vec{x}'} - \ensuremath{\vec{x}}_k) \\
&= ({\distmat}_{j}\ensuremath{\vec{x}}_k)^\top\ensuremath{\vec{x}'} + \frac{1}{2}\ensuremath{\vec{x}}_k^\top{\distmat}_{j}\ensuremath{\vec{x}}_k -({\distmat}_{j}\ensuremath{\vec{x}}_k)^\top\ensuremath{\vec{x}}_k \\
&= \vec{\rho}_{jk}^\top\ensuremath{\vec{x}'} + \tilde{r}_{jk}
\end{split}
\end{equation}
where
\begin{equation}
\vec{\rho}_{jk} = {\distmat}_{j}\ensuremath{\vec{x}}_k
\end{equation}
\begin{equation}
\tilde{r}_{jk} = -\frac{1}{2}\ensuremath{\vec{x}}_k^\top{\distmat}_{j}\ensuremath{\vec{x}}_k
\end{equation}

\begin{footnotesize}

\bibliographystyle{unsrt}


\section*{Introduction}
Single-cell RNA sequencing (scRNA-seq) has emerged as a powerful tool for unbiased identification of previously uncharacterized
molecular heterogeneity at the cellular level \cite{shapiro2013single}. This is in contrast to standard bulk RNA-seq techniques \cite{nagalakshmi2008transcriptional}, which measures average gene expression levels within a cell population, and thus ignore tissue heterogeneity. Consideration of cell-level variability of gene expressions is essential for extracting signals from complex heterogeneous tissues \cite{macosko2015highly}, and also for understanding dynamic biological processes, such as embryo development \cite{deng2014single} and cancer \cite{patel2014single}.

A large body of statistical tools developed for scRNA-seq data analysis include a dimensionality reduction step. This leads to more tractable data, from both statistical and computational point of views. Moreover, the noise in the data can be decreased, while retaining the often intrinsically low-dimensional signal of interest. Dimensionality reduction of scRNA-seq data is challenging. In addition to high gene expression variability due to cell heterogeneity, the excessive amount of zeros in scRNA-seq hinders the application of classical dimensionality reduction techniques such as principal component analysis (PCA). For instance, in real-world datasets, it has been reported that the first or second principal components often depend more on the proportion of detected genes per cell (i.e., genes with at least one read) than on the actual biological signal \cite{finak2015mast}.

Several existing computational tools adopt explicit zero-inflation modeling to infer the latent representation of scRNA-seq data. Zero-inflated factor analysis (ZIFA) \cite{pierson2015zifa} extends the framework of probabilistic PCA \cite{tipping1999probabilistic} to the zero-inflated setting, by modeling the excessive zeros using Bernoulli distributed random variables which indicate the dropout event. Zero-inflated negative binomial-based wanted variation extraction (ZINB-WaVE) \cite{risso2018general} directly models the scRNA-seq counts using a zero-inflated negative binomial distribution, while accounting for both gene- and cell-level covariates. It infers the model parameters using a penalized maximum likelihood procedure.

Despite its popularity, using an explicit zero-inflation term may place unnecessary emphasis on the zero counts, leading to complication in discovering the latent representation of scRNA-seq data. In this paper, we propose a hierarchical gamma-negative binomial (hGNB) model to both perform dimensionality reduction and adjust for the effects of the gene- and cell-level confounding factors simultaneously. Exploiting the hierarchical structure, the proposed hGNB model is capable of capturing the high over-dispersion present in the scRNA-seq data. More precisely, we factorize the logit of the negative-binomial (NB) distribution probability parameter to identify latent representation of the data. In addition to factorization, linear regression terms are also included in that logit function to adjust for the impact of covariates.

In hGNB, a gamma distribution with varying rate parameter is used to model the cell dependent dispersion parameter of the NB distribution. The cell-level dispersion serves as a means of representing the prevalence of the dropout events. For instance, cells that are sequenced deeply will naturally include less dropped-out genes with zero counts, and thus this will be reflected in the cell specific dispersion parameter of NB distribution.

We follow a Bayesian framework, similar to bulk RNA-seq setting \cite{dadaneh2017bnp,dadaneh2018bayesian, zamani2018covariate}, and derive closed-form Gibbs sampling update equations for the model parameters of hGNB, by exploiting sophisticated data augmentation techniques. More specifically, we apply the data augmentation technique of \cite{zhou2015negative} (2015)  for the NB distribution, and the Polya-Gamma distributed auxiliary variable technique of \cite{polson2013bayesian} (2013) for the closed-form inference of regression coefficients and also latent factor parameters, removing the need for non-trivial Metropolis-Hastings correction steps \cite{chib1995understanding}. Experimental results on several real-world scRNA-seq datasets demonstrate the superior performance of hGNB to identify cell clusters, especially in complex settings, and also its potential application in cell lineages inference.

\section*{Methods}

\subsection*{hGNB model}

In this section we present the hierarchical gamma-negative binomial (hGNB) model for factor analysis of scRNA-seq data. The graphical representation of hGNB is shown in Figure~\ref{fig:hgnb}. The parameters of the hGNB model with their interpretations in the context of scRNA-seq experiments are presented in Table~\ref{tab:model}. Let $n_{vj}$ denote the number of sequencing reads mapped to gene $v \in \{1,...,V\}$ in the cell $j \in \{1,...,J\}$. Under the hGNB model, gene counts are distributed according to a negative binomial (NB) distribution:
\begin{equation}
	n_{vj} \sim \mbox{NB}(r_j,p_{vj}),
\end{equation}
where $r_j$ and $p_{vj}$ are dispersion and probability parameters of NB distribution, respectively. The probability mass function (PMF) of this distribution can be expressed as $f_{N}(n_{vj}) = \frac{\Gamma(n_{vj}+r_j)}{n_{vj}!\Gamma(r_j)} p_{vj}^{n_{vj}} (1-p_{vj})^{r_j}$, where $\Gamma(\cdot)$ is the gamma function.

Data from scRNA-seq experiments exhibit high variability between different cells, even for genes with medium or high levels of expression. To capture this variability, we impose a gamma prior on the cell-level dispersion parameters as
\begin{equation}
	r_j \sim \mbox{Gamma}(e_0,1/h),
	\label{eq:rj}
\end{equation}
where for simplification, the hyper-parameter $e_0$ is set to 0.01 in our experiments, and the rate $h$ is learned during the Gibbs sampling inference, presented in the following section. This hierarchical prior on the dispersion parameter, enhances the flexibility of NB distribution to capture the high over-dispersion of scRNA-seq counts, without the need for explicit zero-inflation modeling.

To account for various technical and biological effects common in scRNA-seq technologies, we impose a regression model on the logit of NB probability parameter as
\begin{equation}
	\psi_{vj} = \mbox{logit}(p_{vj}) = \boldsymbol{\beta}_v^T\boldsymbol{x}_j + \boldsymbol{\delta}_j^T\boldsymbol{z}_v + \boldsymbol{\phi}_v^T\boldsymbol{\theta}_j.
	\label{eq:p}
\end{equation}
The three terms in the summation are described below.

In the first term,
$\boldsymbol{x}_j$ is a known vector of $P$ covariates for cell $j$ and $\boldsymbol{\beta}_v$ is the regression-coefficient vector adjusting the effect of covariates on gene $v$. The covariate vector $\boldsymbol{x}_j$ can represent variations of interest, such as cell types, or unwanted variations, such as batch effects or quality control measures. An intercept term can also be included in these cell-level covariates to account for gene dependent baseline expressions.

\begin{figure*}[t]
	\centering
	\begin{tikzpicture}[scale=0.95\textwidth/1cm]

	\node[obs] (n) {$n_{vj}$};
	\node[latent,right=of n] (r) {$r_j$};

	\node[latent,above=of n,yshift=1.2cm] (phi) {$\boldsymbol{\phi}_v$};
	\node[obs,above=of n,yshift=1.2cm,xshift=1cm] (z) {$\boldsymbol{z}_v$};
	\node[latent,above=of n,yshift=1.2cm,xshift=-1cm] (beta) {$\boldsymbol{\beta}_v$};
	\node[latent,left=of n,yshift=1cm,xshift=-1cm] (theta) {$\boldsymbol{\theta}_j$};
	\node[latent,left=of n,yshift=-1cm,xshift=-1cm] (delta) {$\boldsymbol{\delta}_j$};
	\node[obs,below=of n] (x) {$\boldsymbol{x}_j$};
	\node[latent,left=of theta,xshift=-0.2cm] (gamma) {$\boldsymbol{\gamma}$}; %
	\node[latent,left=of delta,xshift=-0.2cm] (eta) {$\boldsymbol{\eta}$};
	\node[latent,left=of beta,xshift=-0.2cm] (alpha) {$\boldsymbol{\alpha}$};
	\node[latent,right=of r] (h) {$h$};

	\plate [inner sep=0.3cm] {plate1} {(n)} {$v=1,...,V$};
	\plate [inner sep=0.3cm] {plate2} {(n)(r)(x)(theta)(delta)} {$j=1,...,J$};
	\plate [inner sep=0.3cm] {plate3} {(phi)(beta)(z)} {$V$};

	\edge {phi,theta,beta,x,z,delta,r} {n};

	\edge {gamma} {theta};
	\edge {h} {r};
	\edge {alpha} {beta};
	\edge {eta} {delta};

	\end{tikzpicture}
	\caption{Graphical representation of the hierarchical gamma-negative binomial (hGNB) model.}
	\label{fig:hgnb}
\end{figure*}

In the second term, $\boldsymbol{z}_v$ is a vector of $Q$ covariates for gene $v$, representing gene length or GC-content for example \cite{risso2011gc}, and $\boldsymbol{\delta}_j$ is its associated regression-coefficient vector. We also include a fixed intercept element in $\boldsymbol{z}_v$ to account for cell-specific expressions, such as the size factors representing differences in sequencing depth.

In the third term, $\boldsymbol{\phi}_v^T\boldsymbol{\theta}_j$ corresponds to the latent factor representation of the count $n_{vj}$, after accounting for the effects of gene- and cell-level covariates. More precisely, the unknown $K \times 1$ vector $\boldsymbol{\phi}_v$ contains the factor loading parameters which determine the association between genes and latent factors. Moreover, the unknown $K \times 1$ vector $\boldsymbol{\theta}_j$ encodes the popularity of the $K$ factors in the expression of cell $j$.

We place independent zero-mean normal distributions on the components of the regression coefficient parameters $\boldsymbol{\beta}_v$ and $\boldsymbol{\delta}_j$ as
\begin{eqnarray}
	\boldsymbol{\beta}_v &\sim& \prod_{p=1}^{P} \mbox{N}(\beta_{vp};0,\alpha_p^{-1}),\nonumber\\
	\boldsymbol{\delta}_j &\sim& \prod_{q=1}^{Q} \mbox{N}(\delta_{jq};0,\eta_q^{-1}),
	\label{eq:reg}
\end{eqnarray}
where $\alpha_p$ and $\eta_q$ are precision parameters of the normal distributions and gamma priors are imposed on them. These priors are known as automatic relevance determination (ARD), which are effective tools for pruning large numbers of irrelevant covariates \cite{wipf2008new,tipping2001sparse}. In addition, by assuming identical precision for components of the regression coefficients
across all genes or samples, hGNB borrows statistical strengths to infer these precision parameters.

We impose independent normal priors on latent factor loading and score parameters $\boldsymbol{\phi}_v$ and $\boldsymbol{\theta}_j$:
\begin{eqnarray}
\boldsymbol{\phi}_v &\sim& \mbox{N}(\boldsymbol{\phi}_v;0,I_K),\nonumber\\
\boldsymbol{\theta}_j &\sim& \prod_{k=1}^{K} \mbox{N}(\theta_{jk};0,\gamma_k^{-1}).
\label{eq:fac}
\end{eqnarray}
Note that the posterior for these terms is not generally independent or normal, but accounts for the statistical
dependence as reflected in the data.

We complete the model by imposing a gamma prior on the precision parameters of normal distributions, and also the rate parameter of gamma distributions. Specifically, throughout the experiments, we set both the shape and rate of these gamma priors to 0.01.

\begin{table*}[h]
	\centering
	\caption{Parameters of the hierarchical gamma-negative binomial (hGNB) model and their interpretations in the context of scRNA-seq data. The inputs of hGNB are gene counts $n_{vj}$ and vector of cell- and gene-level covariates $\mathbf{x}_j$ and $\mathbf{z}_v$.}
	\small
	\begin{tabular}{ c| c | c}
		Parameter & Constraint & Interpretation \\
		\hline \hline
		$r_j$ & $r_j > 0$ & expression heterogeneity of genes in sample $j$\\
		\hline
		$\phi_{vk}$ & $\sum_{v=1}^{V} \phi_{vk}=1$, $\phi_{vk}>0$ & gene-latent factor association \\
		\hline
		$\theta_{jk}$ & $\theta_{kj} > 0$ & popularity of factor $k$ in sample $j$ \\
		\hline
		$\beta_{vp}$ & $\beta_{vp} \in \mathbb{R}$ & impact of cell covariate $p$ on expression
		of gene $v$\\
		\hline
		$\delta_{jq}$ & $\beta_{vp} \in \mathbb{R}$ & impact of gene covariate $q$ on expression
		of cell $j$

	\end{tabular}
	\label{tab:model}
\end{table*}

\subsection*{Inference via Gibbs sampling}
In this section, we provide an efficient inference algorithm that adopts data augmentation techniques tailored to our hGNB model. Algorithm~1 summarizes all the steps in the Gibbs sampling algorithm.

\paragraph{Sample dispersion parameter} We start with the data augmentation technique developed for inferring the NB dispersion parameter \cite{zhou2015negative}. More precisely, the negative binomial random variable $n \sim \text{NB}(r,p)$ can be generated from a compound Poisson distribution as
\begin{equation}
n = \sum_{t=1}^{\ell} u_t, \;\; u_t \sim \text{Log}(p), \;\; \ell \sim \text{Pois}(-r\ln (1-p)), \nonumber
\end{equation}
where $u \sim \text{Log}(p)$ corresponds to the logarithmic random variable \cite{johnson2005univariate}, with the PMF $f_U(u) = -\frac{p^u}{u\ln(1-p)}$, $u \in \{1,2,...\}$. As shown in \cite{zhou2015negative}, given $n$ and $r$, the distribution of $\ell$ is a Chinese Restaurant Table (CRT) distribution, $(\ell | n,r) \sim \text{CRT}(n,r)$, which can be generated  as $\ell = \sum_{t=1}^{n} b_t, b_t\sim \text{Bernoulli}(\frac{r}{r+t-1})$.

Utilizing this augmentation technique, for each observed count $n_{vj}$, an auxiliary count is sampled as
\begin{equation} \label{eq:crt}
(\ell_{vj}|-) \sim \mbox{CRT}(n_{vj},r_j).
\end{equation}
Using gamma-Poisson conjugacy, the cell-dependent dispersion parameters are updated as
\begin{equation} \label{eq:rji}
(r_j | -) \sim \mbox{Gamma} \Big( e_0+ \sum_v \ell_{vj}, \frac{1}{h - \sum_v \ln(1-p_{vj})} \Big).
\end{equation}

\paragraph{Sample regression coefficients} For the regression coefficients modeling potential covariate effects, the lack of conditional conjugacy precludes immediate closed-form inference. Therefore we adopt another data augmentation technique, specifically designed for hGNB, to infer the regression coefficients $\boldsymbol{\beta}_v$ and $\boldsymbol{\delta}_j$, relying on the Polya-Gamma (PG) data augmentation~\cite{LGNB_ICML2012,polson2013bayesian}.

Denote $\omega_{vj}$ as a random variable drawn from the PG distribution as $\omega_{vj} \sim \text{PG}(n_{vj}+r_j,0).$
Since  $\mathbb{E}_{\omega_{vj}}[\exp(- \omega_{vj} \psi_{vj}^2/2)]=\cosh^{(n_{vj}+r_j)}(\psi_{vj}^2/2)$,
the likelihood of $\psi_{vj}$ in~(\ref{eq:p}) can be expressed as
\begin{align}
\mathcal{L}(\psi_{vj}) &\propto \frac{(e^{\psi_{vj}})^{n_{vj}}}{(1+e^{\psi_{vj}})^{n_{vj}+r_j}} \nonumber\\
& \propto \exp\Big(\frac{n_{vj}-r_j}{2} \psi_{vj}\Big) \mathbb{E}_{\omega_{vj}}[\exp(- \omega_{vj} \psi_{vj}^2/2)].
\label{eq:lik}
\end{align}
Exploiting the exponential tilting of the PG distribution in \cite{polson2013bayesian}, we draw
$\omega_{vj}$ as
\begin{equation}\label{eq:omega}
(\omega_{vj}|-) \sim \text{PG}(n_{vj}+r_j,\psi_{vj}).
\end{equation}
Given the values of the auxiliary variables $\omega_{vj}$ for $j=1,...,J$ and the prior in~(\ref{eq:reg}), the conditional posterior of $\boldsymbol{\beta}_v$ can be updated as
\begin{equation}
(\boldsymbol{\beta}_v|-) \sim \mbox{N}(\mu_v^{(\beta)},\Sigma_v^{(\beta)}),
\label{eq:beta}
\end{equation}
where $\Sigma_v^{(\beta)} = \Big( \mbox{diag}(\alpha_1,...,\alpha_P)+\sum_{j} \omega_{vj} \boldsymbol{x}_j \boldsymbol{x}_j^T \Big)^{-1}$ and $\mu_v^{(\beta)} = \Sigma_v^{(\beta)} \Big[ \sum_{j} \big(\frac{n_{vj}-r_j}{2} - \omega_{vj} (\boldsymbol{\delta}_j^T\boldsymbol{z}_v + \boldsymbol{\phi}_v^T\boldsymbol{\theta}_j) \big) \boldsymbol{x}_j \Big]$. \\

A similar procedure can be followed to derive the conditional updates for cell-level regression coefficients as
\begin{equation}
(\boldsymbol{\delta}_j|-) \sim \mbox{N}(\mu_j^{(\delta)},\Sigma_j^{(\delta)}),
\label{eq:delta}
\end{equation}
where $\Sigma_j^{(\delta)} = \Big( \mbox{diag}(\eta_1,...,\eta_Q)+\sum_{v} \omega_{vj} \boldsymbol{z}_v \boldsymbol{z}_v^T \Big)^{-1}$ and $\mu_j^{(\delta)} = \Sigma_j^{(\delta)} \Big[ \sum_{v} \big(\frac{n_{vj}-r_j}{2} - \omega_{vj} (\boldsymbol{\beta}_v^T\boldsymbol{x}_j + \boldsymbol{\phi}_v^T\boldsymbol{\theta}_j) \big) \boldsymbol{z}_v \Big]$.

\paragraph{Sample latent factor parameters} Using the likelihood function in (\ref{eq:lik}) and the priors in (\ref{eq:fac}), we can derive closed-form update steps for factor loading and score parameters. More specifically, the full conditional for factor loading $\boldsymbol{\phi}_v$ is a normal distribution:
\begin{equation}
(\boldsymbol{\phi}_v|-) \sim \mbox{N}(\mu_v^{(\phi)},\Sigma_v^{(\phi)}),
\label{eq:phi}
\end{equation}
where $\Sigma_v^{(\phi)} = \Big( I_K + \sum_{j} \omega_{vj} \boldsymbol{\theta}_j \boldsymbol{\theta}_j^T \Big)^{-1}$ and $\mu_v^{(\phi)} = \Sigma_v^{(\phi)} \Big[ \sum_{j} \big(\frac{n_{vj}-r_j}{2} - \omega_{vj} (\boldsymbol{\beta}_v^T\boldsymbol{x}_j + \boldsymbol{\delta}_j^T\boldsymbol{z}_v) \big) \boldsymbol{\theta}_j \Big]$.

The full conditional for factor score $\boldsymbol{\theta}_j$ is also a normal distribution:
\begin{equation}
(\boldsymbol{\theta}_j|-) \sim \mbox{N}(\mu_j^{(\theta)},\Sigma_j^{(\theta)}),
\label{eq:theta}
\end{equation}
where $\Sigma_j^{(\theta)} = \Big( \mbox{diag}(\gamma_1,...,\gamma_K) + \sum_{v} \omega_{vj} \boldsymbol{\phi}_v \boldsymbol{\phi}_v^T \Big)^{-1}$ and $\mu_j^{(\theta)} = \Sigma_j^{(\theta)} \Big[ \sum_{v} \big(\frac{n_{vj}-r_j}{2} - \omega_{vj} (\boldsymbol{\beta}_v^T\boldsymbol{x}_j + \boldsymbol{\delta}_j^T\boldsymbol{z}_v) \big) \boldsymbol{\phi}_v \Big]$.

\paragraph{Sample precision and rate} The precision parameters of normal distributions in (\ref{eq:reg}) and (\ref{eq:fac}) can be updated using the normal-gamma conjugacy:
\begin{eqnarray}
	\alpha_p &\sim& \mbox{Gamma} \big( e_0+V/2, \frac{1}{f_0 + \sum_{v=1}^{V} \beta_{vp}/2} \big), \nonumber\\
	\eta_q &\sim& \mbox{Gamma} \big( e_0+J/2, \frac{1}{f_0 + \sum_{v=1}^{V} \delta_{jq}/2} \big), \nonumber\\
	\gamma_k &\sim& \mbox{Gamma} \big( e_0+J/2, \frac{1}{f_0 + \sum_{v=1}^{V} \theta_{jk}/2} \big).
	\label{eq:hyp}
\end{eqnarray}

Finally, the rate of gamma distribution in (\ref{eq:rj}) can be updated using the gamma-gamma conjugacy with respect to the rate parameter:
\begin{equation}
	h \sim \mbox{Gamma} \big( e_0(1+J), \frac{1}{f_0 + \sum_{j=1}^{J} r_j} \big).
\end{equation}

\begin{algorithm}[t]
	\caption{hGNB model inference}\label{alg:gibbs}
	\textbf{Inputs}: scRNA-seq counts, design matrix of covariate effects, $N$\\
	\textbf{Outputs}: gene module membership matrix
	\begin{algorithmic}
		\State \textit{Initialize} model parameters
		\State \# Do Gibbs sampling:
		\For {$iter=1$ to $N$}
		\State Sample $\ell_{vj}$ using the CRT distribution (eq. (\ref{eq:crt}))
		\State Update $r_{j}$ using the gamma-Poisson conjugacy (eq. (\ref{eq:rji}))
		\State Sample auxiliary variables $\omega_{vj}$, using the PG distribution (eq. (\ref{eq:omega}))
		\State Update cell- and gene-level regression coefficients (eq. (\ref{eq:delta}),(\ref{eq:beta}))
		\State Update factor loadings and scores (eq. (\ref{eq:phi}),(\ref{eq:theta}))
		\State Update $\alpha_p$, $\eta_q$ and $\gamma_k$ (eq. (\ref{eq:hyp}))
		\EndFor \textbf{end for}
		\State
	\end{algorithmic}
\end{algorithm}

\section*{Results}
We evaluate our hGNB model on four different sets of real-world scRNA-seq data from different platforms, and compare its performance to those of principal component analysis (PCA), ZIFA \cite{pierson2015zifa}, and ZINB-WaVE \cite{risso2018general}. In the following, We briefly describe these scRNA-seq datasets. To pre-process these datasets when needed, we followed the same procedures as in \cite{risso2018general}.

\textbf{V1 dataset}. This dataset characterizes more than 1600 cells from the primary visual cortex (V1) in adult male mice, using a set of established Cre lines \cite{tasic2016adult}. A subset of three Cre lines, including Ntsr1-Cre, Rbp4-Cre, and Scnn1a-Tg3-Cre, that respectively label layer 4, layer 5, and layer 6 excitatory neurons were selected. We only retained 285 cells that passed the authors' quality control (QC) filters. The dimensionality reduction methods were only applied to the 1000 most variable genes.

\textbf{S1/CA1 dataset}. This dataset characterizes 3005 cells from the primary somatosensory cortex (S1) and the hippocampal CA1 region, using the Fluidigm C1 microfluidics cell capture platform followed by Illumina sequencing \cite{zeisel2015cell}. Gene expression is quantified by UMI counts.

\textbf{mESC dataset}. This dataset includes the transcriptome measurement of 704 mouse embryonic stem cells (mESCs), across three culture conditions (serum, 2i, and a2i), using the Fluidigm C1 microfluidics cell capture platform followed by Illumina sequencing \cite{kolodziejczyk2015single}. We excluded the samples that did not pass the authors's QC filters, resulting in a total of 169 serum cells, 141 2i cells, and 159 a2i cells. The dimensionality reduction methods were only applied to the 1000 most variable genes.

\textbf{OE dataset}. This data characterizes 849 FACS-purified cells from the mouse OE, using the Fluidigm C1 microfluidics cell capture platform followed by Illumina sequencing \cite{fletcher2017deconstructing}. We followed the filtering procedure of \cite{perraudeau2017bioconductor}, and filtered the cells that exhibited poor sample quality, retaining a total of 747 cells.

For all datasets, hGNB was run using 2000 MCMC iterations, where after the first 1000 burn-in iterations, the posterior samples with the highest likelihood were collected as the point estimates of model parameters corresponding to latent factors. In the dimensionality reduction analysis below, following \cite{risso2018general}, for S1/CA1 dataset we set the number of latent factors $K=3$, and for V1 and mESC we set $K=2$.

\subsection*{Goodness-of-fit of hGNB model}
We have examined the goodness-of-fit of hGNB model on V1, S1/CA1 and mESC datasets, using the mean-difference (MD) plots. Figure~\ref{fig:md} shows the MD plot for the S1/CA1 dataset, where the y-axis is the difference between observed counts and the expected counts under hGNB, and x-axis is the average of these two sets of counts. The solid red line in this figure, which represents the local regression fit \cite{shyu2017local} to the data, resides near zero for various average levels. This supports the good fit of hGNB model to the highly over-dispersed scRNA-seq data. Similar trends are observed for V1 and mESC datasets (supplementary materials).

\begin{figure}[t]
	\centering
	[width=0.45\textwidth]{fig1}
	\caption{Mean-difference (MD) plot for S1/CA1 dataset. The solid red line represents the local regression fit to the data}
	\label{fig:md}
\end{figure}

\subsection*{Capturing zero-inflation}
Next we evaluate the performance of hGNB on simulated data based on the zero-inflated NB distribution of \cite{risso2018general} to show that hGNB faithfully captures zero inflation without the need of explicit zero-inflation modeling. Specifically, the capability of hGNB to recover true clustering structure of cells under three zero-count prevalence levels with two different total numbers of cells. The parameters of the simulating zero-inflated model were learned based on the S1/CA1 dataset. Genes that did not have at least five reads in at least five cells were filtered out and 1000 genes were then sampled at random for each dataset. The number of latent factors was set to $K=2$. To simulate cell clustering, a $K$-variate Gaussian mixture distribution with three components was fitted to the inferred factor score parameters, and then for each simulated dataset, factor scores were generated from $K$-variate Gaussian distributions. By adjusting the value of regression coefficients in the zero-inflation term of ZINB-WaVE model, we generated synthetic datasets with three levels of zero-count percentages as 40\%, 60\% and 80\% (for details refer to \cite{risso2018general}). The number of cells were set to $J=100$ and $J=1000$. For each scenario, including cell numbers and zero-count prevalence (sparsity) levels, we simulated 10 datasets.

We evaluate the performance of our method for the clustering task based on the average silhouette width measure. The silhouette width $s_j$ of sample $j$ is defined as
$$ s_j = \frac{b_j-a_j}{\max \{a_j,b_j\}}, $$
where $a_j$ is the average distance between sample $j$ and all samples in the cluster that it belongs to, and $b_j$ is the minimum average distance between sample $j$ and samples in other clusters.

Figure~\ref{fig:sim} shows the clustering average silhouette width based on the above simulation setup, for different zero-count prevalence levels and cell numbers. In the setting with small sample size, for 40\% and 60\% zero fractions, hGNB has the best clustering silhouette width, and for the 80\% zero fraction its performance is identical to that of ZINB-WaVE. In the setting with moderate sample size, hGNB has the best clustering silhouette width for 40\% zero fraction, and for 60\% and 80\% zero fractions it closely follows the performance of ZINB-WaVE. This suggests that the hierarchical structure of hGNB equips it with the capacity to capture highly over-dispersed count data, even though an explicit zero-inflation term is not included in its model. Also, ZINB-WaVE requires large enough samples to have robust inference results due to the introduction of zero-inflation terms in its model. Finally, ZIFA and PCA have the worst performance, as they normalize the data before learning its latent representation.

\begin{figure}[t]
	\centering
	[width=0.45\textwidth]{fig0}
	\caption{(\textbf{a})$J=100$, (\textbf{b})$J=1000$. Performance of different methods based on recovering the true cell clusters in synthetic data based on S1/CA1 dataset. Zero-inflated NB model of ZINB-WaVE is used to simulate scRNA-seq data.}
	\label{fig:sim}
\end{figure}

\subsection*{Dimensionality reduction}
We applied hGNB to the three scRNA-seq datasets, V1, S1/CA1 and mESC, to assess its power to separate cell clusters in the low dimensional space, and compared it to PCA, ZIFA, and ZINB-WaVE methods. Figure~\ref{fig:vis} illustrates the projected scRNA-seq expression of profiled cells in the two-dimensional space for S1/CA1 dataset. The proposed hGNB model provides more biologically meaningful latent representations of scRNA-seq gene expressions for S1/CA1 cells, especially compared to PCA and ZIFA that do not model the counts directly. Furthermore, hGNB leads to more separated clusters of cells in the two-dimensional space, compared to ZINB-WaVE. Specifically, hGNB distinguishes {\tt microglia} from {\tt endothelial−mural} cells, while ZINB-WaVE fails to accomplish this task.

To examine the dimensionality reduction results more carefully, we used the average silhouette width as a measure of goodness for clustering.

\begin{figure*}[t]
	\centering
	[width=0.95\textwidth]{fig2}
	\caption{Low-dimensional representations of the S1/CA1 dataset. Panels correspond to (\textbf{a}) PCA (on total-count normalized data), (\textbf{b}) ZIFA (on total-count normalized data), (\textbf{c}) ZINB-WaVE, and (\textbf{d}) hGNB.}
	\label{fig:vis}
\end{figure*}

\begin{figure*}[h]
	\centering
	[width=0.95\textwidth]{fig3}
	\caption{Average silhouette width in scRNA-seq datasets (\textbf{a}) S1/CA1, (\textbf{b}) mESC, and (\textbf{c}) V1. Silhouette widths were
		computed in the low-dimensional space, using the groupings provided by the authors of the original publications. PCA and ZIFA were applied with both unnormalized (RAW) data and after total count (TC) normalization.}
	\label{fig:sil}
\end{figure*}

Figure~\ref{fig:sil} shows the average silhouette width of different methods on V1, S1/CA1, and mESC datasets. For PCA and ZIFA, the results on both raw counts and normalized counts are included in this figure. For S1/CA1 dataset, which has the highest number of clusters, the proposed hGNB method outperforms all other methods in terms of clustering average silhouette. For mESC dataset, performance of hGNB is comparable to ZINB-WaVE, and it is significantly better than PCA and ZIFA. For V1 dataset, however, we observe that hGNB, besides PCA applied to raw counts, possess the lowest average silhouette. By further examination of the latent representations of cells for this dataset (supplementary materials), we observe that all methods split the {\tt Rbp4−Cre\_KL100} cells into two clusters, one of them located near {\tt Scnn1a−Tg3−Cre} cells, suggesting the presence of batch effects, which have led to confounding of latent representations \cite{risso2018general}.

\subsection*{Identification of developmental lineages}
In addition to characterization of cell types, we further demonstrate the capability of hGNB to derive novel biological insights, by analyzing a set of cells from the mouse olfactory epithelium (OE). The samples were collected to identify the developmental trajectories that generate olfactory neurons (mOSN), sustentacular cells (mSUS), and microvillous cells (MV) \cite{fletcher2017deconstructing}.

We first performed dimensionality reduction on the OE dataset by applying hGNB with $K=50$. Next, we clustered the cells using the low-dimensional factor score parameters $\theta_{kj}$. More specifically, the resampling-based sequential ensemble clustering (RSEC) framework implemented in the {\tt RSEC} function from the Bioconductor R package {\tt clusterExperiment} \cite{purdom2017clusterexperiment} was applied to factor scores, leading to identification of 14 cell clusters. The correspondence between the detected clusters and the underlying biological cell types is presented in Table~\ref{tab:oe}. In addition to these already known cell clusters in OE, hGNB is able to detect new clusters, potentially offering novel biological insights.

\begin{table}[t]
	\centering
	\caption{Correspondence between identified clusters and cell types in OE dataset.}
	\small
	\begin{tabular}{| c| c |  }
		\hline
		Cell Type & Clusters \\
		\hline \hline
		GBC & cl4,cl9 \\
		mSUS & cl2,cl3,cl5,cl11 \\
		mOSN & cl8,cl12,cl3 \\
		Immature Neurons & cl10 \\
		MV & cl14\\
		\hline
	\end{tabular}
	\label{tab:oe}
\end{table}

\begin{figure}[!ht]
	\centering
	[width=0.45\textwidth]{fig4}
	\caption{Lineage inference on the OE dataset. The low dimensional data representation derived by hGNB were used to cluster cells by {\tt RSEC}. The minimum spanning tree (MST) of the derived clusters constructed by {\tt slingshot} is also displayed.}
	\label{fig:lin}
\end{figure}

We further investigated the potential benefit of using the learned latent representation by our proposed hGNB model to infer branching cell lineages and order cells by developmental progression along each lineage. To infer the global lineage structure (i.e., the number of lineages and where they branch), a minimum spanning tree (MST) was constructed on the clusters identified above by RSEC. We used the R package {\tt slingshot} \cite{street2018slingshot}. Figure~\ref{fig:lin} illustrates the inferred lineages for the OE dataset, in a two-dimensional space obtained by applying multi-dimensional scaling (MDS) algorithm to the factor scores learned by hGNB. There are three branches in the inferred lineages, with endpoints located in microvillous (MV), mature olfactory sensory neurons (mOSN), and mature sustentacular (mSUS) cells.

\section*{Conclusions}
We propose a hierarchical Bayesian gamma-negative binomial (hGNB) model for extracting low dimensional representations from single-cell RNA sequencing (scRNA-seq) data. hGNB obviates the need for explicit modeling of the zero-inflation prevalent in scRNA-seq count data. Our hGNB can naturally account for covariate effects at both the gene and cell levels, and does not require the commonly adopted pre-processing steps such as normalization. By taking advantage of sophisticated data augmentation techniques, hGNB possesses efficient closed-form Gibbs sampling update equations.
Our experimental results on real-world scRNA-seq data demonstrates that hGNB is capable of identifying insightful cell clusters, especially in complex settings.

\begin{backmatter}

\section*{Competing interests}
  The authors declare that they have no competing interests.

\section*{Author's contributions}
    Conceived the method: SZD, MZ, XQ. Developed the algorithm: SZD, MZ, XQ. Performed the simulations: SZD. Analyzed the results and wrote the paper: SZD, PdF, SHS, MZ, XQ. All authors read and approved the final manuscript.

\section*{Acknowledgements}
  This work was supported by the National Science Foundation (NSF) Grants 1553281 and 1812641. The authors also would like to thank Texas A\&M High Performance Research Computing
  and Texas Advanced Computing Center for providing computational resources to perform experiments in this paper.

\bibliographystyle{bmc-mathphys}