Dana Lahat
email: dana.lahat@irit.fr

Cédric Févotte
email: cedric.fevotte@irit.fr

Positive Semidefinite Matrix Factorization Based on Truncated Wirtinger Flow

Keywords: Positive semidefinite factorization, nonnegative factorizations, phase retrieval, low-rank approximations, optimization

This paper deals with algorithms for positive semidefinite matrix factorization (PSDMF). PSDMF is a recentlyproposed extension of nonnegative matrix factorization with applications in combinatorial optimization, among others. In this paper, we focus on improving the local convergence of an alternating block gradient (ABG) method for PSDMF in a noisefree setting by replacing the quadratic objective function with the Poisson log-likelihood. This idea is based on truncated Wirtinger flow (TWF), a phase retrieval (PR) method that trims outliers in the gradient and thus regularizes it. Our motivation is a recent result linking PR with PSDMF. Our numerical experiments validate that the numerical benefits of TWF may carry over to PSDMF despite the more challenging setting, when initialized within its region of convergence. We then extend TWF from PR to affine rank minimization (ARM), and show that although the outliers are no longer an issue in the ARM setting, PSDMF with the new objective function may still achieves a smaller error for the same number of iterations. In a broader view, our results indicate that a proper choice of objective function may enhance convergence of matrix (or tensor) factorization methods.

I. INTRODUCTION

This paper deals with improving the local convergence of algorithms for positive semidefinite matrix factorization (PSDMF) [START_REF] Gouveia | Lifts of convex sets and cone factorizations[END_REF], [START_REF] Fiorini | Linear vs. semidefinite extended formulations: exponential separation and strong lower bounds[END_REF]. In PSDMF, the (i, j)th entry x ij of a nonnegative matrix X ∈ R I×J is written as an inner product of two K × K symmetric positive semidefinite (psd) matrices A i and B j , indexed by i = [I], j = [J]:

x ij ∼ = A i , B j tr{A T i B j } = tr{A i B j } (1)
where [I] 1, . . . , I, tr{•} denotes the trace of a matrix, and ∼ = means equal or approximately equal, depending on the context. The minimal number K such that a nonnegative matrix X admits an exact PSDMF is called the psd rank of X [START_REF] Gouveia | Lifts of convex sets and cone factorizations[END_REF]. Each psd matrix A i and B j may have a different rank, denoted R Ai and R Bj , respectively. We shall sometimes refer to the psd rank K as the "outer rank" and to R Ai and R Bj as "inner ranks" [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF]. When A i and B j are constrained to be diagonal matrices ∀i, j, PSDMF becomes equivalent to nonnegative matrix factorization (NMF) (e.g., [START_REF] Paatero | Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values[END_REF]).

PSDMF was proposed as an extension to a fundamental result [START_REF] Yannakakis | Expressing combinatorial optimization problems by linear programs[END_REF] in combinatorial optimization that links NMF with geometry. Since it was first introduced, PSDMF attracted a significant amount of attention and has found applications in numerous fields, including combinatorial optimization [START_REF] Fiorini | Linear vs. semidefinite extended formulations: exponential separation and strong lower bounds[END_REF], [START_REF] Kaibel | Extended formulations in combinatorial optimization[END_REF], quantum communication complexity [START_REF] Fiorini | Linear vs. semidefinite extended formulations: exponential separation and strong lower bounds[END_REF], quantum information theory [START_REF] Jain | Efficient protocols for generating bipartite classical distributions and quantum states[END_REF], and recommender systems [START_REF] Stark | Recommender systems inspired by the structure of quantum theory[END_REF]. Despite the increasing interest in PSDMF, only a few algorithms have been proposed so far, e.g., [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF], [START_REF] Stark | Recommender systems inspired by the structure of quantum theory[END_REF], [START_REF] Glasser | Expressive power of tensor-network factorizations for probabilistic modeling[END_REF]. In this paper, we do not focus on a specific application but on understanding some of the optimization properties of PSDMF, in a more general perspective.

Our work is motivated by a recent result [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF] showing a link between PSDMF and affine rank minimization (ARM) (e.g., [START_REF] Zheng | A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements[END_REF], [START_REF] Tu | Low-rank solutions of linear matrix equations via Procrustes flow[END_REF]), and in particular with phase retrieval (PR) (e.g., [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF]), which can be considered as a special case of ARM. We explain their link to PSDMF briefly in Section II, which is dedicated to background material.

