# Structured Sparse Non-negative Matrix Factorization with $\ell_{2,0}$-Norm
for scRNA-seq Data Analysis

Wenwen Min, Taosheng Xu, Xiang Wan and Tsung-Hui Chang  Wenwen Min is with The
Chinese University of Hong Kong, Shenzhen 518172, China, University of Science
and Technology of China, Hefei 230027, China, and the Shenzhen Research
Institute of Big Data, Shenzhen 518172, China. E-mail<EMAIL_ADDRESS>Taosheng Xu is with Warshel Institute for Computational Biology, The Chinese
University of Hong Kong, Shenzhen 518172, China, University of Science and
Technology of China, Hefei 230027, China. E-mail<EMAIL_ADDRESS>Xiang
Wan is with the Shenzhen Research Institute of Big Data, Shenzhen 518172,
China. E-mail<EMAIL_ADDRESS>Tsung-Hui Chang is with The Chinese
University of Hong Kong, Shenzhen 518172, China and the Shenzhen Research
Institute of Big Data, Shenzhen 518172, China. E-mail:
<EMAIL_ADDRESS>Manuscript received XX, 2021; revised XX, 2021.

###### Abstract

Non-negative matrix factorization (NMF) is a powerful tool for dimensionality
reduction and clustering. Unfortunately, the interpretation of the clustering
results from NMF is difficult, especially for the high-dimensional biological
data without effective feature selection. In this paper, we first introduce a
row-sparse NMF with $\ell_{2,0}$-norm constraint (NMF_$\ell_{20}$), where the
basis matrix $\bm{W}$ is constrained by the $\ell_{2,0}$-norm, such that
$\bm{W}$ has a row-sparsity pattern with feature selection. It is a challenge
to solve the model, because the $\ell_{2,0}$-norm is non-convex and non-
smooth. Fortunately, we prove that the $\ell_{2,0}$-norm satisfies the
Kurdyka-Łojasiewicz property. Based on the finding, we present a proximal
alternating linearized minimization algorithm and its monotone accelerated
version to solve the NMF_$\ell_{20}$ model. In addition, we also present a
orthogonal NMF with $\ell_{2,0}$-norm constraint (ONMF_$\ell_{20}$) to enhance
the clustering performance by using a non-negative orthogonal constraint. We
propose an efficient algorithm to solve ONMF_$\ell_{20}$ by transforming it
into a series of constrained and penalized matrix factorization problems. The
results on numerical and scRNA-seq datasets demonstrate the efficiency of our
methods in comparison with existing methods.

###### Index Terms:

$\ell_{2,0}$-norm, feature selection, row sparse NMF and ONMF, non-convex
optimization, scRNA-seq data clustering

##  I Introduction

With the development of single cell RNA sequencing (scRNA-seq) technology, we
can easily obtain biological profile data at single cell level from thousands
of cells at the same time [1]. Clustering such scRNA-seq data has been
becoming increasingly important for biological and medical applications [2].

Non-negative matrix factorization (NMF) and its variants have been widely used
to solve some computational biological problems [3, 4, 5, 6, 7, 8, 9].
Especially, they have achieved lots of successfully applications in scRNA-seq
data clustering analysis [10, 11, 12]. However, the NMF class algorithms
expose some shortcomings when being applied to the clustering analysis of the
high-dimensional biological data. It is well known that the high-dimensional
biological data contains many noisy and redundant features which often impacts
the performance of clustering algorithms.

To overcome the problems, sparse NMF methods have been proposed by adding
sparseness constraints [13, 14, 15]. At present, the proposed sparseness
constraints, such as $\ell_{1}$ and $\ell_{0}$ norms, cannot identify real
row-sparsity patterns of the basis matrix $\bm{W}$ in NMF,
$\bm{X}\approx\bm{W}\bm{H}$ (see Figure 1A and B), such that these sparse NMF
methods cannot select the important features for clustering analysis.

Feature selection is a effective way that extract the informative features to
improve the interpretability and performance of machine learning models [16,
17]. To enhance model interpretability, a common task is to search for some
most important information features when NMF and its variants are used for
high-dimensional data analysis. Previously, $\ell_{2,1}$-norm constraint has
been used in some supervised learning models to perform feature selection [18,
17]. Moreover, $\ell_{2,0}$-norm is more desirable from the sparsity
perspective, because it can select a specific number of the important
information features [19, 20, 21].

To integrate feature selection in the NMF model, we first present a row-sparse
NMF with $\ell_{2,0}$-norm constraint (NMF_$\ell_{20}$). The basis matrix
$\bm{W}$ is constrained by the $\ell_{2,0}$-norm, such that $\bm{W}$ has a
row-sparsity pattern with feature selection (see Figure 1C). However, it is
difficult to find an effective convergence algorithm to solve the
NMF_$\ell_{20}$ model because the $\ell_{2,0}$-norm constraint is non-convex
and non-smooth. Fortunately, we find that the $\ell_{2,0}$-norm satisfies the
Kurdyka-Łojasiewicz (KŁ) property such that the traditional proximal gradient
method can be used to solve a class of optimization problems with
$\ell_{2,0}$-norm constraint. For instance, the proximal alternating
linearized minimization (PALM) algorithm has been proposed to solve a class of
non-convex and non-smooth problems which satisfy KŁ property [22]. Based on
the above point, we introduce the PALM algorithm and its variant, a monotone
accelerated PALM (maPALM) algorithm, to solve the NMF_$\ell_{2,0}$ model. We
prove that both PALM and maPALM algorithms converge to a critical point when
they are used to solve the NMF_$\ell_{20}$ model.

In addition, we also note that the orthogonal NMF (ONMF) which is a variant of
NMF. It improves the clustering performance by adding the non-negative
orthogonal constraint [23, 24, 25]. Non-negative orthogonal matrix has the
following two properties: (1) an orthonormal matrix forms a basis for a
specific subspace, which facilitates geometric interpretation and signal
reconstruction; (2) two non-negative vectors in the matrix are orthogonal if
and only if their nonzero dimensions do not overlap. This may be the reason
why ONMF is sometime more effective than NMF in clustering.

To integrate feature selection and non-negative orthogonal constraint in the
NMF model, we also present a row-sparse ONMF with $\ell_{2,0}$-norm constraint
(ONMF_$\ell_{20}$). We propose an efficient algorithm for ONMF_$\ell_{20}$
model by using a penalty function method. The algorithm transforms
ONMF_$\ell_{20}$ into a series of subproblems so that the PALM and maPALM
algorithms can be used to solve them. Our contributions of this paper are
summarized as follows:

  1. 1.

We prove that the $\ell_{2,0}$-norm satisfies the KŁ property such that a
class of optimization problems with $\ell_{2,0}$-norm constraint can be solved
by the PALM algorithm.

  2. 2.

An efficient algorithm (PALM) for NMF_$\ell_{20}$, and its convergence
property.

  3. 3.

An accelerated version of PALM (maPALM) for NMF_$\ell_{20}$, and its
convergence property.

  4. 4.

An efficient algorithm for ONMF_$\ell_{20}$ by transforming it into a series
of subproblems, and its convergence property.

  5. 5.

The application of our methods and the comparison with the competing methods
using the simulated and scRNA-seq datasets. The results show that our methods
are more effective in clustering accuracy and feature selection.

Figure 1: Illustration of these basis matrices $\bm{W}$ obtained from three
structured sparse NMF models ($\bm{X}\approx\bm{W}\bm{H}$) with different
constraints. (A) showing estimated $\bm{W}$ by NMF_$\ell_{0}$, (B) showing
estimated $\bm{W}$ by NMF_$\ell_{20}$, and (C) showing estimated $\bm{W}$ by
NMF_$\ell_{c0}$, where dark boxes denote the non-zero coefficients and blank
boxes represent zero coefficients. The constraint conditions of $\bm{W}$ in
NMF_$\ell_{0}$, NMF_$\ell_{20}$ and NMF_$\ell_{c0}$ are defined in Table I.
Table I: Summary of notations.

Notation Meaning Normal font, e.g., x A scalar Bold lowercase, e.g., $\bm{w}$
A vector Bold capital, e.g., $\bm{W}$ A matrix $\bm{X}$ A $p$-by-$n$ matrix
$\bm{W}$ A $p$-by-$r$ matrix $\bm{H}$ An $r$-by-$n$ matrix $\bm{I}_{r}$ An
$r\times r$ identity matrix $\bm{h}^{i}$ The $i$-th row of $\bm{H}$
$\bm{h}_{j}$ The $j$-th column of $\bm{H}$ $\|\cdot\|_{1}$ $\ell_{1}$-norm for
a vector $\|\cdot\|_{0}$ $\ell_{0}$-norm for a vector $\|\cdot\|_{2}$ or
$\|\cdot\|$ $\ell_{2}$-norm for a vector $\|\cdot\|_{F}$ Frobenius norm for a
matrix $\|\cdot\|_{2}$ Spectral norm for a matrix $I()$ Indicator function
Model Constraint condition NMF_$\ell_{0}$ $\|\bm{W}\|_{0}\leq k\cdot r$ (see
Eq. 3) NMF_$\ell_{20}$ $\|\bm{W}\|_{2,0}\leq k$ (see Eq. 4) NMF_$\ell_{c0}$
Column-wise sparsity of $\bm{W}$ (see Eq. 8)

##  II Notations and definitions

Given a matrix $\bm{W}\in\mathbb{R}^{p\times r}$, let $\bm{w}^{i}$ and
$\bm{w}_{j}$ denote its $i$-th row and $j$-th column, respectively. The
Frobenius norm of $\bm{W}\in R^{p\times n}$ is defined as:

$\|\bm{W}\|_{F}=\sqrt{\sum_{i=1}^{p}\sum_{j=1}^{r}{w_{ij}^{2}}}=\sqrt{tr(\bm{W}^{T}\bm{W})}.$
(1)

The spectral norm of $\bm{W}$ is the largest singular value of $\bm{W}$ and it
is defined as:

$\|\bm{W}\|_{2}=\sigma_{\max}(\bm{W}),$ (2)

The $\ell_{0}$-norm of $\bm{W}$ is defined as:

$\|\bm{W}\|_{0}=\sum_{i=1}^{p}\sum_{j=1}^{n}{I(w_{ij}\neq
0)}=\sum_{i=1}^{p}\|\bm{w}^{i}\|_{0}=\sum_{j=1}^{r}\|\bm{w}_{j}\|_{0},$ (3)

where $I(x)=1$ if $x\neq 0$, $I(x)=0$ if $x=0$. The $\ell_{2,0}$-norm of
$\bm{W}$ is defined as:

$\|\bm{W}\|_{2,0}=\sum_{i=1}^{p}{I(\|\bm{w}^{i}\|\neq
0)}=\|(\|\bm{w}^{1}\|,\cdots,\|\bm{w}^{p}\|)\|_{0},$ (4)

where $\|\cdot\|$ is the $\ell_{2}$-norm and $\|\bm{x}\|=\sum{x_{i}^{2}}$.
Briefly, $\|\bm{W}\|_{2,0}$ denotes the number of non-zero rows in $\bm{W}$.
More notations are summarized in Table I and the mathematical definitions for
nonconvex optimization are summarized into the appendix A-A. For simplicity,
$\ell_{2,0}$ and $\|\bm{W}\|_{2,0}$ sometimes are abbreviated as $\ell_{20}$
and $\|\bm{W}\|_{20}$, respectively.

##  III Proposed framework

###  III-A NMF and ONMF

Given a data $\bm{X}\in\mathbb{R}^{p\times n}$ with $p$ features and $n$
samples, NMF model [26] can be written as follows:

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$
$\displaystyle\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}$ (5) subject to
$\displaystyle\bm{W}\in\mathbb{R}_{+}^{p\times
r},\bm{H}\in\mathbb{R}_{+}^{r\times n}.$

