A Non-commutative Extension of  Lee-Seung's Algorithm for Positive Semidefinite Factorizations
Authors: Yong Sheng Soh
Antonios Varvitsiotis
Venue: NeurIPS 2021 Poster (2021)
Keywords: PSD factorizations
majorization minimization
Lee-Seung's algorithm
Semidefinite optimization
Lieb's concavity theorem
PSD tensor factorization

Abstract: Given a data matrix $X\in \mathbb{R}_+^{m\times n}$ with non-negative entries, a Positive Semidefinite (PSD) factorization of $X$ is a collection of $r \times r$-dimensional PSD matrices $\{A_i\}$ and $\{B_j\}$ satisfying the condition $X_{ij}= \mathrm{tr}(A_i B_j)$ for all $\ i\in [m],\ j\in [n]$.  PSD factorizations are fundamentally linked to understanding the expressiveness of semidefinite programs as well as the power and limitations of quantum resources in information theory.  The PSD factorization task generalizes the Non-negative Matrix Factorization (NMF) problem in which we seek a collection of $r$-dimensional non-negative vectors $\{a_i\}$ and $\{b_j\}$ satisfying $X_{ij}= a_i^T b_j$,  for all $i\in [m],\ j\in [n]$ -- one can recover the latter problem by choosing matrices in the PSD factorization to be diagonal.  The most widely used algorithm for computing NMFs of a matrix is the Multiplicative Update algorithm developed by Lee and Seung, in which non-negativity of the updates is preserved by scaling with positive diagonal matrices.  In this paper, we describe a non-commutative extension of Lee-Seung's algorithm, which we call the Matrix Multiplicative Update (MMU) algorithm, for computing PSD factorizations.  The MMU algorithm ensures that updates remain PSD by congruence scaling with the matrix geometric mean of appropriate PSD matrices, and it retains the simplicity of implementation that the multiplicative update algorithm for NMF enjoys.  Building on the Majorization-Minimization framework, we show that under our update scheme the squared loss objective is non-increasing and fixed points correspond to critical points.  The analysis relies on a Lieb's Concavity Theorem.  Beyond PSD factorizations, we show that the MMU algorithm can be also used as a primitive to calculate block-diagonal PSD factorizations and tensor PSD factorizations.  We demonstrate the utility of our method with experiments on real and synthetic data. 

