Department of Policy and Planning Sciences
Discussion Paper Series

No.1377

Completely Positive Factorization by Riemannian
Smoothing Method

byZhijian LAI and Akiko YOSHISE

July 2021
Revised September 2022

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANCompletely Positive Factorization by a Riemannian
Smoothing Method
ZHIJIAN LAI∗AKIKO YOSHISE†

Abstract
Copositive optimization is a special case of convex conic programming, and it consistsof optimizing a linear function over the cone of all completely positive matrices under linear constraints. Copositive optimization provides powerful relaxations of NP-hard quadratic
problems or combinatorial problems, but there are still many open problems regarding copositive or completely positive matrices. In this paper, we focus on one such problem; finding acompletely positive (CP) factorization for a given completely positive matrix. We treat it asa nonsmooth Riemannian optimization problem, i.e., a minimization problem of a nonsmoothfunction over a Riemannian manifold. To solve this problem, we present a general smoothing
framework for solving nonsmooth Riemannian optimization problems and show convergence toa stationary point of the original problem. An advantage is that we can implement it quickly
with minimal effort by directly using the existing standard smooth Riemannian solvers, such
as Manopt. Numerical experiments show the efficiency of our method especially for large-scale
CP factorizations.
Keywords: completely positive factorization, smoothing method, nonsmooth Riemannian
optimization problem
AMS: 15A23, 15B48, 90C48, 90C59

1

Introduction

The space of n × n real symmetric matrices Sn is endowed with the trace inner product hA, Bi :=trace(AB). A matrix A ∈ Sn is called completely positive if for some r ∈ N there exists an entrywise
nonnegative matrix B ∈ Rn×r such that A = BB > , and we call B a CP factorization of A. We
define CP n as the set of n × n completely positive matrices, equivalently characterized as
CP n := {BB > ∈ Sn | B is a nonnegative matrix } = conv{xx> | x ∈ Rn+ },where conv(S) denotes the convex hull of a given set S. We denote the set of n × n copositive
matrices by COP n := {A ∈ Sn | x> Ax ≥ 0 for all x ∈ Rn+ }. It is known that COP n and CP n
are duals of each other under the trace inner product; moreover, both CP n and COP n are proper
convex cones [3, Section 2.2]. For any positive integer n, we have the following inclusion relationshipamong other important cones in conic optimization:
CP n ⊆ Sn+ ∩ Nn ⊆ Sn+ ⊆ Sn+ + Nn ⊆ COP n ,where Sn+ is the cone of n × n symmetric positive semidefinite matrices and Nn is the cone of n × nsymmetric nonnegative matrices. See the monograph [3] for a comprehensive description of CP n
and COP n .
Conic optimization is a subfield of convex optimization that studies minimization of linear
functions over proper cones. Here, if the proper cone is CP n or its dual cone COP n , we call theconic optimization problem a copositive programming problem. Copositive programming is closely∗ Graduate School of Systems and Information Engineering, University of Tsukuba, Tsukuba, Ibaraki 305-8573,
Japan. E-mail:s2130117@s.tsukuba.ac.jp† Corresponding author. Faculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba,
Ibaraki 305-8573, Japan. E-mail:yoshise@sk.tsukuba.ac.jp1

CP Factorization by a Riemannian Smoothing Method2

related to many nonconvex, NP-hard quadratic and combinatorial optimizations [29]. For example,
consider the so-called standard quadratic optimization problem,
min{x> M x | e> x = 1, x ∈ Rn+ },(1)

where M ∈ Sn is possibly not positive semidefinite and e is the all-ones vector. Bomze et al. [9]showed that the following completely positive reformulation,min{hM, Xi | hE, Xi = 1, X ∈ CPn },where E is the all-ones matrix, is equivalent to (1). Burer [15] reported a more general result, where
any quadratic problem with binary and continuous variables can be rewritten as a linear program
over CP n . As an application to combinatorial problems, consider the problem of computing theindependence number α(G) of a graph G with n nodes. De Klerk and Pasechnik [24] showed thatα(G) = max{hE, Xi | hA + I, Xi = 1, X ∈ CPn },where A is the adjacency matrix of G. For surveys on applications of copositive programming, see[6, 10, 16, 28, 29].
The difficulty of the above problems lies entirely in the completely positive conic constraint.
Note that because neither COP n nor CP n is self-dual, the primal-dual interior point method forconic optimization does not work as is. Besides this, there are many open problems related tocompletely positive cones. One is checking membership in CP n , which was shown to be NP-hardby [27]. Computing or estimating the cp-rank, as defined later in (3), is also an open problem. We
refer the reader to [4, 28] for a detailed discussion of those unresolved issues.
In this paper, we focus on finding a CP factorization for a given A ∈ CP n , i.e., the CP
factorization problem:
Find B ∈ Rn×r s.t. A = BB > and B ≥ 0,
(CPfact)
which seems to be closely related to the membership problem A ∈ CP n . Sometimes, a matrix isshown to be completely positive through duality, or rather, hA, Xi ≥ 0 for all X ∈ COP n , but inthis case, a CP factorization will not necessarily be obtained.

1.1

Related work on CP factorization

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 solvea particular second-order cone problem. Dickinson and Dür [26] dealt with complete positivityof 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 A-truncated K-moment problem, for which the author developed
an algorithm that solves a series of semidefinite optimization problems. Sponsel and Dür [45]considered the problem of projecting a matrix onto CP n and COP n by using polyhedral approximations of these cones. With the help of these projections, they devised a method to computea CP factorization for any matrix in the interior of CP n . Bomze [7] showed how to construct a
CP factorization of an n × n matrix based on a given CP factorization of an (n − 1) × (n − 1)principal submatrix. Dutour Sikirić 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 Dü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, Boţ and Nguyen [12] proposed a
projected gradient method with relaxation and inertia parameters for the CP factorization problem,
aimed at solvingp(2)
min{kA − XX > k2 | X ∈ Rn×r
∩ B(0, trace(A))},+
X

where B(0, ε) := {X ∈ Rn×r | kXk ≤ ε} is the closed ball centered at 0. The authors argued thatits optimal value is zero if and only if A ∈ CP n .

CP Factorization by a Riemannian Smoothing Method

1.2

3

Our contributions and organization of the paper

Inspired by the idea of Groetzner and Dür [31], wherein (CPfact) is shown to be equivalent toa feasibility problem called (FeasCP), we treat the problem (FeasCP) as a nonsmooth Riemannian optimization problem and solve it through a general Riemannian smoothing method. Ourcontributions 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. Numericalexperiments 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; inparticular, 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 acollection 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.

2

Preliminaries

2.1

cp-rank and cp-plus-rank

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 ∈ Sn , denoted by cp(A), is defined as
cp(A) := min{r ∈ N | A = BB > , B ∈ Rn×r , B ≥ 0},(3)

where cp(A) = ∞ if A ∈/ CP n . Similarly, we can define the cp-plus-rank as
cp+ (A) := min{r ∈ N | A = BB > , B ∈ Rn×r , B > 0}.
Immediately, for all A ∈ Sn , we haverank(A) ≤ cp(A) ≤ cp+ (A).

(4)>

Every CP factorization B of A is of the same rank as A since rank(XX ) = rank(X) holds forany matrix X. The first inequality of (4) comes from the fact that for any CP factorization B,
rank(A) = rank(B) ≤ the number of columns of B.
The second is trivial by definition.
Note that computing or estimating the cp-rank of any given A ∈ CP n is still an open problem.
The following result gives a tight upper bound of the cp-rank for A ∈ CP n in terms of the order n.
Theorem 2.1 (Bomze, Dickinson, and Still [8, Theorem 4.1]). For all A ∈ CPn , we havenfor n ∈ {2, 3, 4}
cp(A) ≤ cpn :=
1
n(n
+
1)
−
4
for n ≥ 5.
2

(5)

The following result is useful for distinguishing completely positive matrices in either the interior
or on the boundary of CP n .
Theorem 2.2 (Dickinson [25, Theorem 3.8]). We haveint(CPn ) = {A ∈ Sn | rank(A) = n, cp+ (A) < ∞}= {A ∈ Sn | rank(A) = n, A = BB > , B ∈ Rn×r , B ≥ 0,
bj > 0 for at least one column bj of B}.

CP Factorization by a Riemannian Smoothing Method

2.2

4

CP factorization as a feasibility problem

Groetzner and Dür [31] reformulated the CP factorization problem as an equivalent feasibility
problem containing an orthogonality constraint.
b with r0 columns for every integer
Given A ∈ CP n , we can easily get another CP factorization B
r0 ≥ r, if we also have a CP factorization B with r columns. The simplest way to construct suchb is to append k := r0 − r zero columns to B, i.e., B
b := [B, 0n×k ] ≥ 0. Anotheran n × r0 matrix B
way is called column replication, i.e.,b := [b1 , . . . , bn−1 , √1 bn , . . . , √1 bn ],
B
m
m
{z
}
|

(6)