[23] has reported that ONMF can improve the clustering performance by adding
the non-negative orthogonal constraint on $\bm{H}$ . Thus, we introduce the
following ONMF model:

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$
$\displaystyle\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}$ (6) subject to
$\displaystyle\bm{W}\in\mathbb{R}_{+}^{p\times
r},\bm{H}\in\mathbb{R}_{+}^{r\times n},\bm{H}\bm{H}^{T}=\bm{I}_{r}.$

Theorem 1 in reference [23] has shown that ONMF is equivalent to k-means
clustering, because the nonnegative orthogonal matrix $\bm{H}$ has a good
property (see remark 1).

######  Remark 1.

For the solution $\bm{H}$ of Eq. (6), it has at most one non-zero entry in
each column, because $\bm{H}$ is non-negative and satisfies orthogonality
$\bm{H}\bm{H}^{T}=\bm{I}_{r}$.

###  III-B Structured sparse NMF (SSNMF)

To integrate feature selection and non-negative orthogonal constraint in NMF
model, we introduce a row-sparse ONMF with $\ell_{2,0}$-norm constraint
(ONMF_$\ell_{20}$):

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$
$\displaystyle\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}$ (7) subject to
$\displaystyle\bm{W}\in\Omega_{w}:=\\{\bm{W}\in\mathbb{R}_{+}^{p\times
r}:\|\bm{W}\|_{2,0}\leq k\\},$
$\displaystyle\bm{H}\in\Omega_{h}:=\\{\bm{H}\in\mathbb{R}_{+}^{r\times
n}:\bm{H}\bm{H}^{T}=\bm{I}_{r}\\},$

where $\|\bm{W}\|_{2,0}\leq k$ encourages $\bm{W}$ to be row sparse and select
some most important features. We can also make $\bm{W}$ column sparse using
the following constraint, named $\ell_{c0}$-norm, which uses $\ell_{0}$-norm
to each column of $\bm{W}$:

$\Omega_{w}:=\\{\bm{W}\in\mathbb{R}_{+}^{p\times r}:\|\bm{w}_{j}\|_{0}\leq
s,\forall j\\}.$ (8)

The key to solve problem (7) is to remove its orthogonal constraint. Based on
the conclusion of Eq. (8) from reference [25], Eq. (7) is equivalent to:

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$
$\displaystyle\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}$ (9) subject to
$\displaystyle\bm{W}\in\Omega_{w}:=\\{\bm{W}\in\mathbb{R}_{+}^{p\times
r}:\|\bm{W}\|_{2,0}\leq k\\},$
$\displaystyle\bm{H}\in\\{\bm{H}\in\mathbb{R}_{+}^{r\times
n}:\|\bm{h}_{j}\|_{1}^{2}=\|\bm{h}_{j}\|_{2}^{2},\forall j\\}.$

We consider its penalized formulation and present a SSNMF framework as
follows:

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$
$\displaystyle\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}+\frac{\rho}{2}\sum_{j=1}^{n}\Big{(}(\bm{1}^{T}\bm{h}_{j})^{2}-\|\bm{h}_{j}\|_{2}^{2}\Big{)}$
(10) subject to
$\displaystyle\bm{W}\in\Omega_{w}:=\\{\bm{W}\in\mathbb{R}_{+}^{p\times
r}:\|\bm{W}\|_{2,0}\leq k\\},$
$\displaystyle\bm{H}\in\Omega_{h}:=\mathbb{R}_{+}^{r\times n}.$

Based on the SSNMF framework (LABEL:equ-10), we introduce the following four
SSNMF models as follows:

  * •

Row-sparse NMF with $\ell_{2,0}$-norm constraint (NMF_$\ell_{20}$). When
$\rho=0$, the penalty term has no effect and Eq. (LABEL:equ-10) reduces to
NMF_$\ell_{20}$.

  * •

Row-sparse ONMF with $\ell_{2,0}$-norm constraint (ONMF_$\ell_{20}$). When
$\rho$ is large enough, Eq. (LABEL:equ-10) reduces to ONMF_$\ell_{20}$. As
$\rho\rightarrow\infty$ in Eq. (LABEL:equ-10), the impact of the penalty
grows, and the estimated $\bm{H}$ will approach a nonnegative orthogonal
matrix.

  * •

Column-sparse NMF with $\ell_{c,0}$-norm constraint (NMF_$\ell_{c0}$). When
$\rho=0$ and $\Omega_{w}$ is defined in Eq. (8), Eq. (LABEL:equ-10) reduces to
NMF_$\ell_{c0}$.

  * •

Column-sparse ONMF with $\ell_{c,0}$-norm constraint (ONMF_$\ell_{c0}$). When
$\rho$ is large enough and $\Omega_{w}$ is defined in Eq. (8), Eq.
(LABEL:equ-10) reduces to ONMF_$\ell_{c0}$.

We first present Proposition 1 to clarify the relationship between the
solutions of Eq. (LABEL:equ-10) and Eq. (7), whose proof is shown in the
appendix A-B. Proposition 1 implies that we can solve Eq. (7) by repeatedly
solving Eq. (LABEL:equ-10) with a gradually increasing $\rho$. Therefore, the
key to solve the above four models (NMF_$\ell_{20}$, NMF_$\ell_{c0}$,
ONMF_$\ell_{20}$, and ONMF_$\ell_{c0}$) is to solve the SSNMF framework
(LABEL:equ-10).

######  Proposition 1.

Let $(\bm{W}^{*},\bm{H}^{*})$ be a local minimizer of Eq. (LABEL:equ-10). Then
for any $\rho>0$, $(\bm{W}^{*},\bm{H}^{*})$ is also the feasible and local
minimizer of Eq. (7).

Alternating minimization is a popular strategy to solve the constrained and
penalized matrix factorization problem in Eq. (LABEL:equ-10). Recently, a PALM
algorithm has been proposed to solve a class of constrained and penalized
matrix factorization problems which satisfies KŁ property [22]. We present
Theorem 1 to show that Eq. (LABEL:equ-10) satisfies KŁ property, such that the
PALM algorithm can be used to solve it.

######  Theorem 1.

$\Phi=F(\bm{W},\bm{H})+\delta_{\bm{W}\geq 0}+\delta_{\bm{H}\geq
0}+\delta_{\|\bm{W}\|_{2,0}\leq k}$ is a semi-algebraic function and it
satisfies the KŁ property, where
$F=\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}+\frac{\rho}{2}\sum_{j=1}^{n}\Big{(}(\bm{1}^{T}\bm{h}_{j})^{2}-\|\bm{h}_{j}\|_{2}^{2}\Big{)}$,
and $\delta_{\|\bm{W}\|_{2,0}\leq k}$ is zero if $\|\bm{W}\|_{2,0}\leq k$,
otherwise $+\infty$.

###### Proof.

Remark 8 in [22] shows the following conclusions: (1) $\|\cdot\|_{0}$ and
$\|\cdot\|_{2}$ are semi-algebraic functions; (2) The indicator function in a
semi-algebraic set is semi-algebraic; (3) Any composition of semi-algebraic
function remains to be semi-algebraic; (4) The finite sums of semi-algebraic
functions remain semi-algebraic. Based the above conclusions, we know that (1)
$F$, $\delta_{\bm{W}\geq 0}$ and $\delta_{\bm{H}\geq 0}$ are semi-algebraic
functions, respectively; (2) $\delta_{\|\bm{W}\|_{2,0}\leq k}$ is a semi-
algebraic function because
$\|\bm{W}\|_{2,0}:=\|(\|\bm{w}^{1}\|_{2},\cdots,\|\bm{w}^{p}\|_{2})\|_{0}$ is
a composition of semi-algebraic function $\|\cdot\|_{0}$ and $\|\cdot\|_{2}$.
Thus, we prove that $\Phi$ is a semi-algebraic function because it is a sum of
four semi-algebraic functions. In addition, we observe that $\Phi$ is a proper
and lower semicontinuous function. Based on the Theorem 3 in [22], a proper,
lower semicontinuous and semi-algebraic function $\Phi$ satisfies the KŁ
property. The relevant mathematical definitions are shown in Appendix A-A. ∎

By the way, the proof of Theorem 1 also implies that the $\ell_{2,0}$-norm
satisfies the KŁ property such that the PALM can be used to solve a class of
optimization problems with $\ell_{2,0}$-norm constraint.

##  IV Optimization method

###  IV-A PALM and maPALM

We first introduce a general constrained and penalized matrix factorization
model as follows:

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$ $\displaystyle
F(\bm{W},\bm{H})=\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}+\phi(\bm{W})+\varphi(\bm{H})$
(11) subject to $\displaystyle\bm{W}\in\Omega_{w},\bm{H}\in\Omega_{h}.$

Obviously, the SSNMF framework in Eq. (LABEL:equ-10) is a special case of Eq.
(LABEL:equ:11). To solve Eq. (LABEL:equ:11) using the PALM algorithm, we need
to perform a projected gradient descent step with respect to $\bm{H}$ and
$\bm{W}$ for $t=1,2,\cdots$:

$\displaystyle\bm{H}^{t+1}$
$\displaystyle=\mathcal{P}_{H\in\Omega_{H}}\Big{\\{}\bm{H}^{t}-\frac{1}{d_{H}^{t}}\nabla_{H}F(\bm{W}^{t},\bm{H}^{t})\Big{\\}},$
$\displaystyle\bm{W}^{t+1}$
$\displaystyle=\mathcal{P}_{W\in\Omega_{W}}\Big{\\{}\bm{W}^{t}-\frac{1}{d_{W}^{t}}\nabla_{W}F(\bm{W}^{t},\bm{H}^{t+1})\Big{\\}},$

where $d_{H}^{t}$ and $d_{W}^{t}$ are two step-size parameters,
$\mathcal{P}_{H\in\Omega_{H}}\\{\cdot\\}$ and
$\mathcal{P}_{W\in\Omega_{W}}\\{\cdot\\}$ are two projection operations onto
$\Omega_{H}$ and $\Omega_{W}$, respectively. A basic algorithmic framework for
solving Eq. (LABEL:equ:11) is shown in Algorithm 1.

Because of the linearization of PALM, its convergence speed may be slow.
Accelerated proximal gradient method has been widely used to solve convex
optimization problems. Unfortunately, if the non-monotone accelerated proximal
gradient uses a bad extrapolation for some non-convex problems, then it may
not converge to a critical point [27]. Fortunately, the monotone accelerated
method guarantees convergence for a non-convex problem by ensuring that its
objective function value decreases every iteration [28, 27].

Algorithm 1 PALM for solving Eq. (LABEL:equ:11).

0: $\bm{X}\in\mathbb{R}^{p\times n}$ and $\epsilon>0$.

0: $\bm{W}\in\mathbb{R}^{p\times r}$ and $\bm{H}\in\mathbb{R}^{r\times n}$.

1: Initialize $(\bm{W}^{0},\bm{H}^{0})$ and set $t=0$

2: repeat

3: Compute $d_{H}^{t}$

4:
$\bm{H}^{t+1}=\mathcal{P}_{H\in\Omega_{H}}\Big{\\{}\bm{H}^{t}-\frac{1}{d_{H}^{t}}\nabla_{H}F(\bm{W}^{t},\bm{H}^{t})\Big{\\}}$

5: Compute $d_{W}^{t}$

6:
$\bm{W}^{t+1}=\mathcal{P}_{W\in\Omega_{W}}\Big{\\{}\bm{W}^{t}-\frac{1}{d_{W}^{t}}\nabla_{W}F(\bm{W}^{t},\bm{H}^{t+1})\Big{\\}}$

7: $t=t+1$

