Completely Positive Factorization by a Riemannian Smoothing Method

Various methods of solving CP factorization problems have been studied. Jarre and Schmallowsky [33] stated a criterion for complete positivity, based on the augmented primal dual method to solve a particular second-order cone problem. Dickinson and Dür [26] dealt with complete positivity of matrices that possess a specific sparsity pattern and proposed a method for finding CP factorizations of these special matrices that can be performed in linear time. Nie [36] formulated the CP factorization problem as an $\mathcal{A}$-truncated $K$-moment problem, for which the author developed an algorithm that solves a series of semidefinite optimization problems. Sponsel and Dür [45] considered the problem of projecting a matrix onto $\mathcal{CP}_{n}$ and $\mathcal{COP}_{n}$ by using polyhedral approximations of these cones. With the help of these projections, they devised a method to compute a CP factorization for any matrix in the interior of $\mathcal{CP}_{n}$. Bomze [7] showed how to construct a CP factorization of an $n\times n$ matrix based on a given CP factorization of an $(n-1)\times(n-1)$ principal submatrix. Dutour Sikirić et al. [43] developed a simplex-like method for a rational CP factorization that works if the input matrix allows a rational CP factorization.

In 2020, Groetzner and Dür [31] applied the alternating projection method to the CP factorization problem by posing it as an equivalent feasibility problem (see (FeasCP)). Shortly afterwards, Chen et al. [19] reformulated the split feasibility problem as a difference-of-convex optimization problem and solved (FeasCP) as a specific application. In fact, we will solve this equivalent feasibility problem (FeasCP) by other means in this paper. In 2021, Boţ and Nguyen [12] proposed a projected gradient method with relaxation and inertia parameters for the CP factorization problem, aimed at solving

where $\mathcal{B}(0,\varepsilon):=\{X\in\mathbb{R}^{n\times r}\mid\|X\|\leq\varepsilon\}$ is the closed ball centered at 0. The authors argued that its optimal value is zero if and only if $A\in\mathcal{CP}_{n}$.

Inspired by the idea of Groetzner and Dür [31], wherein (CPfact) is shown to be equivalent to a feasibility problem called (FeasCP), we treat the problem (FeasCP) as a nonsmooth Riemannian optimization problem and solve it through a general Riemannian smoothing method. Our contributions are summarized as follows:

1. Although it is not explicitly stated in [31], (FeasCP) is actually a Riemannian optimization formulation. We propose a new Riemannian optimization technique and apply it to the problem.

2. In particular, we present a general framework of Riemannian smoothing for the nonsmooth Riemannian optimization problem and show convergence to a stationary point of the original problem.

3. We apply the general framework of Riemannian smoothing to CP factorization. Numerical experiments clarify that our method is competitive with other efficient CP factorization methods, especially for large-scale matrices.

In Section 2, we review the process to reconstruct (CPfact) into another feasibility problem; in particular, we take a different approach to this problem from those in other studies. In Section 3, we describe the general framework of smoothing methods for Riemannian optimization. To apply it to the CP factorization problem, we employ a smoothing function named LogSumExp. Section 4 is a collection of numerical experiments for CP factorization. As a meaningful supplement, in Section 5, we conduct further experiments (FSV problem and robust low-rank matrix completion) to explore the numerical performance of various sub-algorithms and smoothing functions on different applications.

First, let us recall some basic properties of completely positive matrices. Generally, many CP factorizations of a given $A$ may exist, and they may vary in their numbers of columns. This gives rise to the following definitions: the cp-rank of $A\in\mathcal{S}_{n}$, denoted by $\operatorname{cp}(A)$, is defined as

where $\operatorname{cp}(A)=\infty$ if $A\notin\mathcal{CP}_{n}.$ Similarly, we can define the cp-plus-rank as

Immediately, for all $A\in\mathcal{S}_{n}$, we have

Every CP factorization $B$ of $A$ is of the same rank as $A$ since $\operatorname{rank}(XX^{\top})=\operatorname{rank}(X)$ holds for any matrix $X$. The first inequality of (4) comes from the fact that for any CP factorization $B$,