A Non-commutative Extension of Lee-Seung's Algorithm for Positive Semidefinite Factorizations
Yong Sheng Soh Department of Mathematics
National University of Singapore 10 Lower Kent Ridge Road, Singapore 119076 and Institute of High Performance Computing Fusionopolis Way, #16-16 Connexis, Singapore 138632 matsysěnus. edu.sg
Antonios Varvitsiotis
Engineering and Systems Design Pillar Singapore University of Technology and Design 8 Somapah Road, Singapore 487372 avarvits@gmail.cm
Abstract
Given a data matrix X RT with non- non-negatvve entries, a Positive Semidefinite (PSD) factorization of X is collection of r x r-dimensional PSD matrices {A,) and (B,) satisfying the condition X,j B,) for all PSD factorizations are fundamentally linked to understanding the expressiveness semidefinite programs as well as the power and limitations of quantum re- sources in information theory. The PSD factorization task generalizes the Non- negative Matrix Factorization (NMF) problem in which we seek collection r-dimensional non-negative vectors {a,) and {by} satisfying X,₃ a. by. for all E [m]. j E [n] one can recover the latter problem by choosing matrices in the PSD factorization to be diagonal. The most widely used algorithm for computing NMFs of a matrix is the Multiplicative Update algorithm developed by Lee and Seung, in which non-negativity of the updates is preserved by scaling with posi- tive diagonal matrices. In this paper, we describe a non-commutative extension Lee-Seung's algorithm, which we call the Matrix Multiplicative Update (MMU) algorithm, for computing PSD factorizations. The MMU algorithm ensures that updates remain PSD by congruence scaling with the matrix geometric mean of appropriate PSD matrices, and it retains the simplicity of implementation that the multiplicative update algorithm for NMF enjoys. Building on the Majorization- Minimization framework, we show that under our update scheme the squared loss objective is non-increasing and fixed points correspond to critical points. The analysis relies on Lieb's Concavity Theorem. Beyond PSD factorizations, we show that the MMU algorithm can be also used as a primitive to calculate block- diagonal PSD factorizations and tensor PSD factorizations. We demonstrate the utility of our method with experiments on real and synthetic data.
[m]] j E [n].
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Introduction
Let X E " be a X n dimensional matrix with non-negative entries and r E N a user-specified parameter. An r-dimensional positive semidefinite ( factorization of X is given by two families of X r PSD matrices Ai...., Am and B1. ..... Bn satisfying
,, tr(A, i [m], j E [n]]
Every non-negative matrix admits an r-dimensional PSD factorization for an appropriate value of r N-we may, for instance, take A, = (X.:) and B = diag(e;), choice that corresponds to an n-dimensional PSD factorization. The smallest r E N for which X admits an r-dimensional PSD factorization is called the PSD-runk [6].
PSD factorizations are of fundamental importance to a wide range of areas, most notably towards understanding the expressive power of linear optimization over the cone of positive semidefinite matrices, [9] 12]. studying the power and limitations of quantum resources within the framework of
information theory p 15], and as natural non-commutative generalization of the extremely popular dimensionality reduction technique of non-negative matrix factorizations (NMFs) [3] 30.. We
elaborate further on the relevance of PSD factorizations to each of these areas below.
Links to Semidefinite Programming. A Semidefinite Program (SDP) is a convex optimization problem in which we minimize linear function over the of PSD matrices intersected with an affine subspace. SDPs are a powerful generalization of Linear Programs with extensive modeling power and tractable algorithms for solving them, e.g., see (36 and references therein. SDPs are fre- quently used as convex relaxations to combinatorial problems, and have many important applications including optimal power flow computation (21). robustness certification to adversarial examples in neural networks 31]. and inference in graphical models (5)
Given bounded polytope P {or E Rdd cf d,, i (Q) conv("1.... ... ...) basic question concerning the expressive power of SDPs is to find the smallest possible SDP description P, i.e., the minimum N for which we can express P as the projection of an the affine slice of the cone o r x r PSD matrices. Concretely, the goal in this setting is to express P as:
P xS{ nL),
where s is the cone of r x r PSD matrices, L is an affine subspace of the space of r X r sym- metric matrices and π linear projection from the space of r X r symmetric matrices to R A representation of the form 2 is called an extended formulation or PSD-lift of P and extremely useful for optimization purposes. Indeed, the existence of PSD-lift immediately implies that min((c, (c, x) x E P} = min (N (c), y) : y E S. n L), and consequently, linear optimization over P (which can be hard) corresponds to an SDP (which can be solved efficiently).
To describe the connection of PSD factorizations with SDP-lifts, let Sp be the slack matrix of P. namely, Sp is a rectangular matrix whose rows are indexed by the facets of P. its columns indexed by extreme points uy and the ij-entry of Sp corresponds to the slack between the i-th facet and the j-th vertex, S.y d, was shown independently in [9] and [12] that if Sp admits an r-dimensional PSD factorization, the polytope / admits a PSD-lift over the cone of X r PSD matrices. The proof also constructive- given a PSD factorization of S, there is an explicit description of L and π that gives rise to P.
Vy. Generalizing seminal result by Yannakakis for LPs (38),
An important special case of the PSD factorization problem is when the PSD factors are block- diagonal PSD matrices, where both the number of blocks and the size of each block is fixed, i.e.. A,, By E ($ For fixed and user-specified r E N. the least k E N for which X E R Rmxx admits PSD factorization with PSD-factors in (3 )* is called the r-block diagonal PSD-rank of X. In terms of the geometric interpretation of PSD factorizations. black-diagonal PSD factorizations of the slack matrix Sp correspond to extended formulations of over Cartesian product of PSD cones, i.., P = (( ) n L). In terms of relevance to optimization, extended formulations over (S), ( allow to perform linear optimization over by solving block-diagonal SDPs, which can be solved numerically much faster compared to dense SDPs. In fact, most interior-point algorithms for SDPs are designed to exploit block-diagonal structure if it is present in the problem. The first systematic study of block-diagonal PSD-lifts was given in 7 where the focus was mainly on lower bounds.

---

Links to quantum information theory. Consider two parties, Alice and Bob, that try to generate samples (i,j) following some joint distribution P(i.j). For this section, it is crucial think of the distribution P(i, j) as being arranged in an entrywise non-negative matrix, and we use P to simultaneously refer to both the distribution and its matrix representation. Clearly, if P not a product distribution, Alice and Bob should either communicate or share some common information to be able to generate samples according to P. In the correlation generation problem, the goal is to find the least amount of shared resources that are needed to achieve this task. The considered resources can be either classical (shared randomness), quantum (shared entangled state) or hybrid.
In the quantum case, correlation generation boils down to finding a quantum state p E S with tr(p) = 1 and quantum measurements E,, F3 (i.c., E, Fy E S and >. E. F3 I) such that
(3) The least r E N for which a factorization of the form ( is possible is given by the (logarithm) the psd-rank of the matrix P [5]. Moreover, the proof of [1] is constructive, in the sense that, given rr-dimensional PSD factorization of P. there is explicit description of quantum state p E st and measurement operators acting on C" that satisfy 3 c.g., see 66 Proposition 3.8].
P(i,)) = tr(( & Fj)p), i E [m]] j E [n].
Moving beyond purely quantum protocols, there has been recent interest in hybrid classical-quantee protocols, motivated by the fact that near-term quantum devices can only operate reliably on lim- limited to
ited number of qubits (26). Specifically, assuming that their quantum capabilities manipulating 8 qubits, hybrid daassical-quantung protocols that allow to generate samples from a joint distribution P(i.j), correspond to PSD factorizations of P where the PSD factors are block- diagonal, with block-size at most 2 Moreover, the minimum amount of classical resources required in laasiiall-quatuee protocol with qubits, is given by the 2b-block diagonal PSD-rank of P.
Links to nonnegative matrix factorizations. An r-dimensional nonnegative matrix factorizaiioioooonnn (NMF) of X E 1TTx"
30 is specified by two families of r-dimensional entrywise nonnegative am E R4 and .,....,
[23
vectors a1.
......,,, bn E RT satisfying
(4) NMF is widely used dimensionality reduction tool that gives parts- -bsed representation of the input data, as it only allows for additive, not subtractive, combinations. To make this point clear, note that an equivalent reformulation of an r-dimensional NMF (4 is X = AB where RRxx is the matrix whose rows are the a,'s 's and the matrix B E R has as columns the by's, or equivalently,
X1 = (a), E [m]] j [m].
X.j E cone (A.....,, j [n]]
(5)
The equivalent viewpoint for NMFs given in (5) more amenable to interpretation, as it gives a representation of each column of X (i.e., each data point) as nonnegative (and thus additive) combi- nation of the r columns of A. and the columns of B give the coefficients of the conic combination. NMF factorizations have applications in many areas, notable examples including document cluster- ing (37, music analysis [8] speech- separation 32) and cancer-class identification [1O]. For comprehensive discussion on NMFs the reader is referred to the survey [1]] and references therein.
NMF factorizations are special case of PSD factorizations where the X r PSD matrices A, and Bj are diagonal, i.e., we have that A; = diag(a;) and B, diag(b;) for some vectors a;, by E RC (recall that a diagonal matrix is PSD iff its diagonal entries are nonnegative). Moreover, PSI) factorization of X, X1g tr( A, By), for which all the PSD) factors A,, Bj commute corresponds to an NMF factorization. In this sense, PSD factorizations are a non-commutative generalization of NMF factorizations.
Interpretabiliiy of PSD factorizations. An equivalent way to define an r-dimensional NMF for a data matrix X (cf. (4h) is through the existence of a liner mapping A R" - RRm satisfying
X,, E A(RT) for all E [r]] and A C R..
(6)
Consequently, the mapping A (or rather, the image of the extreme rays of the cone RC under A), describe latent space that can generate all data points X.j via nonnegative combinations.
Analogously, in the setting of PSD factorizations, the existence of an r-dimensional PSD factoriza- tion of X (cf. (I) is equivalent to the existence of a linear mapping A: S" + R satisfying
Xy A(S)) for all j [m] and A(S C R
3

---

Comparing (6) and (7). the difference between NMF and PSD factorizations is immediately appar- ent. In the setting of PSD factorizations the latent space is infiiiteddmennnn dimenssioil,,, and specifically, is the image of the extreme rays of the cone of r x PSD matrices (i.e., all matrices uu where E IR") under A. In this latent space, each data point x is represented by a PSD matrix B3 E E and using its spectral decomposition B, = Ei=1 A, leads to the representation X = Σ, A, A(u,u,)) Additional details and explicit examples demonstrating the qualitative difference in expressive power between NMF and PSD factorizations are given in Section
2 Prior Works on PSD Factorizations and Summary of Results
canonical starting point for finding an (approximate) r-dimensional PSD factorization of a given matrix X E I is to solve the non-convex optimization problem
inf EX - trA,
s.t. A A Am, B1, B, S4,
aiming to find an approximate r-dimensional PSD factorization that minimizes the square loss over
all entries of X. Fixing of the two families of matrix variables, say the A, ss problem (8) is
separable with respect to Bi..... Bm- Consequently, reasonable solution approach for (8) is to alternate between updating the A,'s 's's and B, 's by solving the sub-problems:
A; arg inf
(xg,BB((, s.t. AI...., Am E ESt
(9)
3
By f arg inf
Exxy-
s.t. B Bn E S
(10)
The two sub-problems in each update step are symmetric in the variables A, and Bj. with the small modification where replace X with its transpose. As such, for the remainder of this discussion, we only focus on the sub-problem (TO) corresponding to fixing the A,'s 's and updating the B, ss Moreover, (TOp is separable with respect to each variable B,, so it suffices focus on
inf Σ(x, tr(,,) s.t. ET.
(11)
Lastly, to simplify notation we omit subscripts. and specifically, we denote by the j-th column of X and by the PSD matrix variable By. Defining A S - RT to be the linear map A(2) = ((A1, Z). ... (Am, Z)), problem 11] can be then equivalently written
inf | A(B)II s.t. BES,
(12)
The optimization problem 112 is convex, and in fact, falls within the well-studied class of convex quadratic SDPs. Nevertheless, there no closed-form solution for this family of optimization prob- lems, and consequently, typical solution strategies rely on numerical optimization, c.g., see [33].
Summary of results. In this paper we introduce and study an iterative algorithm (Algorithm call the Matrix Multiplicative Update (MMU) algorithm for computing PSD factorizations. The MMU algorithm builds on the Antioriattio-MMninizaan vious section, the main workhorse is iterative algorithm for the convex quadratic SDP 12)
framework, and as discussed in the pre-
From a computational perspective, the iterates of the MMU algorithm are updated via conjugation with appropriately defined matrices, and morcover, the PSDness of the iterates is automatically guaranteed. From a theoretical per- spective, the squared loss objective is non-increasing along the algorithms'' trajectories (Theorem and moreover, its fixed points satisfy the first-order optimality conditions (Theorem 2}. The anal- ysis of the MMU algorithm relies the use of several operator trace inequalities (including Von Neumann's trace inequality and Lieb's Concavity Theorem).
method has the advantage of being simple to implement
An important feature of the MMU algorithm is that if it is initialized with block-diagonal PSD ma- trices, the same block-diagonal structure is preserved throughoul its execution, which leads to an algorithm for calculating block-diagonal PSD factorizations. In particular, if the MMU algorithm

---

is initialized with diagonal PSD matrices, the iterates remain diagonal PSD throughout, and as it turns out, our algorithm this case reduces Lee-Seung's seminal Multiplicative Update algo- rithm for computing NMFs 224] Moreover, we show how the MMU algorithm can be used as primitive to calculate PSD factorizations of nonnegative tensors. In terms of numerical experiments, we demonstrate the utility of our method for both synthetic and real data (CBCL image dataset).
Existing work. All existing algorithms for computing PSD factorizations employ the alternating minimization approach described in the previous section, where we fix one set of variables and min- imize over the other, and essentially boil down into finding algorithms for the convex problem (12)
Projected Gradient Method (PGM). The first approach for computing PSD factorizations is based on applying PGM to (12), alternating between gradient step to minimize the objective and a projection step onto the set of PSD matrices [35]. The latter projection step uses the following useful fact: of
Given the spectral decomposition C U',iag(A,)U of a matrix C E S", the projection onto the PSD cone Udiagg (max((, (0, UU [33. The vanilla PGM has slow convergence rate, so the authors in [35 also propose accelerated variant that incorporates momentum term.
Coordinate Descent. The authors in [35] also propose different algorithm combining the ideas of coordinate descent and change of variables that allows them to also control the rank of the PSD factors, which was popularized by the seminal work of Burer and Monteiro for solving rank- constrained SDPs B. Concretely, the authors use the parameterization A, a,a," and B, = bybj,
where a; E T and by E R for some fixed rA,, "B, E N, and optimize using coordinate descent scheme over the entries of the matrices a; and bj. In this setting, problem (12) is quartic polynomial in the entries of Thus, its gradient is cubic polynomial, and its roots can be found using Cardano's 's method and careful book-keeping (for similar approach see also (27).
Connections Affine Rank Minimization and Phase Retrieval. A different set of algorithms devel- oped in [17 18 19] is based on the connections between computing PSD factorizations with the affine rank minimization (ARM) and the phase retrieval (PR) in signal procesing. First, recall that the PSD-ARM problem focuses on recovering low-rank matrix from affine measurements:
min rank(B) s.t. A(B) = B BE
Here, A is a known linear map representing measurements while x is known vector of observations. Due to the non-convexity of the rank function, a useful heuristic initially popularized in the control community is to replace the rank by the trace function, e.g., see (28) and (29). in which case the resulting problem is an instance of an SDP. A different heuristic for PSD-ARM is to find a PSD matrix of rank at most & that minimizes the squared loss function, i.e.,
inf | A(B) s.t. B E rank(B) Sk
(13)
where alternatively, the rank constraint can be enforced by parametrizing the PSD matrix variable B E S as B = 66 with e R'** The point of departure for the works [[7 [8] 19 is that
problem (13) corresponds exactly to the sub-problem (2) encountered in any alternate minimization strategy for computing PSD) factorizations, albeit with an additional rank constraint. In view of this, any algorithm from the signal processing literature developed for ARM can be applied to (12).
The main algorithms considered in [7 [8] 19 are Singular Value Projection (SVP) [4,, Procrustes Flow (34), and variants thereof. In terms of convergence guarantees, for affine maps A obeying the Restricted Isometry Property 4 both algorithms converge to an optimal solution. Nevertheless, is unclear whether these guarantees carry over when applied to the PSD factorization problem.
Roadmap. In Section 3 we derive our MMU algorithm for computing PSD factorizations and in Section 4 we show that its fixed points correspond to KKT points. In Section 5 we give various theoretical applications of the MMU algorithm and in Section 6 we go from theory to practise and apply the MMU algorithm synthetic and real datasets.
3 A Matrix Multiplicative Update Algorithm for PSD Factorizations
In this section we describe our algorithm for computing (approximate) PSD factorizations of a ma- trix X. As we discussed, our method is an alternating minimization approach in which we alternate
5

---

between optimizing over the variables {A,} and {B,}. The sub- sub-problem in each update step is sym- metric in the variables {A,] and [B,}. with the small modification whereby we replace X with its transpose. As such, in the remainder of this discussion, we assume that the variables {A:} are fixed
and we perform the update on the variables (B,). The resulting sub-problem is given by (12)
(MM) Framework. Our algorithm is an instance of the (MM) framework, e.g. see [20] and references therein. To briefly describe this approach, suppose we need to solve the optimization problem min(F(x) : E X).. The MM framework relies on the existence of parametrized family of auxilliary functions Uy + R, one for each E X, where:
iejiorizttinn-Mnnnizzattii
Fly)= S u,(y), for all y E X' and Fxz) (r)
(14)
Based on these two properties, F is nonincresasing under the update rule:
I new argmin oz (y) y E }
(15)
as can be easily seen by: F(x"ow S Ugaad (x S (x)ld) = F (xood).
We conclude with two important remarks concerning the MM framework. First, note that although the iterates generated by the MM update rule (23) are nonincreasing objective function value, there is in general no guarantee that they converge to minimizer. Secondly, for the MM approach be of any use, the auxilliary functions employed at each iteration need to be easy to optimize.
Matrix Geometric Mean. Our choice of auxilliary functions relies on the well-studied notion of
geometric mean between pair positive definite matrices, whose definition we recall next. For additional details and omitted proof's the reader is referred to [2 222. The matrix geometric mean of two positive definite matrices and D given by
C#D= (D
(16)
or equivalently, it is the unique positive definite solution of the Riccati equation
xc-1x=10
(17)
in the matrix variable X. The matrix geometric mean also has a nice geometric interpretation in terms of the Riemannian geometry of the manifold of positive definite matrices, and specifically, C## D is the midpoint of the unique geodesic joining C and D. Finally, the matrix geometric mean symmetric in its arguments C# D#C and also satisfies (C#D)-11 -)= C-1#D-1
The MMU Algorithm for PSD Factorizations. The main step for deriving our algorithm for approximately computing PSD factorizations is apply the MM framework, with meticulously chosen auxilliary function, to the convex quadratic SDP (12). Our main result is the following:
Theorem 1. Consider a fixed vector x Rmm and let S" + Rm be the linear map defined by Z + A(Z) = (tr(A, ... tr(AAm Z)) for some fixed r x r positive definite matrices A Am. Then, the objective function || A(B) is non-increasing under the update rule
Buew W(AT))W, where w ([AT -1
and moreover, if initialized with a positive definite matrix, the iterates remain positive definite.
Proof. First, note that if the A,'s and Bota are all positive definite, the update rule well-defined. Indeed, we have [A" Al(Baia) = Σ Boia) A& is also positive definite, and thus invertible.
Set F(B) = || A(B) and define the function
B (B) F F(Bata) + (V F(aaa), B - Bota) + (B - Bota, T(B Bata))
(18)
where T + S" the operator given by
T(2) W zw-1 and
W - W Al(aa))) #(Boaa).
The claim of the theorem will follow as an immediate consequence the MM framework, as long as establish that un, (B) is an auxilliary function, i.e., satisfies the two properties given in (14).
Clearly, we have that u (Bola) F(Boa), so it only remains to show the domination property, that is, Baa (B) FBB), for all B St. In fact, show slightly stronger result, namely that

---

BBaaa (B) I F(B) holds for all symmetric matrices B E S' To see this we use the second order Taylor expansion of F Bola, which as quadratic B, is given by
F(B) = ((Boa)) + (VF(Bol)), B Buia) + ||A( Bota)
(19)
Comparing the expressions (18) and (99, to show that F(B) S U.B.ia (B) for all B E S" it suffices to check that the operator - AA' A is positive; i.e., (Z, [T A (Z)) Z 0 for any matrix Z E S".
This claim is the main technical part of the proof, deferred to Lemma ?? in the Appendix.
Furthermore, the fact that T - A A is a positive operator, also implies that T is itself a positive operator. Consequently, the MM update (23) obtained by using the auxilliary function (18) can be
calculated just by setting the gradient equal to zero, and is given by
Buew Bold T (LA) A] A Ax))
(20)
Moreover, as (Z) WWW and W (aaaa) (Boia) it follows that
Bold = W([A" () W = ([ATTA
(21)
where for the first equality we use the unicity property of the matrix geometric mean (recall 17 Subsequently, using (2T). the MM update rule in (20 simplifies to the following:
Bnew = T (Axx x) ((A"x( x)W.
(22)
Lastly, since the A,'s are PSD, it follows that A x E, xgA, a conic combination of PSD matrices (recall that x E RT), ). and thus, it is itself PSD. Consequently, Bnew is PSD. In fact, if the matrices A, and Bold are positive definite, the updated matrix Bnew is also positive definite.
Having established iterative method for problem (12 that is non-increasing in value and retains PSDness, we can incorporate this as a sub-routine in our alternating optimization scheme for com- puting PSD factorizations. The pscudocode of the resulting method is given in Algorithm I
Algorithm Matrix Multiplicatiie Update algorithm for computing PSD factorizations
Input: A matrix X E k parameter E N Output: (Ai.... Am ), (Bi...., B,} ST XiJ so t,), for for i,
while stopping criterion not satisfied:
A Vi(BT V(BT::))V where
Vn = B](A))
4;= Xx: W3 - (LAT A](B,)) xj ===x,
(23)
By (V((AT;;))) where
4 Fixed Points of the MMU Algorithm
In this section, we show that the fixed points of the MMU algorithm satisfy the Karush- Kuhn-Tucker (KKT) optimality conditions for problem (8) Letting (A])emm) (M, hejm] and {Bj [Aj be pairs of primal-dual optimal solutions of (8) with zero duality gap, it is straightforward to verify that the KKT conditions are
tr(A, = tr (BJAA) = O,
i [m]. E [n]
8 (X..) [B BJ(7)) = M
i [m]
A (X..x)) [AT AT(B) = AA.. j
Furthermore, assuming that the primal optimal solutions [stmm)) and (Ej)) are all positive definite, it follows immediately from the complementary slackness conditions that M = A3 for all E [m],j E [n]. Consequently, in the special case of positive definite optimal solutions {AA Ji][m)] and { Je(m] the KKT conditions reduce to
BT (X++) B B]{A),, E [m] and AT(X))= A(B)), j [m).
(24)
Based on the preceding discussion, in the next result (whose proof follows by Lemma ?? in the Appendix) shows that we can interpret our MMU algorithm as a fixed-point method for satisfying the KKT optimality conditions corresponding problem (8).
Theorem 2. If (Am) and (By)se sejn are positive definite fixed points of the update rule of the MWU algorithm given in 23 then they also satisfy the KKT conditions 24
7

---

5 Applications of the MMU algorithm
Block-diagonal (BD) PSD factorizations. the MMU algorithm is initialized with BD posi- tive definite matrices with the same block structure, the BD structure is preserved at each up- date. Indeed, as [AT A[(oaa))
_k=1 tr(As As we see that [A A[(oia) -1 and thus, share the same block structure. Lastly, by definition of the MMU algo- rithm (23). Bnew is also block-diagonal with the same structure. Thus, if initialized with BD-PSD matrices, the MMU algorithm gives a method for computing a BD-PSD factorization.
([A (
Recovering Lee-Seung's algorithm for NME. Diagonal matrices be considered block- diagonal in a trivial manner. Nevertheless, by the preceding discussion, if initialized with diagonal PSD matrices, the iterates of the MMU algorithm remain diagonal PSD throughout, In this special case, our MMU algorithm reduces to Lee-Seung's (I.S) seminal Multiplicative Update algorithm for computing NMFs [2]]. I.S's algorithm is perhaps the most widely used method for computing NMFs as it has succeeded to identify meaningful features in diverse collection of real-life data sets and is extremely simple to implement. Specifically, LS's updates are
XBT
A and B Bo AT AB"
A Ao
(25)
ABBT
where X ool Y, X/Y denote the componentwise multiplication, division of two matrices respec- tively.
Setting A; diag(a;) and Bj diag(b;). the MMU algorithm updates By as
Σ -==l (az, (az,b;) A1) * Em Xij which is also diagonal PSD) matrix. Setting m and B (b1 ... b,), this coincides with LS's update rule 25
Bj f- Bj
A
[a]
PSD factorizations for nonnegative tensors. Motivated by PSD factorizations of nonnegative matrices, [16] define r-dimensional PSD factorization of nonnegative tensor T (with indices of dimension d) as a collection of PSD matrices C, ",,""",,"""""""" E s4, for all is E [d] such that
c!.")
....).".)) Ty = summC (C()) ....o C(")), for all ig E [d], where denotes the Schur product of matrices and sum(X) = Eij Xij. The motivation for studying tensor PSD factorizations comes from the fact that they characterize the quantum correlation complexity for generating multipartite classical distributions [16]. We now show how our MMU algorithm can be used as a primitive to calculate tensor PSD factorizations. For simplicity of presentation we restrict to n = and consider
inf Σ (Thiisss sum( (1) cl clo 11.ta
)
2
subject to c C
(2)
Sr.
(26)
As in the case of PSD factorizations we employ a block coordinate descent approach. Specifi- cally, fixing all matrices C!( cf22 the optimization problem (26) is separable wrt each c(8) for all
is E [d]. Thus, defining the map A : ST + Rdd X + ((X, c (22) Der we need to solve
arg inf (T 1 subject to C CES CeS.
(27)
for all iz E [4]. Note that AA o - Sr where x = (Mz, H Σ En.ze The i1 is
((1) c(() is a
PSD matrix, as the Schur product of PSD matrices is PSD. Thus, Theorem gives an update rule that preserves PSDness and for which the objective function 27 is nonincreasing.
Numerical experiments
Damping. The implementation of the MMU algorithm requires to compute inverses and square- roots of certain positive definite matrices. These operations become ill-conditioned whenever the in- put matrices contain cigenvalues that are close to zero. To mitigate such issues, apply damping operation. First, when computing ([B ()))) -1 in (23) we instead compute ((3 (A,)) Second, when computing the matrix square root / in the process of computing the matrix geo- metric mean 23) using the expression (T6), we instead compute (X + eT)1/2, In our implementa- tions, we apply choice of = 10-8
8

---

Distance matrices. Let E R" be vector and let M be a x n matrix whose entries are M, = (c. - -y) M is known as a distance matrix and admits the following 2-dimensional PSD factorization
where A, = (1) B, -,)))))))))))))
Mi = tr(,,,,),),
We generate a random E IR" with 7 = 20 where each entry is drawn from the standard normal distribution. We apply our algorithm compute a 2-dimensional factorization whereby we perform 500 iterations over 50 random initializations. We compute the normalized squared error loss of the factorization from the data matrix, and we plot the error over each iteration Err (tr(tr(A, M,y²// M2 in Figure Our experiments suggest that, with sufficient random initializations, our algorithm finds a PSD factorization that is close to being exact.
04
-33
61-
00 :02
330
4000
500
400
500
Marataa
Figure 1: Performance of the MWU algorithm for computing a PSD factorization of a distance matrix. Different curves correspond different random initializations.
CBCL Face Image Dataset. In our second experiment, we apply the MMU method to compute PSD factorization of matrix comprising face images from the CBCL Face Database [.] The objective of this experiment is to illustrate how computing a PSD factorization can be viewed as a representation learning algorithm that generalizes NMF. The CBCL dataset comprises 2429 images of faces, each of size 19 x 19 pixels. We process the images so that the pixel intensity has mean 0.5 and standard deviation 0.25, with values subsequently clipped at [0, 1]. The resulting data matrix has size 361 x 2429.
We note that r-dimensional PSD factorization of X specified by {A.), [B;} gives rise to decomposition X = E.= As A(u,u,) ), where B₁ = E Aruuu, is the spectral decomposition. Subsequently, one can view the collection {A(u,u, ||u. = 1} as basic constituents from which all face images are expressible via non-negative linear combinations. This is analogous to NMF in which we express every data point non-negative linear combinations from collection of non- negative basis vectors these are precisely the linear image of a matrix factor B on standard basis vectors. In this experiment, given specific face image X.j. show its decomposition sums of constituents identified by the MMU algorithm these are the images A(u,u,) ), where u, are the eigenvectors of Bj-
As our baseline, we compute 27- -i-dimensional NMF of the CBCL data matrix over 500 iterations. In Figure we illustrate the decomposition of one images from the dataset into these 27 constituents. Next, we calculate 7-dimensional PSD factorization of the CBCL data matrix over 500 iterations and illustrate the decomposition of the same image the new basis in Figure 3 We compute a 7-dimensional PSD factorization because symmetric matrix of dimension 7 x 7 has 28 degrees of freedom, which is comparable with our NMF example. We note that the constituents learned from the PSD factorization appear to capture global features, a phenomenon that has been been also observed for NMF applied to datasets beyond CBCL (25)
Last, we apply the MMU algorithm to compute block-diagonal PSD factorization with 9 blocks of size 2 2 over 500 iterations (the number of degrees of freedom 27), and we illustrate the decom- position in Figure In this instance, constituents contain more localized features. An advantage of learning continuum of basic building blocks that one can express certain geometries in the data that is otherwise not possible using a finite number of building blocks. As an illustration of this intuition, in Figure 5 we show a continuum of atoms corresponding to a single 2 x 2 block which
captures transition between the nose and the nostrils.