8: until
$\frac{\|(\bm{W}^{t},\bm{H}^{t})-(\bm{W}^{t-1},\bm{H}^{t-1})\|}{\|(\bm{W}^{t-1},\bm{H}^{t-1})\|}<\epsilon$

9: return $\bm{W}:=\bm{W}^{t}$ and $\bm{H}:=\bm{H}^{t}$.

Algorithm 2 maPALM for solving Eq. (LABEL:equ:11).

0: $\bm{X}\in\mathbb{R}^{p\times n}$, $\epsilon>0$ and $\omega_{0}$.

0: $\bm{W}\in\mathbb{R}^{p\times r}$ and $\bm{H}\in\mathbb{R}^{r\times n}$.

1: Initialize $(\bm{W}^{-1},\bm{H}^{-1})=(\bm{W}^{0},\bm{H}^{0})$ and $t=0$

2: repeat

3: $\widetilde{\bm{H}}^{t}=\bm{H}^{t}+\omega_{t}(\bm{H}^{t}-\bm{H}^{t-1})$

4: Compute $d_{H}^{t}$

5:
$\widetilde{\bm{H}}^{t+1}=\mathcal{P}_{H\in\Omega_{H}}\Big{\\{}\widetilde{\bm{H}}^{t}-\frac{1}{d_{H}^{t}}\nabla_{H}F(\bm{W}^{t},\widetilde{\bm{H}}^{t})\Big{\\}}$

6: $\widetilde{\bm{W}}^{t}=\bm{W}^{t}+\omega_{t}(\bm{W}^{t}-\bm{W}^{t-1})$

7: Compute $d_{W}^{t}$

8:
$\widetilde{\bm{W}}^{t+1}=\mathcal{P}_{W\in\Omega_{W}}\Big{\\{}\widetilde{\bm{W}}^{t}-\frac{1}{d_{W}^{t}}\nabla_{W}F(\widetilde{\bm{W}}^{t},\widetilde{\bm{H}}^{t+1})\Big{\\}}$

9: if $F(\widetilde{\bm{W}}^{t+1},\widetilde{\bm{H}}^{t+1})\leq
F(\bm{W}^{t},\bm{H}^{t})$ then

10:
$(\bm{W}^{t+1},\bm{H}^{t+1}):=(\widetilde{\bm{W}}^{t+1},\widetilde{\bm{H}}^{t+1})$

11: else

12: update $(\bm{W}^{t+1},\bm{H}^{t+1})$ using Eq. (12)

13: end if

14: Compute $\omega_{t}$ and set $t=t+1$

15: until
$\frac{\|(\bm{W}^{t},\bm{H}^{t})-(\bm{W}^{t-1},\bm{H}^{t-1})\|}{\|(\bm{W}^{t-1},\bm{H}^{t-1})\|}<\epsilon$

16: return $\bm{W}:=\bm{W}^{t}$ and $\bm{H}:=\bm{H}^{t}$.

To this end, we develop a monotone accelerated PALM framework and its details
is shown in Algorithm 2, which can be regarded as a special case of the block
prox-linear method [29]. Similarly, maPALM ensures the objective function
value decreases. Reference [27, 29], $\omega_{k}$ ($k=0,1,2,\cdots$) in
Algorithm 2 is dynamically updated by

$\omega_{k}=\frac{\tau_{k}-1}{\tau_{k+1}},$ (13)

where $\tau_{0}=1$ and $\tau_{k+1}=\Big{(}1+\sqrt{1+4\tau_{k}^{2}}\Big{)}/2$.

We note that maPALM reduces to PALM when $\omega_{t}=0$ for all $t$.
Especially, NMF_$\ell_{20}$ and NMF_$\ell_{c0}$ can be solved by Algorithm 1
and 2 with $\phi(\bm{W})=\varphi(\bm{H})=0$. In addition, the conclusion of
Proposition 1 implies that we can solve Eq. (7) by repeatedly solving Eq.
(LABEL:equ-10) with a gradually increasing $\rho$. To this end, we can develop
an efficient algorithm for ONMF_$\ell_{20}$ or ONMF_$\ell_{c0}$ by turning
ONMF_$\ell_{20}$ or ONMF_$\ell_{c0}$ into a sequence subproblems.

In summary, the key to solving NMF_$\ell_{20}$, NMF_$\ell_{c0}$,
ONMF_$\ell_{20}$ and ONMF_$\ell_{c0}$ is to solve the SSNMF framework in Eq.
(LABEL:equ-10). Below we show the details of using the maPALM algorithm to
solve it.

###  IV-B Solve SSNMF with $\ell_{2,0}$-norm constraint

We use maPALM to solve a row-sparse SSNMF in Eq. (LABEL:equ-10) with
$\ell_{2,0}$-norm constraint:

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$
$\displaystyle\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}+\frac{\rho}{2}\sum_{j=1}^{n}\Big{(}(\bm{1}^{T}\bm{h}_{j})^{2}-\|\bm{h}_{j}\|_{2}^{2}\Big{)}$
(14) subject to $\displaystyle\bm{W}\in\mathbb{R}_{+}^{p\times
r},\|\bm{W}\|_{2,0}\leq k,\bm{H}\in\mathbb{R}_{+}^{r\times n}.$

Let
$F=\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}+\frac{\rho}{2}\sum_{j=1}^{n}\Big{(}(\bm{1}^{T}\bm{h}_{j})^{2}-\|\bm{h}_{j}\|_{2}^{2}\Big{)}$,
then

$\displaystyle\nabla_{H}F=\bm{W}^{T}\bm{W}\bm{H}-\bm{W}^{T}\bm{X}+\rho\bm{1}_{r\times
r}\bm{H}-\rho\bm{H},$ (15a)
$\displaystyle\nabla_{W}F=\bm{W}\bm{H}\bm{H}^{T}-\bm{X}\bm{H}^{T}.$ (15b)

And the Hessian matrices of $F$ with respect to $\bm{W}$ and $\bm{H}$ are

$\displaystyle\nabla_{W}^{2}F=(\bm{H}\bm{H}^{T})\otimes\bm{I}_{p},$ (16a)
$\displaystyle\nabla_{H}^{2}F=\bm{I}_{p}\otimes(\bm{W}^{T}\bm{W}+\rho\bm{1}_{r\times
r}-\rho\bm{I}_{r}).$ (16b)

where $\otimes$ is Kronecker product and $\bm{I}_{p}\in\mathbb{R}^{p\times p}$
is an identity matrix. To use maPALM to solve Eq. (14), we need to calculate
Lipschitz constant to determine the step size. The Lemma 2 in reference [30]
shows that $\nabla_{W}F$ and $\nabla_{H}F$ are Lipschitz continuous, the
Lipschitz constant of $\nabla_{W}F$ is the largest singular value of
$\nabla_{W}^{2}F$, _i.e._ , $L_{W}=\|\bm{H}\bm{H}^{T}\|_{2}$, and the
Lipschitz constant of $\nabla_{H}F$ is the largest singular value of
$\nabla_{W}^{2}F$, _i.e._ , $L_{H}=\|\bm{W}^{T}\bm{W}+\rho\bm{1}_{r\times
r}-\rho\bm{I}_{r}\|_{2}$. Thus, we can set $d_{W}^{t}=L_{W}$ and
$d_{H}^{t}=L_{H}$ in Algorithms 1 to 2.

1) Optimize W. Specifically, to obtain the update of $\bm{W}$ for Eq. (14), we
need to solve a proximal map
$\bm{W}^{t+1}:=\mathcal{P}_{W\in\Omega_{W}}\\{\overline{\bm{W}}\\}$ as
follows:

$\displaystyle\underset{\bm{W}}{\text{minimize}}$
$\displaystyle\|\bm{W}-\overline{\bm{W}}\|_{F}^{2}$ (17) subject to
$\displaystyle\bm{W}\in\mathbb{R}_{+}^{p\times r},\|\bm{W}\|_{2,0}\leq k,$

where
$\overline{\bm{W}}:=\bm{W}^{t}-\frac{1}{d_{W}^{t}}\nabla_{W}F(\bm{W}^{t},\bm{H}^{t})$.
We propose Proposition 2 to solve the proximal map. To this end, we introduce
the following mathematical definitions.

######  Definition 1.

$\mathrm{SupportNorm}(\bm{W},k)$ is a set of indices of $\bm{z}$ with the
largest $k$ values where $\bm{z}=(\|\bm{w}^{1}\|,\cdots,\|\bm{w}^{p}\|)$ and
$\bm{w}^{i}$ denotes $i$-th row of $\bm{W}$.

######  Definition 2.

For a given matrix $\bm{W}\in\mathbb{R}^{p\times r}$,
$\mathrm{RS}_{k}(\bm{W})$ is also a $p\times r$ matrix which is defined as:

$[\mathrm{RS}_{k}(\bm{W})]_{ij}=\begin{cases}W_{ij},&\text{if
}i\in\mathrm{SupportNorm}(\bm{W},k),\\\ 0,&\text{otherwise},\end{cases}$ (18)

where $\mathrm{SupportNorm}(\bm{W},k)$ is defined in Definition 1.
$\mathrm{RS}_{k}(\bm{W})$ only keeps $k$ non-zero rows with the largest
$\ell_{2}$-norm values in $\bm{W}$.

######  Definition 3.

For a given matrix $\bm{W}\in\mathbb{R}^{p\times r}$, $P_{+}(\bm{W})$ is
defined as follows:

$[P_{+}(\bm{W})]_{ij}=\begin{cases}W_{ij},&\text{if}\leavevmode\nobreak\
W_{ij}>0\leavevmode\nobreak\ \forall i\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ j,\\\ 0,&\text{otherwise}.\end{cases}$ (19)

######  Proposition 2.

(Proximal map formula of $\bm{W}$ for ONMF_$\ell_{20}$) Let
$\overline{\bm{W}}\in\mathbb{R}^{p\times r}$, then Eq. (17) has a closed-form
solution

$\mathcal{P}_{W\in\Omega_{W}}\\{\overline{\bm{W}}\\}:=\mathrm{RS}_{k}(P_{+}(\overline{\bm{W}})),$
(20)

where $\mathrm{RS}_{k}(\cdot)$ and $P_{+}(\cdot)$ are defined in Definition 2
and 3, respectively.

###### Proof.

Suppose that the optimal solution of Eq. (17) is $\widehat{\bm{W}}$, then we
can easily observe that $\widehat{W}_{ij}$ must be zero when
$\overline{W}_{ij}<0$. So, Eq. (17) is equivalent to

$\underset{\bm{W}}{\text{minimize}}\leavevmode\nobreak\
\|\bm{W}-P_{+}(\overline{\bm{W}})\|_{F}^{2},\leavevmode\nobreak\ \text{subject
to}\leavevmode\nobreak\ \|\bm{W}\|_{2,0}\leq k.$

Due to the constraint $\|\bm{W}\|_{2,0}\leq k$, the optimal solution only
keeps up to $k$ non-zero rows. So, Eq. (17) has a closed-form solution
$\mathrm{RS}_{k}(P_{+}(\overline{\bm{W}}))$. ∎

2) Optimize H. To obtain the update of $\bm{H}$ for Eq. (14) , we need to
solve a proximal map
$\bm{H}^{t+1}:=\mathcal{P}_{H\in\Omega_{H}}\\{\overline{\bm{H}}\\}$ as
follows:

$\underset{\bm{H}}{\text{minimize}}\leavevmode\nobreak\
\|\bm{H}-\overline{\bm{H}}\|_{F}^{2},\leavevmode\nobreak\ \text{subject
to}\leavevmode\nobreak\ \bm{H}\in\mathbb{R}_{+}^{r\times n},$ (21)

where
$\overline{\bm{H}}:=\bm{H}^{t}-\frac{1}{d_{H}^{t}}\nabla_{H}F(\bm{W}^{t+1},\bm{H}^{t})$.
We propose the following Proposition 3 to solve the above problem.