The second is trivial by definition.

Note that computing or estimating the cp-rank of any given $A\in\mathcal{CP}_{n}$ is still an open problem. The following result gives a tight upper bound of the cp-rank for $A\in\mathcal{CP}_{n}$ in terms of the order $n$.

The following result is useful for distinguishing completely positive matrices in either the interior or on the boundary of $\mathcal{CP}_{n}$.

Groetzner and Dür [31] reformulated the CP factorization problem as an equivalent feasibility problem containing an orthogonality constraint.

Given $A\in\mathcal{CP}_{n}$, we can easily get another CP factorization $\widehat{B}$ with $r^{\prime}$ columns for every integer $r^{\prime}\geq r$, if we also have a CP factorization $B$ with $r$ columns. The simplest way to construct such an $n\times r^{\prime}$ matrix $\widehat{B}$ is to append $k:=r^{\prime}-r$ zero columns to $B,$ i.e., $\widehat{B}:=\left[B,0_{n\times k}\right]\geq 0.$ Another way is called column replication, i.e.,

where $b_{i}$ denotes the $i$-th column of $B$. It is easy to see that $\widehat{B}\widehat{B}^{\top}=BB^{\top}=A.$ The next lemma is easily derived from the previous discussion, and it implies that there always exists an $n\times\mathrm{cp}_{n}$ CP factorization for any $A\in\mathcal{CP}_{n}$. Recall that definition of $\mathrm{cp}_{n}$ is given in (5).

Let $\mathcal{O}(r)$ denote the orthogonal group of order $r$, i.e., the set of $r\times r$ orthogonal matrices. The following lemma is essential to our study. (Note that many authors have proved the existence of such an orthogonal matrix $X$ (see, e.g., [17, Lemma 2.1] and [31, Lemma 2.6]).

The next proposition puts the previous two lemmas together.

This proposition tells us that one can find an orthogonal matrix $X$ which can turn a “bad” factorization $\bar{B}$ into a “good” factorization $\bar{B}X$. Let $r\geq\operatorname{cp}(A)$ and $\bar{B}\in\mathbb{R}^{n\times r}$ be an arbitrary (possibly not nonnegative) initial factorization $A=\bar{B}\bar{B}^{\top}$. The task of finding a CP factorization of $A$ can then be formulated as the following feasibility problem,

We should notice that the condition $r\geq\operatorname{cp}(A)$ is necessary; otherwise, (FeasCP) has no solution even if $A\in\mathcal{CP}_{n}.$ Regardless of the exact $\operatorname{cp}(A)$ which is often unknown, one can use $r=\mathrm{cp}_{n}$ in (5). Note that finding an initial matrix $\bar{B}$ is not difficult. Since a completely positive matrix is necessarily positive semidefinite, one can use Cholesky decomposition or spectral decomposition and then extend it to $r$ columns by using (6). The following corollary shows that the feasibility of (FeasCP) is precisely a criterion for complete positivity.

In this study, solving (FeasCP) is the key to finding a CP factorization, but it is still a hard problem because $\mathcal{O}(r)$ is nonconvex.

Groetzner and Dür [31] applied the so-called alternating projections method to (FeasCP). They defined the polyhedral cone, $\mathcal{P}:=\{X\in\mathbb{R}^{r\times r}:\bar{B}X\geq 0\},$ and rewrote (FeasCP) as

The alternating projections method is as follows: choose a starting point $X_{0}\in\mathcal{O}(r)$; then compute $P_{0}=\operatorname{proj}_{\mathcal{P}}(X_{0})$ and $X_{1}=\operatorname{proj}_{\mathcal{O}(r)}(P_{0})$, and iterate this process. Computing the projection onto $\mathcal{P}$ amounts to solving a second-order cone problem (SOCP), while computing the projection onto $\mathcal{O}(r)$ amounts to a singular value decomposition. Note that we need to solve an SOCP alternately at every iteration, which is still expensive in practice. A modified version without convergence involves calculating an approximation of $\operatorname{proj}_{\mathcal{P}}(X_{k})$ by using the Moore-Penrose inverse of $\bar{B}$; for details, see [31, Algorithm 2].

Our way is to use the optimization form. Here, we denote by $\max(\cdot)$ (resp. $\min(\cdot)$) the max-function (resp. min-function)) that selects the largest (resp. smallest) entry of a vector or matrix. Notice that $-\min(\cdot)=\max(-(\cdot)).$ We associate (FeasCP) with the following optimization problem:

For consistency of notation, we turn the maximization problem into a minimization problem:

The feasible set $\mathcal{O}(r)$ is known to be compact [32, Observation 2.1.7]. In accordance with the extreme value theorem [40, Theorem 4.16], (OptCP) attains the global minimum, say $t$. Summarizing these observations together with Corollary 2.6 yields the following proposition.

Let us briefly review some concepts in Riemannian optimization, following the notation of [13]. Throughout this paper, $\mathcal{M}$ will refer to a complete Riemannian submanifold of Euclidean space $\mathbb{R}^{n}$. Thus, $\mathcal{M}$ is endowed with a Riemannian metric induced from the Euclidean inner product, i.e., $\langle\xi,\eta\rangle_{x}:=\xi^{\top}\eta$ for any $\xi,\eta\in{\rm T}_{x}\mathcal{M}$, where ${\rm T}_{x}\mathcal{M}\subseteq\mathbb{R}^{n}$ is the tangent space to $\mathcal{M}$ at $x$. The Riemannian metric induces the usual Euclidean norm $\|\xi\|_{x}:=\|\xi\|=\sqrt{\langle\xi,\xi\rangle_{x}}$ for $\xi\in{\rm T}_{x}\mathcal{M}$. The tangent bundle ${\rm T}\mathcal{M}:=\bigsqcup_{x\in\mathcal{M}}{\rm T}_{x}\mathcal{M}$ is a disjoint union of the tangent spaces of $\mathcal{M}$. Let $f:\mathcal{M}\rightarrow\mathbb{R}$ be a smooth function on $\mathcal{M}$. The Riemannian gradient of $f$ is a vector field $\operatorname{grad}f$ on $\mathcal{M}$ that is uniquely defined by the identities: for all $(x,v)\in\rm{T}\mathcal{M}$,

where $\mathrm{D}f(x):{\rm T}_{x}\mathcal{M}\to{\rm T}_{f(x)}\mathbb{R}\cong\mathbb{R}$ is the differential of $f$ at $x\in\mathcal{M}$. Since $\mathcal{M}$ is an embedded submanifold of $\mathbb{R}^{n}$, we have a simpler statement for $f$ that is also well defined on the whole $\mathbb{R}^{n}$:

where $\nabla f(x)$ is the usual gradient in $\mathbb{R}^{n}$ and $\operatorname{Proj}_{x}$ denotes the orthogonal projector from $\mathbb{R}^{n}$ to ${\rm T}_{x}\mathcal{M}$. For a subset $D\subseteq\mathbb{R}^{n}$, the function $h\in C^{1}(D)$ is smooth, i.e., continuously differentiable on $D$. Given a point $x\in\mathbb{R}^{n}$ and $\delta>0$, $\mathcal{B}(x,\delta)$ denotes a closed ball of radius $\delta$ centered at $x$. $\mathbb{R}_{++}$ denotes the set of positive real numbers. We use subscript notation $x_{i}$ to select the $i$th entry of a vector and superscript notation $x^{k}$ to designate an element in a sequence $\{x^{k}\}$.

Now let us consider the nonsmooth Riemannian optimization problem (NROP):

where $\mathcal{M}\subseteq\mathbb{R}^{n}$ and $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a proper lower semi-continuous function (maybe nonsmooth or even non-Lipschitzian) on $\mathbb{R}^{n}$. For convenience, the term smooth Riemannian optimization problem (SROP) refers to (NROP) when $f(\cdot)$ is continuously differentiable on $\mathbb{R}^{n}$. To avoid confusion in this case, we use $g$ instead of $f$,

Throughout this subsection, we will refer to many of the concepts in [47].