In [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF], we presented a new type of local optimization scheme for PSDMF that is based on alternating block gradient (ABG). In this ABG scheme, each subproblem is based on Wirtinger flow (WF) [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF] or ARM [START_REF] Zheng | A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements[END_REF], [START_REF] Tu | Low-rank solutions of linear matrix equations via Procrustes flow[END_REF], depending on the values of the inner ranks, and the objective function is based on quadratic loss. This paper is motivated by truncated Wirtinger flow (TWF) [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], a variant of WF obtained by replacing the quadratic objective function with a non-quadratic one, together with special regularization based on truncating excessive values from the gradient and the objective function. These excessive values do not appear in the quadratic objective. TWF improves the convergence rate also in the noise-free case, and requires appropriate initial conditions. In this paper, starting from Section III, we build an ABG algorithm in which each subproblem is based on the principles of TWF, and study to which extent the ideas of TWF can be applied to PSDMF.

Our main result is that in the case where all inner ranks are equal to 1, a setting equivalent to locally minimizing a TWF criterion in each subproblem, we observe phenomena analogous to those in TWF as predicted by our geometrical analysis in Section IV, provided we are in the basin of attraction of a global minimum. This special case of PSDMF, known as Hadamard square root or phaseless rank, is of particular interest (e.g., [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF], [START_REF] Goucha | The phaseless rank of a matrix[END_REF], [START_REF] Fawzi | Positive semidefinite rank[END_REF]). When some of the inner ranks are larger than one, there is no need for truncation, as we show in Section IV. Yet, we still observe faster convergence to the global minimum with our non-quadratic objective, a matter than needs to be further looked into.

Although we focus on the Poisson log-likelihood, as in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF], we keep in mind that the same principles are likely to hold also for other objective functions; see further discussion in Section VI. To the best of our knowledge, this paper is the first to suggest that changing the objective function within the same algorithmic framework in a matrix factorization problem can result in improved numerical properties, in the noise-free case, provided proper initialization and regularization (e.g., truncation) when necessary.

In this paper, we use real-valued notations. However, all our results extend straightforwardly to the complex-valued case.

II. A LINK BETWEEN PSDMF AND PHASE RETRIEVAL

In PSDMF applications, the psd matrices are often of low rank, i.e., R Ai , R Bj < K (e.g., [START_REF] Fawzi | Positive semidefinite rank[END_REF]). Hence, the psd matrices can be written as

A i U i U T i and B j V j V T j , (2)
where U i ∈ R K×R A i and V j ∈ R K×R B j are full-rank factor matrices (factors, for short). With this change of variables, the model in (1) can be written as

x ij ∼ = U i U T i , V j V T j = tr{U i U T i V j V T j } = U T i V j 2 F .( 3
)
In the special case where J = 1, (3) can be written as

y i ∼ = U T i V 2 F , (4)
where we replaced V j with V and x ij with y i . When U i are given, finding V from y i , i = [I], in (4) is equivalent to ARM (e.g., [START_REF] Zheng | A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements[END_REF], [START_REF] Tu | Low-rank solutions of linear matrix equations via Procrustes flow[END_REF]). In this case, U i can be regarded as sensing matrices and V as the signal of interest. If all factors are vectors: v ∈ R K and u i ∈ R K ∀i, (4) can be written as

y i ∼ = u i u T i , vv T = |u T i v| 2 = | u i , v | 2 , i = [I] . (5)
When all vectors u i are given, the problem of estimating v from the phaseless observations y i in ( 5) is known as (generalized) PR [START_REF] Gerchberg | A practical algorithm for the determination of phase from image and diffraction plane pictures[END_REF]. In a PR context, u i are sometimes referred to as sensing vectors and v as a signal of interest. We conclude that PR [START_REF] Yannakakis | Expressing combinatorial optimization problems by linear programs[END_REF] and ARM (4) are special cases of PSDMF [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF].

Back from PR or ARM to PSDMF: Consider a random noise model acting independently on each observation in (4). Then, for a broad range of noise models, the log-likelihood of the model parameters given an observation can be written as

d(y i , U T i V 2 F )
, where d(ζ, η) is a divergence between any two nonnegative numbers ζ, η, and where 2 leads to the quadratic objective function. Then, a maximum likelihood estimate (MLE) to the ARM problem in (4) can be found by minimizing, with respect to (w.r.t.) V, the objective function