---

Computational specifications. Our experiments were conducted in Python on an Intel 7-th Gen 17 processor 2.8GHz.
Code repository. Our codes are available on the following online repository
http::///ghuuu..ccm/ssaa//////
-:
Original
Reconstruction
Figure 2: Image decomposition into building blocks learned from 27-dimensional NME
Original
Reconstruction
Figure 3: Image decomposition into building blocks learned from 7-dimensional PSD factorization.
Original
Reconstruction
Figure 4:
Image decomposition into building blocks learned using 2 x 2-block PSD factorization.
Figure 5: Visualization of continuum of building blocks learned using 2 2-block PSD factorization.
Funding transparency statement. Yong Sheng Soh gratefully acknowledges Ministry of Education (Singapore) Academic Research Fund (Tier 1) R- R-146-000-329-133- Antonios Varvitsiotis gratefully acknowledges Ministry of Education (Singapore) Start-Up Research Grant SRG ESD 2020 154 and NRF2019- NRR--ANR099 ALIAS grant. Antonios Varvitsiotis acknowledges initial discussions with Cédric Févotte, Sandra Tan, and Vincent Tan.
References
[1] MIT-CBCL Face Database, Center for Biological and Computational Leamning Massachus- setts Institute of Technology. html
mit.
[2] Rajendra Bhatia. Positive Definite Matrices. Princeton University Press, 2007.
[3] Samuel Burer and Renato D. C. Monteiro. A Nonlinear Programming Algorithm for Solving Semidefinite Programs via Low-rank Factorization. Mathematical Programming, 95(2):329- 357, 2003.
[4] Emmanuel J. Candes and Terence Tao. Decoding by Linear Programming. IEEE Transactions on Information Theory, 11:4203-4215 2004.
[5] Murat A. Erdogdu, Yash Deshpande, and Andrea Montanari. Inference in graphical models via semidefinite programming hierarchies. In Advances Neural Information Processing Systems, 2017.
10