m:=r 0 −n+1 columnsbBb > = BB > = A. The next lemma
where bi denotes the i-th column of B. It is easy to see that B
is easily derived from the previous discussion, and it implies that there always exists an n × cpn
CP factorization for any A ∈ CP n . Recall that definition of cpn is given in (5).
Lemma 2.3. Suppose that A ∈ Sn , r ∈ N. Then r ≥ cp(A) if and only if A has a CP factorization
B with r columns.
Let O(r) denote the orthogonal group of order r, i.e., the set of r × r orthogonal matrices. The
following lemma is essential to our study. (Note that many authors have proved the existence ofsuch an orthogonal matrix X (see, e.g., [17, Lemma 2.1] and [31, Lemma 2.6]).
Lemma 2.4. Let B, C ∈ Rn×r . BB > = CC > if and only if there exists X ∈ O(r) with BX = C.
The next proposition puts the previous two lemmas together.
Proposition 2.5. Let A ∈ CP n , r ≥ cp(A), A = B̄ B̄ > , where B̄ ∈ Rn×r may possibly be notnonnegative. Then there exists an orthogonal matrix X ∈ O(r) such that B̄X ≥ 0 and A =(B̄X)(B̄X)> .
This proposition tells us that one can find an orthogonal matrix X which can turn a “bad”
factorization B̄ into a “good” factorization B̄X. Let r ≥ cp(A) and B̄ ∈ Rn×r be an arbitrary(possibly not nonnegative) initial factorization A = B̄ B̄ > . The task of finding a CP factorizationof A can then be formulated as the following feasibility problem,
Find X s.t. B̄X ≥ 0 and X ∈ O(r).

(FeasCP)We should notice that the condition r ≥ cp(A) is necessary; otherwise, (FeasCP) has no solution
even if A ∈ CP n . Regardless of the exact cp(A) which is often unknown, one can use r = cpn in(5). Note that finding an initial matrix B̄ is not difficult. Since a completely positive matrix isnecessarily 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.
Corollary 2.6. Set r ≥ cp(A), B̄ ∈ Rn×r an arbitrary initial factorization of A. Then A ∈ CP n if
and only if (F easCP ) is feasible. In this case, for any feasible solution X, B̄X is a CP factorizationof A.
In this study, solving (FeasCP) is the key to finding a CP factorization, but it is still a hard
problem because O(r) is nonconvex.

2.3

Approaches to solving (FeasCP)

Groetzner and Dür [31] applied the so-called alternating projections method to (FeasCP). They
defined the polyhedral cone, P := {X ∈ Rr×r : B̄X ≥ 0}, and rewrote (FeasCP) as
Find X s.t. X ∈ P ∩ O(r).

CP Factorization by a Riemannian Smoothing Method

5

The alternating projections method is as follows: choose a starting point X0 ∈ O(r); then compute
P0 = projP (X0 ) and X1 = projO(r) (P0 ), and iterate this process. Computing the projection onto
P amounts to solving a second-order cone problem (SOCP), while computing the projection ontoO(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 projP (Xk ) by using the Moore-Penrose inverse of B̄; fordetails, see [31, Algorithm 2].
Our way is to use the optimization form. Here, we denote by max(·) (resp. min(·)) the maxfunction (resp. min-function)) that selects the largest (resp. smallest) entry of a vector or matrix.
Notice that − min(·) = max(−(·)). We associate (FeasCP) with the following optimization problem:max {min (B̄X)}.
X∈O(r)For consistency of notation, we turn the maximization problem into a minimization problem:min {max (−B̄X)}.

(OptCP)

X∈O(r)

The feasible set 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. Summarizingthese observations together with Corollary 2.6 yields the following proposition.
Proposition 2.7. Set r ≥ cp(A), and let B̄ ∈ Rn×r be an arbitrary initial factorization of A.
Then the following statements are equivalent:
1. A ∈ CP n .
2. (FeasCP) is feasible.
3. In (OptCP), there exists a feasible solution X such that max(−B̄X) ≤ 0; alternatively,min(B̄X) ≥ 0.
4. In (OptCP), the global minimum t ≤ 0.

3

Riemannian smoothing method

The problem of minimizing a real-valued function over a Riemannian manifold M, which is called
Riemannian optimization, has been actively studied during the last few decades. In particular, theStiefel manifold,St(n, p) = {X ∈ Rn×p | X > X = I},(when n = p, it reduces to the orthogonal group) is an important case and is our main interest here.
We treat the CP factorization problem, i.e., (OptCP) as a problem of minimizing a nonsmoothfunction over a Riemannian manifold, for which variants of subgradient methods [11], proximal
gradient methods [23], and the alternating direction method of multipliers (ADMM) [34] have beenstudied.
Smoothing methods [21], which use a parameterized smoothing function to approximate theobjective function, are effective on a class of nonsmooth optimizations in Euclidean space. Recently, Zhang, Chen and Ma [47] extended a smoothing steepest descent method to the case of
Riemannian submanifolds in Rn . This is not the first time that smoothing methods have beenstudied on manifolds. Liu and Boumal [35] extended the augmented Lagrangian method and exact
penalty method to the Riemannian case. The latter leads to a nonsmooth Riemannian optimization problem to which they applied smoothing techniques. Cambier and Absil [18] dealt with the
problem of robust low-rank matrix completion by solving a Riemannian optimization problem,
wherein they applied a smoothing conjugate gradient method.
In this section, we propose a general Riemannian smoothing method and apply it to the CP
factorization problem.

CP Factorization by a Riemannian Smoothing Method

3.1

6

Notation and terminology of Riemannian optimization

Let us briefly review some concepts in Riemannian optimization, following the notation of [13].
Throughout this paper, M will refer to a complete Riemannian submanifold of Euclidean space
Rn . Thus, M is endowed with a Riemannian metric induced from the Euclidean inner product,i.e., hξ, ηix := ξ > η for any ξ, η ∈ Tx M, where Tx M ⊆ Rn is the tangent
p space to M at x. Thehξ, ξix for ξ ∈ Tx M. The
Riemannian metric induces
the
usual
Euclidean
norm
kξk
:=
kξk
=
x
F
tangent bundle TM := x∈M Tx M is a disjoint union of the tangent spaces of M. Let f : M → R
be a smooth function on M. The Riemannian gradient of f is a vector field grad f on M that is
uniquely defined by the identities: for all (x, v) ∈ TM,
Df (x)[v] = hv, grad f (x)ix
where Df (x) : Tx M → Tf (x) R ∼= R is the differential of f at x ∈ M. Since M is an embedded
submanifold of Rn , we have a simpler statement for f that is also well defined on the whole Rn :grad f (x) = Projx (∇f (x)),
where ∇f (x) is the usual gradient in Rn and Projx denotes the orthogonal projector from Rn toTx M. For a subset D ⊆ Rn , the function h ∈ C 1 (D) is smooth, i.e., continuously differentiable on
D. Given a point x ∈ Rn and δ > 0, B(x, δ) denotes a closed ball of radius δ centered at x. R++
denotes the set of positive real numbers. We use subscript notation xi to select the ith entry of avector and superscript notation xk to designate an element in a sequence {xk }.

3.2

Ingredients

Now let us consider the nonsmooth Riemannian optimization problem (NROP):
min f (x),x∈M

(NROP)

where M ⊆ Rn and f : Rn → R is a proper lower semi-continuous function (maybe nonsmooth oreven non-Lipschitzian) on Rn . For convenience, the term smooth Riemannian optimization problem(SROP) refers to (NROP) when f (·) is continuously differentiable on Rn . To avoid confusion inthis case, we use g instead of f ,min g(x).
(SROP)
x∈M

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
Rn . For a proper lower semi-continuous function f : Rn → R, the Fréchet subdifferential and thelimiting subdifferential of f at x ∈ Rn are defined asˆ (x) := {∇h(x) | ∃δ > 0 such that h ∈ C 1 (B(x, δ)) and∂ff − h attains a local minimum at x on Rn },ˆ x` , x` , f x` → (x, f (x))}.
∂f (x) := { lim v ` | v ` ∈ ∂f`→∞

ˆ (x) above is not the standard one: the standard definition follows [39, 8.3
The definition of ∂f
Definition]. But these definitions are equivalent by [39, 8.5 Proposition]. For locally Lipschitz
functions, the Clarke subdifferential at x ∈ Rn , ∂ ◦ f (x), is the convex hull of the limiting subdifferential. Their relationship is as follows:ˆ (x) ⊆ ∂f (x) ⊆ ∂ ◦ f (x).
∂fNotice that if f is convex, ∂f (x) and ∂ ◦ f (x) coincide with the classical subdifferential in convex
analysis [39, 8.12 Proposition].
Example 1 (Bagirov, Karmitsa, and Mäkelä [2, Theorem 3.23]). From a result on the pointwise
max-function in convex analysis, we have∂ max(x) = conv{ei | i ∈ I(x)},where ei ’s are the standard bases of Rn and I(x) = {i | xi = max(x)}.

CP Factorization by a Riemannian Smoothing Method

7

Next, we extend our discussion to include generalized subdifferentials of a nonsmooth function on submanifolds M. The Riemannian Fréchet subdifferential and the Riemannian limiting
subdifferential of f at x ∈ M (see, e.g., [47, Definition 3.1]) are defined as∂ˆR f (x) := {grad h(x) | ∃δ > 0 such that h ∈ C 1 (B(x, δ)) andf − h attains a local minimum at x on M},∂R f (x) := { lim v ` | v ` ∈ ∂ˆR f x` , x` , f x` → (x, f (x))}.
`→∞n

If M = R , the above definitions coincide with the usual Fréchet and limiting subdifferentials in
Rn . Moreover, it follows directly that, for all x ∈ M, one has ∂ˆR f (x) ⊆ ∂R f (x). According to[47, Proposition 3.2], if x is a local minimizer of f on M, then 0 ∈ ∂ˆR f (x). Thus, we call a pointx ∈ M a Riemannian limiting stationary point of (NROP) if
0 ∈ ∂R f (x).

(7)

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.
Definition 3.1 (Zhang and Chen [46, Definition 3.1]). A function f˜(·, ·) : Rn × R++ → R is calleda smoothing function of f : Rn → R, if f˜(·, µ) is continuously differentiable in Rn for any µ > 0,lim

z→x,µ↓0f˜(z, µ) = f (x)and there exist a constant κ > 0 and a function ω : R++ → R++ such that|f˜(x, µ) − f (x)| ≤ κω(µ)

withlim ω(µ) = 0.
µ↓0

(8)

Example 2 (Chen, Wets, and Zhang [22, Lemma 4.4]). The LogSumExp function, lse(x, µ) :Rn × R++ → R, given byPnlse(x, µ) := µ log( i=1 exp(xi /µ)),is the smoothing function of max(x) because we can see that(i) lse(·, µ) is smooth on Rn for any µ > 0. Its gradient ∇x lse(x, µ) is given by σ(·, µ) : Rn →
n−1
∆
,
1
[ exp(x1 /µ), · · · , exp(xn /µ) ]> ,(9)∇x lse(x, µ) = σ(x, µ) := Pn`=1 exp(x` /µ)Pn
where ∆n−1 := {x ∈ Rn | i=1 xi = 1, xi ≥ 0} is the unit simplex.
(ii) For all x ∈ Rn and µ > 0, we have max(x) < lse(x, µ) ≤ max(x) + µ log(n). Then, theconstant κ = log(n) and ω(µ) = µ. The above inequalities imply that limz→x,µ↓0 lse(z, µ) =max(x).
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).
Definition 3.2 (Zhang, Chen and Ma [47, Definition 3.4 & 3.9]). A smoothing function f˜ of f issaid to satisfy gradient sub-consistency on Rn if, for any x ∈ Rn ,Gf˜(x) ⊆ ∂f (x),(10)