d(ζ, η) ≥ 0 with d(ζ, η) = 0 if and only if ζ = η. For example, d(ζ, η) = 1 2 (ζ -η)
I i=1 d(y i , U T i V 2 F ) . (6)
Similarly, based on (3), an MLE to PSDMF can be found by minimizing, w.r.t. U i and V j , ∀i, j, the objective function

I i=1 J j=1 d(x ij , U T i V j 2 F ) . (7)
Looking at ( 6) and ( 7), we see how any ARM method for (6) can be used to optimize [START_REF] Jain | Efficient protocols for generating bipartite classical distributions and quantum states[END_REF] : given {U i } I i=1 , apply ARM subsequently to update each V j in [START_REF] Jain | Efficient protocols for generating bipartite classical distributions and quantum states[END_REF]. Then change roles: fix {V j } J j=1 and apply the method to update U i subsequently for each i. Repeat until a stopping criterion is achieved. The alternating optimization strategy for PSDMF has already been proposed, e.g., in [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF], [START_REF] Stark | Recommender systems inspired by the structure of quantum theory[END_REF], [START_REF] Glasser | Expressive power of tensor-network factorizations for probabilistic modeling[END_REF], [START_REF] Basu | Computing approximate PSD factorizations[END_REF] . The link between these methods and PR is pointed out in [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF]. The ABG method in [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF] is based on decreasing the objective function in ( 6) by gradient descent with a quadratic objective function, as in [START_REF] Zheng | A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements[END_REF], [START_REF] Tu | Low-rank solutions of linear matrix equations via Procrustes flow[END_REF]. Addressing the PR problem associated with ( 5) via gradient descent with a quadratic objective function is termed WF [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF].

III. EXTENDING TWF TO ARM

Starting from this section, we gradually build an ABG method for PSDMF based on TWF [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]. For the sake of simplicity, we focus on the Poisson log-likelihood, as in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]. In this case, the objective function in ( 7) is based on the Kullback-Leibler divergence (KLD):

d KL (ζ, η) ζ log( ζ η ) - ζ + η.
PSDMF with a Poisson log-likelihood thus consists of decreasing the objective function

f = I i=1 J j=1 U T i V j 2 F -x ij log( U T i V j 2 F ) + x ij log(x ij ) -x ij (8)
whose gradient w.r.t.

V j ∈ C K×R B j is ∇ Vj f = 2 I i=1 U T i V j 2 F -x ij U T i V j 2 F U i U T i V j . (9)
In the complex-valued case, change • T to • H and divide ∇ Vj f by 2. For a specific j, (8) and ( 9) are the low-rank generalizations of the Poisson objective and gradient in TWF [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]. One can thus obtain ARM and PSDMF algorithms based on the Poisson log-likelihood by replacing the quadratic objective function and gradient in [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF]- [START_REF] Tu | Low-rank solutions of linear matrix equations via Procrustes flow[END_REF] with those in ( 8) and ( 9). Our first step, of extending TWF from rank-1 signal recovery to ARM, and of stating PSDMF with a Poisson loglikelihood, is thus complete. We point out that [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]Sec. 8] mention extending TWF to ARM, but only when the sensing vectors are rank-1. However, due to the symmetric nature of the subproblems, in the general PSDMF case we do not restrict the inner ranks in any of the factors.

IV. A GEOMETRIC EXPLANATION

The second step in extending TWF to PSDMF is revisiting the truncation idea. For this purpose, we follow in the steps of [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF] and inspect the locus of the gradient in [START_REF] Glasser | Expressive power of tensor-network factorizations for probabilistic modeling[END_REF] at the proximity of the solution.

We begin with the case where all factors are vectors, which is equivalent to the setting in TWF. The locus of the gradient in (9) for a specific value of j is illustrated in Fig. 1 (left). This figure is a reconstruction of [15, Figure 2.1] using the same parameters. It shows that the desired values of the gradient, around the origin, are rather small (see inset) compared with the outliers, which are few but very large. These large and undesired gradient values occur only in specific directions. A similar examination will reveal that the quadratic objective function does not have outliers (the gradient for the quadratic case is similar to [START_REF] Glasser | Expressive power of tensor-network factorizations for probabilistic modeling[END_REF], but without the denominator).

The idea of TWF [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF] is to avoid these outliers and thus better focus on the desired values. Chen and Candès provide coefficients that allow to detect elements in the gradient and the objective function that are associated with undesired contributions. They show that this procedure achieves faster convergence and better resilience to noise than without truncation, provided an appropriate initialization.