---

[6] Hamza Fawzi, João Gouveia, Pablo A. Parrilo, Richard Z. Robinson, and Rekha R. Thomas. Positive Semidefinite Rank. Mathematical Programming, 153(1):133-177, 2015.
[7] Hamza Fawzi and Pablo A. Parrilo. Exponential lower bounds on fixed-size psd rank and semidefinite extension complexity.
2013,
[8] Cédric Févotte, Naney Bertin, and Jean-Louis Durrieu. Nonnegative Matrix Factorization with the Itakura-Saii Saito Divergence: With Application to Music Analysis. Neural Computation, 21(3):793-830, 2009.
[9] Samuel Fiorini, Serge Massar, Sebastian Pokutta, Hans Raj Tiwary, and Ronald de Wolf. Ex- ponential Lower Bounds for Polytopes in Combinatorial Optimization. Journal of the ACM, (17), 2015.
[10] Yuan Gao and George Church. Improving Molecular Cancer Class Discovery Through Sparse Non-negative Matrix Factorization. Bioinformatics, 21(11)3370-33755 2005.
[11] Nicolas Gillis. The Why and How of Nonnegative Matrix Factorization, chapter Regularization, Optimization, Kernels, and Support Vector Machines, pages 257-291. Number Chapman & Hall/CRC in Machine Learning and Pattern Recognition Series. 2014.
[12] João Gouveia, Pablo A. Parrilo, and Rekha Thomas. Lifts of Convex Sets and Cone Factoriza- tions. Mathematics of Operations Research, 38(2):248-264, 2013.
[13] Nicholas J. Higham. Computing a Nearest Symmetric Positive Semidefinite Matrix. Linear Algebra and its Applications, 103:03-118, 1988.
[14] Prateek Jain, Raghu Meka, and Inderjit Dhillon. Guaranteed Rank Minimization via Singular Value Projection. In Advances in Neural Information Processing Systems, 2010.
[15] Rahul Jain, Yaoyun Shi, Zhaohui Wei, and Shengyu Zhang. Efficient Protocols for Generat- ing Bipartite Classical Distributions and Quantum States. IEEE Transctions on Information Theory, 59:5171 171-5788, 2013.
[16] Rahul Jain, Zhaohui Wei, Penghui Yao, and Shengyu Zhang. Multipartite quantum correlation and communication complexities. computational complexity volume, 26: 199-228, 2017.
[17] Dana Lahat and Cédric Févotte. Positive Semidefinite Matrix Factorization: A Link to Phase Retrieval and Block Gradient Algorithm. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2020.
[18] Dana Lahat and Cédric Févotte. Positive Semidefinite Matrix Factorization Based Trun- cated Wirtinger Flow. In 28th European Signal Processing Conference (EUSIPCO) 2020.
[19] Dana Lahat, Yanbin Lang, Vincent Y. F. Tan, and Cédric Févotte. Positive Semidefinite Ma- trix Factorization: A Connection with Phase Retrieval and Affine Rank Minimization. IEEE Transactions on Signal Processing, in press, 2021.
[20] Kenneth Lange. MM Optimization Algorithms. SIAM, 2016.
[21] Javad Lavaei and Steven H. Low. Zero Duality Gap in Optimal Power Flow Problem. IEEE Transactions on Power Systems, 27(1), 2012.
[22] Jimmie D. Lawson and Yongdo Lim. The Geometric Mean, Metrics, and More. The Amer- ican Mathematical Monthly, volume 108, pages 797-812. 2001.
[23] Daniel D. Lee and H. Sebastian Seung. Learning the Parts of Objects by Non-negative Matrix Factorization, Nature, 401, 1999.
[24] Daniel D. Lee and H. Sebastian Seung. Algorithms for Non Neg-teget Matrix Factorization. In Advances in Neural Information Processing Systems 13, 2000.
[25] Stan Z. Li, Xin Wen Hou, Hong Jiang Zhang, and Qian Sheng Cheng. Learning Spatially Localized, basts-based Representation. Proceedings of the 2001 IEEE Computer Society Conference Computer Vision and Pattern Recognition., 2001.
11

