A polynomial-time algorithm for computing low CP-rank decompositions

Information Processing Letters 118 (2017) 10–14

Contents lists available at ScienceDirect

Information Processing Letters www.elsevier.com/locate/ipl

A polynomial-time algorithm for computing low CP-rank decompositions Khaled Elbassioni, Trung Thanh Nguyen ∗ Masdar Institute of Science and Technology, PO Box 54224, Abu Dhabi, United Arab Emirates

a r t i c l e

i n f o

Article history: Received 1 November 2015 Received in revised form 3 September 2016 Accepted 11 September 2016 Available online 16 September 2016 Communicated by A. Tarlecki Keywords: Completely positive matrix CP-rank CP-decomposition Polynomial-time algorithm Computational complexity

a b s t r a c t This paper investigates computing completely positive (cp) decompositions of positive semi-definite (PSD) matrices, a problem which arises in many applications. We propose the first polynomial-time algorithm for checking if a given PSD matrix has cp-rank of at most r, when r is a given constant. In addition, we present a polynomial-time algorithm for computing a (rational) cp-decomposition of such a matrix if it exists, within any desired accuracy. © 2016 Elsevier B.V. All rights reserved.

1. Introduction A (nonegative) real positive semi-definite (PSD) matrix Q of size n × n is said to be completely positive (cp) if and only if it can be decomposed as Q = U U T , for some nonnegative matrix U ∈ Rn+×r . We refer to r as the inner dimension of the product U U T . The smallest inner dimension r such that such a matrix U exists is called the cprank of Q , denoted cp-rank( Q ) = r, and the corresponding decomposition is called cp-decomposition (a.k.a. symmetric nonnegative factorization). If no such matrix U exists then cp-rank( Q ) = +∞ (we refer to the textbook [1] for a comprehensive overview of this subject). The problem of finding the cp-decomposition (and the cp-rank) of cp-matrices has received considerable attention in the last few decades due to its wide range of applications. One such application was found in studying block design arising in combinatorial analysis [2] where a connection between the existence of block design and cp-decomposition was shown. In

*

Corresponding author. E-mail addresses: [email protected] (K. Elbassioni), [email protected] (T.T. Nguyen). http://dx.doi.org/10.1016/j.ipl.2016.09.004 0020-0190/© 2016 Elsevier B.V. All rights reserved.

data analysis, the cp-decomposition is used in clustering as a useful tool for organizing observed data into relevant clusters (or groups) which maximize the similarity of elements within each cluster while ensuring the dissimilarity between these clusters is maximized [3–6]. Other applications appear in finance and economics [7,8], biology [9], and quadratic programming [10]. The complexity of deciding the existence of a cpdecomposition for a given matrix remained open for a long time. Very recently, Dickinson and Gijben [11] proved that this problem is NP-hard but left open the question of whether or not it belongs to NP. On the other hand, there exist efficient algorithms for computing cpdecompositions of matrices having special structures such as diagonally dominant matrices [12]; matrices of rank k which contain a diagonal principal submatrix of size k × k [13,14]; and matrices whose underlying graphs are bipartite [15], acyclic or circular [16]. Several attempts have been made to solve the problem in the general case: heuristic approach [17], multiplicative update rule [3–6], Newton-like algorithm and ANLS-based algorithm [18], and an algorithm using a so-called Procrustes projections [6]. These approaches, however, look for a good approximate

K. Elbassioni, T.T. Nguyen / Information Processing Letters 118 (2017) 10–14