where the subdifferential of f associated with f˜ at x ∈ Rn is given by
Gf˜(x) := {u ∈ Rn | ∇x f˜ (zk , µk ) → u for some zk → x, µk ↓ 0}.
Similarly, f˜ is said to satisfy Riemannian gradient sub-consistency on M if, for any x ∈ M,
Gf˜,R (x) ⊆ ∂R f (x),where the Riemannian subdifferential of f associated with f˜ at x ∈ M is given by
Gf˜,R (x) = {v ∈ Rn | grad f˜ (zk , µk ) → v for some zk ∈ M, zk → x, µk ↓ 0}.

(11)

CP Factorization by a Riemannian Smoothing Method8

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

Proposition 3.3 (Zhang, Chen and Ma [47, Proposition 3.10]). Let f be a locally Lipschitz function
and f˜ a smoothing function of f. For f˜, if gradient sub-consistency holds on Rn , then Riemannian
gradient sub-consistency holds on M as well.
The next example illustrates Riemannian gradient sub-consistency on M for lse(x, µ) in Example 2, since any convex function is locally Lipschitz continuous.
Example 3 (Chen, Wets, and Zhang [22, Lemma 4.4]). The smoothing function lse(x, µ) of max(x)satisfies gradient consistency on Rn . That is, for any x ∈ Rn ,∂ max(x) = Glse (x) = {lim

xk →x,µk ↓0

σ(xk , µk )}.

Note that the original assertion of [22, Lemma 4.4] is gradient consistency in the Clarke sense, i.e.,∂ ◦ max(x) = Glse (x).

3.3

Riemannian smoothing method

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 thisgeneral method.
Algorithm 1: Basic Riemannian smoothing method for (NROP)
Initialization: Given θ ∈ (0, 1), µ0 > 0, and x−1 ∈ M, select a smoothing function f˜ and a
Riemannian algorithm (called sub-algorithm here) for (SROP).
for k = 0, 1, 2, . . .
Solve
xk = arg min f˜(x, µk )
(12)
x∈M

approximately by using the chosen sub-algorithm, starting at xk−1
if final convergence test is satisfied
stop with approximate solution xk
end ifSet µk+1 = θµk
end forNow let us describe the convergence properties of the basic method. First, let us assume thatthe function f˜(x, µk ) has a minimizer on M for each value of µk .
Theorem 3.4. Suppose that each xk is an exact global minimizer of (12) in Algorithm 1. Thenevery limit point x∗ of the sequence {xk } is a global minimizer of the problem (NROP).
Proof. Let x̄ be a global solution of (NROP), that is,f (x̄) ≤ f (x)for all x ∈ M.

From the Definition 3.1 of the smoothing function, there exist a constant κ > 0 and a function
ω : R++ → R++ such that, for all x ∈ M,
− κω(µ) ≤ f˜(x, µ) − f (x) ≤ κω(µ)(13)

with limµ↓0 ω(µ) = 0. Substituting xk and combining with the global solution x̄, we have thatf˜(xk , µk ) ≥ f (xk ) − κω(µk ) ≥ f (x̄) − κω(µk ).

CP Factorization by a Riemannian Smoothing Method

9

By rearranging this expression, we obtain− κω(µk ) ≤ f˜(xk , µk ) − f (x̄).

(14)

Since xk minimizes f˜(x, µk ) on M for each µk , we have that f˜(xk , µk ) ≤ f˜(x̄, µk ), which leadstof˜(xk , µk ) − f (x̄) ≤ f˜(x̄, µk ) − f (x̄) ≤ κω(µk ).

(15)

The second inequality above follows from (13). Combining (14) and (15), we obtain|f˜(xk , µk ) − f (x̄)| ≤ κω(µk ).

(16)

Now, suppose that x∗ is a limit point of {xk }, so that there is an infinite subsequence K
such that limk∈K xk = x∗ . Note that x∗ ∈ M because M is complete. By taking the limit as
k → ∞, k ∈ K, on both sides of (16), again by the definition of the smoothing function, we obtain|f (x∗ ) − f (x̄)| = lim |f˜(xk , µk ) − f (x̄)| ≤ lim κω(µk ) = 0.
k∈K

k∈K

Thus, it follows that f (x∗ ) = f (x̄). Since x∗ ∈ M is a point whose objective value is equal to thatof the global solution x̄, we conclude that x∗ , too, is a global solution.
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 sequencef˜(xk , µk ) under the condition that f˜ has the following additional property:
0 < µ2 < µ1 =⇒ f˜(x, µ2 ) < f˜(x, µ1 ) for all x ∈ Rn .

(17)

Example 4. The above property holds for lse(x, µ) in Example 2; i.e., we have lse(x, µ2 ) <lse(x, µ1 ) on Rn , provided that 0 < µ2 < µ1 . Note that under the equality,nX

exp(xl /µ) = exp{lse(x, µ)/µ},l=1

the i-th component of σ(x, µ) can be rewritten as
σi (x, µ) = exp{(xi − lse(x, µ))/µ}.
For any fixed x ∈ Rn , consider the derivative of a real function µ → lse(x, ·) : R++ → R. Then we
have
Pn
Pn
xi exp(xi /µ)∇µ lse(x, µ) = lse /µ − i=1=(lse − i=1 xi exp{(xi − lse)/µ})/µµ exp (lse /µ)Pn=(lse − i=1 xi σi )/µ ≤ 0,
where “lse, σ” are shorthand forµ) and σ(x, µ). For the last inequality above, we observe
Plse(x,nn−1that σ∈∆;hence,
theterm
xi=1 i σi is a convex combination of all entries of x, which implies
Pnthat i=1 xi σi ≤ max(x) < lse . This completes the proofs of our claims.
In [18],pthe authors considered a special case of Algorithm 1, wherein the smoothing function˜
f (x, µ) = µ2 + x2 of |x| also satisfies (17) and a Riemannian conjugate gradient method is usedfor (12).
Theorem 3.5. Suppose that f ∗ := inf x∈M f (x) exists and the smoothing function f˜ has property(17). Let f k := f˜(xk , µk ). Then the sequence {f k } generated by Algorithm 1 is strictly decreasing
and bounded below by f ∗ ; hence,
lim |f k − f k−1 | = 0.
k→∞

Proof. For each k ≥ 1, xk is obtained by approximately solving
min f˜(x, µk ),x∈MCP Factorization by a Riemannian Smoothing Method

10starting at xk−1 . Then at least, we havef˜(xk−1 , µk ) ≥ f˜(xk , µk ) = f k .
Since µk = θµk−1 < µk−1 , property (17) ensuresf k−1 = f˜(xk−1 , µk−1 ) > f˜(xk−1 , µk ).
The claim that sequence {f k } is strictly decreasing follows from these two inequalities.
Suppose that, for all µ > 0 and for all x ∈ Rn ,f˜(x, µ) ≥ f (x).

(18)

Then for each k,
f k = f˜(xk , µk ) ≥ f (xk ) ≥ inf f (x) = f ∗ ,x∈Mwhich proves our claims.
Now, we show (18) is true if the smoothing function has property (17). Fix any x ∈ Rn ; (17)implies that f˜(x, ·) is strictly decreasing as µ → 0. Actually, f˜(x, ·) is monotonically increasing onµ > 0. On the other hand, from the definition of the smoothing function, we have thatlim f˜(x, µ) = f (x).
µ↓0Hence, we have inf µ>0 f˜(x, µ) = f (x), as claimed.
Note that the above weak result does not ensure that {f k } → 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).
Algorithm 2: Enhanced Riemannian smoothing method for (NROP)
Initialization: Given θ ∈ (0, 1), µ0 > 0 and a nonnegative sequence {δk } with δk → 0, andx−1 ∈ M, select a smoothing function f˜ and a Riemannian algorithm (called sub-algorithmhere) for (SROP).
for k = 0, 1, 2, . . .
Solve
xk = arg min f˜(x, µk )
(19)
x∈M

approximately by using the chosen sub-algorithm, starting at xk−1 , such thatk grad f˜(xk , µk )k < δk(20)