Next, we study the outlier problem in the general low-rank case. As an example, we look at the locus of the gradient in ( 9) when U i are now 2 × 2 "sensing matrices" instead of 2 × 1 vectors. The "signal" vector is the same as in the previous plot. Figure 1 (right) depicts the locus. We observe that the outliers are significantly smaller and fewer. As in the one-dimensional example in Fig. 1 (left), the outliers occur only in specific directions. Otherwise, the shape of the locus did not change much. We observed similar trends with other values of ranks and matrix sizes. This indicates that truncation will have little to no effect as one gets further away from the case where all inner ranks are equal to one. Our numerical results in Section VI validate this. Comparison of size of outliers in gradient for different dimensions of U i . Each plot depicts the locus of the gradient

T i z 2 F -U T i v 2 F U T i z -2 F U i U T i v with v = 2.7 8 T , z = 3 6
T . The entries of U i are drawn independently from N (0, 1) and then normalized s.t. U i F = 1. Each plot depicts 10 3 independent draws of U i . In each plot, the axes are fit to the data, no zoom and no cropping (except inset in left plot). U i is a 2 × 1 vector as in [15, Fig. 2.1] (left) and a 2 × 2 matrix (right).

V. A PSDMF ALGORITHM BASED ON TWF Algorithm 1 describes one subproblem in our TWF-based ABG scheme for updating {V j } J j=1 given {U i } I i=1 . In the special case where all inner ranks are equal to one, Algorithm 1 is essentially equal to TWF [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF].

The outliers discussed in Section IV are avoided by using a regularized version of the gradient, ∇ Vj f tr , in Algorithm 1. The regularization

I E ij 1 ∩E ij 2
applied to the ith element of the gradient sets it to zero (i.e., does not add it to the sum) if at least one of the following events occurs:

E ij 1 α lb ≤ U T i V j F V j F ≤ α ub (10) E ij 2 U T i V j 2 F -y i ≤ α h U T i V j F • q j -x j 1 I V j F
where x j is the jth column of X and q j is a vector whose ith entry is equal to

U T i V j 2
F . The truncation rules in (10) generalize the truncation rules in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF] to the low-rank case.

If backtracking line search is used, [15, Sec. 2.3] suggest a dedicated truncation procedure. Algorithm 2 extends it to the general low-rank case. In Algorithm 2, in analogy to [15, Sec. 2,3], P j ∇ Vj f tr , and

f (V j ) i∈ T (Vj ) U T i V j 2 F -x ij log( U T i V j 2 F ) (11a)
where

T (V j ) i | U T i V j F ≥ α lb V j F and U T i P j F ≤ α p P j F . (11b)
Equation ( 11) describes the truncation rules for the objective function in the backtracking line search procedure in Algorithm 2, which generalize those in [15, Sec. 2.3] to the case where both "sensing vectors" and "signal" may be matrices.

Algorithm 1 Updating {V j } J j=1 based on TWF Input: X ∈ R I×J + , U 1 , . . . , U I , V 1 , . . . , V J , α ub , α lb , α h . Output: V 1 , . . . , V J .
1: τ ← Choose initial step size 2: for j = 1 : J do 3: τ = β bck τ 3: end while

∇ Vj f tr ← 2 I i=1 ( U T i Vj 2 F -xij ) U T i Vj 2 F U i U T i V j I E ij 1 ∩E ij

VI. NUMERICAL EXPERIMENTS

We compare our proposed Poisson-based ABG, with and without truncation, to ABG [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF] and to the two variants of coordinate descent (CD) [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF], cyclic and Gauss-Southwell (GS), the latter with parameter α GS = 0.5. The competing methods: ABG and CD, use a Gaussian objective function. We set α bck = 0.1, β bck = 0.35 for backtracking line search (Algorithm 2) for ABG [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF] and for our proposed method. The initial step size for the backtracking line search is chosen by evaluating numerically the Lipschitz constant of the gradient, as described in [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF]. The truncation parameters are identical to those suggested as default in [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]Sec. 2.3]:

α lb = 0.1, α ub = 5 = α p , α h = 6.
In our first numerical experiment, whose results are shown in Fig. 2, the input is a 20 × 20 nonnegative matrix (I = 20 = J) generated from factors whose entries are drawn independently from N (0, 1). The factors have K = 5 rows. We test two settings with different values of inner ranks:

(i) R A = 1 = R B (ii) R A = 2, R B = 1, where R Ai = R A ∀i and R Bj = R B ∀j.
In each of the 200 Monte Carlo (MC) trials, the input matrix X and initial factors are drawn anew, and all algorithms use the same initialization. Our input matrix provides 400 observations, and the corresponding numbers of free model variables for each setting are 175, and 255, respectively. Hence, these are low-rank scenarios. The algorithms are given the true values of the ranks. Hence, achieving a global minimum is equivalent to exact factorization. The stopping criterion is 60 iterations. We define as one iteration t a pair of subproblems that update both sets of factors, indexed by i and j. Our algorithms are implemented in Matlab. We use the CD code given in https://sites.google.com/site/exactnmf/. We scale the initial factors, before running the algorithms, using the procedure in [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF].

Initialization: As noted by [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF], the loss function for PSDMF [START_REF] Jain | Efficient protocols for generating bipartite classical distributions and quantum states[END_REF] is in general highly nonconvex. As a result, the chance of achieving a global minimum with random initialization is very small in low-rank settings and in problems whose size is not sufficiently small. Spectral initialization (e.g., [START_REF] Zheng | A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements[END_REF]- [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF], [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]), which is essential for TWF, is not applicable to PSDMF because the roles of "sensing matrix" and "signal" change at every subproblem, and both are unknowns that have to be estimated from the data. For these reasons, and for the sake of our preliminary study of our proposed methods, we chose another strategy, which initializes the algorithms within a sufficiently small region around the true solution. Each entry in an initialization factor is the sum of (1 -p) times the true entry and p times an independently drawn number from N (0, 1), and the overall is normalized such that the variance of this new random variable (since we draw the entries of the true factors from N (0, 1), too) is one. This normalization is chosen to guarantee that we work in a regime that allows for exact recovery in WF/TWF setting (e.g., [START_REF] Candès | Phase retrieval via Wirtinger flow: Theory and algorithms[END_REF], [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]). We set p = 0.1, as this value turned out to guarantee convergence of (most of) our tests of the truncated Poisson-based objective.

The purpose of the first experiment is to compare the performance of the two objective functions when the algorithm is already within its (assumed) region of convergence (ROC). Due to lack of space, other experiments, such as determining the size of the ROC of each method and the rate of success with respect to the various model parameters (matrix size, inner and outer ranks), as well as more realistic initialization strategies, are beyond the scope of this paper.

Results: Figure 2 shows our results. Each boxplot represents the logarithm of the normalized model fit error log 10 ( X-Xt F Xt F ), for each method. The plot on the left in Fig. 2 R A = 2, R B = 1 corresponds to our ARM extension of TWF.

The results for R A = 1 = R B are generally in accordance with TWF [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF]. The smallest error is achieved by our TWFbased method. Non-truncated Poisson stopped with a relatively large error, indicating poor convergence due to outliers. ABG, which is the counterpart of our Poisson-based methods but with a Gaussian objective function, performs better than the Poisson-based objective with no truncation, but the error is still relatively large. The competing CD methods achieved the second-best results. We point out that the computational complexity of the Poisson log-likelihood variant of ABG is similar to that of quadratic ABG; it is identical in the case of all inner ranks equal to 1. The per-subproblem computational complexity of ABG [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF] and its Poisson log-likelihood variant is smaller than that of CD [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF].

The outliers (red crosses) indicating relatively large errors in the truncated Poisson case probably imply that the initialization is close to the boundary of the ROC, and thus not all trials converged towards a global minimum. In general, the results in Fig. 2 indicate that the outlier phenomenon observed in TWF [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF] indeed appears also in PSDMF, despite the evident differences between TWF and PSDMF.

We now turn to the low-rank case, addressed in the right plot of Fig. 2. In this setting, we do not observe an advantage of truncated versus non-truncated Poisson, in accordance with our geometric analysis in Section IV and Fig. 1. Still, the Poisson objective managed to achieve a smaller error per fixed number of iterations than its Gaussian counterpart. In this setting, the CD methods achieved a smaller error than our gradientbased methods for the same number of iterations. However, one should keep in mind that CD has a larger computational complexity per iteration, see [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF], [START_REF] Lahat | Positive semidefinite matrix factorization: A link to phase retrieval and a block gradient algorithm[END_REF]. Similar trends were observed when we increased the values of the inner ranks.