######  Proposition 3.

(Proximal map formula of $\bm{H}$ for ONMF_$\ell_{20}$) Let
$\overline{\bm{H}}\in\mathbb{R}^{r\times n}$, then (21) has a closed-form
solution:

$\mathcal{P}_{H\in\Omega_{H}}\\{\overline{\bm{H}}\\}:=P_{+}(\overline{\bm{H}}),$
(22)

where $P_{+}(\cdot)$ is defined in Definition 3.

3) Algorithm for solving Eq. (14). Based on the Propositions 2 and 3, we
develop a maPALM algorithm to solve Eq. (14) and the detailed algorithm is
given in Algorithm 3. To maintain monotonicity, maPALM needs to use a suitable
$\omega_{t}$ by checking the objective function value. If the objective
function value becomes larger, then we obtain the update rule based on the
traditional projected gradient descent method. Monotonicity can ensure that
Algorithm 3 converges to a critical point for any initial point. The following
Theorem 2 gives the details on the convergence analysis of Algorithm 3.

Algorithm 3 maPALM for solving Eq. (14).

0: $\bm{X}\in\mathbb{R}^{p\times n}$, $\rho\geq 0$, $\epsilon>0$, $k$ (nonzero
rows).

0: $\bm{W}\in\mathbb{R}^{p\times r}$ and $\bm{H}\in\mathbb{R}^{r\times n}$.

1: Initialize $(\bm{W}^{-1},\bm{H}^{-1})=(\bm{W}^{0},\bm{H}^{0})$,
$\tau_{0}=1$ and $t=0$

2: repeat

3: Compute $\omega_{t}=\frac{\tau_{t}-1}{\tau_{t+1}}$ where
$\tau_{t+1}=\frac{1+\sqrt{1+4\tau_{t}^{2}}}{2}$

4: $\widetilde{\bm{H}}^{t}=\bm{H}^{t}+\omega_{t}(\bm{H}^{t}-\bm{H}^{t-1})$

5: $d_{H}^{t}=\|(\bm{W}^{t})^{T}(\bm{W}^{t})+\rho\bm{1}_{r\times
r}-\rho\bm{I}_{r}\|_{2}$

6:
$\bm{H}^{t+1}=P_{+}\big{(}\widetilde{\bm{H}}^{t}-\frac{1}{d_{H}^{t}}\nabla_{H}F(\bm{W}^{t},\widetilde{\bm{H}}^{t})\big{)}$

7: $\widetilde{\bm{W}}^{t}=\bm{W}^{t}+\omega_{t}(\bm{W}^{t}-\bm{W}^{t-1})$

8: $d_{W}^{t}=\|(\bm{H}^{t+1})(\bm{H}^{t+1})^{T}\|_{2}$

9:
$\bm{W}^{t+1}=\mathrm{RS}_{k}\big{(}\widetilde{\bm{W}}^{t}-\frac{1}{d_{W}^{t}}\nabla_{W}F(\widetilde{\bm{W}}^{t},\bm{H}^{t+1})\big{)}$

10: if $F(\widetilde{\bm{W}}^{t+1},\widetilde{\bm{H}}^{t+1})\leq
F(\bm{W}^{t},\bm{H}^{t})$ then

11:
$(\bm{W}^{t+1},\bm{H}^{t+1}):=(\widetilde{\bm{W}}^{t+1},\widetilde{\bm{H}}^{t+1})$

12: else

13: $d_{H}^{t}=\|(\bm{W}^{t})^{T}(\bm{W}^{t})+\rho\bm{1}_{r\times
r}-\rho\bm{I}_{r}\|_{2}$

14:
$\bm{H}^{t+1}=P_{+}\big{(}\bm{H}^{t}-\frac{1}{d_{H}^{t}}\nabla_{H}F(\bm{W}^{t},\bm{H}^{t})\big{)}$

15: $d_{W}^{t}=\|(\bm{H}^{t+1})(\bm{H}^{t+1})^{T}\|_{2}$

16:
$\bm{W}^{t+1}=\mathrm{RS}_{k}\big{(}\bm{W}^{t}-\frac{1}{d_{W}^{t}}\nabla_{W}F(\bm{W}^{t},\bm{H}^{t+1})\big{)}$

17: end if

18: $t=t+1$

19: until
$\frac{\|(\bm{W}^{t},\bm{H}^{t})-(\bm{W}^{t-1},\bm{H}^{t-1})\|}{\|(\bm{W}^{t-1},\bm{H}^{t-1})\|}<\epsilon$

20: return $\bm{W}:=\bm{W}^{t}$ and $\bm{H}:=\bm{H}^{t}$.

Algorithm for NMF_$\ell_{20}$. Algorithm 3 with $\rho=0$ can effectively solve
NMF_$\ell_{20}$.

Initialization. We can adopt two ways to generate the initial point of
Algorithm 3. One is to use the random vectors from a standard normal
distribution to initialize $\bm{W}$ and $\bm{H}$. The other is to use the
solution $\bm{W}$ and $\bm{H}$ derived by the traditional NMF to initialize
them. The second way is a good guess. So, if not specified, we use the second
way to initialize by default.

Step-size. We can use a fixed value of the step-size $1/d_{H}^{t}$ and
$1/d_{W}^{t}$ in Algorithm 3, and also try to perform an approximate
backtracking line search from $t\in(0,0.5)$. Specifically, we use the fixed
$d_{W}^{t}=L_{W}$ and $d_{H}^{t}=L_{H}$ in this paper.

4) Computation cost. The computational complexity of Algorithm 3 depends on
the number of iterations. At each iteration, only two simple closed-form
solutions need to be computed with respect to $\bm{H}$ and $\bm{W}$ in the
steps 6 and 9, respectively. For each update of $\bm{H}$, the most costly step
is the calculation of $\nabla_{H}F$, which requires a computation cost of
$O(pr^{2}+nr^{2})$. For each update of $\bm{W}$, the most costly steps is the
calculation of $\nabla_{W}F$, which requires a computation cost of
$O(pr^{2}+nr^{2}+pnr)$. In addition, the calculation of objective function
value requires a computation cost of $O(pnr+pn+rn)$. Thus, each iteration of
Algorithm 3 requires a computation cost of $O(pr^{2}+nr^{2}+pnr)$.

5) Convergence analysis. Based on some results from references [29, 22], we
propose the following Theorem 2 to show that Algorithm 3 globally converges to
a critical point.

######  Theorem 2.

(Global convergence of Algorithm 3) Let
$\\{(\bm{W}^{(i)},\bm{H}^{(i)})\\}_{i=1}^{\infty}$ be a sequence generated
from any starting point $(\bm{W}^{(0)},\bm{H}^{(0)})$ by Algorithm 3. If
$\\{(\bm{W}^{(i)},\bm{H}^{(i)})\\}$ are bounded, then the objective is non-
increasing and the sequence has finite length and converges to a critical
point.

###### Proof.

Based on Theorem 1, the objective function $\Phi$ is a semi-algebraic function
and satisfies the KL property. Checking the assumptions of Theorem 2 in
reference [29], we observe that all assumptions required in Algorithm 3 are
clearly satisfied. So, we have Theorem 2. ∎

###  IV-C Solve SSNMF with $\ell_{c,0}$-norm constraint

Let $\Omega_{w}:=\\{\bm{W}\in\mathbb{R}_{+}^{p\times r}:\|\bm{w}_{j}\|_{0}\leq
k,\forall j\\}$ in Eq. (LABEL:equ-10), we consider a column-wise SSNMF with
$\ell_{c,0}$-norm constraint by using $\ell_{0}$-norm for each column of
$\bm{W}$ as follows:

$\displaystyle\underset{\bm{W},\bm{H}}{\text{minimize}}$
$\displaystyle\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}+\frac{\rho}{2}\sum_{j=1}^{n}\Big{(}(\bm{1}^{T}\bm{h}_{j})^{2}-\|\bm{h}_{j}\|_{2}^{2}\Big{)}$
(23) subject to $\displaystyle\bm{W}\in\mathbb{R}_{+}^{p\times
r},\|\bm{w}_{j}\|_{0}\leq k,\forall j,\bm{H}\in\mathbb{R}_{+}^{r\times n}.$

To obtain the update of $\bm{W}$, we need to solve a proximal map
$\bm{W}^{t+1}:=\mathcal{P}_{W\in\Omega_{W}}\\{\overline{\bm{W}}\\}$:

$\displaystyle\underset{\bm{W}}{\text{minimize}}$
$\displaystyle\|\bm{W}-\overline{\bm{W}}\|_{F}^{2}$ (24) subject to
$\displaystyle\bm{W}\in\mathbb{R}_{+}^{p\times r},\|\bm{w}_{j}\|_{0}\leq
k,\forall j.$

######  Definition 4.

For a given matrix $\bm{W}\in\mathbb{R}^{p\times r}$,
$\mathrm{CS}_{k}(\bm{W})\in\mathbb{R}^{p\times r}$ ($k\leq p$) is defined as:

$[\mathrm{CS}_{k}(\bm{W})]_{ij}=\begin{cases}W_{ij},&j\in\mathrm{SupportNorm}(\bm{w}_{j},k),\\\
0,&\text{otherwise},\end{cases}$ (25)

where $\bm{w}_{j}$ denotes $j$-th column of $\bm{W}$ and
$\mathrm{SupportNorm}(\cdot,k)$ is defined in Definition 1.

######  Proposition 4.

(Proximal map formula of $\bm{W}$ for ONMF_$\ell_{c0}$) Let
$\overline{\bm{W}}\in\mathbb{R}^{p\times r}$, then Eq. (24) has a closed-form
solution:

$\mathcal{P}_{W\in\Omega_{W}}\\{\overline{\bm{W}}\\}:=\mathrm{CS}_{k}(P_{+}(\bm{W})),$
(26)

where $\mathrm{CS}_{k}(\cdot)$ and $P_{+}(\cdot)$ are defined in Definition 4
and 3, respectively.

###### Proof.

Suppose that the optimal solution of Eq. (24) is $\widehat{\bm{W}}$, then we
can observe that $\widehat{W}_{ij}$ must be zero when $\overline{W}_{ij}<0$.
So, Eq. (24) is equivalent to

$\underset{\bm{W}}{\text{minimize}}\leavevmode\nobreak\
\|\bm{W}-P_{+}(\overline{\bm{W}})\|_{F}^{2}\leavevmode\nobreak\ \text{subject
to}\leavevmode\nobreak\ \|\bm{w}_{j}\|_{0}\leq k,\forall j.$

Due to the constraint $\|\bm{w}_{j}\|_{0}\leq k$, the optimal solution only
keeps $k$ non-zero elements with the largest absolute values for each column
of $\bm{W}$. So, Eq. (24) has a closed-form solution
$\mathrm{CS}_{k}(P_{+}(\overline{\bm{W}}))$. ∎

Algorithm for solving Eq. (23). Based on Proposition 4, we can solve Eq. (23)
by replacing $\mathrm{CS}_{k}(\cdot)$ with $\mathrm{RS}_{k}(\cdot)$ in
Algorithm 3. Similar to Theorem 2, the variant of Algorithm 3 has theoretical
convergence guarantees.

###  IV-D Solve ONMF_$\ell_{20}$ and ONMF_$\ell_{c0}$

We develop an efficient algorithm for ONMF_$\ell_{20}$ by turning it into a
series of constrained and penalized matrix factorization problems, i.e., Eq.
(14) with different $\rho$, which can be solved by the PALM or maPALM
algorithm. The detailed algorithm is given in Algorithm 4.

Algorithm 4 ONMF_$\ell_{20}$ Algorithm