if final convergence test is satisfied
stop with approximate solution xk
end ifSet µk+1 = θµk
end for
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 f˜ on M.
Theorem 3.6. In Algorithm 2, suppose that the chosen sub-algorithm has the following generalconvergence property for (SROP):
lim inf k grad g(x` )k = 0.
`→∞

(21)

Moreover, suppose that, for all µk , the function f˜(·, µk ) satisfies the convergence assumptions of
the sub-algorithm needed for g above and f˜ satisfies the Riemannian gradient sub-consistency on
M. ThenCP Factorization by a Riemannian Smoothing Method

11

1. For each k, there exists an xk satisfying (20); hence, Algorithm 2 is well defined.
2. Every limit point x∗ of the sequence {xk } generated by Algorithm 2 is a Riemannian limiting
stationary point of (NROP) (see (7)).
Proof. Fix any µk . By (21), we have lim inf `→∞ k grad f˜(x` , µk )k = 0. Hence, there is a convergent
subsequence of k grad f˜(x` , µk )k whose limit is 0. This means that, for any  > 0, there exists an
integer ` such that k grad f˜(x` , µk )k < . If  = δk , we get xk = x` . Thus, statement (1) holds.
Next, suppose that x∗ is a limit point of {xk } generated by Algorithm 2, so that there is an
infinite subsequence K such that limk∈K xk = x∗ . From (1), we have
lim k grad f˜(xk , µk )k ≤ lim δk = 0,

k∈K

k∈K

and we find that grad f˜(xk , µk ) → 0 for k ∈ K, xk ∈ M, xk → x∗ , µk ↓ 0. Hence,
0 ∈ Gf˜,R (x∗ ) ⊆ ∂R f (x∗ ).

Now let us consider the selection strategy of the nonnegative sequence {δk } with δk → 0. In[47], when µk+1 = θµk shrinks, the authors set
δk+1 := ρδk(22)with an initial value of δ0 and constant ρ ∈ (0, 1). In the spirit of the usual smoothing methods
described in [21], one can set
δk := γµk
(23)with a constant γ > 0. The latter is an adaptive rule, because µk determines subproblem (19) andits stopping criterion at the same time. The merits and drawbacks of the two rules require morediscussion, but the latter seems to be more reasonable.
We conclude this section by discussing the connections with [18] and [47]. Our work is basedon an efficient unification of them. [18] focused on a specific algorithm and did not discuss theunderlying generalities, whereas we studied a general framework for Riemannian smoothing. Recall
that the “smoothing function” is the core tool of the smoothing method. In addition to what are
required by its definition (see Definition 3.1), it needs to have the following “additional properties”(AP) in order for the algorithms to converge:(AP1) Approximate from above, i.e., (17). (Needed in Algorithm 1)(AP2) (Riemannian) gradient sub-consistency, i.e., Definition 3.2. (Needed in Algorithm 2)
We find that not all smoothing functions satisfy (AP1) and for some functions it is hard toprove whether (AP2) holds. For example, all the functions in Table 1 are smoothing functions of|x|, but only the first three meet (AP1); the last two do not. In [21], the authors showed that thefirst one in Table 1, f˜1 (x, µ), has property (AP2). The others remain to be verified, but doing so
will not be a trivial exercise. To a certain extent, Algorithm 1 as well as Theorem 3.5 guarantee
a fundamental convergence result even if one has difficulty in showing whether one’s smoothing
function satisfies (AP2). Therefore, it makes sense to consider Algorithms 1 and 2 together for thesake of the completeness of the general framework.
Algorithm 2 expands on the results of [47]. It allows us to use any standard method of (SROP),
not just steepest descent, to solve the smoothed problem (19). Various standard Riemannianalgorithms for (SROP), such as the Riemannian conjugate gradient method [41] (which often
performs better than Riemannian steepest descent), the Riemannian Newton method [1, Chapter6], and the Riemannian trust region method [1, Chapter 7], have extended the concepts andtechniques used in Euclidean space to Riemannian manifolds. As shown by Theorem 3.6, nomatter what kind of sub-algorithm is implemented for (19), it does not affect the final convergenceas long as the chosen sub-algorithm has property (21). On the other hand, we advocate that the
sub-algorithm should be viewed as a “Black Box” and the user should not have to care about thecode details of the sub-algorithm at all. We can directly use an existing solver, e.g., Manopt [14],which includes the standard Riemannian algorithms mentioned above. Hence, we can choose themost suitable sub-algorithm for the application and quickly implement it with minimal effort.

CP Factorization by a Riemannian Smoothing Method

12

Table 1. List of smoothing functions of the absolute value function |x| with κ and ω(µ) in (8)f˜1 (x, µ) =


 |x|
2

 x +µµ4p22˜f2 (x, µ) = µ + x

if|x| > µ2if|x| ≤ µ2xµf˜3 (x, µ) = 2µ log(1 + e ) − x
f˜4 (x, µ) = x tanh( µx ), where tanh(z) is the hyperbolic tangent function.
Rz
2f˜5 (x, µ) = x erf( µx ), where erf(z) := √2π 0 e−t dt is the Gauss error function.

4

κ

ω(µ)

1
4

µ

1

µ

2 log(2)

µ

1

µ

2
√

µ

e π

Numerical experiments on CP factorization

The numerical experiments in Section 4 and 5 were performed on a computer equipped with anIntel Core i7-10700 at 2.90GHz with 16GB of RAM using Matlab R2022a. Our Algorithm 2 is
implemented in the Manopt framework [14] (version 7.0). The number of iterations to solve thesmoothed problem (19) with the sub-algorithm is recorded in the total number of iterations. We
refer readers to the supplementary material of this paper for the available codes.1
In this section, we describe numerical experiments that we conducted on CP factorization inwhich we solved (OptCP) using Algorithm 2, where different Riemannian algorithms were employed
as sub-algorithms and lse(−B̄X, µ) was used as the smoothing function. To be specific, we used
three built-in Riemannian solvers of Manopt 7.0 — steepest descent (SD), conjugate gradient (CG),and trust regions (RTR), denoted by SM SD, SM CG and SM RTR, respectively. We compared
our algorithms with the following non-Riemannian numerical algorithms for CP factorization thatwere mentioned in subsection 1.1. We followed the settings used by the authors in their papers.
• SpFeasDC ls [19]: A difference-of-convex functions approach for solving the split feasibility
problem, it can be applied to (FeasCP). The implementation details regarding the parameters
we used are the same as in the numerical experiments reported in [19, Section 6.1].
• RIPG mod [12]: This is a projected gradient method with relaxation and inertia parameters
for solving (2). As shown in [12, Section 4.2], RIPG mod is the best among the many
strategies of choosing parameters.
• APM mod [31]: A modified alternating projection method for CP factorization; it is describedin Section 2.3.
We have shown that lse(x, µ) is a smoothing function of max(x) with gradient consistency. lse(·, µ)of the matrix argument can be simply derived from entrywise operations. Then from the propertiesof compositions of smoothing functions [5, Proposition 1 (3)], we have that lse(−B̄X, µ) is asmoothing function of max(−B̄X) with gradient consistency. In practice, it is important to avoid
numerical overflow and underflow when evaluating lse(x, µ). Overflow occurs when any xi is large
and underflow occurs when all xi are small. To avoid these problems, we can shift each component
xi by max(x) and use the following formula:Pnlse(x, µ) = µ log( i=1 exp((xi − max(x))/µ)) + max(x),whose validity is easy to show.
The details of the experiments are as follows. If A ∈ CP n was of full rank, for accuracy reasons,we obtained an initial B̄ by using Cholesky decomposition. Otherwise, B̄ was obtained by using
spectral decomposition. Then we extended B̄ to r columns by column replication (6). We set
r = cp(A) if cp(A) was known or r was sufficiently large. We used RandOrthMat.m [42] to generate
a random starting point X 0 on the basis of the Gram-Schmidt process.
1 Alternatively, https://github.com/GALVINLAI/General-Riemannian-Smoothing-Method.

CP Factorization by a Riemannian Smoothing Method

13

For our three algorithms, we set µ0 = 100, θ = 0.8 and used an adaptive rule (23) of δk := γµkwith γ = 0.5. Except for RIPG mod, all the algorithms terminated successfully at iteration k,where min(B̄X k ) ≥ −10−15 was attained before the maximum number of iterations (5,000) was
reached. In addition, SpFeasDC ls failed when L̄k > 1010 . Regarding RIPG mod, it terminatedsuccessfully when kA − Xk Xk> k2 /kAk2 < 10−15 was attained before at most 10,000 iterations forn < 100, and before at most 50,000 iterations in all other cases. In the tables of this section, we
report the rounded success rate (Rate) over the total number of trials, although the definitions of“Rate” in the different experiments (described in Sections 4.1-4.4) vary slightly from one experiment
to the other. We will describe them later.

4.1

Randomly generated instances

We examined the case of randomly generated matrices to see how the methods were affected by theorder n or r. The instances were generated in the same way as in [31, Section 7.7]. We computed
C by setting Cij := |Bij | for all i, j, where B is a random n × 2n matrix based on the Matlab
command randn, and we took A = CC > to be factorized. In Table 2, we set r = 1.5n and r = 3n
for the values n ∈ {20, 30, 40, 100, 200, 400, 600, 800}. For each pair of n and r, we generated 50
instances if n ≤ 100 and 10 instances otherwise. For each instance, we initialized all the algorithms
at the same random starting point X 0 and initial decomposition B̄, except for RIPG mod. Notethat each instance A was assigned only one starting point.
Table 2 lists the average time in seconds (Times ) and the average number of iterations (Iters )among the successful instances. For our three Riemannian algorithms, Iters contains the numberof iterations of the sub-algorithm. Table 2 also lists the rounded success rate (Rate) over the total
number (50 or 10) of instances for each pair of n and r. Boldface highlights the two best resultsin each row.
As shown in Table 2, except for APM mod, each method had a success rate of 1 for all pairs ofn and r. Our three algorithms outperformed the other methods on the large-scale matrices withn ≥ 100. In particular, SM CG with the conjugate-gradient method gave the best results.

4.2

A specifically structured instanceLet en denote the all-ones vector in Rn and consider the matrix [31, Example 7.1],

An =

0
en−1

e>
n−1
In−1

> 

0
en−1

e>
n−1
In−1


∈ CP n .

Theorem 2.2 shows that An ∈ int(CP n ) for every n ≥ 2. By construction, it is obvious thatcp(An ) = n. We tried to factorize An for the values n ∈ {10, 20, 50, 75, 100, 150} in Table 3. For
each An , using r = cp(An ) = n and the same initial decomposition B̄, we tested all the algorithms
on the same 50 randomly generated starting points, except for RIPG mod. Note that each instance
was assigned 50 starting points.
Table 3 lists the average time in seconds (Times ) and the average number of iterations (Iters )among the successful starting points. It also lists the rounded success rate (Rate) over the total
number (50) of starting points for each n. Boldface highlights the two best results for each n.
We can see from Table 3 that the success rates of our three algorithms were always 1, whereas
the success rates of the other methods decreased as n increased. Likewise, SM CG with the
conjugate-gradient method gave the best results.

Methodsn(r = n)
10
20
50
75
100
150

Methodsn(r = 1.5n)
20
30
40
100
200
400
600
800
Methodsn(r = 3n)
20
30
40
100
200
400
600
800

SM SD
Rate Times
1
0.03991
0.04861
0.25991
0.384310.745911.8076SM SD
Rate Times
1
0.04091
0.05491
0.073510.23121
1.0723
1
14.6
1
50.6
1
133.3
SM SD
Rate Times
1
0.05971
0.079310.10351
0.5632
1
4.5548
1
46.5
1
209.0
1
609.7

Iters
71
85
295
329
458
647

Iters
52
59
67
103
163
296
446
628

Iters
49
58
65
104
167
314
474
643
Iters
42
45
48
57
68
89
99
114

Iters
41
44
46
56
69
86
105
109

SM RTR
Rate Times
1
0.05141
0.06901
0.085910.40611
1.9855
1
22.1
1
46.4
1
93.5
SM RTR
Rate Times
1
0.0842
10.1161
1
0.19611
1.4128
1
14.3
1
77.1
1
294.5
1
650.7
Iters
37
39
41
50
65
80
76
83

Iters
35
36
37
45
53
69
80
89

SpFeasDC ls
Rate Times
1
0.00271
0.00751
0.02161
0.283112.2504
1
36.9
1
140.1
1
413.3
SpFeasDC ls
Rate Times
1
0.00571
0.01281
0.02561
0.81151
8.1517
1
124.3
1
981.8
1
4027.4
Iters
15
17
19
86
184
453
795
1070

Iters
24
24
46
109
212
636
882
1225

RIPG modRate Times
1
0.0081
1
0.0231
10.05741
0.816915.2908
1
90.6
1
344.7
1
891.1
RIPG modRate Times
1
0.01051
0.03361
0.08221
1.19091
9.2248
1
156.6
1
616.7
1
1289.4SM CG
Rate Times
1
0.03131
0.04081
0.1073
1
0.19231
0.32891
0.7837
Iters
49
63
101
135
168
241

SM RTR
Rate Times
1
0.04241
0.063710.21041
0.429310.90741
2.6030

Iters
45
55
76
93
108
145

SpFeasDC ls
Rate Times
1
0.00430.98
0.01390.98
0.33890.98
1.07060.80
1.66530.70
3.7652

Iters
149
201
770
1186
1083
1170

RIPG modRate Times
1
0.00740.74
0.021200
00-

Table 3. CP factorization of a family of specifically structured instances.

SM CG
Rate Times
1
0.03941
0.04771
0.06061
0.15201
0.54851
4.1453
1
14.7
1
30.0
SM CG
Rate Times
1
0.0551
10.06731
0.0882
1
0.339512.4116
1
19.2
1
65.7
1
160.4

Table 2. CP factorization of random completely positive matrices.

Iters
2085
3478
-

Iters
1062
1127
1460
4753
9402
17563
25336
26820

Iters
1229
1481
1990
4912
9616
17987
26146
34022

APM mod
Rate Times
0.80
0.01740.90
0.05910.76
0.69480.64
1.48090.60
2.81500.35
9.9930

APM mod
Rate Times
0.32
0.35020.04
1.007500
0
0
00
APM mod
Rate Times
0.30
0.7267
0
0
0
0
0
0
0
-

Iters
616
864
1416
1510
1690
2959

Iters
2198
-

Iters
2318
2467
-

CP Factorization by a Riemannian Smoothing Method
14

CP Factorization by a Riemannian Smoothing Method

4.3

15

An easy instance on the boundary of CP n

Consider the following matrix from [44, Example 2.7]:
41 43 80 56
 43 62 89 78
A= 80 89 162 120
 56 78 120 104
50 51 93 62


50
51 

93 
.
62 
65

The sufficient condition from [44, Theorem 2.5] ensures that this matrix is completely positive andcp(A) = rank(A) = 3. Theorem 2.2 tells us that A ∈ bd(CP 5 ), since rank(A) 6= 5.
We found that all the algorithms could easily factorize this matrix. However, our three algorithms returned a CP factorization B whose smallest entry was as large as possible. In fact, theyalso maximized the smallest entry in the n × r symmetric factorization of A, since (OptCP) isequivalent to{min (X)}.
max
A=XX > ,X∈Rn×r

When we did not terminate as soon as min(B̄X k ) ≥ −10−15 , for example, after 1000 iterations,
our algorithms gave the following CP factorization whose the smallest entry is around 2.8573 −10−15 :
3.5771 4.4766 2.8573
 2.8574 3.0682 6.6650 >
A = BB , where B ≈ 
 8.3822 7.0001 6.5374  .
 5.7515 2.8574 7.9219 
2.8574 6.7741 3.3085

4.4

A hard instance on the boundary of CP n

Next, we examined how well these methods worked on a hard matrix on the boundary of CP n .
Consider the following matrix on the boundary taken from [30]:
8 5 1 1 5
 5 8 5 1 1 



A=
 1 5 8 5 1  ∈ bd(CP 5 ).
 1 1 5 8 5 
5 1 1 5 8
Since A ∈ bd(CP 5 ) and A is of full rank, it follows from Theorem 2.2 that cp+ (A) = ∞; i.e.,there is no strictly positive CP factorization for A. Hence, the global minimum of (OptCP), t = 0,
is clear. None of the algorithms could decompose this matrix under our tolerance, 10−15 , in the
stopping criteria. As was done in [31, Example 7.3], we investigated slight perturbations of thismatrix. Given
1 1 0 0 0 0
 1 0 1 0 0 0 
M M > =: C ∈ int(CP 5 ) with M =  1 0 0 1 0 0 ,
 1 0 0 0 1 0 
1 0 0 0 0 1
we factorized Aλ := λA + (1− λ)C for different values of λ ∈ [0, 1) using r = 12 > cp5 = 11. Notethat Aλ ∈ int(CP5 ) provided 0 ≤ λ < 1 and Aλ approached the boundary as λ → 1. We chosethe largest λ = 0.9999. For each Aλ , we tested all of the algorithms on 50 randomly generated
starting points and computed the success rate over the total number of starting points.
Table 4 shows how the success rate of each algorithm changes as Aλ approaches the boundary.
The table sorts the results from left to right according to overall performance. Except for SM RTR,
whose success rate was always 1, the success rates of all the other algorithms significantly decreased
as λ increased to 0.9999. Surprisingly, the method of SM CG, which performed well in the previous
examples, seemed unable to handle instances close to the boundary.

CP Factorization by a Riemannian Smoothing Method

16

Table 4. Success rate of CP factorization of Aλ for values of λ from 0.6 to 0.9999.
λ
0.6
0.65
0.7
0.75
0.8
0.82
0.84
0.86
0.88
0.9
0.91
0.92
0.93
0.94
0.95
0.96
0.97
0.98
0.99
0.999
0.99995SM RTR
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

SM SD
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0.68
0
0

RIPG mod
1
1
1
1
1
1
1
1
1
1
1
1
0.98
0.98
1
0.96
0.86
0.76
0.42
0.14
0

SM CG
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0.98
0.82
0.28
0
0
0

SpFeasDC ls
1
1
1
1
0.96
0.98
0.86
0.82
0.58
0.48
0.4
0.2
0.22
0.1
0.12
0.06
0.06
0.02
0
0
0

APM mod
0.42
0.44
0.48
0.52
0.46
0.4
0.24
0.1
0.18
0.18
0.14
0.18
0.22
0.2
0.32
0.34
0.14
0
0
0
0

Further numerical experiments: comparison with [18, 47]

As described at the end of Section 3, the algorithms in [47] and [18] are both special cases of ouralgorithm. In this section, we compare them to show whether it performs better when we use othersub-algorithms or other smoothing functions. We applied Algorithm 2 to two problems: findinga sparse vector (FSV) in a subspace and robust low-rank matrix completion, which are problems
implemented in [47] and [18], respectively. Since they both involve approximations to the `1 norm,
we applied the smoothing functions listed in Table 1.
We used the six solvers built into Manopt 7.0, namely, steepest descent; Barzilai-Borwein (i.e.,
gradient-descent with BB step size); Conjugate gradient; trust regions; BFGS (a limited-memory
version); ARC (i.e., adaptive regularization by cubics).

5.1

FSV problem

The FSV problem is to find the sparsest vector in an n-dimensional linear subspace W ⊆ Rm ;it has applications in robust subspace recovery, dictionary learning, and many other problems inmachine learning and signal processing [37, 38]. Let Q ∈ Rm×n denote a matrix whose columns
form an orthonormal basis of W : this problem can be formulated as
min kQxk0 ,x∈S n−1

where S n−1 := {x ∈ Rn | kxk = 1} is the sphere manifold, and kzk0 counts the number of nonzero
components of z. Because this discontinuous objective function is unwieldy, in the literature, one
instead focuses on solving the `1 norm relaxation given below:
min kQxk1 ,x∈S n−1(24)

P
where kzk1 := i |zi | is the `1 norm of the vector z.
Our synthetic problems of the `1 minimization model (24) were generated in the same way as
in [47]: i.e., we chose m ∈ {4n, 6n, 8n, 10n} for n = 5 and m ∈ {6n, 8n, 10n, 12n} for n = 10. We
defined a sparse vector en := (1, . . . , 1, 0, . . . , 0)> ∈ Rm , whose first n components are 1 and theremaining components are 0. Let the subspace W be the span of en and some n − 1 random vectors
in Rm . By mgson.m [20], we generated an orthonormal basis of W to form a matrix Q ∈ Rm×n .

CP Factorization by a Riemannian Smoothing Method

17

Table 5. Number of successes from 50 pairs of random instances and random initial points for the`1 minimization model (24) and n = 5.
(n, m)

τ(5, 20)−5(5, 30)(5, 40)(5, 50)

10
10−6
10−7
10−8
10−5
10−6
10−7
10−8
10−5
10−6
10−7
10−8
10−5
10−6
10−7
10−8

(n, m)

τ(5, 20)−9(5, 30)(5, 40)(5, 50)

10
10−10
10−11
10−12
10−9
10−10
10−11
10−12
10−9
10−10
10−11
10−12
10−9
10−10
10−11
10−12

Steepestdescent
21
21
21
16
36
36
36
34
44
44
44
43
47
47
47
46
Steepestdescent
0
0
0
0
8
1
0
0
3
2
1
0
5
2
0
0

BarzilaiBorwein
19
19
19
19
42
42
42
42
47
47
47
47
47
47
47
47
BarzilaiBorwein
19
19
19
18
42
42
42
42
47
47
47
46
47
47
47
47

Conjugategradient
0
0
0
0
0
0
0
0
1
0
0
0
2
2
0
0
Conjugategradient
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

Trustregions
22
22
22
22
34
34
34
34
44
44
44
44
45
45
45
45
Trustregions
22
22
22
22
34
34
34
34
44
44
44
44
45
45
45
45

BFGS

ARC

23
23
23
23
36
36
36
36
47
47
47
47
45
45
45
45

23
23
23
23
35
35
35
35
45
45
45
45
45
45
45
45

BFGS

ARC

23
23
23
22
36
36
36
34
47
47
47
44
45
45
45
45

23
23
19
17
35
35
33
29
45
45
44
36
45
45
45
37

With this construction, the minimum value of kQxk0 should be equal to n. We chose the initialpoints by using the M.rand() tool of Manopt 7.0 that returns a random point on the manifold M and
set x0 = abs(M.rand()). The nonnegative initial point seemed to be better in the experiment.
Regarding the the settings of our Algorithm 2, we chose the same smoothing function f˜1 (x, µ) in
Table 1 and the same gradient tolerance strategy (22) as in [47]: µ0 = 1, θ = 0.5, δ0 = 0.1, ρ = 0.5.
We compared the numerical performances when using different sub-algorithms. Note that withthe choice of the steepest-descent method, our Algorithm 2 is exactly the same as the one in [47].
For each (n, m), we generated 50 pairs of random instances and random initial points. We
claim that an algorithm successfully terminates if kQxk k0 = n, where xk is the k-th iteration.
Here, when we count the number of nonzeros of Qxk , we truncated the entries as(Qxk )i = 0

if |(Qxk )i | < τ,

(25)

where τ > 0 is a tolerance related to the precision of the solution, taking values from 10−5 to10−12 . Tables 5 and 6 report the number of successful cases out of 50 cases. Boldface highlights
the best result for each row.
As shown in Table 5 and 6, surprisingly, the conjugate-gradient method, which performed beston the CP factorization problem in Section 4, performed worst on the FSV problem. In fact, it
was almost useless. Moreover, although the steepest-descent method employed in [47] was not bad

CP Factorization by a Riemannian Smoothing Method

18

Table 6. Number of successes from 50 pairs of random instances and random initial points for the`1 minimization model (24) and n = 10.
(n, m)

τ(10, 60)−5(10, 80)(10, 100)(10, 120)

10
10−6
10−7
10−8
10−5
10−6
10−7
10−8
10−5
10−6
10−7
10−8
10−5
10−6
10−7
10−8

(n, m)

τ(10, 60)−9(10, 80)(10, 100)(10, 120)

10
10−10
10−11
10−12
10−9
10−10
10−11
10−12
10−9
10−10
10−11
10−12
10−9
10−10
10−11
10−12

Steepestdescent
24
24
24
23
39
39
39
39
45
45
45
45
44
44
44
44
Steepestdescent
3
0
0
0
5
0
0
0
13
2
0
0
14
0
0
0

BarzilaiBorwein
28
28
28
28
37
37
37
37
48
48
48
48
46
46
46
46
BarzilaiBorwein
28
28
28
28
37
37
37
37
48
48
48
48
46
46
46
46

Conjugategradient
0
0
0
0
1
0
0
0
3
0
0
0
1
0
0
0
Conjugategradient
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

Trustregions
28
28
28
28
40
40
40
40
45
45
45
45
44
44
44
44
Trustregions
28
28
28
28
40
40
40
40
45
45
45
45
44
44
44
44

BFGS

ARC

28
28
28
28
39
39
39
39
43
43
43
43
44
44
44
44

25
25
25
25
40
40
40
40
41
41
41
41
43
43
43
43

BFGS

ARC

28
28
28
27
39
39
39
37
43
43
43
43
44
44
44
43

25
25
22
12
40
40
39
30
41
41
40
37
43
43
43
40

at obtaining low-precision solutions with τ ∈ {10−5 , 10−6 , 10−7 , 10−8 }, it had difficulty obtaining
high-precision solutions with τ ∈ {10−9 , 10−10 , 10−11 , 10−12 }. The remaining four sub-algorithms
easily obtained high-precision solutions, with the Barzilai-Borwein method performing the best in
most occasions. Combined with the results in Section 4, this shows that in practice, the choice of
sub-algorithm in the Riemannian smoothing method (Algorithm 2) is highly problem-dependent.
For the other smoothing functions in Table 1, we obtained similar results as in Table 5 and 6.

5.2

Robust low-rank matrix completion

Low-rank matrix completion consists of recovering a rank-r matrix M of size m × n from only afraction of its entries with r  min(m, n). The situation in robust low-rank matrix completion is
one where only a few observed entries, called outliers, are perturbed, i.e.,
M = M0 + S,
where M0 is the unperturbed original data matrix of rank r and S is a sparse matrix. This is acase of adding non-Gaussian noise for which the traditional `2 minimization model,
min kPΩ (X − M )k2

X∈MrCP Factorization by a Riemannian Smoothing Method

-9

10

0

20

40

60

80

-9

10 -12

100

0

20

40

80

10 -12

100

40

60

40

80

10 -6

10 -12100050(e) f˜5

100

10 -3

10 -12

150

10

0

50

100

10 -9

0

60

80

100

10 0

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -3
10 -6

150

10 -12

0

50

100

150

Time [s](h) f˜3steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -3

RMSE

RMSE

40

10 -9

(g) f˜2

10 -6

10 -1220(d) f˜4Time [s]steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc-30Iteration

10 -6(f) f˜1
0

10 -12100

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

Time [s]

10

80

10 -9

Iteration

1060(c) f˜3

10 -9

10 -9

20

20

10 0

RMSE

RMSE

10 -6

0

0

Iteration

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -3

10 -6
10 -9

(b) f˜2
10 0steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -3

RMSE

60

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -3

-9

Iteration(a) f˜1
0

10 -1210 -6

10Iteration

10

10 -3

RMSE

10 -12

10 -6

10 0

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

RMSE

10 -6

10 0

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -3

RMSE

RMSE

10 -3

10

10 0

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

RMSE

10 0

19

10 -6
10 -9

0

50

100

150

10 -12050Time [s](i) f˜4100

150Time [s](j) f˜5

Figure 1. Perfect low-rank matrix completion of a rank-10 5000 × 5000 matrix without any outliers
using different smoothing functions in Table 1. (a)–(e) comprise the running iteration comparison;(f)–(j) comprise the time comparison.
is not well suited to recovery of M0 . Here, Mr := {X ∈ Rm×n | rank(X) = r} is a fixed rank
manifold, Ω denotes the set of indices of observed entries, and PΩ : Rm×n → Rm×n is the projection
onto Ω, defined as(
Zij if (i, j) ∈ ΩPΩ
Zij 7−→0if (i, j) ∈/ Ω.
In [18], the authors try to solve
min kPΩ (X − M )k1 ,X∈Mr

because the sparsity-inducing property of the `1 norm leads one to expect exact recovery when thenoise consists of just a few outliers.
In all of the experiments, the problems were generated in the same way as in [18]. In particular,
after picking the values of m, n, r, we generated the ground truth U ∈ Rm×r , V ∈ Rn×r withindependent and identically distributed (i.i.d.) Gaussian entries of zero mean and unit varianceand M := U V > . We then sampled k := ρr(m + n − r) observed entries uniformly at random,
where ρ is the oversampling factor. In our experiments, we set ρ = 5 throughout. We chose theinitial points X0 by using the rank-r truncated SVD decomposition of PΩ (M ).
Regarding the setting of our Algorithm 2, we tested all combinations of the five smoothing
functions in Table 1 and six sub-algorithms mentioned before (30 cases in total). We set µ0 = 100
and chose an aggressive value of θ = 0.05 for reducing µ, as in [18]. The stopping criterion of theloop of the sub-algorithm was set to a maximum of 40 iterations or the gradient tolerance (23),whichever was reached first. We monitored the iterations Xk through the root mean square error(RMSE), which is defined as the error on all the entries between Xk and the original low-rankmatrix M0 , i.e.,
s
Pm Pn
2
i=1
j=1 (Xk,ij − M0,ij ).
RMSE (Xk , M0 ) :=mn

CP Factorization by a Riemannian Smoothing Method

5.2.1

20

Perfect low-rank matrix completion

As in [18], we first tested all the methods on a simple perfect matrix M (without any outliers) of
size 5000 × 5000 and rank 10. The results are shown in Figure 1. We can see that the choice ofsmoothing function does not have much effect on numerical performance. In terms of the numberof iterations ((a)–(e)), our Algorithm 2 inherits the convergence of its sub-algorithm at least Qsuperlinearly when trust regions or ARC are used. But the single iteration cost of trust regions and
ARC is high; they are not efficient in terms of time. Specifically, the conjugate-gradient method

-4

10 -6

10 -6

150

10 -8

200

0

50

(a) f˜1
10 -2

-6

50

100

150

10 -8

200

5

10

10 -4

10 -4

10 -6
-8

10

15

10

20

0

5

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -6

50

100

150

10 -8

200

-6

0

50

0

5

10

Time [s](q) f˜2

100

15

20

10 -2

150

10

Time [s](r) f˜3

15

20

15

20

(h) f˜3
10steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -2
10 -4

150

10 -8

200

0

50

100

150

200

Iteration(l) f˜2steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0
10 -2
10 -4
10 -6

0

50

100

150

10 -8

200

0

5

10 -2

10

15

20

Time [s](p) f˜1steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

10 -4
-6

10 -8

10

0

(o) f˜510

55Iteration

-6

0

0

Time [s]

10 -4

10 -8

200

10 -4

10 -8

10 -8

20

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

10 -2

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

10

100

10 0

RMSE

10 -4

10 -8

50

(n) f˜4

RMSE

10 -2

0

Iteration

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

15

10 -6

(m) f˜3
10 0

10(k) f˜1

10 -6

0

10 -4

Iteration

10 -4

Iteration

10

10

20

RMSE

10 -2

RMSE

RMSE

10 -4

10 -8

15

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

10 -2

10 -6

(j) f˜5(i) f˜4
10 -2

10

200

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

10 0

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

Time [s]Time [s]

10 0

5

10 -4

-8

5

0

10 -2

10 -6

150

(d) f˜4(g) f˜2
10

100

10 -6

0

-8

0

50

Time [s]

10 -6
10

RMSE

10 -8

20

RMSE

10 -2

RMSE

RMSE

10 -2

15

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

00Iteration-6(f) f˜1steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -8

200

10 -4

Time [s](e) f˜5
10

10

150

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -2

10

0

100

10 0

-6

Iteration

050(c) f˜3

10 -4
10

0

0

Iteration

RMSE

10 -4

10 -8

10 -8

200

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

RMSE

RMSE

10 -2

150

(b) f˜2steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

100

Iteration

RMSE

100

10 -4
10 -6

RMSE

50

10 -2

RMSE

0

10

-4

10 -6

Iteration

10

10 -2

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

10 0

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0
10 -2

RMSE

10 -8

10

-4

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

RMSE

10

10 -2

RMSE

RMSE

10 -2

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

RMSE

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

10 -4
10 -6

0

5

10

Time [s](s) f˜4

15

20

10 -8

0

5

10

15

20

Time [s](t) f˜5

Figure 2. Low-rank matrix completion with outliers for two rank-10 500 × 500 matrices by using
different smoothing functions in Table 1. (a)–(j) corresponds to one matrix with outliers createdby using µN = σN = 0.1, while (k)–(t) corresponds to the other with outliers created by usingµN = σN = 1. (a)–(e) and (k)–(o) comprise the running iteration comparison; (f)–(j) and (p)–(t)comprise the time comparison.

CP Factorization by a Riemannian Smoothing Method

21

employed in [18] stagnates at lower precision. Overall, Barzilai-Borwein performed best in termsof time and accuracy.
5.2.2

Low-rank matrix completion with outliers

Given a 500 × 500 matrix for which we observed the entries uniformly at random with an oversampling ρ of 5, we perturbed 5% of the observed entries by adding noise to them in order to create
2
outliers. The added item was a random variable defined as O = S±1 · N (µN , σN
) where S±1 is a2
random variable with equal probability of being equal to +1 or −1, while N (µN , σN
) is a Gaussian

10 -4
10 -6

10 -6

150

10 -8

200

0

50

(a) f˜1

10 -4

10 -4

10 -8

-8

10

200

0

50

(e) f˜5steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10

-4

10 -6
10 -8

10

50

100

150

200

10 -8

250

0

50

-6

100

150

10 -8

200

50

-6

100

150

150

200

250

10 -8

250

100

50

150

200

250

150

200

250

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0
10 -2
10 -4

10 -8

200

0

50

100

150

200

Iteration(l) f˜2steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

10steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

0

10 -2
10 -4
10 -6

0

50

100

150

10 -8

200

0

50

10 -2

100

150

200

250

Time [s](p) f˜1steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

0

10 -4
-6

10 -8

150

(h) f˜3

-6

10

100

Time [s]

10 -4

10

50

0

(o) f˜5

-6

0

10 -4

Iteration

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

0

10 -8

100

0

10 -8

200

10 -4
10

50

50

10 -2

RMSE

10 -4

0

100

10 -2

RMSE

10 -2

10 -8

0

(n) f˜41010 -2(k) f˜1

100steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

Iteration

10

200

10 -6

Iteration

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

0

200

-4

10 -8

250

-6

(m) f˜3

10200

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

0

Iteration

10

150

10 -4
10

50

100

150

10 0

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 -2

RMSE

10 -4

150

10 0

10

100

(d) f˜4

10 -6

10 -2

RMSE

10 -2

10

100

(g) f˜2Time [s]steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

050(j) f˜5

0

10 -8

0

(f) f˜1-4(i) f˜4

10

10

250

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

Time [s]

10

200

10 -6

0

50

10 -6

Time [s]

10 -2

RMSE

RMSE

10 -2

150

10 0

0

Iteration

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgsarc

Time [s]

RMSE

10 0

100

10 -8

200

10 -4
10 -6

150

150

10 -2

-8

100

100

10 0

10 -6

50

50

(c) f˜3

10 -6

0

0

Iteration

RMSE

10 -2

Iteration

RMSE

10 -8

200

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

RMSE

RMSE

10 -2

150

(b) f˜2steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

100

Iteration

RMSE

100

10 -4
10 -6

RMSE

50

10 -4

RMSE

0

10 -2

10 -6

Iteration

RMSE

10 -2

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

10

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

0

10 -2

RMSE

10 -8

10 -4

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

RMSE

10 -2

RMSE

RMSE

10 -2

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

RMSE

steepestdescentbarzilaiborwein
conjugategradienttrustregions
rlbfgs
arc

10 0

10 -4
10 -6

0

50

100

150

200

250

10 -8

0

50

100

150

Time [s]Time [s]Time [s]Time [s](q) f˜2(r) f˜3(s) f˜4(t) f˜5200

250

Figure 3. Low-rank matrix completion with outliers for two rank-10 5000 × 5000 matrices by using
different smoothing functions in Table 1. (a)–(j) corresponds to one matrix with outliers createdby using µN = σN = 0.1, while (k)–(t) corresponds to the other with outliers created by usingµN = σN = 1. (a)–(e) and (k)–(o) comprise the running iteration comparison; (f)–(j) and (p)–(t)comprise the time comparison.

CP Factorization by a Riemannian Smoothing Method222random variable of mean µN and variance σN.
Figure 2 reports the results of two 500 × 500 instances with outliers generated using µN =σN = 0.1 and µN = σN = 1. Again, we can see that the choice of smoothing function does nothave much effect. In most cases, BFGS and trust regions were better than the other methods interms of number of iterations, and BFGS was the fastest. Furthermore, the conjugate-gradientmethod employed in [18] still stagnated at solutions with lower precision, approximately 10−6 , while
steepest descent, BFGS, and trust regions always obtained solutions with at least 10−8 precision.
Next, we ran the same experiment on larger 5000 × 5000 matrices, with 5% outliers. Figure 3
illustrates the results of these experiments, with µN = σN = 0.1 and µN = σN = 1. In most cases,
trust regions still outperformed the other methods in terms of number of iterations, while BFGS
performed poorly. Barzilai-Borwein and the conjugate-gradient method were almost as good interms of time.

6

Concluding remarks

We examined the problem of finding a CP factorization of a given completely positive matrix andtreated it as a nonsmooth Riemannian optimization problem. To this end, we studied a generalframework of Riemannian smoothing for Riemannian optimization. The numerical experiments
clarified that our method can compete with other efficient CP factorization methods, in particularon large-scale matrices.
Let us we summarize the relation of our approach to the existing CP factorization methods.
Groetzner and Dür [31] and Chen et al. [19] proposed different methods to solve (FeasCP). Boţand Nguyen [12] tried to solve another model (2). However, the methods they used do not belong
to the Riemannian optimization techniques, but are rather Euclidean ones, since they treated theset O(r) := {X ∈ Rr×r : X > X = I} as a usual constraint in Euclidean space. By comparison, we
recognize the existence of manifolds, namely, the Stiefel manifold M = O(r), and use optimization techniques specific to them. This change in perspective suggests the possibility of using therich variety of Riemannian optimization techniques. As the experiments in Section 4 show, ourRiemannian approach is faster and more reliable than the Euclidean methods.
In the future, we plan to extend Algorithm 2 to the case of general manifolds and, particularly,
to quotient manifolds. This application is believed to be possible, although effort should be put
into establishing analogous convergence results. In fact, convergence has been verified in a builtin example in Manopt 7.0 [14]: robust pca.m computes a robust version of PCA on data andoptimizes a nonsmooth function over a Grassmann manifold. The nonsmooth term consists of thel2 norm, which is not squared, for robustness. In robust pca.m, Riemannian smoothing with apseudo-Huber loss function is used in place of the l2 norm.
As in the other numerical methods, there is no guarantee that Algorithm 2 will find a CP
factorization for every A ∈ CP n . It follows from Proposition 2.7 that A ∈ CPn if and only if the
global minimum of (OptCP), say t, is such that t 6 0. Since our methods only converge to astationary point, Algorithm 2 provides us with a local minimizer at best. We are looking forward
to finding a global minimizer of (OptCP) in our future work.

Acknowledgments
This work was supported by the Research Institute for Mathematical Sciences, an International
Joint Usage/Research Center, at Kyoto University, JSPS KAKENHI Grant, number (B)19H02373,
and JST SPRING Grant, number JPMJSP2124. The authors would like to sincerely thank theanonymous reviewers and the coordinating editor for their thoughtful and valuable comments whichhave significantly improved the paper.

Data availability
The data that support the findings of this study are available from the corresponding author uponrequest.

CP Factorization by a Riemannian Smoothing Method

23

Conflict of interest
All authors declare that they have no conflicts of interest.

References[1] Absil, P.-A., Mahony, R., Sepulchre, R.: Optimization Algorithms on Matrix Manifolds.
Princeton University Press, Princeton (2009)[2] Bagirov, A., Karmitsa, N., Mäkelä, MM.: Introduction to Nonsmooth Optimization: Theory,
Practice and Software. Springer International Publishing, Berlin (2014)[3] Berman, A., Shaked-Monderer, N.: Completely Positive Matrices. World Scientific, Singapore(2003)[4] Berman, A., Dür, M., Shaked-Monderer, N.: Open problems in the theory of completely
positive and copositive matrices. Electron. J. Linear Algebra. 29, 46–58 (2015)[5] Bian, W., Chen, X.: Neural network for nonsmooth, nonconvex constrained minimization via
smooth approximation. IEEE Trans. Neural Netw. Learn. Syst. 25, 545–556 (2013)[6] Bomze, I.M.: Copositive optimization–recent developments and applications. Eur. J. Oper.
Res. 216, 509–520 (2012)[7] Bomze, I.M.: Building a completely positive factorization. Cent. Eur. J. Oper. Res. 26, 287–
305 (2018)[8] Bomze, I.M., Dickinson, P.J., Still, G.: The structure of completely positive matrices according
to their cp-rank and cp-plus-rank. Linear Algebra Appl. 482, 191–206 (2015)[9] Bomze, I.M., Dür, M., De Klerk, E., Roos, C., Quist, A.J., Terlaky, T.: On copositive programming and standard quadratic optimization problems. J. Glob. Optim. 18, 301–320 (2000)[10] Bomze, I.M., Schachinger, W., Uchida, G.: Think co(mpletely) positive! Matrix properties,
examples and a clustered bibliography on copositive optimization. J. Glob. Optim. 52, 423–445(2012)[11] Borckmans, P.B., Selvan, S.E., Boumal, N., Absil, P.-A.: A Riemannian subgradient algorithm
for economic dispatch with valve-point effect, J. Comput. Appl. Math. 255, 848–866 (2014)[12] Boţ, R. I., Nguyen, D.-K.: Factorization of completely positive matrices using iterative projected gradient steps. Numer. Linear Algebra Appl. 28, e2391 (2021)[13] Boumal, N.: An Introduction to Optimization on Smooth Manifolds. To appear with Cambridge University Press, Mar 2022.
[14] Boumal, N., Mishra, B., Absil, P.-A., Sepulchre, R.: Manopt, a Matlab toolbox for optimization on manifolds. J. Mach. Learn. Res. 15, 1455–1459 (2014)[15] Burer, S.: On the copositive representation of binary and continuous nonconvex quadratic
programs. Math. Program. 120, 479–495 (2009)[16] Burer, S.: A gentle, geometric introduction to copositive optimization. Math. Program. 151,
89–116 (2015)[17] Burer, S., Monteiro, R.D.C.: Local minima and convergence in low-rank semidefinite programming. Math. Program. 103, 427–444 (2005)[18] Cambier, L., Absil, P.-A.: Robust low-rank matrix completion by Riemannian optimization.
SIAM J. Sci. Comput. 38, S440–S460 (2016)

CP Factorization by a Riemannian Smoothing Method

24[19] Chen, C., Pong, T.K., Tan, L., Zeng, L.: A difference-of-convex approach for split feasibility
with applications to matrix factorizations and outlier detection. J. Glob. Optim. 78, 107–136(2020)[20] Chen, M.: Gram-Schmidt orthogonalization. MATLAB Central File Exchange. https://www.
mathworks.com/matlabcentral/fileexchange/55881-gram-schmidt-orthogonalization(2022). Accessed 24 June 2022
[21] Chen, X.: Smoothing methods for nonsmooth, nonconvex minimization. Math. Program. 134,
71–99 (2012)[22] Chen, X., Wets, R.J.B., Zhang, Y.: Stochastic variational inequalities: residual minimizationsmoothing sample average approximations. SIAM J. Optim. 22, 649–673 (2012)[23] De Carvalho Bento, G., da Cruz Neto, J.X., Oliveira, P.R.: A new approach to the proximal
point method: convergence on general Riemannian manifolds. J. Optim. Theory Appl. 168,
743–755 (2016)[24] De Klerk, E., Pasechnik, D.V.: Approximation of the stability number of a graph via copositive
programming. SIAM J. Optim. 12, 875–892 (2002),[25] Dickinson, P.J.: An improved characterization of the interior of the completely positive cone.
Electron. J. Linear Algebra. 20, 723–729 (2010)[26] Dickinson, P.J., Dür, M.: Linear-time complete positivity detection and decomposition ofsparse matrices. SIAM J. Matrix Anal. Appl. 33, 701–720 (2012)[27] Dickinson, P.J., Gijben, L.: On the computational complexity of membership problems forthe completely positive cone and its dual. Comput. Optim. Appl. 57, 403–415 (2014)[28] Dür, M.: Copositive programming — a survey. In: Diehl, M., Glineur, F., Jarlebring, E.,
Michiels W. (eds.) Recent advances in optimization and its applications in engineering, pp. 3–
20. Springer, Berlin, Heidelberg (2010).
[29] Dür, M., Rendl, F.: Conic optimization: A survey with special focus on copositive optimization
and binary quadratic problems. EURO J. Comput. Optim. 9, 100021 (2021)[30] Dür, M., Still, G.: Interior points of the completely positive cone. Electron. J. Linear Algebra.
17, 48-53 (2008).
[31] Groetzner, P., Dür, M.: A factorization method for completely positive matrices. Linear
Algebra Appl. 591, 1–24 (2020)[32] Horn, R.A., Johnson, C.R.: Matrix Analysis. Cambridge University Press, Cambridge (2012)[33] Jarre, F., Schmallowsky, K.: On the computation of C ∗ certificates. J. Glob. Optim. 45,
281–296 (2009)[34] Kovnatsky, A., Glashoff, K., Bronstein, M.M.: MADMM: a generic algorithm for non-smooth
optimization on manifolds. In: European Conference on Computer Vision, pp. 680–696.
Springer, Cham (2016)[35] Liu, C., Boumal, N.: Simple algorithms for optimization on Riemannian manifolds with constraints. Appl. Math. Optim. 82, 949–981 (2020)[36] Nie, J.: The A-truncated K-moment problem. Found. Comput. Math. 14, 1243–1276 (2014)[37] Qu, Q., Sun, J., Wright, J.: Finding a sparse vector in a subspace: Linear sparsity using
alternating directions. Adv. Neural. Inf. Process. Syst. 27 (2014)[38] Qu, Q., Zhu, Z., Li, X., Tsakiris, M.C., Wright, J., Vidal, R.: Finding the sparsest vectors ina subspace: Theory, algorithms, and applications. arXiv preprint arXiv:2001.06970. (2020)

CP Factorization by a Riemannian Smoothing Method

25[39] Rockafellar, R.T., Wets, R.J.B.: Variational Analysis. Springer Science & Business Media,
Berlin (2009)[40] Rudin, W.: Principles of Mathematical Analysis. McGraw-Hill, New York (1964)[41] Sato, H., Iwai, T.: A new, globally convergent Riemannian conjugate gradient method. Optim.
64, 1011–1031 (2015)[42] Shilon O.: Randorthmat. MATLAB Central File Exchange. https://www.mathworks.com/
matlabcentral/fileexchange/11783-randorthmat (2022). Accessed 7 April 2022
[43] Dutour Sikirić, M., Schürmann, A., Vallentin, F.: A simplex algorithm for rational cpfactorization. Math. Program. 187, 25–45 (2020)[44] So, W., Xu, C.: A simple sufficient condition for complete positivity. Oper. Matrices. 9,
233–239 (2015)[45] Sponsel, J., Dür, M.: Factorization and cutting planes for completely positive matrices bycopositive projection. Math. Program. 143, 211–229 (2014)[46] Zhang, C., Chen, X.: A smoothing active set method for linearly constrained non-Lipschitz
nonconvex optimization, SIAM J. Optim. 30, 1–30 (2020)[47] Zhang, C., Chen, X., Ma, S.: A Riemannian smoothing steepest descent method for nonLipschitz optimization on submanifolds. arXiv preprint arXiv:2104.04199. (2021)