---

[26] Xiaodie Lin, Zhaohui Wei, and Penghui Yao. Quantum and Classical Hybrid Generations for Classical Correlations. https://arxi..rr/aa//200.....
2020.
[27] Jakub Marecek and Martin Takac. A low-rank coordinate-descent algorithm for semidefi- nite programming relaxations of optimal power flow. Optimisation Methods and Software, 52(4):849-87)) 2017.
[28] Mehran Mesbahi and George P. Papavassilopoulos. On the Rank Minimization Problem Over Positive Semidefinite Linear Matrix Inequality. IEEE Transactions Automatic Control, 42(2), 1997.
[29] Karthik Mohan and Maryam Fazel. New Restricted Isometry Results for Noisy Low-rank Recovery. In IEEE International Symposium on Information Theory (ISIT), 2010.
[30] Pentti Paatero and Unto Tapper. Positive matrix factorization: A non-negative factor model with optimal utilization of error. Envirommentrics, 5:111-126 1994.
[31] Aditi Raghumathan, Jacob Steinhardt, and Percy Liang. Semidefinite relaxations for certifying robustness to adversarial examples. In Advances in Neural Information Processing Systems 31, 2018.
[32] Mikkel N. Schmidt and Rasmus K. Olsson. Single-channel Speech Separation Using Sparse Non-negative Matrix Factorization. In Ninth International Conference on Spoken Language Processing, 2006.
[33] Kim-Chuan Toh. An inexact primal dual path following algorithm for convex quadratic sdp. Mathematical Programming. 112:221-254 2008
[34] Stephen Tu, Ross Boczar, Max Simchowitz, Mahdi Soltanolkotabi, and Benjamin Recht. Low- rank Solutions of Linear Matrix Equations via Procrustes Flow. In Proceedings of The 33rd International Conference Machine Learning, pages 964-973, 2016.
[35] Arnaud Vandaele, François Glineur, and Nicolas Gillis. Algorithms for Positive Semidefinite Factorization. Computational Optimization and Applications, 71(1):193-219 2018.
[36] Lieven Vandenberghe and Stephen P. Boyd. 48(1):49-95, 1996.
Semidefinite Programming. SIAM Review,
[37] Wei Xu, Xin Liu, and Yihong Gong. Document Clustering Based on Non-negative Matrix Factorization. In the 26th Annual International ACM SIGIR Conference on Research and Development Informaion Retrieval, pages 267-273, 2003.
[38] Mihalis Yannakakis. Expressing Combinatorial Optimization Problems by Linear Programs. Journal of Computer and System Sciences, 43:441-466., 1991.
12