0: $\bm{X}\in\mathbb{R}^{p\times n}$, $\rho\geq 0$ and $\gamma>1$.

0: $\bm{W}\in\mathbb{R}^{p\times r}$ and $\bm{H}\in\mathbb{R}^{r\times n}$.

1: Initialize $(\bm{W}^{0},\bm{H}^{0})$

2: for $k=1,2,\cdots,K$ do

3: Obtain $(\bm{W},\bm{H})$ by solving Eq. (15) with $\rho$ using PALM or
maPALM with the initialization $(\bm{W}^{0},\bm{H}^{0})$.

4: $\rho:=\gamma\rho$ and $(\bm{W}^{0},\bm{H}^{0}):=(\bm{W},\bm{H})$

5: end for

6: return $\bm{W}$ and $\bm{H}$.

In Algorithm 4, we set the default parameters $\rho=0.1$, $\gamma=1.5$, and
$K=10$. As $\rho\rightarrow\infty$ in the step 4 of Algorithm 4, the impact of
the penalty grows, and the estimated $\bm{H}$ will approach a non-negative
orthogonal matrix. Finally, we can also use a method similar to Algorithm 4 to
solve ONMF_$\ell_{c0}$ by turning it into a series of constrained and
penalized matrix factorization problems, i.e., Eq. (23) with different $\rho$.

##  V Experiments

We evaluate the effectiveness of these proposed SSNMF methods for clustering
task and compare them with state-of-the-art matrix factorization, k-means and
two sparse k-means methods on the synthetic and scRNA-seq data. All competing
methods are listed as follows:

  * •

ONMF_$\ell_{20}$: ONMF with $\ell_{2,0}$-norm constraint (Eq. 7).

  * •

ONMF_$\ell_{20}\\_\rho$: ONMF_$\ell_{20}$ with a fixed $\rho$ (Eq. 14).

  * •

ONMF_$\ell_{c0}$: column-wise sparse ONMF (Eq. 23).

  * •

ONMF_$\ell_{0}$: ONMF with $\ell_{0}$-norm constraint.

  * •

NMF_$\ell_{20}$: NMF with $\ell_{2,0}$-norm constraint.

  * •

ONMF: Orthogonal NMF [25].

  * •

NMF: Non-negative Matrix factorization [31].

  * •

NMF_$\ell_{0}$: NMF with $\ell_{0}$-norm constraint [22].

  * •

NMF_$\ell_{c0}$: Column-wise sparse NMF [32].

  * •

Kmeans: A baseline unsupervised method.

  * •

Kmeans_$\ell_{1}$: Sparse $k$ means with $\ell_{1}$ penalty in [33].

  * •

Kmeans_$\ell_{0}$: Sparse $k$ means with $\ell_{\infty}$/$\ell_{0}$ penalty in
[34].

Wherein the proposed SSNMF methods in this study include ONMF_$\ell_{20}$,
ONMF_$\ell_{20}\\_\rho$, ONMF_$\ell_{c0}$, ONMF_$\ell_{0}$ and
NMF_$\ell_{20}$.

###  V-A Evaluation metrics

To evaluate the clustering performance, we use three metrics including NMI
(Normalized Mutual Information) [35], Purity and Entropy [14]. Suppose $c$ is
the number of clustering clusters and $d$ is the number of true categories.
Then NMI is defined as

$\textrm{NMI}=\frac{\sum_{l=1}^{c}\sum_{h=1}^{d}t_{l,h}\log(\frac{n\cdot
t_{l,h}}{t_{l}\hat{t}_{h}})}{\sqrt{(\sum_{l=1}^{c}t_{l}\log\frac{t_{l}}{n})(\sum_{h=1}^{d}\hat{t}_{h}\log\frac{\hat{t}_{h}}{n})}},$
(27)

Purity is defined as

$\textrm{Purity}=\frac{1}{n}\sum_{l=1}^{c}\underset{1\leq h\leq
d}{\max}t_{l,h},$ (28)

Entropy is defined as

$\textrm{Entropy}=\frac{1}{n\log_{2}d}\sum_{l=1}^{c}\sum_{h=1}^{d}t_{l,h}\log_{2}\frac{t_{l,h}}{t_{l}},$
(29)

where $n$ is the number of considered samples, $t_{l}$ is the number of
samples from the $l$-th cluster $C_{l}$, which is obtained by clustering
method and $\hat{t}_{h}$ is the number of samples from the $h$-th ground truth
class $G_{h}$. $t_{l,h}$ denotes the number of overlapping samples between
$C_{l}$ and $G_{h}$. The larger the values of NMI and purity or the smaller
the value of entropy, the better the clustering performance.

###  V-B Application to synthetic data

We generate a synthetic data $\bm{X}\in\mathbb{R}^{p\times n}$ where $p=500$
features and $n=60$ samples from three true classes. Firstly, the elements in
$\bm{X}$ satisfy $X_{ij}\sim N(0,1)$ ($1\leq i\leq 60$, $1\leq j\leq 20$),
$X_{ij}\sim N(0,1)$ ($31\leq i\leq 90$, $21\leq j\leq 40$) and $X_{ij}\sim
N(0,1)$ ($61\leq i\leq 120$, $41\leq j\leq 60$) and $X_{ij}\sim 0.9*N(0,1)$
for the other $i$ and $j$, where $N(0,1)$ denotes standard normal
distribution. Secondly, we set $X_{ij}:=|X_{ij}|$ for any $i$ and $j$ to
ensure every element in the final synthetic data matrix $\bm{X}$ is positive
(see Figure 2A).

First of all, we show the convergence performance of PALM and maPALM when they
are used to solve ONMF_$\ell_{20}$ model in Eq. (14) on the synthetic data
(Figure 2B). We find that the convergence speed of maPALM is significantly
faster than that of PALM.

Secondly, to validate the effectiveness of our proposed methods to perform
feature selection, we compare them with state-of-the-art methods on the
synthetic data. The parameters of these compared methods are carefully
adjusted to give their best performances and all methods are repeated 10 times
using different initial points for comparison. The corresponding NMI scores
are recorded for each method and the methods with higher NMI averages are
regarded as more accurate ones. We find that the SSNMF methods including
ONMF_$\ell_{20}$ and NMF_$\ell_{20}$ outperform other algorithms in terms of
NMI (Figure 3). Interestingly, we also find that as $\rho$ becomes larger in
the Algorithm 3, the orthogonality level of estimated $\bm{H}$ is better, such
that the clustering performance is better (Figure 4). This result implies that
it is possible to improve the clustering accuracy by adding the non-negative
and orthogonal constraint in SSNMF model.

Figure 2:  Results on the synthetic data. (A) Heatmap showing the synthetic
data. (B) Convergence performance of PALM and maPALM for ONMF_$\ell_{20}$ with
a fixed $\rho=0.5$.  Figure 3: Comparison of 12 unsupervised clustering
methods in terms of NMI on the synthetic data. Note that we set $\rho=1$ for
ONMF_$\ell_{20}\\_\rho$.  Figure 4:  Scatter plots showing the change of (A)
$\rho$ in Algorithm 3, (B) orthogonality score of estimated $\bm{H}$ and (C)
NMI score. For the definition of orthogonality level, please see Eq. 23 of
[25].

Finally, we evaluate whether ONMF can detect outliers by sorting the values of
estimated $\bm{H}$. To this end, we generate a new synthetic data
$\bm{X}\in\mathbb{R}^{p\times n}$ where $p=500$ features and $n=60$ samples
from three classes. Firstly, the elements in $\bm{X}$ satisfy $X_{ij}\sim
N(0,1)$ ($1\leq i\leq 60$, $1\leq j\leq 20$), $X_{ij}\sim N(0,1)$ ($31\leq
i\leq 90$, $21\leq j\leq 40$), and $X_{ij}\sim 0.9*N(0,1)$ for other $i$ and
$j$. Secondly, we set $X_{ij}:=|X_{ij}|$ and obtain the final synthetic data.
Note that the columns from 41 to 60 in $\bm{X}$ correspond to outliers. We
apply ONMF_$\ell_{20}$ with parameter $k=90$ to the synthetic data. We find
that the outliers (_i.e._ , noise samples) are those with relatively small
values in the estimated $\bm{H}$ (Figure 5). This result implies
ONMF_$\ell_{20}$ can detect these outliers by checking the values of columns
in the estimated $\bm{H}$.

Figure 5:  (A) Heatmap corresponding to the estimated $\bm{H}$ by
ONMF_$\ell_{20}$ on the second synthetic data. (B) Scatter plot showing the
largest values from each column of the estimated $\bm{H}$.

###  V-C Application to biological data

In this study, the proposed methods and other comparison methods are evaluated
on three scRNA-seq datasets as follows:

  * •

Pollen dataset [36] contains 301 single cells from 11 cell populations which
are divided into 4 classes including Blood cells, Neural cells, Dermal or
epidermal cells and Pluripotent cells.

  * •

Camp1 dataset [37] contains 425 single cells from human liver cells which are
divided into 5 populations, named as iPS_day_0, De_day_6, IH_day_14,
MH_day_21, and HE_day_8.

  * •

Lake dataset [38] contains 3042 single cells from human liver cells which are
divided into 16 populations named as Ex1, Ex3, Ex4, In6, In1, Ex5, In5, Ex7,
In8, Ex2, Ex6, In7, In4, Ex8, In3, In2.

For each scRNA-seq dataset, we first filter out these genes which are not
expressed in more than 70% of cells. We then use a logarithmic transformation
$x=\log 2(x)$ to transform raw expression data. After preprocessing, 8747,
8058 and 5000 genes are retained for Pollen, Camp1 and Lake datasets,
respectively.

Table II: Comparison in terms of (NMI % $\pm$ std), (Purity % $\pm$ std) and
(Entropy % $\pm$ std) on the Pollen, Lake and Camp1 datasets.

