A Survey on Surrogate Approaches to Non-negative Matrix Factorization

Motivated by applications in hyperspectral imaging we investigate methods for approximating a high-dimensional non-negative matrix $\mathbf{\mathit{Y}}$ by a product of two lower-dimensional, non-negative matrices $\mathbf{\mathit{K}}$ and $\mathbf{\mathit{X}}.$ This so-called non-negative matrix factorization is based on defining suitable Tikhonov functionals, which combine a discrepancy measure for $\mathbf{\mathit{Y}}\approx\mathbf{\mathit{KX}}$ with penalty terms for enforcing additional properties of $\mathbf{\mathit{K}}$ and $\mathbf{\mathit{X}}$. The minimization is based on alternating minimization with respect to $\mathbf{\mathit{K}}$ or $\mathbf{\mathit{X}}$, where in each iteration step one replaces the original Tikhonov functional by a locally defined surrogate functional. The choice of surrogate functionals is crucial: It should allow a comparatively simple minimization and simultaneously its first order optimality condition should lead to multiplicative update rules, which automatically preserve non-negativity of the iterates. We review the most standard construction principles for surrogate functionals for Frobenius-norm and Kullback-Leibler discrepancy measures. We extend the known surrogate constructions by a general framework, which allows to add a large variety of penalty terms. The paper finishes by deriving the corresponding alternating minimization schemes explicitely and by applying these methods to MALDI imaging data.


Introduction
Matrix factorization methods for large scale data sets have seen increasing scientific interest recently due to their central role for a large variety of machine learning tasks. The main aim of such approaches is to obtain a low-rank approximation of a typically large data matrix by factorizing it into two smaller matrices. One of the most widely used matrix factorization method is the principal component analysis (PCA), which uses the singular value decomposition (SVD) of the given data matrix. In this work, we review the particular case of non-negative matrix factorization (NMF), which is favorable for a range of applications where the data under investigation naturally satisfies a non-negativity constraint. These include dimension reduction, data compression, basis learning, feature extraction as well as higher level tasks such as classification or clustering [11,26,27,30]. PCA based approaches without any non-negativity constraints would not lead to satisfactory results in this case since possible negative entries of the computed matrices cannot be easily interpreted for naturally non-negative datasets. Typically, the NMF problem is formulated as a minimization problem. The corresponding cost function includes a suitable discrepancy term, which measures the difference between the data matrix and the calculated factorization, as well as penalty terms to tackle the non-uniqueness of the NMF, to deal with numerical instabilities but also to provide the matrices with desirable properties depending on the application task. The NMF cost functions are commonly non-convex and require tailored minimization techniques to ensure the minimization but also the non-negativity of the matrix iterates. This leads us to the so called surrogate minimization approaches, which are also known as majorize-minimization algorithms [18,23,36]. Such surrogate methods have been investigated intensively for some of the most interesting discrepancy measures and penalty terms [11,13,16,23,25,33,36]. The idea is to replace the original cost function by a so called surrogate functional, such that its minimization induces a monotonic decrease of the objective function. It should be constructed in such a way that it is easier to minimize and that the deduced update rules should preserve the non-negativity of the iterates, which typically leads to alternating, multiplicative update rules. It appears, that these constructions are obtained case-by-case employing different analytical approaches and different motivations for their derivation. The purpose of this paper, first of all, is to give a unified approach to surrogate constructions for NMF discrepancy functionals. This general construction principle is then applied to a wide class of functionals obtained by different combinations of divergence measures and penalty terms, thus extending the present state of the art for surrogate based NMF-constructions. Secondly, one needs to develop minimization schemes for these functionals. Here we develop concepts for obtaining multiplicative minimization schemes, which automatically preserve non-negativity without the need for further projections. Finally, we exemplify some characteristic properties of the different functionals with MALDI imaging data, which are particularly high-dimensional and challenging hyperspectral data sets. The paper is organized as follows. Section 2 introduces the basic definition of the considered NMF problems. Section 3 gives an overview about the theory of surrogate functionals as well as the construction principles. This is then exemplified in Section 4 for the most important cases of discrepancy terms, namely the Frobenius norm and the Kullback-Leibler divergence as well as for a variety of penalty terms. Section 5 discusses alternating minimization schemes for these general functionals with the aim to obatin non-negativity-preserving, multiplicative iterations. Finally, Section 6 contains numerical results for MALDI imaging data.