Our second experiment shows the effect of the Poisson log-likelihood on a different data set, a slack matrix of the regular 10-gon (a ten-sided polygon), denoted as S 10 . S 10 is a 10 × 10 matrix, whose explicit form can be found, e.g., in [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF]. PSDMF of geometric data, and slack matrices in particular, is a central application of PSDMF [START_REF] Gouveia | Lifts of convex sets and cone factorizations[END_REF], [START_REF] Vandaele | Algorithms for positive semidefinite factorization[END_REF]. Our input matrix X is S 10 normalized such that X F = 1. Here, we use random initialization: initial factors with dimensions K = 5, R A = 1, R B = 3, whose entries are drawn independently from N (0, 1).

Results: Figure 3 (left) shows our results in terms of the number of iterations each algorithm needed to achieve the same normalized model fit error

1 2 X-Xt 2 F Xt F < 10 -7 .
We run 50 MC trials, each with a new initialization that is the same for all algorithms. Figure 3 (right) illustrates the evolution of the model fit error of each algorithm for one run in this setting; here, all algorithms are implemented in Matlab. Figure 3 does not show the results for truncated Poisson because this method did not converge properly, in this setting. The reason is probably not only the lack of outliers in this low-rank setting, as explained in Section IV, but also due to the fact that the truncation constants proposed by [START_REF] Chen | Solving random quadratic systems of equations is nearly as easy as solving linear systems[END_REF] were designed for factors whose entries are drawn from the standard normal distribution N (0, 1), an assumption that does not hold in the case of geometric data. In Fig. 3 (left), we observe that our proposed method with the Poisson objective (and no truncation) often achieves the desired model fit error with the fewest number of iterations; note that the per-iteration complexity must be taken into account. The CPU time comparison in Fig. 3 (right) is only illustrative because our Matlab implementation is not optimized for speed. Still, it indicates that the proposed method can converge much faster than its competitors, consistently with Fig. 3

(left).

Conclusion:

Our results provide proof of concept to the idea of improving the convergence of exact PSDMF by changing the objective function. However, our results depend heavily on the type of data and on initialization within the ROC, in certain cases. There is still need for further study of the ROC of each PSDMF method and its dependence on the model parameters, ranks and dimensions. Our results call for further theoretical validation. Our results provide further motivation to design more efficient PSDMF methods with better global convergence. A possible direction is more recent techniques used in PR, see, e.g., [START_REF] Yang | Parallel coordinate descent algorithms for sparse phase retrieval[END_REF] and references therein.

We mention that Gaussian and Poisson log-likelihoods are special cases of the β-divergence, a popular tool in NMF (e.g., [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF]). Extending our algorithm to β-divergence is straightforward: the only change in the gradient in [START_REF] Glasser | Expressive power of tensor-network factorizations for probabilistic modeling[END_REF] is the power of the denominator, which becomes (2-β) instead of 1. One can readily verify, by plotting diagrams similar to those in Fig. 1 for other values of β, that the size of the outliers is continuous in the value of β. It is thus worth exploring the benefit of truncation and use of other objectives also for other values of β in all concerned applications, namely PR, ARM, and PSDMF. This paper addressed only the noise-free case: study of the usefulness of PSDMF with non-Gaussian objectives in the presence of noise is left for follow-up work.

In a broader view, the idea of accelerating factorizations by changing the objective function and truncating outliers might allow to accelerate other methods, and thus serve purposes beyond providing better resilience to non-Gaussian types of noise.

Fig. 1 .

 1 Fig. 1.Comparison of size of outliers in gradient for different dimensions of U i . Each plot depicts the locus of the gradientT i z 2 F -U T i v 2

2 4 : 5 : 2 F do 2 :

 24522 µ j ← Choose step size (e.g., τ or Algorithm 2)V j ← V j -µj I ∇ Vj f tr 6: end for Algorithm 2 Backtracking line search Input: τ , α bck , β bck , α lb , α p 1: while f (V j -τ P j ) > f (V j ) -α bck τ P j

Fig. 2 .

 2 Fig. 2. Comparison of model fit error for a fixed number of iterations. PSDMF of a 20 × 20 matrix with randomly-generated factors. Initialization with perturbed true factors. Left: R A = 1 = R B , right: R A = 2, R B = 1.

Fig. 3 .

 3 Fig. 3. PSDMF of S 10 with K = 5, R A = 1, R B = 3. Random initialization. Left: Comparison of number of iterations for the same model fit error. Right: Evolution of error versus CPU time in one run.

This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No. 681839 (project FACTORY).

ACKNOWLEDGMENT

The authors would like to thank Vincent Y. F. Tan for bringing the work on TWF to our attention.