Pollen data #gene NMI$\pm$sd Purity$\pm$sd Entropy$\pm$sd NMF_$\ell_{20}$ 2000
85.87 $\pm$ 1.29 91.86 $\pm$ 0.23 19.10 $\pm$ 0.97 ONMF_$\ell_{20}\\_\rho$
2000 84.96 $\pm$ 0.86 91.69 $\pm$ 0.16 19.79 $\pm$ 0.61 ONMF_$\ell_{20}$ 2000
80.71 $\pm$ 3.56 89.30 $\pm$ 2.76 24.34 $\pm$ 4.30 ONMF_$\ell_{c0}$ * 81.62
$\pm$ 0.32 90.73 $\pm$ 0.11 22.28 $\pm$ 0.22 ONMF_$\ell_{0}$ * 51.67 $\pm$
4.58 67.08 $\pm$ 5.65 27.56 $\pm$ 1.92 NMF_$\ell_{0}$ * 54.46 $\pm$ 6.86 67.34
$\pm$ 5.78 24.97 $\pm$ 5.59 NMF_$\ell_{c0}$ * 81.45 $\pm$ 0.88 91.89 $\pm$
0.17 21.28 $\pm$ 0.69 ONMF all 83.70 $\pm$ 2.81 91.23 $\pm$ 2.52 20.52 $\pm$
4.07 Kmeans all 71.39 $\pm$ 7.89 82.52 $\pm$ 5.27 31.37 $\pm$ 5.17
Kmeans$\\_\ell_{0}$ 2000 72.36 $\pm$ 7.30 81.40 $\pm$ 8.25 29.48 $\pm$ 7.92
Kmeans$\\_\ell_{1}$ 2039 78.09 $\pm$ 9.60 88.54 $\pm$ 7.01 21.45 $\pm$ 5.43
Camp1 data #gene NMI$\pm$sd Purity$\pm$sd Entropy$\pm$sd NMF_$\ell_{20}$ 2000
89.94 $\pm$ 1.30 96.07 $\pm$ 0.39 9.88 $\pm$ 0.88 ONMF_$\ell_{20}\\_\rho$ 2000
89.32 $\pm$ 3.16 92.40 $\pm$ 8.13 10.96 $\pm$ 3.68 ONMF_$\ell_{20}$ 2000 71.00
$\pm$ 9.25 72.12 $\pm$ 5.54 29.04 $\pm$ 6.15 ONMF_$\ell_{c0}$ * 47.31 $\pm$
8.06 44.35 $\pm$ 5.13 10.25 $\pm$ 3.95 ONMF_$\ell_{0}$ * 18.57 $\pm$ 18.69
34.14 $\pm$ 7.56 5.85 $\pm$ 4.99 NMF_$\ell_{0}$ * 12.41 $\pm$ 17.04 31.25
$\pm$ 6.72 5.94 $\pm$ 6.38 NMF_$\ell_{c0}$ * 85.54 $\pm$ 5.20 91.67 $\pm$ 2.74
17.27 $\pm$ 4.10 ONMF all 78.98 $\pm$ 9.88 78.31 $\pm$ 9.89 23.19 $\pm$ 8.71
Kmeans all 78.06 $\pm$ 1.18 79.72 $\pm$ 5.76 24.09 $\pm$ 4.16
Kmeans$\\_\ell_{0}$ 2000 73.13 $\pm$ 5.36 71.62 $\pm$ 7.18 28.67 $\pm$ 4.58
Kmeans$\\_\ell_{1}$ 2005 71.96 $\pm$ 1.32 79.11 $\pm$ 1.10 28.07 $\pm$ 1.33
Lake data #gene NMI$\pm$sd Purity$\pm$sd Entropy$\pm$sd NMF_$\ell_{20}$ 2000
72.43 $\pm$ 2.15 76.73 $\pm$ 2.25 28.88 $\pm$ 2.86 ONMF_$\ell_{20}$ 2000 57.19
$\pm$ 4.16 63.92 $\pm$ 4.30 42.67 $\pm$ 5.01 ONMF_$\ell_{20}\\_\rho$ 2000
71.99 $\pm$ 2.19 76.23 $\pm$ 2.49 29.60 $\pm$ 2.60 ONMF_$\ell_{c0}$ * 10.23
$\pm$ 10.54 36.97 $\pm$ 2.42 8.33 $\pm$ 8.70 ONMF_$\ell_{0}$ * 0.92 $\pm$ 0.64
35.01 $\pm$ 0.21 1.53 $\pm$ 1.42 NMF_$\ell_{0}$ * 17.08 $\pm$ 14.82 41.30
$\pm$ 6.59 11.05 $\pm$ 11.93 NMF_$\ell_{c0}$ * 66.35 $\pm$ 2.20 72.00 $\pm$
1.95 34.61 $\pm$ 3.22 ONMF all 61.44 $\pm$ 3.84 67.08 $\pm$ 4.00 40.20 $\pm$
4.43 Kmeans all 59.96 $\pm$ 1.31 68.43 $\pm$ 1.68 42.10 $\pm$ 1.59
Kmeans$\\_\ell_{0}$ 2000 62.05 $\pm$ 2.67 68.93 $\pm$ 2.62 40.01 $\pm$ 2.91
Kmeans$\\_\ell_{1}$ 2270 70.73 $\pm$ 1.34 77.62 $\pm$ 2.12 32.32 $\pm$ 1.28

Figure 6: Results are shown in (A) to (D) when NMF_$\ell_{20}$ is applied to
the Pollen and Camp1 scRNA-seq datasets. (A) Heatmap showing $\bm{W}$ and
$\bm{H}$ obtained from NMF_$\ell_{20}$ with $k=2000$ on the Pollen dataset
where 2000 genes are selected across 301 cells. (B) A Venn diagram showing
overlap level for the selected genes from four biclusters (M1 to M4) on the
Pollen dataset and total of 329 genes are selected. (C) Heatmap showing
$\bm{W}$ and $\bm{H}$ obtained from NMF_$\ell_{20}$ with $k=2000$ on the Camp1
dataset where 2000 genes are selected across 425 cells. (D) A Venn diagram
showing overlap level for the selected genes from five biclusters (M1 to M5)
on the Camp1 dataset and total of 443 genes are selected.

For ONMF_$\ell_{20}$, we set $k=2000$ (to extract 2000 genes) which is to
ensure that the number of selected genes is about 2000 for further analysis of
biological function; rank $r$ equal to the number of true classes of scRNA-seq
datasets (herein $r=4$ for the Pollen dataset, $r=5$ for the Camp1 dataset and
$r=16$ for the Lake dataset); $\rho=0.1$ and $\gamma=1.5$ which are for
updating next $\rho:=\gamma\rho$ (herein $\rho$ is updated up to 10 times);
and $\epsilon=1e-3$ in Algorithm 3. For fairness of comparison, we ensure
these estimated $\bm{W}$ of all sparse learning methods including
NMF_$\ell_{20}$, ONMF_$\ell_{c0}$, NMF_$\ell_{c0}$, ONMF_$\ell_{0}$,
NMF_$\ell_{0}$ and Kmeans$\\_\ell_{0}$ have the same sparsity level, and all
methods are repeated 10 times using different initial points for comparison.

We evaluate the clustering performance of all methods in terms of NMI, Purity
and Entropy. The detailed results on the pollen, camp1 and lake datasets are
summarized in Table II. These results show that the proposed SSNMF methods,
especially NMF_$\ell_{20}$, outperform other methods. The use of
$\ell_{2,0}$-norm enables some SSNMF methods to select some important features
by checking the non-zero rows of $\bm{W}$. Finally, we also discuss the
influence of $k$-choice on the clustering performance for these proposed SSNMF
with $\ell_{2,0}$-norm methods and the results show that the proposed
NMF_$\ell_{20}$ outperforms Kmeans_$\ell_{0}$ in different situations (Table
III).

Table III: Comparison of (NMI %) averages with different number of genes on the Pollen, Camp1 and Lake datasets. Pollen data ($k$=) | 500 | 1000 | 2000 | 3000 | 4000 | 5000
---|---|---|---|---|---|---
NMF_$\ell_{20}$ | 81.77 | 84.46 | 85.87 | 85.70 | 85.04 | 85.04
ONMF_$\ell_{20}$ | 69.36 | 73.06 | 80.71 | 84.27 | 84.65 | 84.34
ONMF_$\ell_{20}\\_\rho$ | 82.35 | 84.09 | 84.96 | 86.02 | 85.21 | 85.04
Kmeans$\\_\ell_{0}$ | 80.79 | 76.98 | 72.36 | 71.04 | 71.57 | 69.49
Camp1 data ($k$=) | 500 | 1000 | 2000 | 3000 | 4000 | 5000
NMF_$\ell_{20}$ | 91.56 | 91.84 | 89.94 | 90.11 | 90.28 | 90.28
ONMF_$\ell_{20}$ | 62.57 | 63.93 | 71.00 | 72.54 | 77.34 | 80.77
ONMF_$\ell_{20}\\_\rho$ | 90.59 | 90.80 | 89.32 | 89.21 | 89.24 | 89.03
Kmeans$\\_\ell_{0}$ | 74.26 | 73.41 | 73.13 | 76.28 | 75.31 | 74.23
Lake data ($k$=) | 500 | 1000 | 2000 | 3000 | 4000 | 5000
NMF_$\ell_{20}$ | 72.96 | 72.88 | 72.43 | 72.34 | 72.29 | 72.23
ONMF_$\ell_{20}$ | 52.94 | 55.67 | 57.19 | 58.64 | 59.95 | 61.44
ONMF_$\ell_{20}\\_\rho$ | 72.48 | 72.37 | 71.99 | 71.67 | 71.55 | 71.52
Kmeans$\\_\ell_{0}$ | 70.83 | 69.25 | 62.05 | 58.83 | 56.17 | 55.65

###  V-D Biological analysis

In this section, we show that the SSNMF with $\ell_{2,0}$-norm methods can be
used for the identification of subpopulation and gene selection for scRNA-seq
data. Based on the experiment results in the previous section, the
NMF_$\ell_{20}$ achieved the best performance. Therefore, we select the
computing results from NMF_$\ell_{20}$ ($k=2000$) on Pollen and Camp1 datasets
as a example for further biological analysis. Some of the significant gene and
sample expression patterns can be identified based on the output $\bm{W}\in
R^{p\times r}$ and $\bm{H}\in R^{r\times n}$ from NMF_$\ell_{20}$.

To be simplify, a biological functional bicluster is defined as a gene subset
with a sample subset (also as cell subset). For each pair of $\bm{w}_{i}$ (in
$\bm{W}$) and $\bm{h}^{i}$ (in $\bm{H}$), a bicluster is extracted based on
the following computational steps.

  * •

Step 1: For the $i$-th column of $\bm{W}$ ($\bm{w}_{i}$), the higher numerical
values, the more important the corresponding genes are. We do a z-score
normalization for the $\bm{w}_{i}$ using the formula
$\bm{z}_{i}=(\bm{w}_{i}-mean(\bm{w}_{i}))/sd(\bm{w}_{i})$. The genes with the
corresponding coefficient values of $\bm{w}_{i}$, whose z-scores are larger
than a given threshold $T$, are extracted as the bicluster gene subset.

  * •

Step 2: For $i$-th row of $\bm{H}$ ($\bm{h}^{i}$), we screen the cells with
the largest coefficient in their corresponding column, i.e., $\\{j|w_{ij}\geq
w_{tj}\leavevmode\nobreak\ \mbox{for}\leavevmode\nobreak\
j=1,\cdots,n,\leavevmode\nobreak\ \forall\leavevmode\nobreak\ t\\}$. Thus, we
obtain a single cell set for the bicluster.

For the output $\bm{W}$ and $\bm{H}$ from Pollen dataset, we first calculate
z-scores normalization for each column of $\bm{W}$ and then rank all genes
according to their z-socre values. These genes with z-score more than the
threshold $T=1.5$ are regarded as the gene set of biclusters. We extract four
biclusters with total 329 genes (Figure 6A and B). Bicluster 1 contains 170
genes and 68 cells. The selected cells are all Blood cells. Bicluster 2
contains 136 genes and 99 cells and all selected cells are Dermal/Epidermal
cells. Bicluster 3 contains 191 genes and 65 cells and all selected cells are
Neural cells. Bicluster 4 contains 71 genes and 69 cells while 45 of 69 cells
are Blood cells and 24 of 69 are Pluripotent cells. To be interesting, some
genes are shared in multiple biclusters. 27 genes are shared on all four
biclusters. It shows that the genes may be hub genes and play a joint role in
multiple biological functions (pathways).

To demonstrate the biological function of the gene sets from these identified
biclusters, we perform the gene function enrichment analysis. We retrieve the
KEGG pathways data information from Molecular Signatures Database (MSigDB,
http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). KEGG pathways are a class
of collection of manually drawn pathway maps representing the biological
knowledge of the molecular interaction, reaction and relation networks. Each
KEGG pathway is consisted by a set of functionally gene set which expresses
complex regulatory mechanism among different genes. The hypergeometric test is
applied for the biological statistical analysis. A class of critical KEGG
signaling pathway is significantly enriched for the identified four bicluster
gene sets with the Benjamini-Hochberg adjusted $p<0.05$.

As expected, most of the enriched KEGG signaling pathways are related on the
developing cerebral cortex which is highly consistent to the Pollen dataset
from the diverse neural cell types [36]. Interestingly, we find that multiple
bicluster gene sets are enriched on the ribosome pathway which is the cell
factories responsible for making proteins and the ribosome pathway has been
reported to play an important role in brain development [39]. In addition,
some brain disease-related pathways have been discovered. For example,
bicluster 1 is enriched in the pathways like Parkinson’s disease ($p=$7.5e-14)
and Alzheimer’s disease ($p=$3.9e-11). Bicluster 2 is enriched in some KEGG
pathways like proteasome ($p=$1.2e-3), leukocyte transendothelial migration
($p=$1.2e-3), regulation of actin cytoskeleton ($p=$1.2e-3), and focal
adhesion ($p=$1.4e-3).