Notation
Throughout this work, we will denote matrices in bold capital Latin or Greek letters (e.g. Y , K , Ψ , Λ) while vectors will be written in small bold Latin or Greek letters (e.g. c, d , β, ζ). The entries of matrices and vectors will be indicated in a non-bold format to distinguish between the i-th entry x i of a vector x and n different vectors x j for j = 1, . . . , n. In doing so, we write for the entry of a matrix M in the i-th row and the j-th column M ij and the i-th entry of a vector x the symbol x i . The same holds for an entry of a matrix product: the ij-th entry of the matrix product MN will be indicated as (M N ) ij . Furthermore, we will use a dot notation to indicate rows and columns of matrices. For a matrix M we will write M •,j for the j-th column and M i,• for the i-th row of the matrix. What is more, we will use · for the usual Euclidean norm, M 1 := ij |M ij | for the 1-norm and M F for the Frobenius norm of a matrix M . Besides that, we will use equivalently the terms function and functional for a mapping into the real numbers. Finally, the dimensions of the matrices in the considered NMF problem are reused in this work and will be introduced in the following section.

Non-negative Matrix Factorization
Before we introduce the basic NMF problem, we give the following definition to clarify the meaning of a non-negative matrix.
The non-negativity of an arbitrary matrix M will be abbreviated for simplicity as M ≥ 0 in the later sections of this work. The basic NMF problem requires to approximately decompose a given non-negative matrix Y ∈ R n×m ≥0 into two smaller non-negative matrix factors K ∈ R n×p ≥0 and X ∈ R p×m ≥0 , such that p min(n, m) and For an interpretation let us assume, that we are given m data vectors Y •,j ∈ R n for j = 1, . . . , m, which are stored column-wise in the matrix Y . Similarly for k = 1, . . . , p we denote by K •,k , respectively X k,• , the column vectors of K , respectively the row vectors of X . We then obtain the following approximation for Furthermore, we define accordingly The corresponding matrix divergences are defined componentwise, i.e. β = 2 yields the Frobenius norm and β = 1 the Kullback-Leibler divergence. These discrepancy measures are typically amended by so called penalty terms for stabilization and for enforcing additional properties such as sparsity or orthogonality. This yields the following general minimization task.
Definition 4 (NMF Minimization Problem) For a data matrix Y ∈ R n×m ≥0 , we consider the following generalized NMF minimization task The functional is called the cost functional. Furthermore, we call (i) D β (Y , KX ) the discrepancy term, (ii) α the regularization parameters or weights (iii) and ϕ (K , X ) the penalty terms.
The functional in (4) is typically non-convex in (K , X ). Hence, algortihms based on alternating minimization with respect to K or X are favourable, i.e.
where the index d denotes the iteration index of the corresponding matrices. This yields simpler, often convex restricted problems with respect to either K or X . Considering for example the minimization of the NMF functional with Frobenius norm and without any penalty term, yields a high dimensional linear system K Y = K K X , which, however, would need to be solved iteratively. Instead, so called surrogate methods for computing NMF decompositions have been proposed recently and are introduced in the next section. They also consider alternating minimization steps for K and X , but they replace the restricted minimization problems in (5) and (6) by simpler minimization tasks, which are obtained by locally replacing F by surrogate functionals for K and X separately.

Surrogate Functionals
In this section, we discuss general surrogate approaches for minimizing general non-convex functionals, which are then exemplified for specific NMF functionals in later sections. Let us consider a general functional F : Ω → R where Ω ⊂ R N and the minimization problem min x ∈Ω F (x ).
We will later add suitable conditions guaranteeing the existence of minimizers or at least the existence of stationary points. Surrogate concepts replace this task by solving a sequence of comparatively simpler and convex surrogate functionals, which can be minimized efficiently. These methods are also commonly referred to as surrogate minimization (or maximization) algorithms (SM) or also as MM algorithms, where the first M stands for majorize and the second M for minimize (see also [18,23,36]). Such approaches have been demonstrated to be very useful in many fields of inverse problems, in particular for hyperspectral imaging [11], medical imaging applications such as transmisson tomography [13,14] as well as MALDI imaging and tumor typing applications [26]. Replacing a non-convex functional by a series of convex problems is the main motivation for such surrogate approaches. However, if constructed appropriately, they can also be used to replace non-differentiable functionals by a series of differentiable problems and they can be tailored such that gradient-descent methods for minimization yield multiplicative update rules which automatically incorporate non-negativity constraints without further projections. From this point on it is important to note that possible zero denominators during the derivation of the NMF algorithms as well as in the multiplicative update rules themselves will not be discussed explicitly throughout this work. Usually, this issue is handled in practice by adding a small positive constant in the denominator during the iteration scheme. In fact, the instability of NMF algorithms due to the convergence of some entries in the matrices to zero is not sufficiently discussed in the literature and still needs proper solution techniques. We will not focus on this problem and turn now to the basic definition and properties of surrogate functionals.

Definitions and Basic Properties
As in [25], we use the following definition of a surrogate functional.

Definition 5 (Surrogate Functional)
Let Ω ⊆ R N denote an open set and F : Ω → R a functional defined on Ω. Then Q F : Ω ×Ω → R is called a surrogate functional or a surrogate for F, if it satisfies the following conditions: This is the most basic definition, which does not require any convexity or differentiability of the functional. However, it already allows to prove that the iteration yields a sequence which monotonically decreases F .
Proof The monotone decrease (9) follows directly from the defining properties of surrogate functionals, see Definition 5: We obtain where ( ) follows from the definition of x [d+1] in (8).

