Efficient algorithms for regularized Poisson Non-negative Matrix Factorization
Nathanaël Perraudin Adrien Teutrie Cécile Hébert Guillaume Obozinski




1]Nathanaël Perraudin
2]Adrien Teutrie
3]Cécile Hébert
1]Guillaume Obozinski
[1]Swiss Data Science Center, EPFL and ETH Zürich, Andreasstrasse 5, Zürich, 8050, Zürich, Switzerland
[2]Unité Matériaux et Transformations, UMR-CNRS 8207, Université de Lille, Cité scientifique, Bâtiment C6, Villeneuve d’Ascq, 59655, Nord, France
[3]Electron Spectrometry and Microscopy Laboratory, Institute of Physics, EPFL, Bâtiment PH, Station 3, Lausanne, 1015, Vaud, Switzerland






Efficient algorithms for regularized Poisson Non-negative Matrix Factorization
    [
    
==============================================================================
































We consider the problem of regularized Poisson Non-negative Matrix Factorization (NMF) problem, encompassing various regularization terms such as Lipschitz and relatively smooth functions, alongside linear constraints. This problem holds significant relevance in numerous Machine Learning applications, particularly within the domain of physical linear unmixing problems. A notable challenge arises from the main loss term in the Poisson NMF problem being a KL divergence, which is non-Lipschitz, rendering traditional gradient descent-based approaches inefficient. In this contribution, we explore the utilization of Block Successive Upper Minimization (BSUM) to overcome this challenge. We build approriate majorizing function for Lipschitz and relatively smooth functions, and show how to introduce linear constraints into the problem. This results in the development of two novel algorithms for regularized Poisson NMF. We conduct numerical simulations to showcase the effectiveness of our approach.




  
*Disclaimer
This document is a technical report and has not undergone peer review. 
The findings and conclusions presented herein are solely based on the authors' research and analysis.
We apologize for any potential errors or shortcomings in the content.



§ INTRODUCTION

The problem of factorizing a matrix ≈ as Non Negative components ≥0,≥0 is central in many Machine Learning (ML) applications <cit.>.
The motivation for performing such a factorization is that  is often associated with a probability distribution density of the form P_(,)=P̃_().
Typically, the optimal decomposition is found by minimizing the negative log-likelihood of that distribution:

    _≥0,≥0-∑_i,jlog(P_(,))


If  is assumed to be perturbed with Normal noise, we obtain a Gaussian distribution, i.e. P_(,)∝ e^--^2, and we end up with the classic non-negative matrix factorization (NMF) problem <cit.>, where the quadratic function-^2 is minimized. 
For other distribution families, and in particular exponential families, a log term often appears. In particular, the Poisson negative log-likelihood model <cit.> leads to a loss of the form

    ℒ_(,):=-⟨,log()⟩ +⟨1,⟩∝-∑_i,jlog(P_(,)),

where the inner product over matrices is the Frobenius inner product defined as ⟨A,B⟩ =∑_ija_ijb_ij. 



  
Regularized Poisson Non Negative Matrix Factorisation
In many problems (see, for example, <cit.>), additional prior information about the matrices , is known. For example, it might be known that the columns of  are smooth, or that the rows of  are sparse. One might also be interested in normalizing to unity the columns of  or the rows of  because they might quantify physical quantities for which normalization is necessary. For example, in the analysis of hyperspectral imaging data, the images  are assumed to be smooth and the components  are summing to the unity <cit.>.
Typically, this information can be encoded via an extra regularization term R(,) and/or additional constraints ∈𝒞_1,∈𝒞_2.
This leads to the general optimization problem we solve in this contribution:

    _ ,     ℒ_(,)+R_W()+R_H(),
    subject to     ≥ϵ,≥ϵ,  e_H^⊤=1 or e_W=1^⊤

where ϵ>0. The different colors emphasize the changes compared to the traditional problem of <cit.>. First, in violet, we slightly simplify the problem by imposing strict non-negtativity. While, this is not strictly necessary[The case ϵ=0 could be handled with an approach similar to <cit.>.], this assumption significantly simplify our analysis. We believe that handling ϵ=0 is an unnecessary complication, as the results with ϵ close to machine precision will be practically identical. Second, in light blue, we consider the case where the constraints are linear, specifically e_H^⊤=1 or e_W=1^⊤.
In general, it is only meaningful to use one of the constraints, as it will fix the ratio between  and . We note that this includes the simplex constraint when e=1. 
Third, in brown, we consider regularizations R_W()+R_H() of the form:

    r(x)=s_L(x)+s_R(x)+∑_j=1^ns_C(x_j),

where x is the vector of a row of  or a colum of H., i.e R_W()=∑_ir_H(w_i) and R_H()=∑_jr_W(h_i).



    
    
    




  * The term s_L is assumed to be σ_L gradient Lipschitz, i.e., there exists a σ_L such that

    ∇ s_L(x)-∇ s_L(y)_2≤σ_Lx-y_2for all x,y∈𝒞.

Alternatively, this condition could be rewriten as 

    s_L(y) ≤ s_L(x) 
    + ⟨∇ s_L (x),y-x⟩
     + σ_L/2x -y_2^2 for all x,y∈𝒞.

In this contribution, we will consider in particular s_L(x)=x^⊤Δx,
where Δ is the Laplacian operator, favoring smoothness in the columns of . In this case σ_L=2 λ_max(Δ).

  * The term s_R(x) is assumed to be σ_R relatively smooth with respect to κ(x)=-1^⊤log(x). Relative smoothness is a generalization of Lipschitz smoothness, and is defined as follows <cit.>:

    f(y)≤ f(x)+⟨∇ f(x),y-x⟩ +σ_Rℬ_κ(y,x)

where ℬ_κ is the Bregman divergence <cit.> associated with κ <cit.>:

    ℬ_κ(y,x):=κ(y)-κ(x)-⟨∇κ(x),y-x⟩for all x,y∈𝒞.

While Lipschitz functions can be upper bounded with quadratic functions, we observe from (<ref>) that relative smoothness allows us to upper bound a function with the Bregman divergence of a function κ. This allows us to use a much wider range of functions to regularize our problem, and in particular non-gradient Lipschitz ones.
As an example, let us consider κ(x)=-1^⊤log(x), with its Bregman divergence:

    ℬ_κ(y,x)=∑_i(y_i/x_i-log(y_i/x_i)-1)

One can observe that the objective function (<ref>) is relatively smooth with respect to κ(x)=-1^⊤log(x).
This term could also be used to introduce soft contraints such as a log-barrier:
s_R(x)=-1^⊤log(x-ϵ1)
for x>ϵ and +∞ otherwise. 

  * Eventually, s_C(x) is a smooth point-wise concave function (i.e. -s_C is convex). A typical example of this regularisation could be s_C(x)=log(x+α^-1)
which favor sparsity in the vector x without penalizing large values too heavily. Its slope starts at α for x=0 and
tends to 0 for x→∞. 

Our approach can handle regularizations of the form R(,). However, for simplicity of notation, we restrict ourselves to separable regularizations.

The fidelity term ℒ_(,) is convex in  and in , but jointly non convex. We note that depending on the regularization and the constraints, multiple equivalent scaled solutions (stationary points) could exist, i.e., =^'^' for ^'=α^-1 and ^'=α. 
However, this is generally not the case with additional constraints or regularizations. 




  
Why is this problem challenging?
In general, Poisson Non Negative Matrix factorization, i.e. minimizing ℒ_(,) is challenging because it is not
gradient Lipschitz[The function ℒ_ is, according to the definition, gradient Lipschitz because of the constraint ,≥ϵ. However, in practice, the constant ϵ is chosen to be so small that its actual Lipschitz constant is too large to be useful.] despite being differentiable for ,>0.
In practice, this implies that there exists no fixed learning rate that ensures convergence of gradient descent and that line search would have to be used.
For that reason, solving the more general problem (<ref>) is a difficult task, and to the best of our knowledge, there are no existing algorithms that can be directly applied to it. 
Although there exist many algorithms to solve the traditional Poisson NMF problem <cit.>, none of them focuses on the regularized case (see the related work Section <ref>). 
Our main contribution is to fill this gap by providing multiple algorithms that minimize (<ref>) for a wide range of regularizationsR(,) described by (<ref>) and some linear constraints.



  
Our approach
A natural approach to minimize (<ref>) is to optimize for each variable , at a time. 
For example, Block Coordinate Descent (BCD) can be expressed as

    ^t+1←_∈𝒞_ℒ(,^t),


    ^t+1←_∈𝒞_ℒ(^t+1,).

This type of iterative scheme ensures that the loss does not increase between iterations and has been successfully used for the L2 case. However, in the Poisson case, there is no closed-form solution for problems (<ref>) and (<ref>), making this approach generally computationally expensive.

Fortunately, in practice, one does not need to find the global minima of (<ref>) and (<ref>) at each iteration. Instead, using a Block Successive Minimization (BSUM) algorithm <cit.>, it is sufficient to minimize approximations of ℒ which are locally tight upper bounds of ℒ.
To use the BSUM efficiently, these approximation functions need to have three properties: (1) to satisfy the hypotheses of the BSUM Theorem <cit.>, (2) to be as tight as possible, and (3) to be easy to optimize, i.e. to lead to a closed-form solution for each subproblem.

Our contributions can be summarized as follows. We show how regularized Poisson NMF can be efficiently solved using BSUM. We derive tight upper bounds for multiple regularizers and compare our approach with traditional algorithms. We also propose a simple way to introduce linear constraints into the problem and suggest using line search to build even tighter upper bounds. Finally, we propose multiple algorithms for regularized Poisson NMF and conduct numerical simulations to demonstrate the effectiveness of our approach.



  
Outline of this contribution

In Section <ref>, we provide a review of the literature. In Section <ref>, we clarify the notation and provide the necessary definitions for the BSUM Theorem <cit.> that will be used for the convergence of our algorithm. In Section <ref>, we develop convenient approximations of the objective and regularization functions leading to sub-problems with closed-form solutions. In Section <ref>, we explore how to modify the optimization scheme to introduce generalized simplex constraints. In Section <ref>, we present our algorithms. Section <ref> provides numerical applications of our algorithm, and Section <ref> concludes this contribution.




§ RELATED WORK











































































































































































































































































































  
Applications of Poisson Distribution likelihood Maximization

The maximization of likelihood in a Poisson distribution finds relevance in various applications, prompting the resolution of the problem outlined in Equation (<ref>).

Many such applications arise in the domain of physical constrained linear unmixing problems <cit.>. Some noteworthy instances encompass:
1. Scanning transmission electron microscopy (STEM) <cit.>,
2. Hyperspectral Raman and optical imaging <cit.>,
3. Tensor SVD applied to denoise atomic-resolution 4D scanning transmission electron microscopy <cit.>,
and 4. Non-local Poisson PCA denoising <cit.>.
It is noteworthy that many of these applications predominantly employ the L2 case, which offers a comparatively simpler solution. As a result, data is often renormalized to convert Poisson distributions into Gaussian distributions <cit.>. Nevertheless, the efficacy of these applications could be significantly enhanced by the development of algorithms tailored explicitly for the Poisson case.

Furthermore, Problem (<ref>) also surfaces in hyperspectral image denoising, where noise is assumed to follow a Poisson distribution <cit.>. In the domain of text mining, the Poisson distribution assumption is frequently utilized for modeling word occurrences based on latent variables such as categories, leading to the problem formulation depicted by (<ref>) <cit.>. Additionally, within the context of recommender systems, several matrix factorization problems can be reformulated into the structure of (<ref>) <cit.>.



  
Other Optimization Approaches

The literature offers a limited number of optimization methods suitable for addressing the problem presented by Equation (<ref>). This constraint arises from the requirement of many optimization techniques to have a continuously differentiable gradient Lipschitz function. Examples of such techniques include gradient descent <cit.>, perturbed gradient descent <cit.>, nonlinear conjugate gradient method <cit.>, various proximal point minimization algorithms <cit.>, and second-order methods like the Newton-CG algorithms <cit.>.

One potentially attractive direction is the utilization of Proximal Alternating (Linearized) Minimization (PALM) <cit.> or Proximal Alternating Minimization (PAM) <cit.>. These algorithms are designed to solve problems of the form:

    _,R_W()+R_H()+ℒ(,)

PALM employs a Gauss-Seidel iteration scheme, consisting of the following sub-problems:

    ^k+1     = _  R_W ()+ℒ( ,^k) + c_W  - ^k _2^2  
    ^k+1    = _  R_H()+ℒ(^k+1,) + c_H  - ^k _2^2
 

Unfortunately, PALM requires the objective function ℒ to possess a Lipschitz gradient, which is not the case in our scenario. Additionally, the Gauss-Seidel iterations generally lack a closed-form solution, resulting in a slow algorithm with sub-iterations.

To overcome the non-Lipschitz gradient issue, Bregman gradient descent (B-GD) <cit.> can be considered. This type of algorithm has been extended to alternating minimization <cit.>. Such an approach can be adapted to our case, as the objective function (<ref>) exhibits relative smoothness for most regularization scenarios (see (<ref>)). However, this optimization scheme involves non-tight majorization functions, leading to slow convergence, as discussed in Section <ref>.

Given the presence of two blocks of variables, Block Coordinate Descent (BCD) algorithms <cit.> naturally emerge as a potential solution. In fact, previous work <cit.> demonstrates that many existing approaches can be viewed as Block Coordinate Descent (BCD) problems. However, a primary challenge with BCD lies in its propensity to necessitate full minimization of the sub-problems, which proves to be challenging in the Poisson case. In the L2 case, the subproblems often have closed-form solutions. To address this concern, we explore BSUM algorithms <cit.> in this study, a generalization of BCD that avoids the requirement for full minimization of the subproblems, instead using upper bounds for the objective function.



  
Non-Negative Matrix Factorization
Non-Negative Matrix Factorization (NMF) algorithms have been extensively studied, and for a comprehensive review, one can refer to <cit.>. Initially formulated as Positive Matrix Factorization for the Gaussian (L2-NMF) case by <cit.>, NMF has seen numerous algorithmic developments. In the L2 case, popular approaches include the Alternating Nonnegative Least Squares (ANLS) framework <cit.> and the Hierarchical Alternating Least Squares (HALS) method <cit.>. As for the Poisson case, known as KL NMF, the first algorithm using Multiplicative Updates (MU) was proposed by <cit.>, with later demonstrations of its convergence provided in <cit.> for both Poisson and Gaussian cases. A more rigorous convergence analysis is presented in <cit.>.

Considering the specific problem of Poisson KL-NMF, there have been a few notable contributions. For instance, <cit.> employed the Alternating Direction Method of Multipliers (ADMM) with the variable change X=WH. The Primal Dual algorithm, based on the framework from Chambolle-Pock, was explored by <cit.>. Moreover, <cit.> conducted a comparative study of various optimization algorithms for KL NMF, including MU <cit.>, ADMM <cit.>, Primal Dual <cit.>, and Cyclic Coordinate Descent Method <cit.>. They also introduced three new algorithms for KL-NMF, namely Block Mirror Descent Method, A Scalar Newton-Type Algorithm, and A Hybrid SN-MU Algorithm.

In the introduction, we mentioned that there are few contributions that address the problem of regularized NMF, with many of them focusing on the L2 case. <cit.> demonstrated how many existing works can be cast as Block Coordinate Descent (BCD) problems, allowing the derivation of MU update rules for different regularizers, such as L1 for sparsity. However, their work is limited to L2 NMF. Xu et al. <cit.> proposed a general optimization scheme for block multiconvex optimization using block coordinate descent, which can accommodate regularization on each block. Although applied to L2-NMF, such an approach may lead to algorithms with sub-iterations.
In the context of L2-NMF with sparsity constraints, <cit.> presented an approach to address this scenario. Additionally, <cit.> introduced a framework for handling L2-NMF with Lipschitz regularizers, akin to our term s_L. Other forms of regularization have also been explored, such as graph-based <cit.> or simplex constraint <cit.>.

Regarding the specific problem of regularized KL loss, <cit.> provided a notable contribution. However, their work focused solely on the subproblems of the NMF problem, rather than addressing the NMF problem itself. Notably, one could potentially employ a similar approach to solve the subproblems of (<ref>) using BCD. Nonetheless, this would result in a less efficient algorithm with sub-iterations.




§ PRELIMINARIES




 §.§ Notation




We reserve capital letters for matrices and vectors, e.g., A,a. We use a_i to refer to the i-th numbered vector, and a_i to denote the i-th element of vector a. 
The j-th element of vector a_i or the element at the i-th row and j-th column of matrix A is denoted as a_ij.

A≥0 and a≥0 indicate that all entries of matrix A or vector a are greater than or equal to 0, i.e., a_ij≥0 for all i,j. 
A^⊤,a^⊤ represent the transpose of A and a, respectively. 
We use A^t to denote A at step t. 
A^t^⊤ denotes the transpose of A^t. 
∘ and ⊘ denote elementwise multiplication (also known as the Hadamard product) and division for matrices, respectively. 
For example, [A∘B]_ij=a_ijb_ij and [A⊘B]_ij=a_ij/b_ij.

As mentioned in the introduction, the matrix to be factorized is generally denoted as ≈, where  and  are its factors. 
We will use w,h as the vectorized versions of ,, while w_i,h_i will denote the i-th column of ,. 
x,x are general variables that can replace either w or h. 
We use ℒ(,) as the general loss function, and f(x) is broadly used to denote a multivariate scalar function.

Finally, we commonly employ calligraphic notation for variable domains. Let 𝒳 serve as a generic domain for the variable x. In practical terms, it is frequently defined by the ≥ϵ constraint, specifically as 𝒳={x∈^m|x≥ϵ}. Additionally, we utilize C_w and C_h to represent the domains of w and h, respectively.



 §.§ Definitions



In this contribution, we consider the loss function ℒ(w,h) with two blocks of variables: w∈𝒞_w and h∈𝒞_h, where 𝒞_w⊂^m_w and 𝒞_h⊂^m_h are both non-empty convex sets. Here, w and h correspond to the vectorized versions of the two matrices  and , respectively. Therefore, 𝒞_w and 𝒞_h often correspond to the sets w≥ϵ and h≥ϵ with ϵ>0. Let us use z=[w,h] to denote all the variables. We have z∈𝒞=𝒞_w×𝒞_h⊂^m=^m_w×^m_h, where the total dimension of the problem is m=m_w+m_h.

Let ℒ:𝒞→ be a scalar function, where 𝒞⊂^m is a convex set. The directional derivative of ℒ at point x in the direction d is defined by

    ℒ^'(z;d)lim_λ↓0ℒ(z+λd)-ℒ(z)/λ


Note that when ℒ is differentiable, ℒ'(z;d) = ∇ℒ(z)d^⊤ since d and z are row vectors. In this contribution, almost all functions are differentiable on the domain of interest.

The point z = [w,h]∈𝒞 is a coordinatewise minimum of a function ℒ if 

    ℒ(w+d_w,h)≥ℒ(w,h)∀d_w∈ℝ^m_wwithw+d_w∈𝒞_w


    ℒ(w,h+d_h)≥ℒ(w,h)∀d_h∈ℝ^m_hwithh+d_h∈𝒞_h


A coordinatewise minimum is a natural termination point for an alternating minimization algorithm. However, it is important to note that a coordinatewise minimum is not equivalent to a local minimum, as it does not guarantee minimality in all directions. Figure <ref> (left) provides a counterexample illustrating this.

Another significant concept is the notion of a stationary point, where the gradient is non-negative in all directions.

Let ℒ:𝒞→ be a scalar function, where 𝒞⊂^m is a convex set. A point x is a stationary point of ℒ if

    ℒ^'(z;d)≥0∀d|z+d∈𝒳



We emphasize that a stationary point is not equivalent to a strict local minimum as there might be directions where the directional derivative equals 0.
For example, in the simple function f([w,h]) = (wh-2)^2, the point [w,h] = [√(2),√(2)] has a zero derivative in the direction [1,-1], which corresponds to rescaling the solution as [α,1/α].
Even worse, a stationary point is not necessarily a local minimum, even if it is a coordinatewise minimum, as shown in Figure <ref> (a). Here, in the diagonal directions, the directional derivative equals 0 but the function is concave in this direction. 


When it comes to the Poisson loss function, there is, at least, a continuous set of local minimas corresponding to rescaling the solution.
This is illustrated in Figure <ref> (b). Note that the introduction of regularization or constraints can lead to strict local minima, as shown in Figure <ref> (c).

In this contribution, we prove convergence to a coordinatewise minimum that is also a stationary point. To accomplish this, we will consider a class of functions that are regular at their coordinatewise minima.


The function ℒ:𝒞→ is said to be regular at the point z∈𝒞 if ℒ'(z;d) ≥ 0 for all d = [d_w,d_h]∈^m such that ℒ'(z;[d_w,0]) ≥ 0 and ℒ'(z;[0,d_h]) ≥ 0.



lemmacontiuousdifferentiableregular

Continuously differentiable functions are regular at their coordinatewise minimums.


Proof is provided in Appendix <ref>. 
Lemma <ref> plays a crucial role, as it ensures that the coordinate-wise minimum we converge to in Theorem <ref> is also a stationary point. In this work, we assume the regularizer to be continuously differentiable on the domain.






 §.§ Approximation functions

In order to facilitate optimization algorithms, it is beneficial to work with approximation functions that majorize or approximate the objective function at a given point. One commonly used class of approximation functions is known as first-order majorization functions. These functions provide a convenient framework for constructing surrogates and facilitating optimization. We adopt the definition of first-order majorization functions from <cit.>.

<cit.>
A function g(x,x^t) is said to be a first-order majorization of f at the point x^t if it satisfies the following properties:

A.1      
g(x,x^t) ≥f(x)  ∀x,x^t∈𝒳,

A.2    
g(x^t,x^t) = f(x^t)  ∀x^t∈𝒳,

A.3     
g^'(x,x^t;d)|_x=x^t = f^'(x^t;d)  ∀d  such that  x^t+d∈𝒳,

A.4     
g(x,x^t)  is continuous in  (x,x^t). 


It is worth noting that for continuously differentiable functions, the third statement can be equivalently expressed as ∇_xg(x^t,x^t)=∇ f(x^t). Although the definition of first-order majorization functions resembles the concept of surrogate functions introduced in <cit.>, the additional requirement for a surrogate function is that g(x,x^t)-f(x) is L gradient Lipschitz as defined in (<ref>). 

Importantly, all majorization functions defined in the following Section <ref> satisfy this condition and can thus serve as majorization functions.

Conveniently, majorization functions can be built term by term, leveraging their additivity property. This property allows us to combine multiple majorization functions to obtain a new majorization function. 


First-order majorization functions are additive. If g_1(x,x^t) and g_2(x,x^t) majorize f_1(x) and f_2(x) at x^t, respectively, then g_1(x,x^t)+g_2(x,x^t) majorizes f(x)=f_1(x)+f_2(x) at x^t.



The additivity property preserves each property of (<ref>).

Lemma <ref> provides a valuable tool for constructing majorization functions by combining simpler majorization functions. Additionally, when proving that a function is majorizing, it is often unnecessary to explicitly demonstrate the equality of partial derivatives or gradients at x^t. Instead, in the case of differentiable functions, it is typically sufficient to establish the first two properties (A.1 and A.2). According to <cit.>, properties A.3 and A.4 follow as a consequence. Intuitively, one can observe that the continuity of the gradient ensures that the majorization function g shares the tangent spaces with f at the point x^t.










    


 §.§ Two Blocks Successive Minimization (TBSUM)

















The TBSUM algorithm is designed to solve the following problem:

    _h,w    ℒ(w,h) 
    such that     h∈𝒞_h,w∈𝒞_w

It relies on two first-order majorizing functions: g_w(w,w^t,h^t) and g_h(h,h^t,w^t), which majorize ℒ(w,h) at (w^t,h^t) for all w^t∈𝒞_w and h^t∈𝒞_h. The construction of these functions will be presented in Section <ref>. The TBSUM algorithm, outlined in Algorithm <ref>, alternates between minimizing g_w and g_h. It is assumed that the subproblem solutions are unique.
Theorem <ref> establishes the convergence of the TBSUM algorithm, which is a variant of the algorithm presented in <cit.> adapted for solving the specific problem at hand.


Given two quasi-convex first order majorizing functions g_h(w,w^t,h^t)
and g_w(h,h^t,w^t+1) of ℒ(w,h)
at (w^t,h^t),∀w^t,h^t∈𝒞_w×𝒞_h.
Furthermore assuming that the two subproblems in the TBSUM Algorithm
<ref> have unique solutions for any points w^t∈ C_w,
h^t∈ C_h. Then, every limit point z=[w,h]
of the iterates generated by the TBSUM Algorithm <ref>
is a coordinatewise minimum of (<ref>). In addition,
if ℒ is regular at any point z∈𝒞,
then z is a stationary point of (<ref>).






§ SUBPROBLEM MINIMIZATION

In this section, we focus on constructing the appropriate majorization functions g_h(w,w^t,h^t) and g_w(h,h^t,w^t+1) for our problem (<ref>). Since we consider the same type of regularization for w and h, both subfunctions have the same form.

Practically, the loss function can be rewritten as

    ℒ_(,)+R(,)=∑_jf_w(w_j)=∑_if_h(h_i)

where w_j and h_i are the i^th row and j^th column of  and . The functions f_w and f_h have the form:

    f(x)=-∑_i=1^m(b_ilog(a_i^⊤x)+a_i^⊤x)+s_L(x)+s_R(x)+∑_j=1^ns_C(x_j).

Therefore, in this section, our objective is to find majorization functions for (<ref>). Once this is done, we will provide closed-form solutions for steps 1 and 2 of Algorithm <ref>. It is worth noting that each term of (<ref>) can be handled separately using the additivity property of majorizing functions (Lemma <ref>).



 §.§ Majorizing functions


The following four lemmas provide majorizing functions for the different term of our objective function. Proofs are provided in Appendix <ref>.

In order to develop an efficient algorithm, our objective is to identify majorizing functions that result in sub-problems with closed-form tractable solutions. Often, this can be accomplished under two conditions: 1. all the majorizing functions are of the same form, and, 2. the majorization function is separable with respect to the variables x, i.e., g(x) = ∑_i g_i (x_i). Within the scope of this contribution, we consider two forms of majorizing functions: quadratic g(x) = a + bx + cx^2 and logarithmic g(x) = a + bx - c log(x). 




First, we propose a majorization scheme for the logarithmic term log(a_i^⊤x) in the objective function (<ref>). We utilize a widely used majorization technique based on the concavity of the logarithm function. This technique has been employed in the original work by Lee and Seung <cit.> as well as in many EM (Expectation-Maximization) schemes.
[Log majorization]lemmalogmajorization
     Assuming a∘x>0, for x∈𝒞,
    let us define q_j=a_jx_j^t/∑_ka_kx_k^t
    for x^t∈𝒞, then g(x,x^t)=-∑_jq_jlog(a_jx_j/q_j)
    is a first order majorizing function of f(x)=-log(a^⊤x)=-log(∑_ja_jx_j).


We now proceed to majorize the different terms of the regularisation function s_L(x), s_R(x), and s_C(x_j). We can majorize any Lipschitz function using the following lemma.
[Lipschitz-majorization]lemmalipschitzmajorization
Given s_L(x) a gradient Lipschitz function with
constant σ_L over the domain x∈𝒞. The
functions

    g_1(x,x^t)     =s_L(x^t)+(x-x^t)^⊤∇ s_L(x^t)+σ_Lx-x^t_2^2

and 

    g_2(x,x^t)=s_L(x^t)+(x-x^t)^⊤∇ s_L(x^t)+2σ_L(max_jx_j^t)(∑_jx_j^tlog(x_j^t/x_j)-x_j^t+x_j)

are first oder majorizing functions at x^t∈𝒞.

We note that (<ref>) (quadratic majorisation) is tighter than (<ref>) (logarithmic majorisation) . However, the looser majorisation function is needed to obtain a close form solution for the MU (see Section <ref>). 

Next, the term that is relatively smooth can be majorized using the following lemma.
[Relative smoothness majorization]lemmarelativesmoothnessmajorization

Assuming s_R(x)
a σ_R relatively smooth function with respect to κ(x)=-1^⊤log(x)
for x∈𝒞⊂ℝ_+^n. Then the function

    g(x,x^t)=s_R(x^t)+⟨∇ s_R(x^t),x-x^t⟩ +σ_R∑_i^n(x_i/x_i^t-log(x_i/x_i^t)-1)

is a first order majorizing function of s_R(x)
for x^t∈𝒞.


[Concave majorisation]lemmaconcavemajorisation
Given s(x) a concave function defined on x∈𝒞⊂ℝ,
it's linear approximation at the point x^t

    g(x,x^t)=s(x^t)+∂ s(x^t)/∂ x(x-x^t)

is a first order majorization function for x^t∈𝒞.





 §.§ Subproblem updates

Now that we have defined majorizing functions for each term of (<ref>), we can apply the additivity property of Lemma <ref> to obtain a general majorizing function for f(x):




    g(x,x^t)     =-∑_i=1^m(b_i∑_jq_ijlog(a_ijx_j/q_ij)+a_i^⊤x)
       +s_L(x^t)+(x-x^t)^⊤∇ s_L(x^t)+2σ_L(max_jx_j^t)(∑_jx_j^tlog(x_j^t/x_j)-x_j^t+x_j)
       +s_R(x^t)+⟨∇ s_R(x^t),x-x^t⟩ +σ_R∑_j(x_j/x_j^t-log(x_j/x_j^t)-1)
       +∑_j=1^ns_C(x_j^t)+∂ s_C(x_j^t)/∂ x_j(x_j-x_j^t)

where q_ij=a_ijx_j^t/∑_ka_ikx_k^t. 
We use the colors green, blue, orange, and purple to denote and keep track of the dependencies of the different terms in (<ref>).
Finding the local optimum of the majorizing function will provide us with an update for Algorithm <ref>.
[Generalized MU for (<ref>)]propositionmumarjorization

    Assuming x^t, x, b, A > 0, the first-order majorizing function defined in (<ref>) is strictly convex, and its global minimum x^t+1 is given by
    
    x_j^t+1 = x_j^tα_j^t/β_j^t

    where
   
    α_j^t    =∑_ib_ia_ij/∑_ka_ikx_k^t+2(max_ix_i^t)σ_L+σ_R/x_j^tand
    β_j^t    =∑_ia_ij+∇_x_js_L(x^t)+2(max_ix_i^t)σ_L+∇_x_js_R(x^t)+σ_R/x_j^t+∂ s_C(x_j^t)/∂ x.

    

The proof is provided in Appendix <ref>.



  
Generalization of the traditional MU Rule
We note that (<ref>) serves as a generalization of the original Multiplicative Update (MU) rule presented in <cit.>. Removing the regularization terms (blue, orange, and purple) results in precisely the MU rule as outlined in <cit.>.



  
Connection with (Block) Mirror Descent <cit.>
Another interesting observation is that the majorization of the relatively smooth term is done similarly to a Bregman proximal method algorithm <cit.>. Since the objective function -∑_i=1^m(b_ilog(a_i^⊤x)+a_i^⊤x) is relatively smooth, one could drop all terms except for s_R and optimize using Block Bregman Proximal Gradient (BBPG) <cit.>. This would result in an algorithm very similar to Block Mirror Descent (BMD), which has recently been proposed for solving Poisson NMF <cit.>. Nevertheless, we advice against this this solution as discussed further in Section <ref>.



  
Alternative majorizing function and Quadratic Update (QU)
As shown experimentally in Section <ref> and illustrated in Figure <ref>, having a majorization function as tight as possible leads to faster convergence of the algorithm. 
In (<ref>), we deliberately choose to use a looser majorizing function for the term s_L in order to recover an algorithm with multiplicative update that generalizes the original approach from <cit.>. However, instead of using (<ref>), one can also use (<ref>) when constructing the majorizing function:

    g(x,x^t)     =-∑_i=1^m(b_i∑_jq_ijlog(a_ijx_j/q_ij)+a_i^⊤x)
       +s_L(x^t)+(x-x^t)^⊤∇ s_L(x^t)+σ_L‖x-x^t‖ _2^2
       +s_R(x^t)+⟨∇ s_R(x^t),x-x^t⟩ +σ_R∑_j(x_j/x_j^t-log(x_j/x_j^t)-1)
       +∑_j=1^ns_C(x_j^t)+∂ s_C(x_j^t)/∂ x_j(x_j-x_j^t)

which is also a strictly convex function. 
[QU for (<ref>)]propositionqumarjorization

Assuming x^t, x, b, A > 0, the first-order majorizing function defined in Equation (<ref>) is strictly convex, and its global minimum x^t+1 is given by


    x_j^t+1=-β_j^t+√((β_j^t)^2+4αζ_j^t)/2α

where

    α=2σ_Lβ_j^t=∑_ia_ij+∇_x_js_L(x^t)-2σ_Lx_j^t+∇_x_js_R(x^t)+σ_R/x_j^t+∂ s_C(x_j^t)/∂ xζ_j^t=∑_ib_ia_ijx_j^t/∑_ka_ikx_k^t+σ_R.




The proof is provided in Appendix <ref>. Both of these propositions lead to the update rule for our MU and QU algorithms detailed in Section <ref>. We note also that, with the appropriate assumptions, the update rule <ref> and <ref> will preserve positivity of the variable x. However, since our desire is also to handle extra constraint, we develop in the next section rigorous approach.



 §.§ Generalized simplex constraint 







We need to handle two constraints: 1. the linear constraint x≥ϵ, and, 2. the scale constraint e^⊤x=1, where e≥0. 
While the first one is used to keep the variable non-negative, typically with a strictly positive small ϵ, the second one can set the scale of one of the variables ( or ) in the factorization problem. Furthermore, in the case e=1, the simplex constraint is recovered. 
It turns out that the update rules of (<ref>) and (<ref>) can simply be updated to handle this constraint. 
The actual optimization problem we want to solve becomes:

    _x≥ϵf(x)such thate^⊤x=1.

where f is given in (<ref>). 

To solve this problem, we used the KKT approach, i.e, we find points that satisfy the KKT (Karush-Kuhn-Tucker) conditions:

    1. Stationarity     
∇_xL(ẋ,ν,μ)=0,     

    2. Primal feasibility 
    
    e^⊤ẋ-1=0, 

    ẋ≥ϵ
,     

    3. Dual feasibility 
    μ≥0,     

    4. Complementary slackness 
    μ^⊤(-x+ϵ1)= 0 ,     













where the Lagrangian is defined as:

    L(x,ν,μ)=f(x)+ν(e^⊤x-1)+μ^⊤(-x+ϵ1).



We follow the same method as developed in Section <ref>, except that we majorize the Lagrangian L(x,ν,μ). The resulting first-order majorizing function is given by:

    g^'(x,x^t,ν,μ)=g(x,x^t)+ν(e^⊤x-1)+μ^⊤(-x+ϵ1)

where g(x,x^t) is given in (<ref>) or (<ref>). We repeat the development of Section <ref> (and the proofs of Appendix <ref>). We end up with and update that is very similar to (<ref>) or (<ref>). In the MU case, we end up with:

    x_j^t+1=x_j^tα_j/β_j+ν e_j-μ_j,

where the only final difference consists of two terms in cyan and violet (α_j and β_j remain identical). For the QU, we stick to the same update rule (<ref>), where only β_j^t is modified:

    β_j^t'=β_j^t+ν e_j-μ_j.

This update rule ensures the first of the KKT conditions (stationarity). We now find ν, μ such that the second KKT condition holds (primal feasibility). It turns out that μ does not need to be computed explicitly. In the MU case, μ_j is selected to be large enough such that


    x_j^t+1    =max(x_j^tα_j/β_j+ν e_j,ϵ)=x_j^tα_j/min(x_j^tα_j/ϵ,β_j+ν e_j).

In the QU case, we obtain

    x_j^t+1    =max(-(β_j^t+ν e_j)+√((β_j^t+ν e_j)^2+4αζ_j^t)/2α,ϵ)
        =-min(ζ_j^t/ϵ-ϵα,β_j^t+ν e_j)+√((min(ζ_j^t/ϵ-ϵα,β_j^t+ν e_j))^2+4αζ_j^t)/2α


Note that dual feasibility and complementary slackness could be verified, but we leave them out for simplicity. We then need to find the value of ν such that e^⊤x=1, which is equivalent to searching for

    h_1(ν)     =∑_je_jx_j^tα_j/min(x_j^tα_j/ϵ,β_j^t+ν e_j)-1=0

Similarly, for the quadratic update of (<ref>), we search for ν that satisfies

    h_2(ν)=∑_je_j-min(ζ_j^t/ϵ-ϵα,β_j^t+ν e_j)+√((min(ζ_j^t/ϵ-ϵα,β_j^t+ν e_j))^2+4αζ_j^t)/2α-1=0.

There is no closed-form solution for ν; however, the value can be found using a simple dichotomy search. Bounds for starting the dichotomy are computed in Appendix <ref>.




  
Case ϵ=0
Most of our reasoning relies on the fact that x_i>0 and, therefore, on the domain constraint ϵ>0. 
We have found that setting ϵ to a small non-zero value works well in practice. 
However, our approach can likely be generalized to the case where ϵ=0, following the approach of <cit.>, which studies the unregularized Poisson NMF case.




§ ALGORITHMS FOR POISSON MATRIX FACTORISATION


Equipped with the update rules developed in the previous Sections <ref> and <ref>, we are ready to tackle the general problem of this contribution[Here we show the problem with the linear constraint on , however, by symmetry, a similar algorithm can be developed with the constraint on .] which consists of minimizing (<ref>):

    ,=_,    -⟨,log()⟩ +⟨1,⟩ +R_W()+R_H()
    such that   ≥ϵ,≥ϵ,e^⊤=1

We first observe that with respect to each variable ,, the problem is separable by column/row. 
For example, given , finding the optimal  can be done for each column:

    ḣ_̇i̇=_h_i-⟨y_i,log(h_i)⟩ +⟨1,h_i⟩ +r_H(h_i)such thath_i≥ϵ,e^⊤h_i=1

We therefore apply the TBSUM Algorithm <ref>, where all lines of  and all columns of  are updated independently, and obtain the two Algorithms <ref> and <ref>. We note here that the function max(·, ϵ) ensures that the solutions are ≥ϵ (non-negativity).









  
Convergence
The two update rules in steps 5 and 10 correspond to minimizing first-order strictly convex majorization functions. As a result, we can apply Theorem <ref> to guarantee convergence towards a coordinate-wise minimum. 
It is important to note that this coordinate-wise minimum is also a stationary point, given that the objective function remains regular for any point ∈𝒞_W,∈𝒞_H.



 §.§ Algorithm complexity

Let's examine the complexity of both Algorithms <ref> and <ref> when considering ∈ℝ^n× k and ∈ℝ^k× m. In each iteration, the following complexities are observed: 

(a) Step 4 has a complexity of 𝒪(nmk). 

(b) Step 5 has a complexity of 𝒪(nk).

(c) Step 8 has a complexity of 𝒪(nmk).

(d) Step 9 has a complexity of 𝒪(c_dkm), where c_d denotes the number of iterations performed by the dichotomy.

(e) Step 10 has a complexity of 𝒪(mk).

Thus, the overall complexity per iteration can be expressed as 𝒪(nmk)+𝒪(c_dkm)=𝒪((n+c_d)mk). This indicates that the computational complexity per iteration is linear with respect to the problem size, i.e., nm, multiplied by the number of components, i.e., k. 


  
Impact of the dichotomy
When n is small, the computational cost of the dichotomy in step 4 becomes dominant. Nevertheless, in general, for larger values of n, the impact of the dichotomy becomes negligible.



 §.§ Tight Majorizing Functions 


While we do not make any theoretical contributions concerning the speed of convergence of the algorithm, we want to emphasize the natural fact that tighter majorizing functions lead to faster convergence. Therefore, when evaluating an algorithm, we believe that the analysis of the underlying majorizing function is as insightful as the experimental evaluation. As an example, we could have used Block Mirror Descent to solve (<ref>), as was done in <cit.>. This algorithm uses a Bregman Difference to create a majorization function for the subproblem. However, this would result in a much looser majorization function, which partly explains the slow convergence of this algorithm observed in <cit.>. This difference between majorization functions is exemplified in Figure <ref>.






  
Linesearch
By tightening the bounds we used to construct the surrogate function, we can develop a more efficient algorithm. Here, we apply a classic "linesearch" method to the functions s_L. However, the same technique can be trivially applied to s_R as well.
First, in (<ref>) or (<ref>), replace the constant σ_L with a parameter γ and initialize it with σ_L.
Second, at each iteration, update the parameter γ according to the following rule:

    γ^t+1=υγ^t   ifs_L(x)≥ g(x,x^t,γ)
    1/τγ^t   otherwise.

Here, υ and τ are two update rates that determine how fast γ is updated. Choosing values that are too small for these parameters leads to an inefficient linesearch, while selecting values that are too large can result in strong oscillation patterns. Typical values for υ and τ range from 1.05 to 1.5.
However, it is important to note that when using linesearch, we are not guaranteed to converge, as we might invalidate the assumptions of Theorem <ref>.



§ NUMERICAL SIMULATION 





  
Problem
In this section, we analyze the speed of convergence of Algorithms <ref> and <ref> through numerical simulations. As a regularizer for , we consider the Laplacian regularization R(,λ) = λ/2tr(^TΔ), where Δ represents the two-dimensional Laplacian for the kth line of ∈ℝ^k × p^2 reshaped as k images of size p × p. Since a straightforward approach to minimize R(,λ) is to reduce the amplitude of , we add the simplex constraint 1^T = 1. This leads to the following optimization problem:

    ,  =    min_,  -⟨, log() ⟩ + ⟨1, ⟩ + λ/2tr(^TΔ)    subject to ≥ϵ, ≥ϵ, 1^T = 1

This particular problem can be applied in various domains, such as Non-Negative Matrix Factorization for hyperspectral images <cit.> and remote sensing <cit.> (See Related Work Section <ref> for more references and applications).
Our algorithms and regularisations were specifically developed for the espm python package <cit.>. All algorithms and experiments can be found in the espm package.



  
Dataset
We construct two datasets consisting of 50 randomly drawn samples. In the first dataset, both matrices  and  are randomly generated from a uniform distribution. In the second dataset, each column of  corresponds to the sum of Gaussian functions that are randomly centered and scaled. The matrix  represents random smooth images. This second dataset is created using the espm package <cit.>, where the toy model is used for  and  is generated using the "laplacian" weight type.

This choice of dataset is selected because it can benefit from the Laplacian regularization on .

Once  and  are generated, the noiseless matrix  is obtained as =. We introduce noise by independently sampling each element ij∼Poisson(λij)/λ, where λ can be regarded as the noise control parameter. For all samples, we set ∈ℝ^n × k and ∈ℝ^k × p^2 with k = 3, p = 64, and n selected from the set 25, 100, 500, 1000. Thus, the images in the dataset have dimensions of 64 × 64.



  
Results
We compare the performance of Algorithm <ref> (MU), Algorithm <ref> (QU), Block Mirror Descent (similar to <cit.>) and the projected gradient algorithm applied to (<ref>). 
Figure <ref> displays the convergence curves for 1000 iterations and n=25. 
Although the overall complexity of all algorithms is the same, the time per iteration differs due to the different operations performed within each iteration and the time spent on dichotomy to compute the dual variable ν. 
Therefore, we provide the time in seconds for each algorithm in Figure <ref> for various value of n. All results are averaged over 50 repetitions.







  
Discussion
Let's discuss the results in more detail: 

 - Number of iterations: Figure <ref> illustrates the convergence behavior of QU and MU algorithms. It is evident that QU converges faster per iteration compared to MU, which aligns with our expectations due to the tighter majorizing function used in QU. However, it is important to note that the introduction of the linesearch technique, while accelerating convergence, can lead to occasional instability, as indicated by occasional increases in the loss function. This observation supports our earlier discussion in Section <ref>. 
The challenge with Projected Gradient Descent is that we need to find an initial learning rate that is not too large, as the algorithm can diverge and not too small as the algorithm can be slow. Overall, since the selected learning rate cannot be selected optimally, the algorithm is slower than QU and MU.
The Block Mirror Descent algorithm is also slower than QU and MU, which is consistent with the results of <cit.> and can be explained by the fact that the majorizing function used in the algorithm is looser.


 - Time per iteration: Figure <ref> presents the total time taken by the algorithms to complete 100 iterations. The results demonstrate that for small values of n, the computation of the p^2 dual variable during the dichotomy process dominates the overall execution time. However, as n increases, the time spent on dichotomy becomes negligible in comparison. These findings align with the complexity per iteration discussed in Section <ref>.





§ CONCLUSION

This contribution is the first to address the Poisson NMF problem with general regularization terms, such as Lipschitz functions, relatively smooth functions, or those expressed as linear constraints. We introduce two new algorithms and demonstrate their convergence to a coordinate-wise minimum, which is also a stationary point. Emphasizing the impact of the majorizing function choice on convergence speed, we validate our findings through numerical simulations. In essence, we believe that this work serves as a helpful guide for developing efficient algorithms suited for regularized Poisson NMF problems.





§ PROOFS

In this Appendix, we provide the different proofs used in the paper.



 §.§ Proof of Lemma <ref>


We note that this lemma and its proof likely exist in the literature, but we were unable to find a reference.



*

If ℒ is continuously differentiable, the directional derivative can be written as ℒ'(z;d)=d^⊤∇ℒ(z). At a coordinatewise minimum z, we have by definition:

    ℒ^'(z;[d_w,0])=∇ℒ(z)[d_w,0]^⊤≥0

and

    ℒ^'(z;[0,d_h])=∇ℒ(z)[0,d_h]^⊤≥0.

Therefore,

    ℒ'(z;d)    = ∇ℒ(z)[d_w,d_h]^⊤
       = ∇ℒ(z)([d_w,0]^⊤ + [0,d_h]^⊤) 
       = ∇ℒ(z)[d_w,0]^⊤ + ∇ℒ(z)[0,d_h]^⊤
       ≥ 0.






 §.§ Proof of majorizing functions



*


First let us observe that 

    g(x^t,x^t) 
        =-∑_ja_jx_j^t/∑_ka_kx_k^tlog(a_jx_j^t/a_jx_j^t/∑_ka_kx_k^t) 
        =-∑_ja_jx_j^t/∑_ka_kx_k^tlog(∑_ka_kx_k^t) 
        =-log(∑_ka_kx_k^t)=f(x^t).

Therefore, we have g(x,x^t)≥ f(x). The inequality follows from the convexity of the -log function:

    -log(∑_ja_jx_j)     =-log(∑_jq_ja_jx_j/u_j)≤-∑_jq_jlog(a_jx_j/u_j)

where we set q_j=a_jx_j^t/∑_ka_kx_k^t.
Finally, by continuity, we obtain the 3^rd and 4^th properties of majorizing functions.


We now proceed to majorize s_L(x), s_R(x), and s_C(x_j).
*

First it can be trivially observed that 

    s_L(x^t)=g_1(x^t,x^t)=g_2(x^t,x^t),

which satisfies the first property. We then take 1^st
order Taylor expension of s_L around x^t and find 


    s_L(x)=s_L(x^t)+(x-x^t)^⊤∇ s_L(x^t)+

where ≤σ_Lx-x^t_2^2
since the function s_L is gradient Lipschitz with constant σ_L.
Therefor we have

    s_L(x)    ≤ s_L(x^t)+(x-x^t)^⊤∇ s_L(x^t)+σ_Lx-x^t_2^2=g_1(x,x^t) 
       ≤ s_L(x^t)+(x-x^t)^⊤∇ s_L(x^t)+2σ_L(max_jx_j^t)(∑_jx_j^tlog(x_j^t/x_j)-x_j^t+x_j) 
        =g_2(x,x^t),

where (<ref>) will be shown later in this proof.
By continuity, we obtain the 3^rd and 4^th
property of majorizing functions. We now need to prove (<ref>) and reformulate it as 




    x-x^t_2^2≤2(max_jx_j^t)(∑_jx_j^tlog(x_j^t/x_j)-x_j^t+x_j)=2(max_jx_j^t)D_GKL(xx^t).

For simplicity, let us define the function

    q(x)=∑_jx_jlog(x_j),

with the gradient ∇_x_jq(x)=log(x_j)+1
and the Hessian 

    H_ij^q(x)=∂ q/∂ x_i∂ x_j(x)=1/x_j   if i=j
    
    0    otherwise.

Note that q is a strictly convex function for x>0.
We expand the generalized KL divergence:

    D_GKL(xx^t)     =∑_jx_j^tlog(x_j^t)-∑_jx_j^tlog(x_j)-∑_jx_j^t+∑_jx_j
        =∑_jx_j^tlog(x_j^t)-∑_jx_jlog(x_j)-∑_j(log(x_j)+1)(x_j^t-x_j)
        =q(x^t)-(q(x)+∇ q(x)^⊤(x^t-x))
        =1/2(x^t-x)^TH^q(x̃)(x^t-x),

where x̃ is selected such that the last equality holds.
Since q is a strictly convex function, we know that x̃=ρx+(1-ρ)x^t
for some give ρ∈[0,1]. Now we bound the Hessian as

    H^q(x)≥1/max_jx_jI

and introducing this inquality in (<ref>), we obtain

    D_GKL(xx^t)≥1/2max_jx_jx-x^t_2^2,

which is equivalent to (<ref>) and completes the proof.


*

The first property g(x^t,x^t)=s_R(x^t)
can be trivially verified. Then using by the definition of relatively
smoot function

    s_R(x)≤ s_R(x^t)+⟨∇ s_R(x^t),x-x^t⟩ +σ_Rℬ_κ(x,x^t)

where

    ℬ_κ(x,x^t):=κ(x)-κ(x^t)-⟨∇κ(x^t),x-x^t⟩ .

Given that κ(x)=-1^⊤log(x),
we compute

    ℬ_κ(x,x^t)=∑_i^n(x_i/x_i^t-log(x_i/x_i^t)-1)

Finally, by continuity, we obtain the 3^rd and 4^th
property of majorizing functions.


*

One can simply observe s(x^t)=g(x^t,x^t).
Then by concavity, we have

    s(x_i)≤ s(x_i^t)+∂ s(x_i^t)/∂ x_i(x_i-x_i^t)

Finally, by continuity, we obtain the 3^rd and 4^th
property of majorizing functions.




 §.§ Proof of subproblem updates


In this subsection, we present the proof of the subproblem updates used in  the MU and QU algorithms. Let us start with the MU updates. 
*

    Assuming a_ij, b_i > 0, (<ref>) is strictly convex, as the green term is strictly convex, and the remaining terms are convex. Consequently, (<ref>) possesses a global minimum. To identify this minimum, we seek the stationary point ∇_xg=0. Due to our meticulous selection of majorizing functions, this subproblem becomes separable. When computing the gradient with respect to the variable x_j, we obtain:
    
    
    ∇_x_jg(x,x^t)     =-1/x_j∑_ib_ia_ijx_j^t/∑_ka_ikx_k^t+∑_ia_ij
       +∇_x_js_L(x^t)+2(max_ix_i^t)σ_L(1-x_j^t/x_j)
       +∇_x_js_R(x^t)+σ_R(1/x_j^t-1/x_j)
       +∂ s_C(x_j^t)/∂ x=0

    where we assume that x_j, x_j^t > 0 since 0 ∉𝒞. Transforming the above expression, we find a multiplicative update rule (<ref>) for x_j.
    
    
The proof of the QU updates is similar to the MU updates, expect that we need to solve a quadratic equation to obtain a closed-form solution.
*


    Assuming a_ij, b_i > 0, Equation (<ref>) is strictly convex. This convexity arises from the strict convexity of the green term, coupled with the convexity of the other terms. Consequently, (<ref>) possesses a global minimum. To identify this minimum, we seek the stationary point by computing ∇_xg=0:
    
    ∇_x_jg(x,x^t)     =-1/x_j∑_ib_ia_ijx_j^t/∑_ka_ikx_k^t+∑_ia_ij
       +∇_x_js_L(x^t)+2σ_L(x_j-x_j^t)
       +∇_x_js_R(x^t)+σ_R(1/x_j^t-1/x_j)
       +∂ s_C(x_j^t)/∂ x=0

    We observe that it is a separable quadratic function, hence the update rule named Quadratic Update (QU). 
    Solving (<ref>) for x_j > 0 can be rewritten as
    
    α x_j^2 + β_j^tx_j - ζ_j^t = 0,

    where α, β_j^t, and ζ_j^t are given in (<ref>). Assuming ζ_j^t≠ 0, we have ζ_j^t > 0 and 4αζ_j^t > 0. 
    Therefore, the previous quadratic equation has two real solutions. 
    Since √((β_j^t)^2 + 4αζ_j^t) > β_j^t, they are of opposite sign. Due to the constraint x_j≥ϵ > 0, we select the positive one, leading to the update rule of (<ref>) for x_j.
    




§ COMPUTATION OF LOWER AND UPPER BOUNDS FOR DICHOTOMY


In Section <ref>, we introduced modifications to the MU and QU algorithms to incorporate the positivity constraint x≥ϵ and the linear constraint e^⊤x = 1. 
However, solving for the dual parameter ν in equations (<ref>) for MU or (<ref>) for QU is intractable. To address this, we propose using the dichotomy method to solve for h(ν) = 0. 
Therefore, this appendix provides the computation of lower bound ν_low and upper bound ν_up such that h(ν_low) < 0 and h(ν_up) > 0. 
These bounds will serve as convenient initializations for the dichotomy algorithm.



 §.§ Case 1: MU


For MU, we aim to solve equation (<ref>) for ν:

    h_1(ν)=∑_je_jx_j^tα_j/min(x_j^tα_j/ϵ,β_j^t+ν e_j)-1=0

Terms where e_j = 0 can be ignored since they do not contribute to the sum. Assuming ϵ > 0, the function h_1 is well-defined for ν∈ℝ. 
Since x_j^t α_j > 0, we have h_1(ν)=e_1/ϵ-1 for ν≤ν_lim=min_j(x_j^tα_j/ϵ-β_j^t/e_j), and h_1 is monotonically decreasing for ν∈[ν_lim,∞[. 
Assuming e_1/ϵ-1≥0 (which ensures the feasibility of the constraints x≥ϵ and e^⊤x = 1), we have h_1(ν_lim) ≥ 0 and lim_ν→∞ h_1(ν) = -1. Thus, there exists exactly one root for the function h_1.



  
Negative bound
First, let's find ν_low such that h_1(ν) < 0 for ν_low≤ν < ∞. We can bound h_1 as follows:

    h_1(ν)     =∑_jx_j^tα_j/min(x_j^tα_j/ϵ e_j,β_j^t/e_j+ν)-1
       ≤ nmax_jx_j^tα_j/min_jmin(x_j^tα_j/ϵ e_j,β_j^t/e_j+ν)-1<0

Therefore one possible bound is

    ν_low=nmax_jx_j^tα_j-min_jβ_j^t/e_j.





  
Positive bound
Similarly, let's find ν_up such that h_1(ν) ≥ 0 for ν_up≥ν≥ν_lim. Note that ν_lim is not a good bound when ϵ is small. We can bound h_1 as follows:

    h_1(ν)     =∑_j=1^nx_j^tα_j/min(x_j^tα_j/ϵ e_j,β_j^t/e_j+ν)-1
       ≥max_jx_j^tα_j/min(x_j^tα_j/ϵ e_j,β_j^t/e_j+ν)-1
       ≥max_jx_j^tα_j/β_j^t/e_j+ν-1>0,

As a result, we have 

    ν_up=max_j(x_j^tα_j-β_j^t/e_j)

To improve numerical stability, one could use ν_low^'=2nmax_jx_j^tα_j-min_jβ_j^t/e_j
and ν_up^'<max_jx_j^tα_j/2-β_j^t/e_j. 



 §.§ Case 2: QU


For QU, we aim to solve equation (<ref>) for ν:

    h_2(ν)=∑_je_jmax(-β_j^t-ν e_j+√((β_j^t+ν e_j)^2+4αζ_j^t)/2α,ϵ)-1=0.

Let's analyze the function h_2. We observe that the function -x + √(x^2 + δ) is strictly decreasing over ℝ for δ > 0, since its derivative -1 + x/√(x^2 + δ) is strictly negative for δ > 0. 
Therefore, each term of the sum is decreasing, and h_2 is a decreasing function.
Note that once at least for one j, we have 
    -β_j^t - ν e_j + √((β_j^t + ν e_j)^2 + 4αζ_j^t)≥ 2αϵ,
 and therefore the function h_2 becomes strictly decreasing. 
In the limit, we have lim_ν→ -∞ h_2(ν) = ∞ and lim_ν→∞ h_2(ν) = ϵe_1 - 1.
Assuming ϵe_1>1 (which ensures the feasibility of the constraints x≥ϵ and e^⊤x = 1), we know that the function h_2 has exactly one root.




  
Negative bound
Let's find ν_low such that h_2(ν) ≤ 0 for ν≥ν_low. We start by bounding the term

    -β_j^t-ν e_j+√((β_j^t+ν e_j)^2+4αζ_j^t)≤δ_j,

where δ_j > ϵ > 0. We move -β_j^t - ν e_j to the left and square the inequality to remove the square root:

    (β_j^t+ν e_j)^2+4αζ_j^t≤(β_j^t+ν e_j+ϵ_j)^2

Eventually, we can extract a bound for ν to ensure that the inequality is satisfied for a chosen δ_j:

    ν≥4αζ_j^t-ϵ_j^2/2δ_je_j-β_j^t/e_j.

Let's set δ_j = 2α/me_j, where m is the number of elements in the sum, and take the maximum over j to obtain the bound:

    ν_low≥max_j(mζ_j^t-α/me_j^2-β_j^t/e_j)

We observe the validity of this bound provided that δ_j = 2α/me_j≥ϵ for all j. Then, it can be verified that

    h_2(ν_low)     =∑_je_jmax(-β_j^t-ν_lowe_j+√((β_j^t+ν_lowe_j)^2+4αζ_j^t)/2α,ϵ)-1
        =∑_je_j-β_j^t-ν_lowe_j+√((β_j^t+ν_lowe_j)^2+4αζ_j^t)/2α-1
       ≤∑_je_jδ_j/2α-1
        =∑_j^me_j2α/me_j/2α-1=0,

which proves that ν_low is a valid negative bound.



  
Positive bound
Similarly, let's find ν_up such that h_2(ν) ≥ 0 for ν≤ν_up. 
This time, we will bound h_2 from below and obtain

    h_2(ν)     =∑_je_jmax(-β_j^t-ν e_j+√((β_j^t+ν e_j)^2+4αζ_j^t)/2α,ϵ)-1
       ≥∑_je_j-β_j^t-ν e_j+√((β_j^t+ν e_j)^2+4αζ_j^t)/2α-1
       ≥1/2α∑_je_j(-β_j^t-ν e_j)-1
        =-1/2α(∑_je_jβ_j^t+ν∑_je_j^2)-1>0.

We can, therefore, define

    ν_up=-2α+∑_je_jβ_j^t/∑_je_j^2.



plain