First, let us review the usual concepts and properties related to generalized subdifferentials in $\mathbb{R}^{n}$. For a proper lower semi-continuous function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, the Fréchet subdifferential and the limiting subdifferential of $f$ at $x\in\mathbb{R}^{n}$ are defined as

The definition of $\hat{\partial}f(x)$ above is not the standard one: the standard definition follows [39, 8.3 Definition]. But these definitions are equivalent by [39, 8.5 Proposition]. For locally Lipschitz functions, the Clarke subdifferential at $x\in\mathbb{R}^{n}$, $\partial^{\circ}f(x)$, is the convex hull of the limiting subdifferential. Their relationship is as follows:

Notice that if $f$ is convex, $\partial f(x)$ and $\partial^{\circ}f(x)$ coincide with the classical subdifferential in convex analysis [39, 8.12 Proposition].

Next, we extend our discussion to include generalized subdifferentials of a nonsmooth function on submanifolds $\mathcal{M}$. The Riemannian Fréchet subdifferential and the Riemannian limiting subdifferential of $f$ at $x\in\mathcal{M}$ (see, e.g., [47, Definition 3.1]) are defined as

If $\mathcal{M}=\mathbb{R}^{n}$, the above definitions coincide with the usual Fréchet and limiting subdifferentials in $\mathbb{R}^{n}$. Moreover, it follows directly that, for all $x\in\mathcal{M}$, one has $\hat{\partial}_{\mathcal{R}}f(x)\subseteq\partial_{\mathcal{R}}f(x)$. According to [47, Proposition 3.2], if ${x}$ is a local minimizer of $f$ on $\mathcal{M}$, then $0\in\hat{\partial}_{\mathcal{R}}f({x})$. Thus, we call a point $x\in\mathcal{M}$ a Riemannian limiting stationary point of (NROP) if

In this paper, we will treat it as a necessary condition for a local solution of (NROP) to exist.

The smoothing function is the most important tool of the smoothing method.

Gradient sub-consistency or consistency is crucial to showing that any limit point of the Riemannian smoothing method is also a limiting stationary point of (NROP).

If one substitutes the inclusion with equality in (10), then $\tilde{f}$ satisfies gradient consistency on $\mathbb{R}^{n}$, and similarly in (11) for $\mathcal{M}$. Thanks to the following useful proposition from [47], we can induce gradient sub-consistency on $\mathcal{M}$ from that on $\mathbb{R}^{n}$ if $f$ is locally Lipschitz.

The next example illustrates Riemannian gradient sub-consistency on $\mathcal{M}$ for $\operatorname{lse}(x,\mu)$ in Example 2, since any convex function is locally Lipschitz continuous.

Motivated by the previous papers [18, 35, 47] on smoothing methods and Riemannian manifolds, we propose a general Riemannian smoothing method. Algorithm 1 is the basic framework of this general method.

Now let us describe the convergence properties of the basic method. First, let us assume that the function $\tilde{f}(x,\mu_{k})$ has a minimizer on $\mathcal{M}$ for each value of $\mu_{k}$.

This strong result requires us to find a global minimizer of each subproblem, which, however, cannot always be done. The next result concerns the convergence properties of the sequence $\tilde{f}(x^{k},\mu_{k})$ under the condition that $\tilde{f}$ has the following additional property:

In [18], the authors considered a special case of Algorithm 1, wherein the smoothing function $\tilde{f}(x,\mu)=\sqrt{\mu^{2}+x^{2}}$ of $\lvert x\rvert$ also satisfies (17) and a Riemannian conjugate gradient method is used for (12).

Note that the above weak result does not ensure that $\{f^{k}\}\rightarrow f^{*}$. Next, for better convergence (compared with Theorem 3.5) and an effortless implementation (compared with Theorem 3.4), we propose an enhanced Riemannian smoothing method: Algorithm 2. This is closer to the version in [47], where the authors use the Riemannian steepest descent method for solving the smoothed problem (19).

The following result is adapted from [47, Proposition 4.2 & Theorem 4.3]. Readers are encouraged to refer to [47] for a discussion on the stationary point associated with $\tilde{f}$ on $\mathcal{M}$.