Remark 1 (Addition of Surrogate Functionals)
Let Ω ⊆ R n be an open set, F, G : Ω → R pointwise defined functionals and Q F , Q G corresponding surrogates. Then Q F + Q G is a surrogate functional for F + G.
For each functional F there typically exist a large variety of surrogate functionals and we can aim at optimizing the structure of surrogate functionals. The following additional property is the key to simple and efficient minimization schemes for surrogate functionals.
Lemma 1 above only ensures the monotonic decrease of the cost functional, which is not sufficient to guarantee convergence of the sequence {x [d] } to a minimizer of F or at least to a stationary point of F . The convergence theory for surrogate functionals is far from being complete (see also the works [23] and [36]). Despite this lack of theoretical foundation, surrogate based minimization yields strictly decreasing sequences for a large variety of applications. In particular, surrogate based methods can be constructed such that first order optimality conditions lead to multiplicative update rules, which -in view of applications to NMF constructions -is a very desirable property. We now turn to discussing three different construction principles for surrogate functionals.

Jensen's Inequality
The starting point is the well known Jensen's inequality for convex functions, see [10].
Lemma 2 (Jensen's Inequality) Let Ω ⊆ R N denote a convex set, F : Ω → R a convex function and λ i ∈ [0, 1] non-negative numbers for i ∈ {1, . . . , k} with In this subsection we consider functionals F which are derived from continuously differentiable and convex functions f : for Ω ⊆ R N ≥0 and some auxiliary variable c ∈ Ω. This also implies, that F is convex, since We now choose λ i ∈ [0, 1] with N i=1 λ i = 1 and α ∈ R N and define for some b ∈ Ω. This implies The functional Q F : Ω × Ω → R defines a surrogate for F , which can be seen by the inequality above and by observing

Low Quadratic Bound Principle
This concept is based on a Taylor expansion of F in combination with a majorization of the quadratic term. This so called low quadratic bound principle (LQBP) has been introduced in [5] and was used in particularly for the computation of maximum-likelihood estimators. These methods do not require that F itself is convex and its construction is based on the following lemma.
We then obtain a quadratic majorization =: and Q f is a surrogate functional for f.
Proof The proof of this classical result is based on the second-order Taylor polynomial of f and shall be left to the reader.
The related update rule for surrogate minimization can be stated explicitly under natural assumptions on the matrix Λ.
Corollary 1 Assume that the assumptions of Lemma 3 hold. In addition, assume that Λ is a positive definite and symmetric matrix. Then, the corresponding surrogate Q f is strictly convex in its first variable and we have from (8) Proof For an arbitrary α ∈ {1, . . . , N }, we have that where ( ) utilizes the symmetry of Λ. Hence, it holds that . The Hessian matrix of Q f then satisfies This implies the positive definiteness of the functional, hence, it has a uniqie minimizer, which is given by This is the update rule above.
The computation of the inverse of Λ is particularly simple if Λ is a diagonal matrix. Furthermore, the diagonal structure ensures the separability of the surrogate functional mentioned in Definition 6. Therefore, we consider matrices of the form where κ i ≥ 0 has to be chosen individually depending on the considered cost function. We will see that an appropriate choice of κ i will lead finally to the desired multiplicative update rules of the NMF algorithm.
The matrix Λ(a) in (18) fulfills the conditions in Corollary 1 as it can be seen by the following lemma. Therefore, if Λ is constructed as in (18), the update rule in (17) can be applied immediately.
denote a symmetric matrix. With a ∈ R N >0 and κ i ≥ 0, we define the diagonal matrix Λ, such that for i = 1, . . . , N. Then Λ and Λ − M are positive semi-definite.
Proof Let ζ ∈ R N denote an arbitrary vector and let δ denote the Kronecker symbol. Then The positive semi-definiteness of Λ follows from its diagonal structure.

Further Construction Principles
So far we have discussed two major construction principles based on either Jensen's inequality or on upper bounds for the quadratic term in Taylor expansions. [23] lists further construction principles, which however will not be used for NMF constructions in the subsequent sections of this paper. For completeness we shortly list their main properties. A relaxation of the approach based on Jensen's inequality is achieved by choosing A typical choice is |c j | p which leads to surrogate functionals for p ≥ 0. This type of surrogate was originally introduced in the context of medical imaging, see [12], for positron emission tomography. Another approach is based on combining arithmetic with geometic means and can be used for constructing surrogates for posynomial functions. For α, v , a ∈ R N >0 , we obtain

Surrogates for NMF Functionals
In this section we apply the general construction principles of Section 3 to the NMF problem as stated in (3). The resulting functional F (K , X ) depends on both factors of the matrix decomposition and minimization is attempted by alternating minimization with respect to K and X as described in (5) and (6). However, we replace the functional F in each iteration by suitable surrogate functionals, which allow an explicit minimization. Hence, we avoid the minimization of F itself, which even for the most simple quadaratic formulation requires to solve a high dimensional linear system. We start by considering the discrepancy terms for β = 2 (Frobenius norm) and β = 1 (Kullback-Leibler divergence) and determine surrogate functionals with respect to X and K . We then add several penalty terms and develop surrogate functionals accordingly. With regard to the construction of surrogates for the case of β = 0, which leads to the so-called Itakura-Saito divergence, we refer to the works [15,16,32].

Frobenius Discrepancy and Low Quadratic Bound Principle
We start by constructing a surrogate for the minimization with respect to X for the Frobenius discrepancy Let Y •,j , resp. X •,j , denote the column vectors of Y , resp. X . The separability of F yields Hence, the minimization separates for the different f Y•,j terms. The Hessian of these terms is given by and the LQBP construction principle of the previous section with κ k = 0 yields leading to the surrogate functionals An appropriate choice of κ k ensures the multiplicativity of the final NMF algorithm. In the case of the Frobenius discrepancy term, we will see that suitable κ k can be chosen dependent on 1 regularization terms in the cost function, which are not included up to now (see Subsection 4.4 and Appendix A.1 for more details on this issue). Due to the absent 1 terms, we set κ k = 0 to get the desired multiplicative update rules. Summing up the contributions of the columns of X yields the final surrogate The equivalent construction for K can be obtained by regarding the rows of K separately, which for We summarize this surrogate construction in the following theorem.

Theorem 1 (Surrogate Functional for the Frobenius Norm with LQBP)
We consider the cost functionals F (X ) : define separable and convex surrogate functionals.

Frobenius Discrepancy and Jensen's Inequality
Again we focus on deriving a surrogate functional for X , the construction for K will be very similar. Expanding the Frobenius discrepancy yields Putting v := X •,j ∈ R p ≥0 and c := K i,• ∈ R p ≥0 allows us to define such that Hence we have separated the Forbenius discrepancy suitably for applying Jensen's inequality. Following the construction principle in Subsection 3.2, we define with the auxiliary variable A ∈ R p×m Inserting this into the decomposition of the Frobenius discrepancy yields the surrogate Q F,2 (X , A) by The construction of a surrogate for K proceeds in the same way. We summarize the results in the following theorem.
Theorem 2 (Surrogate Functional for the Frobenius Norm with Jensen's Inequality) We consider the cost functionals F (X ) : define separable and convex surrogate functionals.
These surrogates are equal to the ones proposed in [11]. We will later use first order necessary conditions of the surrogate functionals for obtaining algorithms for minimization. We already note i.e. despite the rather different derivations, the update rules for the surrogates obtained by LQBP and Jensen's inequality will be identical.

Surrogates for Kullback-Leibler Divergence
The case β = 1 in Definition 3 yields the so-called Kullback-Leibler divergence (KLD). For matrices M , N ∈ R n×m >0 , it is defined as and has been investigated intensively in connection with non-negative matrix factorization methods [11,16,24,25]. In our context, we define the cost functional for the NMF construction by We will focus in this subsection on Jensen's inequality for constructing surrogates for the KLD since they will lead to the known classical NMF algorithms (see also [11,24,25]). However, it is also possible to use the LQBP principle to construct a suitable surrogate functional for the KLD which leads to different, multiplicative update rules (see Appendix B).
We start by deriving a surrogate for the minimization with respect to X , i.e. we consider Using the same λ k and α k as in the section above and applying it to the convex function f (t) := − ln(t), we obtain Multiplication with Y ij ≥ 0 and the addition of appropriate terms yields The condition Q F,1 (X , X ) = F (X ) follows by simple algebraic manipulations, such that Q F,1 is a valid surrogate functional for F. The approach by Jensen's inequality is very flexible and we obtain different surrogate functionals Q F,2 and Q F,3 by using i.e.
Inserting the same λ k and α k as before in Equation (15), we obtain immediately the surrogates It easy to check, that the partial derivatives for all three variants are the same, hence, the update rules obtained in the next section based on first order optimality conditions will be identical. Applying the same approach for obtaining a surrogate for K yields the following theorem. Then define separable and convex surrogate functionals.

Surrogates for 1 -and 2 -Norm Penalties
Computing an NMF is an ill-posed problem, see [11], hence one needs to add stabilizing penalty terms for obtaining reliable matrix decompositions. The most standard penalties are 1 -and 2 -terms for the matrix factors leading to for β ∈ {1, 2}. The 2 -penalty prohibits exploding norms for each matrix factor and the 1 -term promotes sparsity in the minimizing factors, see [21,28] for a general exposition.
Combinations of 1 -and 2 -norms are sometimes called elastic net regularization, [20], due to there importance in medical imaging. These penalty terms are convex and they separate, hence, they can be used as surrogates themselves. For the case of Kullback-Leibler divergences this leads to the following surrogate for minimization with respect to X : where Q KL is the surrogate for the Kullback-Leibler divergence of Theorem 3 for X .
The Frobenius case cannot be treated in the same way. If we use the penalty terms as surrogates themselves and obtain the standard minimization algorithm by first order optimality conditions, then this does not lead to a multiplicative algorithm, which preserves the non-negativity of the iterates. It can be easily seen that the which yields the Hessian ∇ 2 f y (a) = K K + νI . The choice of κ k is done dependent on the 1 regularization term of the cost function f y as already described in Subsection 4.1. It can be shown in the derivation of the NMF algorithm that κ k = λ for all k leads to multiplicative update rules. A more general cost function is considered in Appendix A.1, where the concrete effect of κ k is described in more detail. This yields the following diagonal matrix Λ f y (a): The surrogate for minimization with respect to X is then Similar, for minimization with respect to K we obtain the surrogate by using the diagonal matrix Λ g y (a) kk := (a(XX + µI p×p )) k + ω a k .

Surrogates for Orthogonality Constraints
The observation that a non-negative matrix with pairwise orthogonal rows has at most one non-zero entry per column is the motivation for introducing orthogonality constraints for K or X . This will lead to strictly uncorrelated feature vectors, which is desirable in several applications e.g. for obtaining discriminating biomarkers from mass spectra, see Section 6 on MALDI Imaging. We could add the orthogonality constraint K K = I as an additional penalty term σ K K K − I 2 . However, this would introduce fourth order terms. Hence we introduce additional variables V and W and split the orthogonality condition into two second order terms leading to Surrogates for the terms I − V K 2 F and I − X W 2 F can be calculated via Jensen's inequality (see Subsection 4.2). The other penalties can be used as surrogates themselves and therefore, we obtain Theorem 4 (Surrogate Functionals for Orthogonality Constraints) We consider the cost functionals define separable and convex surrogate functionals.

Surrogates for Total Variation Penalties
Total variation (TV) penalty terms are the second important class of regularization terms besides p -penalty terms. TV-penalties aim at smooth or even piecewise constant minimizers, hence they are defined in terms of first order or higher order derivatives [7]. Originally, they were introduced for denoising applications in image processing [31] but have since been applied to inpainting, deconvolution and other inverse problems, see e.g. [8]. The precise mathematical formulation of the total variation in the continuous case is described in the following definition.

Definition 7 (Total Variation (Continuous))
Let Ω ⊂ R N be open and bounded. The total variation of a function u ∈ L 1 loc (Ω) is defined as There exist several analytic relaxations of TV based on 1 -norms of the gradient, which are more tractable for analytical investigations. For numerical implementations one rather uses the L 1 -norm of the gradient ∇f L 1 as a more computationally tractable substitute. For discretization the gradient is typically replaced by a finite difference approximation [9]. For applying TV-norms to data, we assume that the row index in the data matrix Y refers to spatial locations and the column index to so-called channels. In this case, we consider the most frequently used isotropic TV for applying it to measured, discretized hyperspectral data.

Definition 8 (Total Variation (Discrete))
For fixed ε TV > 0, the total variation of a matrix K ∈ R n×p ≥0 is defined as where ψ k ∈ R ≥0 is a weighting of the k-th data channel and N i ⊆ {1, . . . , n} \ {i} denotes the index set referring to spatially neighboring pixels.
We will use the following short hand notation which can be seen as a finite difference approximation of the gradient magnitude of the image K •,k at pixel K ik for some neighbourhood pixels defined by N i . A typical choice for neighbourhood pixels in two dimensions for the pixel (0, 0) is N (0,0) := {(1, 0), (0, 1)} to get an estimate of the gradient components along both axes. Finally, by introducing the positive constant ε TV > 0, we get a differentiable approximation of the total variation penalty. In Section 6, we will discuss the application of NMF methods to hyperspectral MALDI imaging datasets, which has a natural 'spatial structure' in its columns. In this section we stay with a generic choice of N i as well as of the ψ k and we construct a surrogate following the approach of the groundbreaking works of [13] and [29]. For t ≥ 0 and s > 0 we use the inequality (linear majorization) and apply it in order to compare ∇ ik K with values obtained by an arbitrary non-negative matrix A: Summation with respect to i, multiplication with ψ k and summation with respect to k leads to This yields a candidate for a surrogate Q Oli TV for the TV-penalty term, which is the same as the one used in [29]. However, it is not separable, hence we aim at a second, separable approximation. For arbitrary a, b, c, d ∈ R we have This leads to Therefore, we have the following Theorem 5 (Surrogate Functional for TV Penalty Term) We consider the cost functional F (K ) := TV(K ) with the total variation defined in (34). Then defines a separable and convex surrogate functional.
The separability of the surrogate is not obvious. The proof (see Appendix C) delivers the following notation, which we also need for an description of the algorithms in the next section. First of all we need the definition of the so-called adjoint neighborhood pixelsN i given by One then introduces matrices P (A) ∈ R n×p ≥0 and Z(A) ∈ R n×p ≥0 via Using these notations, it can be shown that the surrogate can be written as such that we obtain the desired separability. Here, C(A) denotes some function depending on A. The description of Q TV with the help of P (A) ik and Z(A) ik will also allow us to compute the partial derivatives in a way more comfortable way (see also Appendix A.2).

Surrogates for Supervised NMF
As a motivation for this section, we consider classification tasks. We view the data matrix Y as a collection of n data vectors, which are stored in the rows of Y . Moreover, we do have an expert annotation u i for i = 1, . . . , n, which assigns a label to each data vector. For a classification problem with two classes we have u i ∈ {0, 1}. As already stated, the rows of the matrix X of an NMF decomposition can be regarded as a basis for approximating the rows of Y . Hence, one assumes that the correlations between a row Y i,• of Y and all row vectors of X , i.e. computing Y i,• X , contains the relevant information of Y i,• . The vector of correlations yields a so-called feature vector of length p. A classical linear regression model, which uses these feature vectors, then asks to compute weights β k for k = 1, . . . , p, such that Y i,• X β ≈ u i (for more details on linear discriminant analysis methods, we refer to Chapter 4 in [4]).
In matrix notation and using least squares, this is equivalent to computing β as a minimizer of u − Y X β 2 .
We now use X and β to define In tumor typing classifications, where the data matrix Y is obtained by MALDI measurements, the vector x * can be interpreted as a characteristic mass spectrum of some specific tumor type and can be directly used for classification tasks in the arising field of digital pathology (see also Section 6 and [26]). The classification of a new data vector y is then simply obtained by computing the scalar product w = x * y and assigning either the class label 0 or 1 by comparing w with a pre-assigned threshold s. This threshold is typically obtained in the training phase of the classification procedure by computing YX β for some given training data Y and choosing s, such that a performance measure of the classifier is optimized. The approach we have described is based on first computing an NMF, i.e. K and X , and then computing the weights β of the classifier. Hence, the computation of the NMF is done independently of the class labels u, which is also referred to as an unsupervised NMF approach. We might expect, that computing the NMF by minimizing a functional which includes the class labels, i.e.
will lead to an improved classifier. In the application field of MALDI imaging, this supervised approach yields an extraction of features from the given training data, which allow a better distinction between spectra obtained from different tissue types such as tumorous and non-tumorous regions (see also [26]). Surrogates for the first term have been determined in the previous section for the case of the Frobenius norm and the Kullback-Leibler divergence. Hence, we need to determine surrogates of the new penalty term for minimization with respect to X and β: Surrogates can be obtained by extending the Jensen principle to the matrix valued case. Here, we consider a convex subset Ω ⊂ R N ×M >0 and definẽ with a convex and continuously differentiable function f and auxiliary variables c ∈ R N >0 and d ∈ R M >0 . We now use the following generalized Jensen's inequality Setting for some i ∈ {1, . . . , n} yields by inserting λ jk and α jk into (45) The computation of a surrogate for minimization with respect to β proceeds analogously. We summarize the results in the following theorem.
A big advantage of linear regression models are their simplicity and manageability. However, they are by far not the optimal approach to approximate the binary output data u with a continuous input. Logistic regression models offer a way more natural method for binary classification tasks. Together with the supervised NMF as a feature extraction method, this overall workflow leads in [26] to excellent classification results and outperformed classical approaches. However, the proposed model is based on a gradient descent approach, such that the non-negativity of the iterates can only be guaranteed by a projection step. Appropriate surrogate functionals for this workflow are still ongoing research and could lead to even better outcomes (see also [35,36]).

Surrogate Based NMF Algorithms
In the previous section we have defined surrogate functionals for various NMF cost functions. Besides the necessary surrogate properties we also expect that the minimization of these surrogates is straightforward and can be computed efficiently. In our case we demand additionally, that the minimization schemes based on solving the first order optimality conditions leads to a separable algorithm and that it only requires multiplicative updates, which automatically preserve the nonnegativity of its iterates. Let us start with denoting the most general functional with Kullback-Leibler divergence, the Frobenius case follows similarly. For constructing non-negative matrix factorizations, we incorporate 2 -, sparsity-, orthogonality-, TV-constraints and also the penalty terms coming from the supervised NMF. Of course, in most applications one only uses a subset of these constraints for stabilization and for enhancing certain properties. These algorithms can readily be obtained from the general case by putting the respective regularization parameters to zero. The corresponding update rules are classical results and can be found in numerous works [11,24,25].
≥0 , β ∈ R p ≥0 and a set of regularization parameters λ, µ, ν, ω, τ, σ K ,1 , σ K ,2 , σ X ,1 , σ X ,2 , ρ ≥ 0, we define the NMF minimization problem by The choice of the various regularization parameters occuring in Definition 9 is often based on heuristic approaches. We will not focus on that issue in this work and refer instead to [19] and the references therein, where two methods are investigated for the general case of multi-parameter Tikhonov regularization. The algorithms studied in this section will start with positive initializations for K , X , V , W and β. These matrices are updated alternatingly, i.e. all matrices except one matrix are kept fixed and only the selected matrix is updated by solving the respective first order optimality condition. We will focus in this section on the derivation of the update rules of K (see also Appendix A.2). The iteration schemes for the other matrices follow analogously. For that, we only have to consider those terms in the general functional which depend on K , i.e. we aim at determining a minimizer for Instead of minimizing this functional we exchange it with the previously constructed surrogate functionals, which leads to with the surrogates Q KL for the Kullback-Leibler divergence in Theorem 3, Q TV for the TV penalty term in Theorem 5 and Q Orth for the orthogonality penalty terms in Theorem 4. The computation of the partial derivatives leads to a system of equations for ξ ∈ {1, . . . , n} and ζ ∈ {1, . . . , p}. This leads to a system of quadratic equations Solving for K ξζ and denoting the Hadamard product by • as well as the matrix division for each entry separately by a fraction line, yields the following update rule for K . (Note that the notation for P (A) and Z(A) was introduced in the section on TV-regularization above.) In the above update rule, 1 n×p denotes an n × p matrix with ones in every entry and Ψ ∈ R n×p ≥0 is defined as Details on the derivation can be found in Appendix A.2.
The partial derivatives with repect to X are computed similarly. Defining leads to the update The updates for V , W are straight forward and we obtain the following theorem.

Theorem 7 (Alternating Algorithm for NMF Problem in Definition 9)
The initializations and the iterative updates lead to a monotonically decrease of the cost functional It is easy to see that the classical, regularized NMF algorithms described in [11,24,25] can be regained by putting the corresponding regularization parameters to zero. In the case of 1 -and 2 -regularized NMF, this leads to the cost function The classical update rule for X is obtained by setting which -in connection with the update rule for X of the previous theorem -leads to = 2Λ [d] which is the update rule described in [11]. By the same approach and with the surrogate functionals derived in Section 4, we obtain the update rules for the Frobenius discrepancy term, i.e. we consider the functional A monotone decrease of this functional is obtained by the following iteration in combination with the update rules for V , W , β as in Theorem 7 (see also Appendix A.1 for more details on the derivation of these algorithms.) .

MALDI Imaging
As a test case we analyse MALDI imaging data (matrix assisted laser desorption/ionization) of a rat brain. MALDI imaging is a comparatively novel modality, which unravels the molecular landscape of tissue slices and allows a subsequent proteomic or metabolic analysis [1,6,22]. Clustering this data reveals for example different metabolic regions of the tissue, which can be used for supporting pathological diagnosis of tumors. The data used in this paper was obtained by a MALDI imaging experiment, see Figure 2 for a schematic experimental setup. In our numerical experiments, we used a classical rat brain dataset which has been used in several data processing papers before [2,3,22]. It constitutes a standard test set for hyperspectral data analysis. The tissue slice was scanned at 20185 positions. At each position a full mass spectrum with 2974 m/z (mass over charge) values was collected. I.e. instead of three color channels, as it is usual in image processing, this data has 2974 channels, each channel containing the spatial distribution of molecules having the same m/z value.
The following numerical examples were obtained with the multiplicative algorithms described in the previous section. We just illustrate the effect of the different penalty terms for some selected functionals. One can either display the columns of K as the pseudo channels of the NMF decomposition or the rows of X as pseudo spectra characterizing the different metabolic processes present in the tissue slice, see the Figures 3-6. Both ways of visualization do have their respective value. Looking at the pseudo spectra in connection with orthogonality constraints leads to a clustering of the spectra and to a subdivision of the tissue slice in regions with potentially different metabolic activities, see [22]. Considering instead the different pseudo spectra, which were constructed in order to have a bases which allows a low dimensional approximaton of the data set, is the basis for subsequent proteomic analysis. E.g. one may target pseudospectra where the related pseudo channels are concentrated in regions, which were annotated by pathological experts. Mass values which are dominant in those spectra may stem from proteins/peptides relating to biomarkers as indicators for certain diseases. Hence, classification schemes based on NMF decompositions have been widely investigated [26,30,34].  Fig. 3: NMF of the rat brain dataset for p = 6. Orthogonality constraints on the channels with σ K ,1 = 1 and σ K ,2 = 1.

Conclusion
In this paper, we investigated methods based on surrogate minimization approaches for the solution of NMF problems. The interest in NMF methods is related to its importance for several machine learning tasks. Application for large data sets require that the resulting algorithms are very efficient and that iteration schemes only need simple matrix-vector multiplications. The state of the art for constructing appropriate surrogates are based on case-bycase studies for the different, considered NMF models. In this paper, we embedded the different approaches in a general framework, which allowed us to analyze several extensions to the NMF cost functional, including 1 -and 2 -regularization, orthogonality constraints, total variation penalties as well as extensions, which leaded to supervised NMF concepts. Secondly, we analyzed surrogates in the context of the related iteration schemes, which are based on first order optimality conditions. The requirement of separability as well as the need of having multiplicative updates, which preserve nonnegativity without additional projections, were analyzed. This resulted in a general description of algorithms for alternating minimization of constrained NMF functionals. The potential of these methods is confirmed by numerical tests using hyperspectral data from a MALDI imaging experiment. Several further directions of research would be of interest. First of all, besides the most widely used penalty terms discussed in this paper further penalty terms, e.g. higher order TV-terms, could be considered. Secondly, construction principles for more general discrepancy terms could be analyzed (see also [16]). Potentially more importantly, this paper contains only very first results for combining NMF constructions directly with subsequent classification tasks. The question of an appropriate surrogate functional for the supervised NMF model with logistic regression used in [26] remains unanswered and also the comparison with algorithmic alternatives such as ADMM methods needs to be explored.  6: NMF of the rat brain dataset for p = 6. Orthogonality constraints on the channels with σ K ,1 = 10 and σ K ,2 = 10 and sparsity penalty term on X with λ = 0.06. The sparsity penalty term has in connection with the orthogonality constraint a comparatively strong influence on the NMF computation: The sparsity in the spectra increases significantly and thus their biological interpretability, whereas the anatomic structure in the pseudo channels diminishes.

Appendix A Details on the Derivation of the Algorithms in Section 5
In this section, we give a more detailed derivation of the algorithms presented in Section 5. We start with the less complex case of the Frobenius norm as discrepancy term and then turn to the Kullback-Leibler divergence. To cover both aspects, we derive the update rules of X for the Frobenius discrepancy term and of K in the case of the KLD. We will also take a closer look at the effect of κ in equation (18) with respect to the LQBP construction principle.

A.1 Frobenius Norm
We consider the general cost function described in Section 5 for the case of the Frobenius norm. To compute the update rules for X , it is enough to examine the function where all terms independent from X are omitted. Based on Remark 1 and following the discussion of Section 4.4, the construction of a surrogate for F can be done separately for F 1 and the remaining penalty terms. The construction of a surrogate for F 1 with the LQBP principle as it has been done similarly in Subsection 4.4 is essential. If we would use instead a surrogate for the discrepancy term 1 /2 Y − KX 2 F from Subsection 4.1 or 4.2 and take the 1 -penalty term λ X 1 as surrogate itself, it is easy to see that this would not lead to multiplicative update rules. It is the 1 -penalty term which causes the difficulty. Computing the first order optimality condition for the corresponding surrogateQ F (X , A) =: λ X 1 +Q F (X , A) with respect to X would lead to where the second term on the right hand side does not depend on λ. Hence, we get a sign in front of λ by solving the equation for X ξζ and we will not obtain multiplicative updates for X . A correct surrogate is obtained by using the LQBP principle to F 1 and leads to and with the diagonal matrix The functionals Q Orth resp. Q LR are the surrogates obtained from Theorem 4 resp. Theorem 6. It will turn out that an appropriate choice of κ k will ensure a multiplicative NMF algorithm. The computation of the first order optimality condition for Q F leads to One can see immediately, that the choice of κ ξ := λ for all ξ ∈ {1, . . . , p} is appropriate to get rid of the problematic term λ. Hence, we obtain Reordering the terms leads to (K Y ) ξζ + (σ X ,1 + σ X ,2 )W ξζ + ρβ ξ (Y u) ζ = X ξζ A ξζ (K KA) ξζ + νA ξζ + λ + ρβ ξ (Y Y A β) ζ + σ X ,1 (AW W ) ξζ + σ X ,2 A ξζ .
By exploiting the surrogate minimization principle as described in Lemma 1, we get finally the update rule for X presented in Section 5.
This equation holds for arbitrary ξ ∈ {1, . . . , n} and ζ ∈ {1, . . . , k}. We therefore can extend this relation to the whole matrix K and obtain which is exactly the described update rule in Section 5.

Appendix B Kullback-Leibler Divergence Discrepancy and LQBP
In this Section, we will use the LQBP construction principle to derive a multiplicative algorithm for the cost function It follows for the partial derivatives of f such thatP st (A) · K 2 st takes all quadratic terms of the matrix entries of K in the surrogate functional into account. The same can be done with the linear terms K st , which leads to the coefficient Therefore, the surrogate Q TV can be written as for some functionC, which only depends on A. This shows the separability of the surrogate.