For the Camp1 dataset, five biclusters are extracted from the $\bm{W}$ and
$\bm{H}$ (Figure 6C and D). Similarly, we calculate z-scores normalization for
each column of $\bm{W}$ and then rank all genes according to their z-socre
values. These genes with z-score more than the threshold $T=1.2$ are regarded
as the gene set of biclusters. Bicluster 1 contains 227 genes and 73 cells
where all selected cells are MH day 21 cells, bicluster 2 contains 65 genes
and 81 cells where 80 of 81 cells are iPS day 0 cells, bicluster 3 contains
166 genes and 67 cells where all selected cells are De day 6 cells, bicluster
4 contains 111 genes and 82 cells where 77 of 82 cells are IH day 14 cells and
5 of 82 are MH day 21 cells, and bicluster 5 contains 205 genes and 122 cells
where 113 of 122 cells are HE day 8 cells, 3 of 122 are De day 6 cells, 4 of
122 are IH day 14 cells and 2 of 122 are MH day 21 cells. Similar to the
results on Pollen dataset, the cells in each identified bicluster are highly
pure of population. Meanwhile, each bicluster possesses different domain genes
while a small number of genes overlap each other (see Figure 6D). Due to these
extracted single cells in the Camp1 dataset are from human liver [37], the
biological function analysis show that the identified bicluster gene sets are
significantly enriched a series of KEGG pathways (Benjamini-Hochberg adjusted
$p<0.05$) which are highly related to biological processes associated with
liver. For example, bicluster 1 is significantly enriched in the oxidative
phosphorylation pathway ($p=7.6e-13$ and there are 26 genes in this pathway)
which have been reported to be related to liver [40]. Bicluster 2 is
significantly enriched in the glycolysis gluconeogenesis pathway ($p=6.2e-5$
and there are 6 genes in this pathway) which have been reported to be related
to liver [41]. All these results show that the proposed method can be used for
single cell type discovery, gene selection and biological process analysis.

##  VI Conclusion

In this paper, we present a class of SSNMF models with $\ell_{2,0}$-norm
constraint. We prove that $\ell_{2,0}$-norm satisfies the KŁ property, such
that the PALM algorithm can be used to solve a class of non-convex and non-
smooth optimization problems with $\ell_{2,0}$-norm constraint. Especially, we
first introduce the NMF_$\ell_{20}$ model which integrates feature selection
in the NMF model. To improve the convergence rate of PALM, we further develop
an accelerated version of PALM (maPALM) to solve NMF_$\ell_{20}$. We also
prove the convergence of proposed algorithms (PALM and maPALM) when they are
used to solve NMF_$\ell_{20}$. To integrate feature selection and non-negative
orthogonal constraint in the NMF model, we furthermore introduce the
ONMF_$\ell_{20}$ model. We develop an efficient algorithm to solve it by using
a penalty function method. Briefly, the algorithm converts ONMF_$\ell_{20}$
into a series of constrained and penalized NMF problems which can be solved by
the PALM and maPALM algorithms. Finally, we compare these proposed SSNMF
methods with other methods for clustering task on the synthetic and scRNA-seq
data. The results show that the proposed SSNMF methods can be used not only
for clustering (single cell type discovery), but also for gene selection and
biological function analysis.

##  Appendix A Definitions and Proofs

###  A-A Mathematical definitions for non-convex analysis

We introduce some mathematical definitions which are used in this study for
non-convex analysis [42, 22, 43].

######  Definition 5.

(Proper) $f(\bm{x})$ is proper if
$\mathbf{dom}(f):=\\{\bm{x}\in\mathbb{R}^{n}:f(\bm{x})<+\infty\\}$ is nonempty
and $f(\bm{x})>-\infty$.

######  Definition 6.

(Lower semi-continuous) $f(\bm{x})$ is lower semi-continuous if
$\liminf\limits_{\bm{x}\rightarrow\bm{x}_{0}}f(\bm{x})\geq f(\bm{x}_{0})$ at
any point $\bm{x}_{0}\in\mathbf{dom}(f)$.

######  Definition 7.

(Coercive Function) $f(\bm{x})$ is called coercive if $f(\bm{x})$ is bounded
from below and $f(\bm{x})\rightarrow\infty$ if $\|\bm{x}\|\rightarrow\infty$.

######  Definition 8.

(Lipschitz smooth) $f(\bm{x})$ is Lipschitz smooth if it is differentiable and
there exists $L>0$ and such that

$\|\nabla f(\bm{x})-\nabla f(\bm{y})\|\leq
L\|\bm{x}-\bm{y}\|,\forall\bm{x},\bm{y}\in\mathbb{R}^{n}.$

Any such $L$ is considered to as a Lipschitz constant for $f(\bm{x})$.

######  Definition 9.

(Subdifferential) Let $f:\mathbb{R}^{n}\rightarrow(-\infty,\infty]$ be a
proper and lower semi-continuous function. Then the Frecht sub-differential of
$f$, denoted as $\widehat{\partial}f$, at point $\bm{x}\in\mathbf{dom}(f)$ is
the set of all vectors $\bm{z}$ which satisfies

$\liminf\limits_{\bm{x}\neq\bm{y},\bm{y}\rightarrow\bm{x}}\frac{f(\bm{y})-f(\bm{x})-\langle\bm{z},\bm{y}-\bm{x}\rangle}{\|\bm{y}-\bm{x}\|}\geq
0,$

where $\langle\cdot,\cdot\rangle$ denotes the inner product. Then the limiting
Frecht sub-differential, or simply the sub-differential, denoted as $\partial
f$, at $\bm{x}\in\mathbf{dom}(f)$ is the following closure of
$\widehat{\partial}f$:

$\\{\bm{z}\in\mathbb{R}^{n}:\exists(\bm{x}^{k},g(\bm{x}^{k}))\rightarrow(\bm{x},f(\bm{x}))\\},$

where $\bm{z}^{k}\in\widehat{\partial}f(\bm{x}^{k})\rightarrow\bm{z}$ when
$k\rightarrow\infty$.

######  Definition 10.

(Critical Point) A point $\bm{x}$ is called a critical point of function $f$
if $0\in\partial f(\bm{x})$.

######  Definition 11.

(Semi-algebraic set and function) A subset $\Omega$ of $\mathbb{R}^{n}$ is a
real semi-algebraic set if there exist a finite number of real polynomial
functions $r_{ij},h_{ij}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ such that

$\Omega=\bigcup_{j=1}^{p}\bigcap_{i=1}^{q}\big{\\{}\bm{x}\in\mathbb{R}^{n}:r_{ij}(\bm{x})=0\leavevmode\nobreak\
\mbox{and}\leavevmode\nobreak\ h_{ij}(\bm{x})<0\big{\\}}.$

Function $f$ is called semi-algebraic function if its graph
$\\{(\bm{x},z)\in\mathbb{R}^{n+1}:f(\bm{x})=z\\}$ is a semi-algebraic subset
of $\mathbb{R}^{n+1}$.

######  Definition 12.

(Kurdyka-Łojasiewicz property and function) Function $f$ has the Kurdyka-
Łojasiewicz (KŁ) property at $\bar{\bm{x}}\in\mathbf{dom}(\partial
f):=\\{\bm{x}\in\mathbb{R}^{n}:\partial f(\bm{x})\neq\emptyset\\}$ if there
exist $\eta\in(0,\infty]$, a neighborhood $U$ of $\bar{\bm{x}}$ and a function
$\phi:[0,\eta)\rightarrow\mathbb{R}_{+}$ which satisfies (1) $\phi$ is
continuous at 0 and $\phi(0)=0$; (2) $\phi$ is concave and $C^{1}$ on
$(0,\eta)$; (3) for all $s\in(0,\eta):\phi^{\prime}(s)>0$, such that for all
$\bm{x}\in U\cap[f(\bar{\bm{x}})<f(\bm{x})<f(\bar{\bm{x}})+\eta],$ the
following inequality holds

$\phi^{\prime}(f(\bm{x})-f(\bar{\bm{x}}))dist(0,\partial f(\bm{x}))\geq 1.$

Function $f$ is called a KŁ function if $f$ satisfies the KŁ property at each
point of $\mathbf{dom}(\partial f)$. Moreover, Theorem 3 in [22] shows that
all proper, lower semicontinuous and semi-algebraic functions satisfy KŁ
property property.

###  A-B Proof of Proposition 1

###### Proof.

Reference [25], we give the following proof. Let
$F_{\rho}(\bm{W},\bm{H})=\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}+\frac{\rho}{2}\sum_{j=1}^{n}\Big{(}(\bm{1}^{T}\bm{h}_{j})^{2}-\|\bm{h}_{j}\|_{2}^{2}\Big{)}$,
$\bm{Z}^{*}:=(\bm{W}^{*},\bm{H}^{*})$ is a local minimizer of (LABEL:equ-10).
If $\Omega_{1}$ is the feasible set of (7), then there exist a neighborhood of
$\bm{Z}^{*}$, $N_{\epsilon}(\bm{Z}^{*})=\\{\bm{Z}|\|\bm{Z}-\bm{Z}^{*}|\\}$
where $\epsilon>0$ and it satisfies $F_{\rho}(\bm{Z}^{*})\leq
F_{\rho}(\bm{Z})$ for any $\bm{Z}\in N_{\epsilon}(\bm{Z}^{*})\cap\Omega$. Let
$\phi(\bm{h}_{j})=(\bm{1}^{T}\bm{h}_{j})^{2}-\|\bm{h}_{j}\|_{2}^{2}$, we note
that $\phi(\bm{h}_{j})>0$ if $\bm{h}_{j}$ has at least two non-zeros entries
and $\phi(\alpha\bm{h}_{j})$ is an increasing function with respect to
$\alpha$. Suppose that there exists an $j^{\prime}$ and
$\bm{h}_{j^{\prime}}^{*}$ is infeasible to (7) and thus
$\phi(\bm{h}_{j^{\prime}}^{*})>0$. Then, for a scalar $\alpha\in(0,1)$ and
$\bm{Z}_{\alpha}:=(\bm{W}^{*}/\alpha,\alpha\bm{H}^{*})$ is a feasible point to
(LABEL:equ-10). Specifically, since
$\|\bm{Z}_{\alpha}-\bm{Z}^{*}\|=(1/\alpha-1)^{2}\|\bm{W}^{*}\|_{F}^{2}+(\alpha-1)^{2}\|\bm{H}^{*}\|_{F}^{2}\leq\max\\{(1/\alpha-1)^{2},(\alpha-1)^{2}\\}\|\bm{Z}^{*}\|_{F}^{2}$.
To have $\bm{Z}_{\alpha}\in N_{\epsilon}(\bm{Z}^{*})$, it is sufficient to let
$\max\\{(1/\beta-1)^{2},(\beta-1)^{2}\\}<\alpha<1$ where
$\beta=\sqrt{\epsilon/(2*\|\bm{Z}^{*}\|_{F}^{2})}$. Thus, we have
$F_{\rho}(\bm{Z}_{\alpha})-F_{\rho}(\bm{Z}^{*})=\phi(\alpha\bm{h}_{j^{\prime}}^{*})-\phi(\bm{h}_{j^{\prime}}^{*})<0.$
This is a contradiction. So, we have proved that $(\bm{W}^{*},\bm{H}^{*})$ is
a feasible solution of (7). Next, will prove that $(\bm{W}^{*},\bm{H}^{*})$ is
also a local minimizer of (7). Let $\Omega_{2}$ as the feasible set of (7),
then we have
$\\{N_{\epsilon}(\bm{Z}^{*})\cap\Omega_{2}\\}\subseteq\\{N_{\epsilon}(\bm{Z}^{*})\cap\Omega_{1}\\}$.
This implies that $G_{\rho}(\bm{Z}^{*})=F_{\rho}(\bm{Z}^{*})\leq
F_{\rho}(\bm{Z})=G_{\rho}(\bm{Z})$ for any $\bm{Z}\in
N_{\epsilon}(\bm{Z}^{*})\cap\Omega_{2}$ where
$G_{\rho}(\bm{W},\bm{H})=\frac{1}{2}\|\bm{X}-\bm{W}\bm{H}\|_{F}^{2}$. So, we
prove that $(\bm{W}^{*},\bm{H}^{*})$ is a local minimizer of (7). ∎