solution rather than an exact one, and they all do not give any theoretical guarantee on the quality of the solution found. The only known exact algorithm was proposed by Berman and Rothblum in [19], which employs quantifierelimination for the first order theory of the reals [20]. Their basic idea is to transform the factorization problem into a problem of finding a solution of a system of polynomial inequalities which is known to be solved by an algorithm (e.g., [20,21]) whose running time is polynomial in the number of equations but is exponential in the number of variables. In fact, the number of variables used in the formulation in [19] is exactly nr, which makes the complexity of the algorithm exponential in both n and r. An interesting question is whether there exists an algorithm whose running time is singly-exponential in r. Of course, such an algorithm will not have reduced running time in case when the cp-rank of the cp-matrix is as large as its dimension.1 Nevertheless, in a number of practical applications such as document clustering, information retrieval, and image processing, cp-ranks of considered matrices are often much smaller than their dimensions. In this case, therefore, a singly-exponential in r algorithm will be very helpful to know. To the best of our knowledge, the answers to the question above are only known for very special cases when2 r = 1 and r = 2 (see [14]). It is worth noting that there is a notion related to the cp-rank, the nonnegative rank of a nonnegative (not necessarily symmetric) matrix, which has been studied intensively in the literature: a positive integer r is called the nonnegative rank of a matrix P ∈ Rn+×m , denoted rank+ ( P ) = r, if and only if it is the smallest number such that there exist two matrices U ∈ Rn+×r , V ∈ Rr+×m such that P = U V . Like cp-rank, computing nonnegative rank is also NP-hard (see [23]). Arora et al. [24] and Moitra [25] proposed exact polynomial-time algorithms for computing a nonnegative factorization, which have running time 2 r 2 (nm)2r 2 and (nm) O (r ) , respectively. Even though they build on the general idea in [26,19], this significantly improves the result of Cohen and Rothblum [26] which used totally mr + nr variables. Here, we would like to emphasize that, in general, the nonnegative rank and the cp-rank are different from each other and, in fact, it is well-known that cp-rank( Q ) ≥ rank+ ( Q ) ≥ rank( Q ) (see [19]) for every PSD matrix Q , where rank( Q ) is the (ordinary) rank of the matrix Q (one can easily construct an example for which the first inequality is strict by a non-constant margin). For a few very specific circumstances a nonnegative factorization will result in a completely positive one [27], but it would be not true in general. The main contribution of this paper is an n O (rs) -time deterministic algorithm that, given a PSD matrix Q of size n × n and of rank s, decides if a cp-decomposition of inner dimension at most r exists. In terms of the rank s, our al-

1 We refer the interested readers to the paper by Bomze, Dickinson, and Still [22] for more details on the generic properties of cp-rank and the distribution of cp-matrices according to their cp-ranks. 2 This is a trivial case because of the fact that if Q is a completely positive matrix of rank 1 then cp-rank( Q ) = rank( Q ) = 1 ([1]). Therefore, checking if Q is completely positive is equivalent to finding a nonnegative vector v such that Q = v v T which can be done in linear time.

11 3

gorithm’s running time is bounded by n O (s ) which could be also of independent interest as it contains only s in the exponent. Furthermore, we will show that a rational approximation to the solution to within an additive accuracy of δ > 0 can be computed in time poly(L, n O (rs) , log 1δ ), where L is the bit length of the input. 2. Completely positive decompositions A naive approach for checking if a given matrix Q = (qi j )n×n has cp-rank( Q ) ≤ r is to determine if there exists a nonnegative matrix U = (u ik )n×r such that Q = U U T . This amounts to checking the solvability of the following system over the reals:

 r

k=1 u ik u jk

u ik ≥ 0,

= q i j , i , j = 1, . . . , n , i = 1, . . . , n ; k = 1, . . . , r .

(1)

Generally, whether or not a system of N polynomial inequalities of degree at most d with T variables has a solution is decidable in ( Nd) O ( T ) time by using, e.g. the quantifier elimination algorithm of Renegar [20]. This, in fact, yields an n O (nr ) -time algorithm when applied to the system (1). In the following, we will show an exact n O (rs) -time algorithm for the decision problem of whether or not cp-rank( Q ) ≤ r, where s = rank( Q ). The idea is to use an equivalent formulation of the system (1) which requires a number of variables depending on rs only. This along with the Renegar’s algorithm yields an algorithm that is exponential in rs, but runs in polynomial time when r is fixed. Theorem 1. There is an n O (rs) -time deterministic algorithm that given a positive semi-definite matrix Q of size n × n and of rank s checks if a nonnegative factorization Q = U U T of inner dimension at most r exists. Proof. We proceed essentially along the same lines as in the proof of Theorem 2.1 in [24] (for finding a simplicial factorization), but using a different approach based on taking the singular value decomposition of the unknown matrix U to prove completeness for the system of polynomial equations. Let Q be a nonnegative positive semi-definite matrix of size n × n. Assume that Q has a decomposition Q = U U T , where U is a nonnegative matrix of size n × r. A basic fact from linear algebra is that Q , U and U T have the same rank: rank( Q ) = rank(U ) = rank(U T ) = s ≤ r (see, e.g. [28]). Fix an arbitrary basis V of the column vectors of Q , let B be an n × s matrix corresponding to this basis (more precisely, B is the matrix containing exactly s linearly independent columns of V ), and let Q V be the matrix of size s × n corresponding to the (unique) representation of Q in the basis V , that is, B Q V = Q . To prove the claim of the theorem, it will suffice to prove the following. Claim 1. Q has a nonnegative factorization Q = U U T of inner dimension at most r if and only if there is an r × s matrix H satisfying the following two conditions:

12

K. Elbassioni, T.T. Nguyen / Information Processing Letters 118 (2017) 10–14

(i) H Q V is a nonnegative matrix, and (ii) ( H Q V ) T ( H Q V ) = Q . Proof. (⇐) The sufficient condition is easy to show. Suppose that the conditions (i) and (ii) are satisfied. Let U = ( H Q V )T . This matrix is nonnegative and has size of n × r and thus, Q has cp-rank at most r. (⇒) Now suppose that there is a nonnegative factorization Q = U U T , where U is an n × r matrix, with r ≤ r. For simplicity, we can assume r = r as one can add r − r columns of zero entries to the matrix U in case r < r. We will prove that there exists a matrix H of size r × s satisfying condition (i) and (ii). It is well known that there exists a factorization of U , called a singular value decomposition of U , of the form U = L S R T (see, e.g. [29]), where L , R are unitary matrices3 of size n × n and r × r, respectively, and S is an n × r diagonal matrix with the first s diagonal entries are non-zero whilst other entries are zero (these first s non-zero diagonal entries of S are exactly the singular values of U ). Since all the entries of the matrix U are real numbers, L and R are orthogonal matrices, i.e., L −1 = L T and R −1 = R T , or equivalently, L T L = I n×n and R T R = I r ×r , where I k×k denotes the identity matrix of size k × k. We then can rewrite Q in the form

Q = ( L S R T )( L S R T ) T = ( L S R T )( R S T L T ).

(2)

R S + L T (an r

Let N = × n matrix), where S + is the pseudoinverse4 of matrix S. Note that such a matrix S + always exists. More precisely, S + is an r × n matrix such that:

S S+ = +





S S=

I s×s 0 0 0 I s×s 0 0 0



and n×n



.

(3)

r ×r

On the other words, S S + and S + S are the diagonal matrices in which only the first s diagonal entries are non-zero and equal to 1, whilst the other entries are zero. Define an r × s matrix H = N B. We now prove that this matrix H satisfies the two conditions (i) and (ii). Indeed, since N = R S + L T and B Q V = Q (see the construction of B and Q V above) we have:

H Q V = N B Q V = (R S+ LT )B Q V = (R S+ LT ) Q

(4)

From (2) and (4), it follows

H Q V = ( R S + L T )( L S R T )( R S T L T )

= R S+ S S T LT = R S T LT = U T . The second equality is due to the orthogonality of the matrices L and R, while the third equality is obtained by using (3). The claim follows. 2 Turning back to the proof the theorem. By the claim, to check if cp-rank( Q ) ≤ r is equivalent to determining if a 3

A matrix is unitary if its conjugate transpose is also its inverse. We refer to the textbook [29] for a formal definition of a pseudoinverse matrix. 4

system of at most nr linear inequalities and n2 quadratic equations on at most rs variables is feasible. This decision problem can be solved in n O (rs) time using quantifier elimination algorithms [21]. This completes the proof. 2 It is worth noting that cp-rank( Q ) ≤ max{1, 12 s(s + 1) − 1} holds for any cp-matrix Q of rank s (see [30]). Therefore, this and the result of Theorem 1 together immediately imply the following result: 3

Corollary 1. There is an n O (s ) -time deterministic algorithm to check if Q , a given positive semi-definite matrix of size n × n and of rank s, is completely positive. In case Q has a cp-rank of at most r, finding exactly matrix U (with finite precision) could be impossible as it might involve irrational entries. In the following we instead propose a polynomial-time algorithm for computing a (rational) cp-decomposition of Q within any desired accuracy. We denote by M ∞ := maxi , j | M i j | the max norm of any real matrix M = ( M i j )n×n . Theorem 2. Given a rational positive semi-definite matrix Q of size n × n and of rank s such that cp-rank( Q ) = r, and  > 0, one can compute a rational nonnegative n × r-matrix  U such that  U U T ≥ Q and Q −  U U T ∞ ≤  , in time poly(L, n O (rs) , log 1 ). Proof. Let s = rank( Q ) and B and Q V be the matrices corresponding to basis V of Q , as in the proof of Theorem 1. Then from this proof, it follows that we can discover in n O (rs) time that there is an r × s-real matrix H such that conditions (i) and (ii) of Claim 1 hold. Furthermore, in time poly(L, n O (rs) , log 1δ ), we can find an r × s-rational matrix  H , such that H −  H ∞ ≤ δ , for any desired accuracy δ > 0 (see [31,20]). Let B be a non-singular s × s-submatrix of B, obtained by selecting s linearly independent rows of B, and let Q be the corresponding s × n submatrix of Q . Then Q V = ( B )−1 Q . Let us choose