Reviews:
[{"review_id": "o8Xbkck_M5Y", "reviewer": "Reviewer_QbiV", "summary": "The submission \"A Non-commutative Extension of Lee-Seung's Algorithm for Positive Semidefinite Factorizations\" develops a new algorithm for the problem of computing a PSD factorization of a nonnegative matrix. This generalizes the nonnegative matrix factorization (NMF) problem and is of interest both as a succinct representation of data, but also due to connections to lifted representations of convex sets and problems in quantum information. The algorithm is a non-obvious generalization of a well-known method for NMF due to Lee and Seung to the matrix setting. The authors establish that the updates are non-increasing for the squared loss objective and the fixed points satisfy the KKT conditions, if the iterates are initialized to be positive definite, and show that the algorithm can be used in practice via numerical experiments.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality: This paper introduces an algorithm for PSD factorization that is new, and generalizes a widely-used algorithm for nonnegative matrix factorization. The work is well situated in the existing literature on algorithms for PSD factorization, and what is different about this method (when compared to, for instance, directly solving the PSD constrained least squares subproblems via IPM or Burer-Monetiro-style factorization methods) is clearly articulated. \n\nQuality: The submission appears to be correct. It motivates the algorithm (via the MM framework and its connection to the Lee-Seung algorithm), describes it clearly, and establishes key results about the algorithm. Some of the proofs are deferred to the supplementary material. It would perhaps have been nice to include the 3 key steps involved in proving Lemma 3 in the main body (normalizing so B=I, the use of Lieb's theorem, and the use of the inequality tr(AB)^2<= tr(A^2)tr(B^2)), to make the body even more self-contained. \n\nIn the experimental section, it would have been nice to discuss two issues a little more:\n-- the addition of small positive quantities for \"conditioning purposes\" is a little mysterious, and may deserve more discussion. Is this required in the LS (commutative case)? Can this damped update also be interpreted naturally in the MM framework and shown to converge, etc?\n-- what is the stopping criterion used and/or suggested, or are the experiments just run for a fixed number of iterations? Are there any interesting issues related to choosing a stopping criterion?\n\nClarity: The paper is clearly written and easy to follow. \n\nI found a few minor typos:\n-- Line 72: \"if it present\" should be \"if it is present\"\n-- Line 111: consider \"Moreover, a PSD factorization...\" rather than \"Moreover, given a PSD factorization...\"\n-- Line 161: \"how to the MMU\" should be \"how the MMU\"\nI also felt that the proof of Lemma 3 could be streamlined a little. The concavity of \n$$f(X)= \\textrm{tr}(ZX^{1/2}ZX^{1/2})$$\n(a very special case of Lieb's concavity theorem) and the fact that $f$ is positively homogeneous of degree one in X, implies that this function satisfies, $f(X_1+X_2 + \\cdots + X_k) >= f(X_1)+f(X_2) + \\cdots + f(X_k)$ which is really the first main inequality used. This saves the argument by induction, etc. \n\nSignificance: The authors devote a significant amount of space to motivating the PSD factorization problem. As a tool in studying semidefinite representations of convex bodies, and as a tool in the correlation generation problem, PSD factorization is a well-known tool. It is, perhaps, less widely used as a concise representation of nonnegative data. The last paragraph of the introduction, discussing the interpretation of a PSD factorization in this context, is welcome, and may help stimulate further use of PSD factorizations in this context.\n\nThe main theoretical results of the paper focus on the case where the computed factors are positive definite. This is appropriate in the context of representing nonnegative data, but is not always appropriate in the context of PSD lifts of polytopes, for instance, where the factors are never positive definite (since the slack matrix always has some zero entries). To what extent it is possible to extend the algorithm and/or its analysis to the case where the matrices are not initialized to be positive definite? Certainly any discussion in this direction would be welcome, and could further extend the significance of the results in this submission.\n\nUpdate: after reading the other reviews and the response of the authors, I have not changed my (positive) assessment of this paper. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have adequately addressed the limitations of the work, by keeping it focused on precisely formulated mathematical statements (the limitations of which are inherent in their precise formulation). ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "maWbe1JgDcM", "reviewer": "Reviewer_k2ys", "summary": "This paper suggest an approach to PSD (positive-semidefinite) matrix factorisation (NMF), which generalizes the algorithm of Lee and Seung’s (Multiplicative Update algorithm, MUa) to what they call Matrix Multiplicative Update (MMU) algorithm.  It is shown that the fixed points of the MMU algorithm correspond to the KKT optimality conditions for the the problem (problem 8 in the main body of the paper). While there are papers showing similar result for the non-negative matrix factorization (https://arxiv.org/abs/2002.11323), the present paper simplified the proof technique. This is accompanied by code, which uses a \"damped variant\" of the algorithm (not covered by the analysis?), and is shown to work on an 361 x 2429 matrix.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper suggest an approach to PSD (positive-semidefinite) matrix factorisation (NMF), which generalizes the algorithm of Lee and Seung’s (Multiplicative Update algorithm, MUa) to what they call Matrix Multiplicative Update (MMU) algorithm. The idea of the MMU algorithm is to choose carefully a surrogate function and apply SM to the convex SDP problem (12) which is equivalent to the original formulation of the PSD factorisation problem (8). This is quite elegant. The main theorem (Theorem 1) shows that what essentially is the vectorisation operation from the space of symmetric rank r matrices to m dimensional real vector space allows one to have an objective function to be minimised that does not increase under the update rule they provide. \n\nOriginality:\n\nWhile there are some recent algorithms for PSD factorization (utilizing SDP hierarchies for NCPOP in https://link.springer.com/article/10.1007/s10208-018-09410-y, or heuristic approaches https://arxiv.org/abs/1707.07953), the proposed algorithm nicely combines some scalability and some guarantees. Notably, authors prove that \n-- if one starts with a feasible solution, the iterations lead to a feasible solution\n-- if one terminates with a feasible solution, then the fixed point of the iteration satisfies the KKT conditions for the the problem (problem 8 in the main body of the paper). \n\nQuality:\n\nThe analysis is of interest. Similar (or stronger) result had been available only for the non-negative matrix factorization (which is a special case of the problem with diagonal PSD matrices; https://arxiv.org/abs/2002.11323), plus the present paper simplified the proof technique. \n\nThis is accompanied by code, which uses a \"damped variant\" of the algorithm (not covered by the analysis?), and is shown to work on an 361 x 2429 matrix. The code is ok -- the key method is pasted below.\n\nClarity:\n\nThe 9-page writeup is mostly describing the motivation (pages 1-3) and related work (pages 4 and 5), in a rather confusing manner. (E.g. Section \"Existing work\" does not mention the \"MM Framework\".) While one may see the lengthy motivation as adding an \"interdisciplinary\" aspect, it may actually detract from the paper's value. \n\nStill, there is a proof of Theorem 1, which is quite elegant. The technical analysis is relatively well written, although the notation can be puzzling occasionally, for example the operator T introduced in (18) is not motivated. The definition of W is repeated several times but this actually eases the readability of the proof of Theorem 1. (Notice that Appendix A contains Lemma 3 which is a technical requirement for the proof of Theorem 1 and Appendix B contains Lemma 3 which is proves Theorem 2.)\n\nSignificance:\n\nFor the \"PSD-Rank community\", this will be an important paper. Considering that even papers by Nesterov on the subject (e.g. https://arxiv.org/abs/1707.07953) have less than 10 citations, though, it's not clear whether this is of wide interest. \n\nMinor comments:\n\nNeither the iteration complexity of the algorithm, nor the computational complexity of a single update is discussed in detail. \n\nThe discussion of related work on lines 173-180 is flawed: the use of Cardano's method for a cubic polynomial in a step of a coordinate descent on has been first proposed by https://arxiv.org/abs/1506.08568, rather than by Nesterov (https://arxiv.org/abs/1707.07953). See the paragraph \"The Closed-Form Step\" in https://arxiv.org/abs/1506.08568. \n\nA detailed table comparing the MMU algorithm with the relevant existing work (part of it described in Section 2) could have been provided.\n\nThe matrix geometric mean is used in the main algorithm but a fully geometric picture from the viewpoint of Riemannian geometry is not presented although briefly motivated around line 216.\n\nThere is a typo in line 161.\n\n\n--\ndef update_blockmultiplicative(A,B,X):\n    \"\"\"\n    Updates the factor B via the multiplicative updates\n    \"\"\"\n    \n    l = len(B)\n    d = B[0].shape[0]\n    \n    AAB = compute_blockAAB(A,B)\n    AtX = compute_block(A,X)\n    \n    # Invert AAB\n    AABinv = []\n    for k in range(l):\n        AABinv.append(np.zeros(B[k].shape))\n        for i in range(d):\n            AABinv[k][i,:,:] = np.linalg.inv(AAB[k][i,:,:])\n    \n    # Compute the geometric mean with a damping term\n    H = []\n    for k in range(l):\n        H.append(np.zeros(B[k].shape))\n        for i in range(d):\n            H[k][i,:,:] = matrix_gm(AABinv[k][i,:,:],B[k][i,:,:])\n            \n    # Update B_new\n    B_new = []\n    for k in range(l):\n        B_new.append(np.zeros(B[k].shape))\n        for i in range(d):\n            M = ( H[k][i,:,:] @ AtX[k][i,:,:] ) @ H[k][i,:,:].T\n            B_new[k][i,:,:] = M\n            \n    return B_new", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "This will not have any societal impact.", "review_text": "", "extra_scores": {"ethical_concerns": "There are no ethical issues.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "iddz4XSDaTb", "reviewer": "Reviewer_8tm7", "summary": "The paper studies algorithms for computing rank r positive semi-definite decompositions of a matrix with non-negative entries. They cast this non-convex problem as an alternating minimization problem and propose an update rule that is guaranteed to be non-increasing. They also show that all the fixed-points of their algorithm satisfy the first order optimality conditions for the objective. ", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Given an mxn matrix $X$ with non-negative entries, there is a rank $r$ positive semi-definite decomposition for the matrix $X$ if there exist $r \\times r$ positive semi-definite matrices $A_1, \\ldots, A_m$ and $B_1, \\ldots, B_n$, one for each row and column respectively such that for all $i,j$ : $X_{ij} = \\text{tr}(A_i B_j)$. This is a generalization of non-negative matrix factorization where all the r x r positive semi-definite matrices $A_1, ..., A_m$ and $B_1, ..., B_n$ have an additional constraint that they should be diagonal matrices. It is clear that relaxing the constraint of diagonality increases the expressibility and therefore the positive semi-definite rank of the matrix may be much smaller than the non-negative rank of a matrix. The problem is cast in an alternating-minimization framework where we fix $A_1,...,A_m$ and optimize over $B_1,...,B_n$ and then fix $B_1,...,B_n$ and optimize over $A_1,...,A_m$ until we reach a fixed point. Fixing matrices $A_1,...,A_m$ makes the objective to be an \"independent\" sum of contributions of the matrices $B_1,...,B_n$ and therefore the problem that needs to be solved is optimizing only over a single positive semi-definite matrix. This is the problem they solve. \n\nThey give an update rule for this problem and show that their updates never increase the objective value. The update rule is simple and falls-back to Lee-Seung's update rule when initialized with diagonal matrices. An advantage of the update rule presented by the authors is that it preserves the block diagonal structure of the initial iterates. The authors also argue that the fixed points of their update rule satisfy KKT conditons for the objective. They also run experiments using their algorithms and show that the solutions converge to an optimal decompositon with sufficient number of random initializations.\n\nStrengths:\n- Paper is written well with enough motivation for why the problem is important to be solved and how it relates to other important problems.\n- Introduces new techniques to the problem where the earlier works were only using Projected Gradient Method or Coordinate Descent.\n- Update rule falling back to the current algorithms for non-negative matrix factorization when initialized with diagonal matrices is good as this means their update rule is a natural extension of the earlier work.\n- Proofs are concise and clear.\n\nWeaknesses:\n- Experiments : It would be good to compare with the earlier coordinate descent and projected gradient methods.\n- Not having any bounds on how quick the algorithm converges to a fixed point. I understand it may be hard to get any analysis at all. Do you have any ideas here?\n\nQuestions:\n1. You mention in the experiments that with sufficient initializations, the algorithm for PSD factorization of a distance matrix is close to being exact. It does seem like all other initializations converge to some fixed point from your graphs. Is there anything interesting at all about the factorizations your algorithm converges to in these other initializations?\n\n2. Does rank of the optimal decomposition play any role at all? Say for example, in the distance matrix decomposition, the optimal PSD-matrices all have rank 1. Does the algorithm converge faster for instances with low rank optimal decompositions.\n\n3. I don't think it is fair to call the decomposition in Figure 3 as building blocks. They are just the reconstructions for the components in the eigen decomposition for that particular image unlike the building blocks that are same for all the images in the case of non-negative matrix factorization. That is why I don't agree the latent space being called \"infinite dimensional\". Am I missing something here? This isn't particularly a weakness -- I am just looking for a clarification.\n\nSuggestions:\n- $ \\text{tr}(XY)^2 \\le \\text{tr}(X^2)\\text{tr}(Y^2)$ is just the Cauchy-Schwarz inequality for the standard inner-product on the space of symmetric matrices. So you need not use von Neumann trace inequality and can obtain an easier proof.\n\nTypos/Mistakes:\n- Line 122 pace -> space\n- Line 161 \"how to the\" -> \"how the\"\n- Line 211 $f$ -> $F$\n- Define $\\mathcal{A}^{\\mathsf{T}}$ before its use\n- Line 440 In the display equation, $^2$ is missing\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "NA", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "NpvvcrbxamV", "reviewer": "Reviewer_wGDo", "summary": "This paper studies the problem fo finding a certain type positive semidefinite (PSD) factorization of a matrix with non-negative entries. The main result is a new algorithm for this problem which has some nice properties. The authors begin by describing various applications the PSD factorization problem. One special case is the problem of non-negative matrix factorization, which corresponds to the constraint that the PSD matrices are diagonal. Then, the authors described their method which falls in the majorization-minimization framework where the majorization is based on a matrix version of the geometric mean. The authors show that that the iterates of this algorithm preserve certain structural properties of the initialization point, and that in the case of diagonally restricted matrices, it reduces to the well-known Lee-Seung algorithm for non-negative matrix factorizations. Some numerical results are also provided. \n\n", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I enjoyed reading this paper. I think the authors did a good job of explaining the problem and some of its applications, and then clearly describing their contribution and the connection with prior work. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": " I think the discussion is appropriate.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The proposed method for computing a positive semidefinite factorization of a non-negative matrix and its analysis are original and of interest to the community. The paper is also clearly written and motivated. The authors are encouraged to take into account the reviewers' comments when preparing the revision.
Decision: Accept (Poster)