## Acknowledgment

This work was supported by Key-Area Research and Development Program of
Guangdong Province [2020B0101350001], and the National Science Foundation of
China [61272274], and Natural Science Foundation of Jiangxi Province of China
[20192BAB217004] and China Postdoctoral Science Foundation [2020M671902].

## References

  * [1] M. D. Luecken and F. J. Theis, “Current best practices in single-cell rna-seq analysis: a tutorial,” _Mol. Syst. Biol._ , vol. 15, no. 6, p. e8746, 2019\.
  * [2] V. Y. Kiselev, T. S. Andrews, and M. Hemberg, “Challenges in unsupervised clustering of single-cell RNA-seq data,” _Nat. Rev. Genet._ , vol. 20, no. 5, pp. 273–282, 2019.
  * [3] J.-P. Brunet, P. Tamayo, T. R. Golub, and J. P. Mesirov, “Metagenes and molecular pattern discovery using matrix factorization,” _Proc. Natl. Acad. Sci._ , vol. 101, no. 12, pp. 4164–4169, 2004.
  * [4] M. Stražar, M. Žitnik, B. Zupan, J. Ule, and T. Curk, “Orthogonal matrix factorization enables integrative analysis of multiple rna binding proteins,” _Bioinformatics_ , vol. 32, no. 10, pp. 1527–1535, 2016.
  * [5] J.-X. Liu, D. Wang, Y.-L. Gao, C.-H. Zheng, Y. Xu, and J. Yu, “Regularized non-negative matrix factorization for identifying differentially expressed genes and clustering samples: a survey,” _IEEE/ACM Trans. Comput. Biol. Bioinform._ , vol. 15, no. 3, pp. 974–987, 2017.
  * [6] S. Zhang, C.-C. Liu, W. Li, H. Shen, P. W. Laird, and X. J. Zhou, “Discovery of multi-dimensional modules by integrative analysis of cancer genomic data,” _Nucleic Acids Res._ , vol. 40, no. 19, pp. 9379–9391, 2012.
  * [7] J. Chen and S. Zhang, “Discovery of two-level modular organization from matched genomic data via joint matrix tri-factorization,” _Nucleic Acids Res._ , vol. 46, no. 12, pp. 5967–5976, 2018.
  * [8] L. Zhang and S. Zhang, “Learning common and specific patterns from data of multiple interrelated biological scenarios with matrix factorization,” _Nucleic Acids Res._ , vol. 47, no. 13, pp. 6606–6617, 2019.
  * [9] X. Fu, K. Huang, N. D. Sidiropoulos, and W.-K. Ma, “Nonnegative matrix factorization for signal and data analytics: Identifiability, algorithms, and applications.” _IEEE Signal Process. Mag._ , vol. 36, no. 2, pp. 59–80, 2019\.
  * [10] C. Shao and T. Höfer, “Robust classification of single-cell transcriptome data by nonnegative matrix factorization,” _Bioinformatics_ , vol. 33, no. 2, pp. 235–242, 2017.
  * [11] Z. Duren, X. Chen, M. Zamanighomi, W. Zeng, A. T. Satpathy, H. Y. Chang, Y. Wang, and W. H. Wong, “Integrative analysis of single-cell genomics data by coupled nonnegative matrix factorizations,” _Proc. Natl. Acad. Sci._ , vol. 115, no. 30, pp. 7723–7728, 2018.
  * [12] J. D. Welch, V. Kozareva, A. Ferreira, C. Vanderburg, C. Martin, and E. Z. Macosko, “Single-cell multi-omic integration compares and contrasts features of brain cell identity,” _Cell_ , vol. 177, no. 7, pp. 1873–1887, 2019.
  * [13] P. O. Hoyer, “Non-negative matrix factorization with sparseness constraints,” _J. Mach. Learn. Res._ , vol. 5, no. Nov, pp. 1457–1469, 2004.
  * [14] H. Kim and H. Park, “Sparse non-negative matrix factorizations via alternating non-negativity-constrained least squares for microarray data analysis,” _Bioinformatics_ , vol. 23, no. 12, pp. 1495–1502, 2007.
  * [15] R. Peharz and F. Pernkopf, “Sparse nonnegative matrix factorization with $\ell_{0}$-constraints,” _Neurocomputing_ , vol. 80, pp. 38–46, 2012.
  * [16] J. Li, K. Cheng, S. Wang, F. Morstatter, R. P. Trevino, J. Tang, and H. Liu, “Feature selection: A data perspective,” _ACM Comput. Surv._ , vol. 50, no. 6, pp. 1–45, 2017.
  * [17] J. Gui, Z. Sun, S. Ji, D. Tao, and T. Tan, “Feature selection based on structured sparsity: A comprehensive study,” _IEEE Trans. Neural Netw. Learn. Syst._ , vol. 28, no. 7, pp. 1490–1507, 2017.
  * [18] F. Nie, H. Huang, X. Cai, and C. H. Ding, “Efficient and robust feature selection via joint $\ell_{21}$-norms minimization,” in _Adv. Neural. Inf. Process. Syst._ , 2010, pp. 1813–1821.
  * [19] H. Huang, C. Ding, and D. Luo, “Towards structural sparsity: An explicit $\ell_{2}$/$\ell_{0}$ approach,” in _IEEE 10th Int. Conf. Data Mining_ , 2010, pp. 344–353.
  * [20] T. Pang, F. Nie, J. Han, and X. Li, “Efficient feature selection via $\ell_{2,0}$-norm constrained sparse regression,” _IEEE Trans. Knowl. Data Eng._ , vol. 31, no. 5, pp. 880–893, 2019.
  * [21] X. Du, F. Nie, W. Wang, Y. Yang, and X. Zhou, “Exploiting combination effect for unsupervised feature selection by $\ell_{20}$ norm,” _IEEE Trans. Neural Netw. Learn. Syst._ , vol. 30, no. 1, pp. 201–214, 2018.
  * [22] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized minimization for nonconvex and nonsmooth problems,” _Math. Program._ , vol. 146, no. 1-2, pp. 459–494, 2014.
  * [23] C. Ding, T. Li, W. Peng, and H. Park, “Orthogonal nonnegative matrix t-factorizations for clustering,” in _Proc. 12th ACM Int. Conf. Knowl. Discovery Data Mining_ , 2006, pp. 126–135.
  * [24] K. Zhang, S. Zhang, J. Liu, J. Wang, and J. Zhang, “Greedy orthogonal pivoting algorithm for non-negative matrix factorization,” in _Int. Conf. Mach. Learn._ , 2019, pp. 7493–7501.
  * [25] S. Wang, T.-H. Chang, Y. Cui, and J.-S. Pang, “Clustering by orthogonal nmf model and non-convex penalty optimization,” _arXiv preprint arXiv:1906.00570_ , 2019.
  * [26] D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factorization,” _Nature_ , vol. 401, no. 6755, pp. 788–791, 1999\.
  * [27] H. Li and Z. Lin, “Accelerated proximal gradient methods for nonconvex programming,” in _Adv. Neural. Inf. Process. Syst._ , 2015, pp. 379–387.
  * [28] Q. Li, Y. Zhou, Y. Liang, and P. K. Varshney, “Convergence analysis of proximal gradient with momentum for nonconvex optimization,” in _Int. Conf. Mach. Learn._ , 2017, pp. 2111–2119.
  * [29] Y. Xu and W. Yin, “A globally convergent algorithm for nonconvex optimization based on block coordinate update,” _J. Sci. Comput._ , vol. 72, no. 2, pp. 700–734, 2017.
  * [30] N. Guan, D. Tao, Z. Luo, and B. Yuan, “NeNMF: An optimal gradient method for nonnegative matrix factorization,” _IEEE Trans. Signal Process._ , vol. 60, no. 6, pp. 2882–2898, 2012.
  * [31] Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” _SIAM J. Imaging Sci._ , vol. 6, no. 3, pp. 1758–1789, 2013.
  * [32] T. Pock and S. Sabach, “Inertial proximal alternating linearized minimization (iPALM) for nonconvex and nonsmooth problems,” _SIAM J. Imaging Sci._ , vol. 9, no. 4, pp. 1756–1787, 2016.
  * [33] D. M. Witten and R. Tibshirani, “A framework for feature selection in clustering,” _J. Am. Stat. Assoc._ , vol. 105, no. 490, pp. 713–726, 2010\.
  * [34] X. Chang, Y. Wang, R. Li, and Z. Xu, “Sparse k-means with $\ell_{\infty}/\ell_{0}$ penalty for high-dimensional data clustering,” _Stat. Sin._ , vol. 28, no. 3, pp. 1265–1284, 2018.
  * [35] Z. Li, J. Liu, Y. Yang, X. Zhou, and H. Lu, “Clustering-guided sparse structural learning for unsupervised feature selection,” _IEEE Trans. Knowl. Data Eng._ , vol. 26, no. 9, pp. 2138–2150, 2014.
  * [36] A. A. Pollen _et al._ , “Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex,” _Nat. Biotechnol._ , vol. 32, no. 10, p. 1053, 2014.
  * [37] J. G. Camp _et al._ , “Multilineage communication regulates human liver bud development from pluripotency,” _Nature_ , vol. 546, no. 7659, p. 533, 2017.
  * [38] B. B. Lake, R. Ai, G. E. Kaeser _et al._ , “Neuronal subtypes and diversity revealed by single-nucleus RNA sequencing of the human brain,” _Science_ , vol. 352, no. 6293, pp. 1586–1590, 2016.
  * [39] K. F. Chau, M. L. Shannon, R. M. Fame, E. Fonseca, H. Mullan, M. B. Johnson, A. K. Sendamarai, M. W. Springel, B. Laurent, and M. K. Lehtinen, “Downregulation of ribosome biogenesis during early forebrain development,” _Elife_ , vol. 7, p. e36998, 2018.
  * [40] F. Santacatterina, L. Sánchez-Cenizo _et al._ , “Down-regulation of oxidative phosphorylation in the liver by expression of the atpase inhibitory factor 1 induces a tumor-promoter metabolic state,” _Oncotarget_ , vol. 7, no. 1, p. 490, 2016.
  * [41] R. Ma, W. Zhang, K. Tang _et al._ , “Switch of glycolysis to gluconeogenesis by dexamethasone for treatment of hepatocarcinoma,” _Nat. Commun._ , vol. 4, no. 1, pp. 1–12, 2013.
  * [42] R. Liu, S. Cheng, Y. He, X. Fan, Z. Lin, and Z. Luo, “On the convergence of learning-based iterative methods for nonconvex inverse problems,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 42, no. 12, pp. 3027–3039, 2020.
  * [43] C. Bao, H. Ji, Y. Quan, and Z. Shen, “Dictionary learning for sparse coding: Algorithms and convergence analysis,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 38, no. 7, pp. 1356–1369, 2015.