δ := min

 

3/2 1/2 4rs2 ( B )−1 ∞ Q ∞ (1+s2 ( B )−1 ∞ Q ∞ )

,

1 2



(5)

and set δ := sδ ( B )−1 ∞ and  H :=  H + δ E r ×s B , where E r ×s is the r × s-matrix with all-ones. Define

U := ( H Q V ) T ,  U := ( H Q V ) T , and  U := ( H Q V )T . Then, it follows that

 U T = U T + [δ E r ×s + ( H − H )( B )−1 ] Q ≥ U T + (δ − sδ ( B )−1 ∞ ) E r ×s Q = UT, where the first inequality follows from the fact that Q ≥ 0, and the second equality follows from our choice of δ . Let us further note that

K. Elbassioni, T.T. Nguyen / Information Processing Letters 118 (2017) 10–14

 U U T − U U T = ( Q )T T U T + U  Q T

T



+ (Q )  Q ,

(6)

+ ( H− To bound Q −  U U T ∞ , we bound the norm of each term in (6). We first note that where  := δ E r ×s



H )( B )−1 .

−1

 ∞ ≤ δ + δ s ( B )

13

very interesting to complement this with a hardness result. In particular, it would be interesting to see whether an no(r ) -algorithm for the problem would imply a subexponential time algorithm for 3-SAT, as in the nonnegative factorization case [24]. Acknowledgements



∞ = 2δ ,

The authors thank the anonymous referees for their valuable comments and suggestions that substantially improved the paper. This work was supported by the MI-MIT Flagship Project 13CAMA1.

and thus

U  Q ∞ ≤ 2δ U E r ×s Q ∞ ≤ 2δ rs U ∞ Q ∞ ≤ 2δ rs Q ∞ , 3/2

(7)

where the last inequality follows from the fact that 1/ 2

U ∞ ≤ Q ∞ . Now we write

T  = (δ )2 E s×r E r ×s + δ ( E s×r X + X T E r ×s ) + X T X,

(8)

where X := ( H − H )( B )−1 . Note that

X ∞ ≤ s  H − H ∞ ( B )−1 ∞ ≤ sδ ( B )−1 ∞ = δ .

(9)

From (8) and (9) follows the inequality   ∞ ≤ 4r (δ )2 , and then the inequality T

( Q )T T  Q ∞ ≤ 4r (δ )2 ( Q )T E s×s Q ∞ ≤ 4r (sδ )2 Q 2∞ .

(10)

Using (7) and (10) in (6), we get

 U U T − Q ∞ ≤ 4δ rs Q ∞ + 4r (sδ )2 Q 2∞ 3/2

≤ 4δrs2 ( B )−1 ∞ Q ∞ (1 + s2 ( B )−1 ∞ Q ∞ ) 3/2

1/2

≤ , by our selection (5) of δ . Finally one can show that log 1δ ≤ poly(L, log 1 ), where L is the size of the input and is

bounded by O (n2 log || Q ||∞ ). To prove the claim it is observed from (5) that log 1δ is polynomially bounded by

log || Q ||∞ , log ||( B )−1 ||∞ , log r, and log 1 . Since r ≤ O (n2 ), it remains to show that log ||( B )−1 ||∞ is bounded by O (L). This is because the fact that B is a (non-singular) submatrix of Q and its invertible matrix can be computed in polynomial time in L (see, e.g. [29]). The proof is completed. 2 3. Conclusion

In this paper, we gave an algorithm for checking whether the cp-rank of a given n × n positive semidefinite matrix is at most r, and provided a rational cpdecomposition within any accuracy, if exists. The running 2 time of the algorithm is bounded by n O (rs) (or n O (r ) since s ≤ r, where s is rank of the matrix) which is polynomial when r is a constant. For the future work, it would be

References [1] A. Berman, N. Shaked-Monderer, Completely Positive Matrices, World Scientific, 2003. [2] M. Hall, Combinatorial Theory, 2nd ed., John Wiley & Sons, Inc., New York, USA, 1986. [3] R. Zass, A. Shashua, A unifying approach to hard and probabilistic clustering, in: ICCV, 2005, pp. 294–301. [4] C. Ding, X. He, On the equivalence of nonnegative matrix factorization and spectral clustering, in: SDM, 2005, pp. 606–610. [5] Z. He, S. Xie, R. Zdunek, G. Zhou, A. Cichocki, Symmetric nonnegative matrix factorization: algorithms and applications to probabilistic clustering, IEEE Trans. Neural Netw. 22 (12) (2011) 2117–2131. [6] K. Huang, N.D. Sidiropoulos, A. Swami, Non-negative matrix factorization revisited: uniqueness and algorithm for symmetric decomposition, IEEE Trans. Signal Process. 62 (1) (2014) 211–224. [7] A. Vandendorpe, D. Ho, S. Vanduffel, P. Dooren, On the parameterization of the creditrisk+ model for estimating credit portfolio risk, Insur. Math. Econ. 42 (2) (2008) 736–745. [8] L.J. Gray, D.G. Wilson, Nonnegative factorization of positive semidefinite nonnegative matrices, Linear Algebra Appl. 31 (1980) 119–127. [9] C. Kelly, A test of the Markovian model of DNA evolution, Biometrics 50 (1994) 653–664. [10] S. Burer, On the copositive representation of binary and continuous nonconvex quadratic programs, Math. Program. 120 (2) (2009) 479–495. [11] P. Dickinson, L. Gijben, On the computational complexity of membership problems for the completely positive cone and its dual, Comput. Optim. Appl. 57 (2) (2014) 403–415. [12] M. Kaykobad, On nonnegative factorization of matrices, Linear Algebra Appl. 96 (1987) 27–33. [13] N. Shaked-Monderer, A note on upper bounds of the cp-rank, Linear Algebra Appl. 431 (2009) 2407–2413. [14] V. Kalofolias, E. Gallopoulos, Computing symmetric nonnegative rank factorizations, Linear Algebra Appl. 436 (2012) 421–435. [15] A. Berman, R. Grone, Bipartite completely positive matrices, Math. Proc. Camb. Philos. Soc. 103 (1988) 269–276. [16] P. Dickinson, M. Dür, Linear-time complete positivity detection and decomposition of sparse matrices, SIAM J. Matrix Anal. Appl. 33 (3) (2012) 701–720. [17] F. Jarre, K. Schmallowsky, On the computation of C∗ certificates, J. Glob. Optim. 45 (2) (2009) 281–296. [18] D. Kuang, S. Yun, H. Park, SymNMF: nonnegative low-rank approximation of a similarity matrix for graph clustering, J. Glob. Optim. 62 (3) (2015) 545–574. [19] A. Berman, U. Rothblum, A note on the computation of the cp-rank, Linear Algebra Appl. 419 (2006) 1–7. [20] J. Renegar, On the computational complexity and geometry of the first-order theory of the reals, J. Symb. Comput. 13 (3) (1992) 255–352. [21] S. Basu, R. Pollack, M. Roy, On the combinatorial and algebraic complexity of quantifier elimination, J. ACM 43 (6) (1996) 1002–1045. [22] I.M. Bomze, P.J. Dickinson, G. Still, The structure of completely positive matrices according to their cp-rank and cp-plus-rank, Linear Algebra Appl. 482 (2015) 191–206. [23] S. Vavasis, On the complexity of nonnegative matrix factorization, SIAM J. Optim. 20 (3) (2009) 1364–1377. [24] S. Arora, R. Ge, R. Kannan, A. Moitra, Computing a nonnegative matrix factorization – provably, in: STOC, 2012, pp. 145–162.

14

K. Elbassioni, T.T. Nguyen / Information Processing Letters 118 (2017) 10–14

[25] A. Moitra, An almost optimal algorithm for computing nonnegative rank, in: SODA, 2013, pp. 1454–1464. [26] J. Cohen, U. Rothblum, Nonnegative ranks, decompositions and factorizations of nonnegative matrices, Linear Algebra Appl. 190 (1993) 149–168. [27] M. Catral, L. Han, M. Neumann, R. Plemmons, On reduced rank for symmetric nonnegative matrices, Linear Algebra Appl. 393 (2004) 107–126.

[28] R. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985. [29] G. Strang, Introduction to Linear Algebra, fifth edition (5th ed.), Wellesley-Cambridge Press, MA, USA, 2016. [30] F. Barioli, A. Berman, The maximal cp-rank of rank k completely positive matrices, Linear Algebra Appl. 363 (2003) 17–33. [31] D. Grigoriev, N. Vorobjov, Solving systems of polynomial inequalities in subexponential time, J. Symb. Comput. 5 (1/2) (1988) 37–64.