Non-negative Matrix Factorization via Archetypal Analysis

TL;DR

This paper presents a novel approach to non-negative matrix factorization (NMF) via archetypal analysis that does not require data separability, offering a generally unique decomposition and robustness under certain conditions.

Contribution

It introduces a new NMF method based on optimizing convex envelopes of data and archetypes, extending archetypal analysis without the separability constraint.

Findings

The proposed method achieves a generally unique decomposition.

Standard optimization methods find good solutions in practice.

The approach is robust under a new uniqueness condition.

Abstract

Given a collection of data points, non-negative matrix factorization (NMF) suggests to express them as convex combinations of a small set of `archetypes' with non-negative entries. This decomposition is unique only if the true archetypes are non-negative and sufficiently sparse (or the weights are sufficiently sparse), a regime that is captured by the separability condition and its generalizations. In this paper, we study an approach to NMF that can be traced back to the work of Cutler and Breiman (1994) and does not require the data to be separable, while providing a generally unique decomposition. We optimize the trade-off between two objectives: we minimize the distance of the data points from the convex envelope of the archetypes (which can be interpreted as an empirical risk), while minimizing the distance of the archetypes from the convex envelope of the data (which can be interpreted as a data-dependent regularization). The archetypal analysis method of (Cutler, Breiman, 1994) is recovered as the limiting case in which the last term is given infinite weight. We introduce a `uniqueness condition' on the data which is necessary for exactly recovering the archetypes from noiseless data. We prove that, under uniqueness (plus additional regularity conditions on the geometry of the archetypes), our estimator is robust. While our approach requires solving a non-convex optimization problem, we find that standard optimization methods succeed in finding good solutions both for real and synthetic data.

Equations

$${\mathbf{x}}_i\approx \sum _{\ell =1}^r{w}_{i,\ell }{\mathbf{h}}_{\ell },w_{i,\ell }\ge 0,\sum _{\ell =1}^r{w}_{i,\ell }=1.$$

$$\Vert \mathbf{X}-\mathbf{W}\mathbf{H}{\Vert }_F^2,$$

$$\mathbf{W}\ge 0,\mathbf{H}\ge 0.$$

$$\text{mathscrsfs}D(\mathbf{u};\mathbf{V})$$

$${\Delta }^m$$

$$\text{mathscrsfs}D(\mathbf{U};\mathbf{V})\equiv \sum _{\ell =1}^k\text{mathscrsfs}D({\mathbf{u}}_{\ell };\mathbf{V}).$$

$$\displaystyle\widehat{\boldsymbol{H}}_{\lambda}\in\arg\min\Big{\{}\mathscrsfs{D}({\boldsymbol{X}};{\boldsymbol{H}})+\lambda\,\mathscrsfs{D}({\boldsymbol{H}};{\boldsymbol{X}}):\;\;{\boldsymbol{H}}\in{\mathcal{Q}}^{r}\Big{\}}\,,$$

$$\text{mathscrsfs}D(\mathbf{H};\mathbf{X}),$$

$$\text{mathscrsfs}D({\mathbf{x}}_i;\mathbf{H})\le {\delta }^2\text{mbox}foralli\in \{1,\dots ,n\}.$$

$$\mathbf{X}={\mathbf{W}}_0{\mathbf{H}}_0+\mathbf{Z},$$

$$\text{mathscrsfs}L({\mathbf{H}}_0,\hat{\mathbf{H}})\equiv \sum _{\ell =1}^r\underset{{\ell }^{\prime }\le r}{\min }D({\mathbf{h}}_{0,\ell },{\hat{\mathbf{h}}}_{{\ell }^{\prime }}).$$

$$\displaystyle\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{X}}_{0})^{1/2}\geq\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{X}}_{0})^{1/2}+\alpha\,\big{\{}\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{H}}_{0})^{1/2}+\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{H}})^{1/2}\big{\}}\,.$$

$$\text{mathscrsfs}D({\mathbf{h}}_i,{\mathbf{X}}_0)\le \text{mathscrsfs}D({\mathbf{h}}_i,{\mathbf{H}}_0).$$

$$\displaystyle\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{X}}_{0})^{1/2}\leq\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{X}}_{0})^{1/2}+\big{\{}\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{H}}_{0})^{1/2}+\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{H}})^{1/2}\big{\}}\,.$$

$$\underset{i\le n}{\max }\Vert {\mathbf{Z}}_{i,\cdot }{\Vert }_2\le \frac{\alpha \mu }{30r^{3/2}},$$

$$\text{mathscrsfs}L({\mathbf{H}}_0,\hat{\mathbf{H}})\le \frac{C_{\ast }^2{r}^5}{{\alpha }^2}\underset{i\le n}{\max }\Vert {\mathbf{Z}}_{i,\cdot }{\Vert }_2^2,$$

$$\text{mathscrsfs}{R}_{\lambda }(\mathbf{H})\equiv \text{mathscrsfs}D(\mathbf{X};\mathbf{H})+\lambda \text{mathscrsfs}D(\mathbf{H};\mathbf{X}),$$

$$\mathrm{minimize}\Psi (\mathbf{x},\mathbf{y})=f(\mathbf{x})+g(\mathbf{y})+h(\mathbf{x},\mathbf{y})$$

$$\text{mathscrsfs}{R}_{\lambda }(\mathbf{H})=\underset{\mathbf{W}}{\min }\Psi (\mathbf{H},\mathbf{W})=f(\mathbf{H})+g(\mathbf{W})+h(\mathbf{H},\mathbf{W}).$$

$$f(\mathbf{H})=\lambda \text{mathscrsfs}D(\mathbf{H},\mathbf{X}),$$

$$g(\mathbf{W})=\sum _{i=1}^n\mathbf{I}\left({\mathbf{w}}_i\in {\Delta }^r\right),$$

$$h(\mathbf{H},\mathbf{W})={\Vert \mathbf{X}-\mathbf{W}\mathbf{H}\Vert }_F^2.$$

$${\tilde{\mathbf{H}}}^k={\mathbf{H}}^k-\frac{1}{{\gamma }_1^k}({\mathbf{W}}^k{)}^{\mathsf{T}}\left({\mathbf{W}}^k{\mathbf{H}}^k-\mathbf{X}\right),$$

$${\mathbf{H}}^{k+1}={\tilde{\mathbf{H}}}^k-\frac{\lambda }{\lambda +{\gamma }_1^k}\left({\tilde{\mathbf{H}}}^k-{\mathbf{\Pi }}_{\mathrm{conv}(\mathbf{X})}\left({\tilde{\mathbf{H}}}^k\right)\right),$$

$${\mathbf{W}}^{k+1}={\mathbf{\Pi }}_{{\Delta }^r}\left({\mathbf{W}}^k-\frac{1}{{\gamma }_2^k}\left({\mathbf{W}}^k{\mathbf{H}}^{k+1}-\mathbf{X}\right)({\mathbf{H}}^{k+1}{)}^{\mathsf{T}}\right),$$

$$\nabla \text{mathscrsfs}{R}_{\lambda }(\mathbf{H})$$

$${\mathbf{\alpha }}_i^{\ast }$$

$${\mathbf{x}}_i=\sum _{\ell =1}^4{w}_{i,\ell }{\mathbf{h}}_{\ell }+{\mathbf{z}}_i,$$

$$\mathrm{conv}(\mathbf{H})=\{\mathbf{x}\in {\mathbb{R}}^d:\mathbf{x}={\mathbf{H}}^{\mathsf{T}}\mathbf{\pi },\mathbf{\pi }\in {\Delta }^r\},$$

$$\mathrm{aff}(\mathbf{H})=\{\mathbf{x}\in {\mathbb{R}}^d:\mathbf{x}={\mathbf{H}}^{\mathsf{T}}\mathbf{\alpha },⟨1,\mathbf{\alpha }⟩=1\}.$$

$$Q_{r,n}=\left\{\mathbf{\Pi }\in {\mathbb{R}}_{\ge 0}^{r\times n}:⟨{\mathbf{\Pi }}_{i,.},1⟩=1\right\}.$$

$$S_r=\left\{\mathbf{\Pi }\in Q_r:{\Pi }_{i,j}\in \{0,1\}\right\}.$$

$$\text{mathscrsfs}D({\mathbf{H}}_1,\mathbf{X})$$

$$\text{mathscrsfs}L({\mathbf{H}}_1,{\mathbf{H}}_2)$$

$${\mathsf{B}}_r(\rho ;\mathbf{H}):=\{\mathbf{x}\in {\mathbb{R}}^d:\text{mathscrsfs}D(\mathbf{x},\mathbf{H})\le {\rho }^2\}.$$

Full text

In this paper, we study an approach to NMF that can be traced back to the work of Cutler and Breiman [CB94] and does not require the data to be separable, while providing a generally unique decomposition. We optimize the trade-off between two objectives: we minimize the distance of the data points from the convex envelope of the archetypes (which can be interpreted as an empirical risk), while minimizing the distance of the archetypes from the convex envelope of the data (which can be interpreted as a data-dependent regularization). The archetypal analysis method of [CB94] is recovered as the limiting case in which the last term is given infinite weight.

We introduce a ‘uniqueness condition’ on the data which is necessary for exactly recovering the archetypes from noiseless data. We prove that, under uniqueness (plus additional regularity conditions on the geometry of the archetypes), our estimator is robust. While our approach requires solving a non-convex optimization problem, we find that standard optimization methods succeed in finding good solutions both for real and synthetic data.

1 Introduction

Given a set of data points ${\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2},\cdots,{\boldsymbol{x}}_{n}\in\mathbb{R}^{d}$, it is often useful to represent them as convex combinations of a small set of vectors (the ‘archetypes’ ${\boldsymbol{h}}_{1},\dots,{\boldsymbol{h}}_{\ell}$):

[TABLE]

Decompositions of this type have wide ranging applications, from chemometrics [PT94] to image processing [LS99] and topic modeling [XLG03]. As an example, Figure 1 displays the infrared reflection spectra111Data were retrieved from the NIST Chemistry WebBook dataset [LM]. of four molecules (caffeine, sucrose, lactose and trioctanoin) for wavenumber between $1186\;{\rm cm}^{-1}$ and $1530\;{\rm cm}^{-1}$. Each spectrum is a vector ${\boldsymbol{h}}_{0,\ell}\in\mathbb{R}^{d}$, with $d=87$ and $\ell\in\{1,\dots,4\}$. If a mixture of these substances is analyzed, the resulting spectrum will be a convex combination of the spectra of the four analytes. This situation arises in hyperspectral imaging [MBDC*+*14], where a main focus is to estimate spatially varying proportions of a certain number of analytes. In order to mimic this setting, we generated $n=250$ synthetic random convex combinations ${\boldsymbol{x}}_{1},\dots,{\boldsymbol{x}}_{n}\in\mathbb{R}^{d}$ of the four spectra ${\boldsymbol{h}}_{0,1}$, …, ${\boldsymbol{h}}_{0,4}$, each containing two or more of these four analytes, and tried to reconstruct the pure spectra from the ${\boldsymbol{x}}_{i}$’s. Each column in Figure 1 displays the reconstruction obtained using a different procedure. We refer to Appendix A for further details.

Without further constraints, the decomposition (1.1) is dramatically underdetermined. Given a set of valid archetypes $\{{\boldsymbol{h}}_{0,\ell}\}_{\ell\leq r}$, any set $\{{\boldsymbol{h}}_{\ell}\}_{\ell\leq r}$ whose convex hull contains the $\{{\boldsymbol{h}}_{0,\ell}\}_{\ell\leq r}$ also satisfies Eq. (1.1). For instance, we can set ${\boldsymbol{h}}_{\ell}={\boldsymbol{h}}_{0,\ell}$ for $\ell\leq r-1$, and ${\boldsymbol{h}}_{r}=(1+s){\boldsymbol{h}}_{0,r}-s{\boldsymbol{h}}_{0,1}$ for any $s\geq 0$, and obtain an equally good representation of the data $\{{\boldsymbol{x}}_{i}\}_{i\leq n}$.

How should we constrain the decomposition (1.1) in such a way that it is generally unique (up to permutations of the $r$ archetypes)? Since the seminal work of Paatero and Tapper [PT94, Paa97], and of Lee and Seung [LS99, LS01], an overwhelming amount of work has addressed this question by making the assumptions that the archetypes are componentwise non-negative ${\boldsymbol{h}}_{\ell}\geq 0$. Among other applications the non-negativity constraint is justified for chemometrics (reflection or absorption spectra are non-negative), and topic modeling (in this case archetypes correspond to topics, which are represented as probability distributions over words). This formulation has become popular as non-negative matrix factorization (NMF).

Under the non-negativity constraint ${\boldsymbol{h}}_{\ell}\geq 0$ the role of weights and archetypes becomes symmetric, and the decomposition (1.1) is unique provided that the archetypes or the weights are sufficiently sparse (without loss of generality one can assume $\sum_{\ell=1}^{r}h_{\ell,i}=1$). This point was clarified by Donoho and Stodden [DS03], introduced a separability condition that ensure uniqueness. The non-negative archetypes ${\boldsymbol{h}}_{1},\cdots,{\boldsymbol{h}}_{r}$ are separable if, for each $\ell\in[r]$ there exists an index $i(\ell)\in[d]$ such that $({\boldsymbol{h}}_{\ell})_{i(\ell)}=1$, and $({\boldsymbol{h}}_{\ell^{\prime}})_{i(\ell)}=0$ for all $\ell^{\prime}\neq\ell$. If we exchange the roles of weights $\{w_{i,\ell}\}$ and archetypes $\{h_{\ell,i}\}$, separability requires that $\ell\in[r]$ there exists an index $i(\ell)\in[n]$ such that $w_{i(\ell),\ell}=1$, and $w_{i(\ell),\ell^{\prime}}=0$ for all $\ell^{\prime}\neq\ell$ This condition has a simple geometric interpretation: the data are separable if for each archetype ${\boldsymbol{h}}_{\ell}$ there is at least one data point ${\boldsymbol{x}}_{i}$ such that ${\boldsymbol{x}}_{i}={\boldsymbol{h}}_{\ell}$. A copious literature has developed algorithms for non-negative matrix factorization under separability condition or its generalizations [DS03, AGKM12, RRTB12, AGH*+*13, GZ15].

Of course this line of work has a drawback: in practice we do not know whether the data are separable. (We refer to the Section 5 for a comparison with [GZ15], which relaxes the separability assumption.) Further, there are many cases in which the archetypes ${\boldsymbol{h}}_{1},\dots,{\boldsymbol{h}}_{\ell}$ are not necessarily non-negative. For instance, in spike sorting, the data are measurements of neural activity ad the archetypes correspond to waveforms associated to different neurons [RCP09]. In other applications the archetypes ${\boldsymbol{h}}_{\ell}$ are non-negative, but –in order to reduce complexity– the data $\{{\boldsymbol{x}}_{i}\}_{i\leq n}$ are replaced by a random low-dimensional projection [KSD08, WL10]. The projected archetypes loose the non-negativity property. Finally, the decomposition (1.1) is generally non-unique, even under the constraint ${\boldsymbol{h}}_{\ell}\geq 0$. This is illustrated, again, in Figure 1: all the spectra are strictly positive, and hence we can find archetypes ${\boldsymbol{h}}_{1},\dots,{\boldsymbol{h}}_{4}$ that are still non-negative and whose convex envelope contains ${\boldsymbol{h}}_{0,1},\dots,{\boldsymbol{h}}_{0,4}$.

Since NMF is underdetermined, standard methods fail in such applications, as illustrated in Figure 1 . We represent the data as a matrix ${\boldsymbol{X}}\in\mathbb{R}^{n\times d}$ whose $i$-th row is the vector ${\boldsymbol{x}}_{i}$, the weights by a matrix ${\boldsymbol{W}}=(w_{i,\ell})_{i\leq n,\ell\leq r}\in\mathbb{R}^{n\times d}$ and the prototypes by a matrix ${\boldsymbol{H}}=(h_{\ell,j})_{\ell\leq r,j\leq d}\in\mathbb{R}^{r\times d}$. The third column of Figure 1 uses a projected gradient algorithm from [Lin07] to solve the problem

[TABLE]

Empirically, projected gradient converges to a point with very small fitting error $\|{\boldsymbol{X}}-{\boldsymbol{W}}{\boldsymbol{H}}\|_{F}^{2}$, but the reconstructed spectra (rows of ${\boldsymbol{H}}$) are inaccurate. The second column in the same figure shows the spectra reconstructed using an algorithm from [AGH*+*13], that assumes separability: as expected, the reconstruction is not accurate.

In a less widely known paper, Cutler and Breiman [CB94] addressed the same problem using what they call ‘archetypal analysis.’ Archetypal analysis presents two important differences with respect to standard NMF: $(1)$ The archetypes ${\boldsymbol{h}}_{\ell}$ are not necessarily required to be non-negative (although this constraint can be easily incorporated); $(2)$ The under-determination of the decomposition (1.1) is addressed by requiring that the archetypes belong to the convex hull of the data points: ${\boldsymbol{h}}_{\ell}\in{\rm{conv}}(\{{\boldsymbol{x}}_{i}\}_{i\leq n})$.

In applications the condition ${\boldsymbol{h}}_{\ell}\in{\rm{conv}}(\{{\boldsymbol{x}}_{i}\}_{i\leq n})$ is too strict. This paper builds on the ideas of [CB94] to propose a formulation of NMF that is uniquely defined (barring degenerate cases) and provides a useful notion of optimality. In particular, we present the following contributions.

Archetypal reconstruction. We propose to reconstruct the archetypes ${\boldsymbol{h}}_{1},\dots,{\boldsymbol{h}}_{r}$ by optimizing a combination of two objectives. On one hand, we minimize the error in the decomposition (1.1). This amounts to minimizing the distance between the data points and the convex hull of the archetypes. On the other hand, we minimize the distance of the archetypes from the convex hull of data points. This relaxes the original condition imposed in [CB94] which required the archetypes to lie in ${\rm{conv}}(\{{\boldsymbol{x}}_{i}\})$, and allows to treat non-separable data.

Robustness guarantee. We next assume that that the decomposition (1.1) approximately hold for some ‘true’ archetypes ${\boldsymbol{h}}^{0}_{\ell}$ and weights $w_{i,\ell}^{0}$, namely ${\boldsymbol{x}}_{i}={\boldsymbol{x}}_{i}^{0}+{\boldsymbol{z}}_{i}$, where ${\boldsymbol{x}}^{0}_{i}=\sum_{\ell=1}^{r}w^{0}_{i,\ell}{\boldsymbol{h}}^{0}_{\ell}$ and ${\boldsymbol{z}}_{i}$ captures unexplained effects. We introduce a ‘uniqueness condition’ on the data $\{{\boldsymbol{x}}^{0}_{i}\}_{i\leq n}$ which is necessary for exactly recovering the archetypes from the noiseless data. We prove that, under uniqueness (plus additional regularity conditions on the geometry of the archetypes), our estimator is robust. Namely it outputs archetypes $\{\hat{\boldsymbol{h}}_{\ell}\}_{\ell\leq r}$ whose distance from the true ones $\{{\boldsymbol{h}}^{0}_{\ell}\}_{\ell\leq r}$ (in a suitable metric) is controlled by $\sup_{i\leq n}\|{\boldsymbol{z}}_{i}\|_{2}$.

Algorithms. Our approach reconstructs the archetypes ${\boldsymbol{h}}_{1},\dots,{\boldsymbol{h}}_{r}$ by minimizing a non-convex risk function $\mathscrsfs{R}_{\lambda}({\boldsymbol{H}})$. We propose three descent algorithms that appear to perform well on realistic instances of the problem. In particular, Section 4 introduces a proximal alternating linearized minimization algorithm (PALM) that is guaranteed to converge to critical points of the risk function. Appendix E discusses two alternative approaches. One possible explanation for the success of such descent algorithms is that reasonably good initializations can be constructed using spectral methods, or approximating the data as separable, cf. Section 4.1. We defer a study of global convergence of this two-stages approach to future work.

2 An archetypal reconstruction approach

Let ${\mathcal{Q}}\subseteq\mathbb{R}^{d}$ be a convex set and $D:{\mathcal{Q}}\times{\mathcal{Q}}\to\mathbb{R}$, $({\boldsymbol{x}},{\boldsymbol{y}})\mapsto D({\boldsymbol{x}};{\boldsymbol{y}})$ a loss function on ${\mathcal{Q}}$. For a point ${\boldsymbol{u}}\in{\mathcal{Q}}$, and a matrix ${\boldsymbol{V}}\in\mathbb{R}^{m\times d}$, with rows ${\boldsymbol{v}}_{1},\dots,{\boldsymbol{v}}_{m}\in{\mathcal{Q}}$, we let

[TABLE]

In other words, denoting by ${\rm{conv}}({\boldsymbol{V}})={\rm{conv}}(\{{\boldsymbol{v}}_{1},\dots,{\boldsymbol{v}}_{m}\})$ the convex hull of the rows of matrix ${\boldsymbol{V}}$, $\mathscrsfs{D}({\boldsymbol{u}};{\boldsymbol{V}})$ is the minimum loss between ${\boldsymbol{x}}$ and any point in ${\rm{conv}}({\boldsymbol{V}})$. If ${\boldsymbol{U}}\in\mathbb{R}^{k\times d}$ is a matrix with rows ${\boldsymbol{u}}_{1},\dots,{\boldsymbol{u}}_{k}\in{\mathcal{Q}}$, we generalize this definition by letting

[TABLE]

While this definition makes sense more generally, we have in mind two specific examples in which $D({\boldsymbol{x}};{\boldsymbol{y}})$ is actually separately convex in its arguments ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$. (Most of our results will concern the first example.)

Example 2.1** (Square loss).**

In this case ${\mathcal{Q}}=\mathbb{R}^{d}$, and $D({\boldsymbol{x}};{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$. This is the case originally studied by Cutler and Breiman [CB94].

Example 2.2** (KL divergence).**

We take ${\mathcal{Q}}=\Delta^{d}$, the $d$-dimensional simplex, and $D({\boldsymbol{x}};{\boldsymbol{y}})$ to be the Kullback-Leibler divergence between probability distributions ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$, namely $D({\boldsymbol{x}};{\boldsymbol{y}})\equiv\sum_{i=1}^{d}x_{i}\log(x_{i}/y_{i})$.

Given data ${\boldsymbol{x}}_{1},\dots,{\boldsymbol{x}}_{n}$ organized in the matrix ${\boldsymbol{X}}\in\mathbb{R}^{n\times d}$, we estimate the archetypes by solving the problem222This problem can have multiple global minima if $\lambda=0$ or in degenerate settings. One minimizer is selected arbitrarily when this happens.

[TABLE]

where we denote by ${\mathcal{Q}}^{r}$ the set of matrices ${\boldsymbol{H}}\in\mathbb{R}^{r\times d}$ with rows ${\boldsymbol{h}}_{1},\dots,{\boldsymbol{h}}_{r}\in{\mathcal{Q}}$. A few values of $\lambda$ are of special significance. If we set $\lambda=0$, and ${\mathcal{Q}}=\Delta^{d}$, we recover the standard NMF objective (1.2), with a more general distance function $D(\,\cdot\,,\,\cdot\,)$. As pointed out above, in general this objective has no unique minimum. If we let $\lambda\to 0+$ after the minimum is evaluated, $\widehat{\boldsymbol{H}}_{\lambda}$ converges to the minimizer of $\mathscrsfs{D}({\boldsymbol{X}};{\boldsymbol{H}})$ which is the ‘closest’ to the convex envelope of the data ${\rm{conv}}({\boldsymbol{X}})$ (in the sense of minimizing $\mathscrsfs{D}({\boldsymbol{H}};{\boldsymbol{X}})$). Finally as $\lambda\to\infty$, the archetypes ${\boldsymbol{h}}_{\ell}$ are forced to lie in ${\rm{conv}}({\boldsymbol{X}})$ and hence we recover the method of [CB94].

Figure 2 illustrates the advantages of the estimator (2.4) on a small synthetic example, with $d=2$, $r=3$, $n=500$: in this case the data are non separable. We first use the successive projections algorithm of [AGH*+*13] (that is designed to deal with separable data) in order to estimate the archetypes. As expected, the reconstruction is not accurate because this algorithm assumes separability and hence estimates the archetypes with a subset of the data points. We then use these estimates as initialization in the alternate minimization algorithm of [CB94], which optimizes the objective (2.4) with $\lambda=\infty$. The estimates improve but not substantially: they are still constrained to lie in ${\rm{conv}}({\boldsymbol{X}})$. A significant improvement is obtained by setting $\lambda$ to a small value. We (approximately) minimize the cost function (2.4) by generalizing the alternate minimization algorithm, cf. Section 4. The optimal archetypes are no longer constrained to ${\rm{conv}}({\boldsymbol{X}})$, and provide a better estimate of the true archetypes. The last column in Figure 1 uses the same estimator, and approximately solves problem (2.4) by gradient descent algorithm.

In our analysis we will consider a slightly different formulation in which the Lagrangian of Eq. (2.4) is replaced by a hard constraint:

[TABLE]

We will use this version in the analysis presented in the next section, and denote the corresponding estimator by $\widehat{\boldsymbol{H}}$.

3 Robustness

In order to analyze the robustness properties of estimator $\widehat{\boldsymbol{H}}$, we assume that there exists an approximate factorization

[TABLE]

where ${\boldsymbol{W}}_{0}\in\mathbb{R}^{n\times r}$ is a matrix of weights (with rows ${\boldsymbol{w}}_{0,i}\in\Delta^{r}$), ${\boldsymbol{H}}_{0}\in\mathbb{R}^{r\times d}$ is a matrix of archetypes (with rows ${\boldsymbol{h}}_{0,\ell}$), and we set ${\boldsymbol{X}}_{0}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}$. The deviation ${\boldsymbol{Z}}$ is arbitrary, with rows ${\boldsymbol{z}}_{i}$ satisfying $\max_{i\leq n}\|{\boldsymbol{z}}_{i}\|_{2}\leq\delta$. We will assume throughout $r$ to be known.

We will quantify estimation error by the sum of distances between the true archetypes and the closest estimated archetypes

[TABLE]

In words, if $\mathscrsfs{L}({\boldsymbol{H}}_{0},\widehat{\boldsymbol{H}})$ is small, then for each true archetype ${\boldsymbol{h}}_{0,\ell}$ there exists an estimated archetype $\hat{\boldsymbol{h}}_{\ell^{\prime}}$ that is close to it in $D$-loss. Unless two or more of the true archetypes are close to each other, this means that there is a one-to-one correspondence between estimated archetypes and true archetypes, with small errors.

Assumption** (Uniqueness).**

We say that the factorization ${\boldsymbol{X}}_{0}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}$ satisfies uniqueness with parameter $\alpha>0$ (equivalently, is $\alpha$-unique) if for all ${\boldsymbol{H}}\in{\mathcal{Q}}^{r}$ with ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{conv}}({\boldsymbol{H}})$, we have

[TABLE]

The rationale for this assumption is quite clear. Assume that the data lie in the convex hull of the true archetypes ${\boldsymbol{H}}_{0}$, and hence Eq. (3.1) holds without error term ${\boldsymbol{Z}}=0$, i.e. ${\boldsymbol{X}}={\boldsymbol{X}}_{0}$. We reconstruct the archetypes by demanding ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{conv}}({\boldsymbol{H}})$: any such ${\boldsymbol{H}}$ is a plausible explanation of the data. In order to make the problem well specified, we define ${\boldsymbol{H}}_{0}$ to be the matrix of archetypes that are the closest to ${\boldsymbol{X}}_{0}$, and hence $\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{X}}_{0})\geq\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{X}}_{0})$ for all ${\boldsymbol{H}}$. In order for the reconstruction to be unique (and hence for the problem to be identifiable) we need to assume $\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{X}}_{0})>\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{X}}_{0})$ strictly for ${\boldsymbol{H}}\neq{\boldsymbol{H}}_{0}$. The uniqueness assumption provides a quantitative version of this condition.

Remark 3.1**.**

Given ${\boldsymbol{X}}_{0}$, ${\boldsymbol{H}}_{0}$, the best constant $\alpha$ such that Eq. (3.3) holds for all ${\boldsymbol{H}}$ is a geometric property that depend on ${\boldsymbol{X}}_{0}$ only through ${\rm{conv}}({\boldsymbol{X}}_{0})$. In particular, if ${\boldsymbol{X}}_{0}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}$ is a separable factorization, then it satisfies uniqueness with parameter $\alpha=1$. Indeed in this case ${\rm{conv}}({\boldsymbol{H}}_{0})={\rm{conv}}({\boldsymbol{X}}_{0})$, whence $\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{X}}_{0})=\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{H}}_{0})$ and $\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{X}}_{0})=\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{H}})=0$.

It is further possible to show that $\alpha\in[0,1]$ for all ${\boldsymbol{H}}_{0},{\boldsymbol{X}}_{0}$. Indeed, we took ${\boldsymbol{H}}_{0}$ to be the matrix of archetypes that are closest to ${\boldsymbol{X}}_{0}$. In other words, $\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{X}}_{0})\geq\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{X}}_{0})$ and hence, $\alpha\geq 0$. In addition, since ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{conv}}({\boldsymbol{H}}_{0})$, for ${\boldsymbol{h}}_{i}$ an arbitrary row of ${\boldsymbol{H}}$ we have

[TABLE]

Hence, $\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{X}}_{0})\leq\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{H}}_{0})$ and therefore

[TABLE]

Thus, $\alpha\leq 1$.

We say that the convex hull ${\rm{conv}}({\boldsymbol{X}}_{0})$ has internal radius (at least) $\mu$ if it contains an $r-1$-dimensional ball of radius $\mu$, i.e. if there exists ${\boldsymbol{z}}_{0}\in\mathbb{R}^{d}$, ${\boldsymbol{U}}\in\mathbb{R}^{d\times(r-1)}$, with ${\boldsymbol{U}}^{\sf T}{\boldsymbol{U}}={\bf I}_{d}$ , such that ${\boldsymbol{z}}_{0}+{\boldsymbol{U}}{\sf B}_{r-1}(\mu)\subseteq{\rm{conv}}({\boldsymbol{X}}_{0})$. We further denote by $\kappa({\boldsymbol{M}})$ the condition number of matrix ${\boldsymbol{M}}$.

Theorem 1**.**

Assume ${\boldsymbol{X}}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}+{\boldsymbol{Z}}$ where the factorization ${\boldsymbol{X}}_{0}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}$ satisfies the uniqueness assumption with parameter $\alpha>0$, and that ${\rm{conv}}({\boldsymbol{X}}_{0})$ has internal radius $\mu>0$. Consider the estimator $\widehat{\boldsymbol{H}}$ defined by Eq. (2.5), with $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$ (square loss) and $\delta=\max_{i\leq n}\|{\boldsymbol{Z}}_{i,\cdot}\|_{2}$. If

[TABLE]

then, we have

[TABLE]

where $C_{*}$ is a coefficient that depends uniquely on the geometry of ${\boldsymbol{H}}_{0}$, ${\boldsymbol{X}}_{0}$, namely $C_{*}=120(\sigma_{\rm max}({\boldsymbol{H}}_{0})/\mu)\cdot\max(1,\kappa({\boldsymbol{H}}_{0})/\sqrt{r})$.

4 Algorithms

While our main focus is on structural properties of non-negative matrix factorization, we provide evidence that the optimization problem we defined can be solved in practical scenarios. A more detailed study is left to future work.

From a computational point of view, the Lagrangian formulation (2.4) is more appealing. For the sake of simplicity, we denote the regularized risk by

[TABLE]

and leave implicit the dependence on the data ${\boldsymbol{X}}$. Notice that this function is non-convex and indeed has multiple global minima: in particular, permuting the rows of a minimizer ${\boldsymbol{H}}$ yields other minimizers. We will describe two greedy optimization algorithms: one based on gradient descent, and one on alternating minimization, which generalizes the algorithm of [CB94]. In both cases it is helpful to use a good initialization: two initialization methods are introduced in the next section.

4.1 Initialization

We experimented with two initialization methods, described below.

*(*1) Spectral initialization. Under the assumption that the archetypes $\{{\boldsymbol{h}}_{0,\ell}\}_{\ell\leq r}$ are linearly independent (and for non-degenerate weights ${\boldsymbol{W}}$), the ‘noiseless’ matrix ${\boldsymbol{X}}_{0}$ has rank exactly $r$. This motivates the following approach. We compute the singular value decomposition ${\boldsymbol{X}}=\sum_{i=1}^{n\wedge d}\sigma_{i}{\boldsymbol{u}}_{i}{\boldsymbol{v}}_{i}^{{\sf T}}$, $\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{n\wedge d}$, and initialize $\widehat{\boldsymbol{H}}$ as the matrix $\widehat{\boldsymbol{H}}^{(0)}$ with rows $\hat{\boldsymbol{h}}^{(0)}_{1}={\boldsymbol{v}}_{1},\dots,\hat{\boldsymbol{h}}^{(0)}_{r}={\boldsymbol{v}}_{r}$.

*(*2) Successive projections initialization. We initialize $\widehat{\boldsymbol{H}}^{(0)}$ by choosing archetypes $\{\hat{\boldsymbol{h}}^{(0)}_{\ell}\}_{1\leq\ell\leq r}$ that are a subset of the data $\{{\boldsymbol{x}}_{i}\}_{1\leq i\leq n}$, selected as follows. The first archetype $\hat{\boldsymbol{h}}^{(0)}_{1}$ is the data point which is farthest from the origin. For each subsequent archetype, we choose the point that is farthest from the affine subspace spanned by the previous ones.

[TABLE]

This coincides with the successive projections algorithm of [ASG*+*01], with the minor difference that $V_{\ell}$ is the affine subspace spanned by the first $\ell$ vectors, instead of the linear subspace333The same modification is also used in [AGH*+*13], but we do not apply the full algorithm of this paper. This method can be proved to return the exact archetypes if data are separable the archetypes are affine independent [AGH*+*13, GV14]. When data are not separable it provides nevertheless a good initial assignment.

4.2 Proximal alternating linearized minimization

The authors of [BST14] develop a proximal alternating linearized minimization algorithm (PALM) to solve the problems of the form

[TABLE]

where $f:\mathbb{R}^{m}\to(-\infty,+\infty]$ and $g:\mathbb{R}^{n}\to(-\infty,\infty]$ are lower semicontinuous and $h\in C^{1}(\mathbb{R}^{m}\times\mathbb{R}^{n})$. PALM is guaranteed to converge to critical points of the function $\Psi$ [BST14].

We apply this algorithm to minimize the cost function (4.1), with $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$ which we write as

[TABLE]

where,

[TABLE]

In above equations ${\boldsymbol{w}}_{i}$ are the rows of ${\boldsymbol{W}}$ and the indicator function ${\rm\bf I}({\boldsymbol{x}}\in\Delta^{r})$ is equal to zero if ${\boldsymbol{x}}\in\Delta^{r}$ and is equal to infinity otherwise.

By using this decomposition, the iterations of the PALM iteration reads

[TABLE]

where $\gamma_{1}^{k}$, $\gamma_{2}^{k}$ are step sizes and, for ${\boldsymbol{M}}\in\mathbb{R}^{m_{1}\times m_{2}}$, and ${\cal S}\subseteq\mathbb{R}^{m_{2}}$ a closed convex set, ${\boldsymbol{\Pi}}_{{\cal S}}({\boldsymbol{M}})$ is the matrix obtained by projecting the rows of ${\boldsymbol{M}}$ onto the simplex ${\cal S}$.

Proposition 4.1**.**

Consider the risk (4.1), with loss $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$, and the corresponding cost function $\Psi({\boldsymbol{H}},{\boldsymbol{W}})$. If the step sizes are chosen such that $\gamma_{1}^{k}>\left\|{\boldsymbol{W}}^{k^{\sf T}}{\boldsymbol{W}}^{k}\right\|_{F}$, $\gamma_{2}^{k}>\max\left\{\left\|{\boldsymbol{H}}^{k+1}{\boldsymbol{H}}^{k+1^{\sf T}}\right\|_{F},{\varepsilon}\right\}$ for some constant ${\varepsilon}>0$, then $({\boldsymbol{H}}^{k},{\boldsymbol{W}}^{k})$ converges to a stationary point of the function $\Psi({\boldsymbol{H}},{\boldsymbol{W}})$.

The proof of this statement is deferred to Appendix D.

It is also useful to notice that the gradient of $\mathscrsfs{R}_{\lambda}({\boldsymbol{H}})$ can be computed explicitly (this can be useful to devise a stopping criterion).

Proposition 4.2**.**

Consider the risk (4.1), with loss $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$, and assume that the rows of ${\boldsymbol{H}}$ are affine independent. Then, $\mathscrsfs{R}_{\lambda}$ is differentiable at ${\boldsymbol{H}}$ with gradient

[TABLE]

where we recall that ${\boldsymbol{\Pi}}_{{\rm{conv}}({\boldsymbol{X}})}({\boldsymbol{H}})$ denotes the matrix with rows ${\boldsymbol{\Pi}}_{{\rm{conv}}({\boldsymbol{X}})}({\boldsymbol{H}}_{1,\cdot}),\dots,{\boldsymbol{\Pi}}_{{\rm{conv}}({\boldsymbol{X}})}({\boldsymbol{H}}_{r,\cdot})$.

The proof of this proposition is given in Appendix C. Appendix E also discusses two alternative algorithms.

4.3 Numerical experiments

We implemented both the PALM algorithm described in the previous section, and the two algorithms described in Appendix E. The outcomes are generally similar.

Figures 3 and 4 repeat the experiment already described in the introduction. We generate $n=250$ convex combinations of $r=4$ spectra ${\boldsymbol{h}}_{0,1},\dots,{\boldsymbol{h}}_{0,4}\in\mathbb{R}^{d}$, $d=87$, this time adding white Gaussian noise with variance $\sigma^{2}$. We minimize the Lagrangian $\mathscrsfs{R}_{\lambda}({\boldsymbol{H}})$, with444These values were chosen as to approximately minimize the estimation error. $\lambda=4$ (for Figure 3) and $\lambda=0.8$ (for Figure 4). The reconstructed spectra of the pure analytes appear to be accurate and robust to noise.

In Figure 5 we repeated the same experiment systematically for $10$ noise realizations for each noise level $\sigma$, and report the resulting average loss. Among various reconstruction methods, the approach described in this paper seem to have good robustness to noise and achieves exact reconstruction as $\sigma\to 0$.

5 Discussion

We introduced a new optimization formulation of the non-negative matrix factorization problem. In its Lagrangian formulation, our approach consists in minimizing the cost function $\mathscrsfs{R}_{\lambda}({\boldsymbol{H}})$ defined in Eq. (4.1). This encompasses applications in which only one of the factors is required to be non-negative. A special case of this formulation ($\lambda\to\infty$) corresponds to the ‘archetypal analysis’ of [CB94]. In this case, the archetype estimates coincide with a subset of the data points, which is appropriate only under the separability assumption of [DS03].

Our main technical result (Theorem 1) is a robustness guarantee for the reconstructed archetypes. This holds under the uniqueness assumption, which appears to hold for generic geometries of the dataset. In particular, while separability implies uniqueness (with optimal constant $\alpha=1$), uniqueness holds for non-separable data as well. To the best of our knowledge, similar robustness results have been obtained under separability [RRTB12, AGH*+*13, GV14, GV15] (albeit these works obtain a better dependence on $r$). The only exception is the recent work of Ge and Zou [GZ15] who prove robustness under a ‘subset separability’ condition, which provides a significant relaxation of separability. Under this condition, [GZ15] develops a polynomial-time algorithm to estimate the archetypes by identifying and intersecting the faces of ${\rm{conv}}({\boldsymbol{H}}_{0})$. However, the algorithm of [GZ15] exploits collinearities to identify the faces, and this requires additional ‘genericity’ assumptions.

Admittedly, the uniqueness constant $\alpha$ is difficult to evaluate analytically, even for simple geometries of the data. However, by definition it does not vanish except in the case of multiple minimizers, and we expect it typically to be of order one. Figure 6 illustrate this point by computing numerically $\alpha$ for a simple one-parameter family of geometries with $r=3$, $d=2$. The parameter $\alpha$ vanishes at a single point, corresponding to a degenerate problem with multiple solutions.

Finally, several earlier works addressed the non-uniqueness problem in classical non-negative matrix factorization. Among others, Miao and Qi [MQ07] penalize a matrix of archetypes ${\boldsymbol{H}}$ by the corresponding volume. Closely related to our work is the approach of Mørup and Hansen [MH12] tha also builds on archetypal analysis. To the best of our knowledge, none of these works establishes robustness of the proposed methods.

We conclude by mentioning three important problems that are not addressed by this paper: $(1)$ Are there natural condition under which the risk function $\mathscrsfs{R}_{\lambda}({\boldsymbol{H}})$ of Eq. (4.1) can be optimized in polynomial time? We only provided an algorithm that is guaranteed to converge to a critical point. $(2)$ We assumed the rank $r$ to be known. In practice it will need to be estimated from the data. $(3)$ Similarly, the regularization parameter $\lambda$ should be chosen from data.

Acknowledgements

This work was partially supported by the NSF grant CCF-1319979 and a Stanford Graduate Fellowship.

Appendix A Further details on numerical experiments

The data in Figures 1, 3, and 4 were generated as follows. We retrieved infrared reflection spectra of caffeine, sucrose, lactose and trioctanoin from the NIST Chemistry WebBook dataset [LM]. We restricted these spectra to the wavenumbers between $1186\;{\rm cm}^{-1}$ and $1530\;{\rm cm}^{-1}$, and denote by ${\boldsymbol{h}}_{0,1},\dots,{\boldsymbol{h}}_{0,4}\in\mathbb{R}^{d}$, $d=87$ the vector representations of these spectra. We then generates data ${\boldsymbol{x}}_{i}\in\mathbb{R}^{d}$, $i\leq n=250$ by letting

[TABLE]

where ${\boldsymbol{z}}_{i}\sim{\sf N}(0,\sigma^{2}{\rm I}_{d})$ are i.i.d. Gaussian noise vectors. The weights ${\boldsymbol{w}}_{i}=(w_{i,\ell})_{\ell\leq 4}$ were generated as follows. The weight vectors $\{{\boldsymbol{w}}_{i}\}_{1\leq i\leq 9}$ are generated such that they have $2$ nonzero entries. In other words, $9$ data points are on one dimensional facets of the polytope generated by ${\boldsymbol{h}}_{0,1},\dots,{\boldsymbol{h}}_{0,4}$. In order to randomly generate these weight vectors, for each $1\leq i\leq 9$, a pair of indices $(\ell_{1},\ell_{2})$ between $1$ and $4$ is chosen uniformly at random. Then $\{\widetilde{\boldsymbol{w}}\}_{1\leq i\leq 9}$, $\widetilde{\boldsymbol{w}}\in\mathbb{R}^{2}$ are generated as independent Dirichlet random vectors with parameter $(5,5)$. Then we let $w_{i,\ell_{1}}=\tilde{w}_{i,1}$ and $w_{i,\ell_{2}}=\tilde{w}_{i,2}$ for $1\leq i\leq 9$. The weight vectors $\{{\boldsymbol{w}}_{i}\}_{10\leq i\leq 20}$ each have $3$ nonzero entries. Similar to above, for each of these weight vectors a $3$-tuple of indices $(\ell_{1},\ell_{2},\ell_{3})$ between $1$ and $4$ is chosen uniformly at random. Then we let $w_{i,\ell_{1}}=\tilde{w}_{i,1}$, $w_{i,\ell_{2}}=\tilde{w}_{i,2}$, $w_{i,\ell_{3}}=\tilde{w}_{i,4}$ for $10\leq i\leq 20$, where $\{\widetilde{\boldsymbol{w}}\}_{10\leq i\leq 20}$, $\widetilde{\boldsymbol{w}}\in\mathbb{R}^{3}$ are i.i.d. Dirichlet random vectors with parameter $(5,5,5)$. The rest of the weight vectors have cardinality equal to $4$. Hence, for $21\leq i\leq 250$, ${\boldsymbol{w}}_{i}$ are generated as i.i.d. Dirichlet random vectors with parameter $(5,5,5,5)$.

Appendix B Proof of Theorem 1

In this appendix we prove Theorem 1. We start by recalling some notations already defined in the main text, and introducing some new ones. We will then state a stronger form of the theorem (with better dependence on the problem geometry in some regimes). Finally, we will present the actual proof.

Throughout this appendix, we assume the square loss $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$.

B.1 Notations and definitions

We use bold capital letters (e.g. ${\boldsymbol{A}}$, ${\boldsymbol{B}}$, ${\boldsymbol{C}}$,…) for matrices, bold lower case for vectors (e.g. ${\boldsymbol{x}}$, ${\boldsymbol{y}}$, …) and plain lower case for scalars ($aa$, $b$, $c$ and so on). In particular, ${\boldsymbol{e}}_{i}\in\mathbb{R}^{d}$ denotes the $i$’th vector in the canonical basis, $E^{r,d}=\{{\boldsymbol{e}}_{1},{\boldsymbol{e}}_{2},\dots,{\boldsymbol{e}}_{r}\}$ and for $r\leq d$, ${\boldsymbol{E}}_{r,d}\in\{0,1\}^{r\times d}$ is the matrix whose $i$’th column is ${\boldsymbol{e}}_{i}$, and whose columns after the $r$-th one are equal to ${\boldsymbol{0}}$. For a matrix ${\boldsymbol{X}}$, ${\boldsymbol{X}}_{i,.}$ and ${\boldsymbol{X}}_{.,i}$ are its $i$’th row and column, respectively.

As in the main text, we denote by $\Delta^{m}$ the $m$-dimensional standard simplex, i.e. $\Delta^{m}=\{{\boldsymbol{x}}\in\mathbb{R}_{\geq 0}^{m},\langle{\boldsymbol{x}},{\bf 1}\rangle=1\}$, where ${\bf 1}\in\mathbb{R}^{m}$ is the all ones vector. For a matrix ${\boldsymbol{H}}\in\mathbb{R}^{r\times d}$, we use $\sigma_{\max}({\boldsymbol{H}})$, $\sigma_{\min}({\boldsymbol{H}})$ to denote its largest and smallest nonzero singular values and $\kappa({\boldsymbol{H}})=\sigma_{\max}({\boldsymbol{H}})/\sigma_{\min}({\boldsymbol{H}})$ to denote its condition number. We denote by ${\rm{conv}}({\boldsymbol{H}}),{\rm{aff}}({\boldsymbol{H}})$ the convex hull and the affine hull of the rows of ${\boldsymbol{H}}$, respectively. In other words,

[TABLE]

We denote by $Q_{r,n}$ is the set of $r$ by $n$ row stochastic matrices. Namely,

[TABLE]

with use $Q_{r}\equiv Q_{r,r}$. Further, $S_{r}$ is defined as

[TABLE]

As a consequence, given ${\boldsymbol{X}}\in\mathbb{R}^{n\times d}$, ${\boldsymbol{H}}_{1},{\boldsymbol{H}}_{2}\in\mathbb{R}^{r\times d}$, the loss functions $\mathscrsfs{D}(\,\cdot\,,\,\cdot\,)$ and $\mathscrsfs{L}(\,\cdot\,,\,\cdot\,)$ take the form

[TABLE]

We use ${\sf B}_{m}(\rho)$ to denote the closed ball with radius $\rho$ in $m$ dimensions, centered at [math]. In addition, for ${\boldsymbol{H}}\in\mathbb{R}^{m\times d}$ we define the $\rho$-neighborhood of ${\rm{conv}}({\boldsymbol{H}})$ as

[TABLE]

For a convex set $\mathcal{C}$ we denote the set of its extremal points by ${\rm{ext}}(\mathcal{C})$ and the projection of a point ${\boldsymbol{x}}\in\mathbb{R}^{d}$ onto $\mathcal{C}$ by ${\boldsymbol{\Pi}}_{\mathcal{C}}({\boldsymbol{x}})$. Namely,

[TABLE]

Also, for a matrix ${\boldsymbol{X}}\in\mathbb{R}^{n\times d}$, and a mapping (not necessarily linear) ${\boldsymbol{P}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, ${\boldsymbol{P}}({\boldsymbol{X}})\in\mathbb{R}^{n\times d}$ is the matrix whose $i$’th row is ${\boldsymbol{P}}({\boldsymbol{X}}_{i,.})$.

B.2 Theorem statement

The statement below provides more detailed result with respect to the one in Theorem 1.

Theorem 2**.**

Assume ${\boldsymbol{X}}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}+{\boldsymbol{Z}}$ where the factorization ${\boldsymbol{X}}_{0}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}$ satisfies the uniqueness assumption with parameter $\alpha>0$, and that ${\rm{conv}}({\boldsymbol{X}}_{0})$ has internal radius $\mu>0$. Consider the estimator $\widehat{\boldsymbol{H}}$ defined by Eq. (2.5), with $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$ (square loss) and $\delta=\max_{i\leq n}\|{\boldsymbol{Z}}_{i,\cdot}\|_{2}$. If

[TABLE]

then, setting $\delta=\max_{i\leq n}\|{\boldsymbol{Z}}_{i,\cdot}\|_{2}$ in the problem (2.5) we get

[TABLE]

where $C_{*}$ is a coefficient that depends uniquely on the geometry of ${\boldsymbol{H}}_{0}$, ${\boldsymbol{X}}_{0}$, namely $C_{*}=120(\sigma_{\rm max}({\boldsymbol{H}}_{0})/\mu)\cdot\max(1,\kappa({\boldsymbol{H}}_{0})/\sqrt{r})$.

Further, if

[TABLE]

then, setting $\delta=\max_{i\leq n}\|{\boldsymbol{Z}}_{i,\cdot}\|_{2}$ in the problem (2.5) we get

[TABLE]

where $C_{**}=120\max(\kappa({\boldsymbol{H}}_{0}),(\sigma_{\max}({\boldsymbol{H}}_{0})/r+\|{\boldsymbol{z}}_{0}\|_{2})/(\mu r^{1/2}))\cdot\max(1,\kappa({\boldsymbol{H}}_{0})/\sqrt{r})$.

B.3 Proof

B.3.1 Lemmas

Lemma B.1**.**

Let $\mathcal{R}$ be a convex set and $\mathcal{C}$ be a convex cone. Define

[TABLE]

We have

[TABLE]

An illustration of this lemma in the case of $\mathcal{R}\subset\mathcal{C}$ is given in Figure 7. Note that, $\gamma_{\mathcal{C}}$ measures the pointedness of the cone $\mathcal{C}$. Geometrically (for ${\mathcal{R}}\subseteq{\mathcal{C}}$) the lemma states that the cosine of the angle between $\arg\min_{{\boldsymbol{x}}\in\mathcal{R}}\|{\boldsymbol{x}}\|_{2}$ and $\arg\min_{{\boldsymbol{x}}\in\mathcal{{\rm{ext}}}(R)}\|{\boldsymbol{x}}\|_{2}$ is smaller than $\gamma_{\mathcal{C}}$.

Proof.

We write

[TABLE]

Replacing

[TABLE]

we get

[TABLE]

Hence, using the definition of $\gamma_{\mathcal{C}}$, we have

[TABLE]

Note that

[TABLE]

Therefore,

[TABLE]

and this completes the proof. ∎

The next lemma is a consequence of Lemma B.1.

Lemma B.2**.**

Let ${\boldsymbol{H}},{\boldsymbol{H}}_{0}\in\mathbb{R}^{r\times d}$, $r\leq d$, be matrices with linearly independent rows. We have

[TABLE]

Proof.

Consider the cone $\mathcal{C}_{1}\subset\mathbb{R}^{d}$, generated by vectors ${\boldsymbol{e}}_{2}-{\boldsymbol{e}}_{1},\dots,{\boldsymbol{e}}_{r}-{\boldsymbol{e}}_{1}\in\mathbb{R}^{d}$, i.e,

[TABLE]

For ${\boldsymbol{v}}\in\mathcal{C}_{1},\|{\boldsymbol{v}}\|_{2}=1$ we have

[TABLE]

where ${\boldsymbol{x}}\in\mathbb{R}_{\geq 0}^{r-1}$ and

[TABLE]

Since, $\langle{\bf 1},{\boldsymbol{x}}\rangle=\|{\boldsymbol{x}}\|_{1}\geq\|{\boldsymbol{x}}\|_{2}$, we get $\langle{\bf 1},{\boldsymbol{x}}\rangle\geq 1/\sqrt{2}$. Thus, for ${\boldsymbol{u}}=-{\boldsymbol{e}}_{1}$, we have $\langle{\boldsymbol{u}},{\boldsymbol{v}}\rangle\geq 1/\sqrt{2}$. Therefore, for $\gamma_{{\mathcal{C}}_{1}}$ defined as in Lemma B.1, we have $\gamma_{\mathcal{C}_{1}}\geq 1/\sqrt{2}$. In addition, by symmetry, for $i\in\{1,2,\dots,r\}$, for the cone $\mathcal{C}_{i}\subset\mathbb{R}^{d}$, generated by vectors ${\boldsymbol{e}}_{1}-{\boldsymbol{e}}_{i},{\boldsymbol{e}}_{2}-{\boldsymbol{e}}_{i},\dots,{\boldsymbol{e}}_{r}-{\boldsymbol{e}}_{i}\in\mathbb{R}^{d}$ we have $\gamma_{\mathcal{C}_{i}}=\gamma\geq 1/\sqrt{2}$. Hence, using Lemma B.1 for ${\boldsymbol{H}}\in\mathbb{R}^{r\times d}$, $\mathcal{R}={\rm{conv}}({\boldsymbol{H}})-{\boldsymbol{e}}_{j}$ (the set obtained by translating ${\rm{conv}}({\boldsymbol{H}})$ by $-{\boldsymbol{e}}_{j}$), $\mathcal{C}=\mathcal{C}_{j}$ we get for $j=1,2,\dots,r$

[TABLE]

Hence,

[TABLE]

Therefore,

[TABLE]

Now consider ${\boldsymbol{H}}_{0}\in\mathbb{R}^{r\times d}$ where ${\boldsymbol{H}}_{0}={\boldsymbol{E}}_{r,d}{\boldsymbol{M}}$, ${\boldsymbol{H}}={\boldsymbol{Y}}{\boldsymbol{M}}$, where ${\boldsymbol{M}}\in\mathbb{R}^{d\times d}$ is invertible. We have

[TABLE]

Thus, using the fact that $\sigma_{\max}({\boldsymbol{M}})/\sigma_{\min}({\boldsymbol{M}})=\kappa({\boldsymbol{M}})=\kappa({\boldsymbol{H}}_{0})$,

[TABLE]

Therefore,

[TABLE]

Finally, note that the function $f(x)=(1+x)/x$ is monotone decreasing over $\mathbb{R}_{>0}$. Hence, for $\gamma\geq 1/\sqrt{2}$, $(1+\gamma)/\gamma\leq 1+\sqrt{2}$. Therefore, we get

[TABLE]

and this completes the proof. ∎

We continue with the following lemmas on the condition number of the matrix ${\boldsymbol{H}}$.

Lemma B.3**.**

Let ${\boldsymbol{H}}_{0},{\boldsymbol{H}}\in\mathbb{R}^{r\times d}$, $r\leq d$,with ${\boldsymbol{H}}$ having full row rank. We have

[TABLE]

In addition, if

[TABLE]

then

[TABLE]

Further, if

[TABLE]

then

[TABLE]

Proof.

For the sake of simplicity, we will write $\mathscrsfs{D}_{1}=\mathscrsfs{D}({\boldsymbol{H}},{\boldsymbol{H}}_{0})^{1/2}$, $\mathscrsfs{D}_{2}=\mathscrsfs{D}({\boldsymbol{H}}_{0},{\boldsymbol{H}})^{1/2}$ Note that using the assumptions of Lemma B.3 we have

[TABLE]

where ${\boldsymbol{P}},{\boldsymbol{R}}\in\mathbb{R}_{\geq 0}^{r\times r}$ are row-stochastic matrices and ${\boldsymbol{A}}_{1},{\boldsymbol{A}}_{2}\in\mathbb{R}^{r\times d}$. Also, $\sigma_{\max}({\boldsymbol{A}}_{1})\leq\|{\boldsymbol{A}}_{1}\|_{F}=\mathscrsfs{D}_{1}$, $\sigma_{\max}({\boldsymbol{A}}_{2})\leq\|{\boldsymbol{A}}_{2}\|_{F}=\mathscrsfs{D}_{2}$. Therefore,

[TABLE]

In addition, note that for a row stochastic matrix ${\boldsymbol{P}}\in Q_{r}$, we have

[TABLE]

Hence, for $\mathscrsfs{D}_{2}\leq\sigma_{\min}({\boldsymbol{H}}_{0})$ we get

[TABLE]

In addition,

[TABLE]

Hence, using (B.48), (B.49), for $\mathscrsfs{D}_{2}\leq\sigma_{\min}({\boldsymbol{H}}_{0})$ we have

[TABLE]

Thus, for $\mathscrsfs{D}_{2}\leq\sigma_{\min}({\boldsymbol{H}}_{0})/2$, we get Eqs. (B.39), (B.41).

Now assume that $\mathscrsfs{D}_{1}+\mathscrsfs{D}_{2}\leq\sigma_{\min}({\boldsymbol{H}}_{0})/(6\sqrt{r})$. In this case, using (B.45) we have

[TABLE]

Therefore,

[TABLE]

hence,

[TABLE]

and

[TABLE]

where ${\boldsymbol{H}}_{0}^{\dagger}$ is the right inverse of matrix ${\boldsymbol{H}}_{0}$. Note that

[TABLE]

By permuting the rows and columns of ${\boldsymbol{H}}_{0}$, without loss of generality, we can assume that $R_{ii}=\ \|{\boldsymbol{R}}_{.,i}\|_{\infty}$. We can write

[TABLE]

Hence, for all $i,j\in[r],i\neq j$, since ${\boldsymbol{R}}$ is row-stochastic,

[TABLE]

Thus,

[TABLE]

Therefore, using (B.60),

[TABLE]

Thus, we can write

[TABLE]

Therefore,

[TABLE]

Hence,

[TABLE]

From (B.45) we have $\sigma_{\min}({\boldsymbol{P}}{\boldsymbol{H}})\geq\sigma_{\min}({\boldsymbol{H}}_{0})-\mathscrsfs{D}_{2}$. Using $\sigma_{\min}({\boldsymbol{P}}{\boldsymbol{H}})\leq\sigma_{\max}({\boldsymbol{P}})\sigma_{\min}({\boldsymbol{H}})$, we get

[TABLE]

Further, from (B.45) we have $\sigma_{\max}({\boldsymbol{P}}{\boldsymbol{H}})\leq\sigma_{\max}({\boldsymbol{H}}_{0})+\mathscrsfs{D}_{2}$. Using $\sigma_{\max}({\boldsymbol{P}}{\boldsymbol{H}})\geq\sigma_{\min}({\boldsymbol{P}})\sigma_{\max}({\boldsymbol{H}})$, we get

[TABLE]

Hence, for $\mathscrsfs{D}_{1}+\mathscrsfs{D}_{2}\leq\sigma_{\min}({\boldsymbol{H}}_{0})/(6\sqrt{r})$, we have $\sigma_{\max}({\boldsymbol{H}})\leq 35\sigma_{\max}({\boldsymbol{H}}_{0})/18<2\sigma_{\max}({\boldsymbol{H}}_{0})$. In addition,

[TABLE]

and this completes the proof. ∎

Lemma B.4**.**

Let ${\boldsymbol{X}}_{0}={\boldsymbol{W}}_{0}{\boldsymbol{H}}_{0}\in\mathbb{R}^{n\times d}$ be such that ${\rm{conv}}({\boldsymbol{X}}_{0})$ has internal radius at least $\mu>0$, and ${\boldsymbol{X}}={\boldsymbol{X}}_{0}+{\boldsymbol{Z}}$ with $\max_{i\leq n}\|{\boldsymbol{Z}}_{i,.}\|_{2}\leq\delta$. If ${\boldsymbol{H}}\in\mathbb{R}^{r\times d},{\boldsymbol{H}}_{i,.}\in{\rm{aff}}({\boldsymbol{H}}_{0})$ is feasible for problem (2.5) and has linearly independent rows, then we have

[TABLE]

Proof.

Let

[TABLE]

Note that since ${\boldsymbol{H}}$ is feasible for problem (2.5) and $\max_{i\leq n}\|{\boldsymbol{Z}}_{i,.}\|_{2}\leq\delta$

[TABLE]

Therefore, for any ${\boldsymbol{x}}_{0}\in{\rm{conv}}({\boldsymbol{X}}_{0})$, writing ${\boldsymbol{x}}_{0}={\boldsymbol{X}}_{0}^{\sf T}{\boldsymbol{a}}_{0}$, ${\boldsymbol{a}}_{0}\in\Delta^{n}$, we have

[TABLE]

Since ${\rm{conv}}({\boldsymbol{X}}_{0})$ has internal radius at least $\mu$, there exists ${\boldsymbol{z}}_{0}\in\mathbb{R}^{d}$, and an orthogonal matrix ${\boldsymbol{U}}\in\mathbb{R}^{d\times r^{\prime}}$, $r^{\prime}=r-1$, such that ${\boldsymbol{z}}_{0}+{\boldsymbol{U}}{\sf B}_{r^{\prime}}(\mu)\subseteq{\rm{conv}}({\boldsymbol{X}}_{0})$. Hence, for every ${\boldsymbol{z}}\in\mathbb{R}^{r^{\prime}}$, $\|{\boldsymbol{z}}\|_{2}=1$ there exists ${\boldsymbol{a}}\in\Delta^{n}$ such that

[TABLE]

Therefore, for any unit vector ${\boldsymbol{u}}$ in column space of ${\boldsymbol{U}}$, for the line segment

[TABLE]

Thus,

[TABLE]

where ${\boldsymbol{P}}_{{\boldsymbol{u}}}$ is the orthogonal projection onto the line containing $l_{{\boldsymbol{u}},\mu}$. Note that using (B.76), for any ${\boldsymbol{x}}_{0}\in{\rm{conv}}({\boldsymbol{X}}_{0})$ we have

[TABLE]

In other words, for any ${\boldsymbol{x}}_{0}\in{\boldsymbol{P}}_{{\boldsymbol{u}}}({\rm{conv}}({\boldsymbol{X}}_{0}))$, $D({\boldsymbol{x}}_{0},{\boldsymbol{P}}_{{\boldsymbol{u}}}({\rm{conv}}({\boldsymbol{X}}^{\prime})))\leq 2\delta$. Therefore, using (B.78) for any ${\boldsymbol{u}}$ in column space of ${\boldsymbol{U}}$, we have

[TABLE]

This implies that

[TABLE]

Hence, for every ${\boldsymbol{z}}\in\mathbb{R}^{r^{\prime}}$, $\|{\boldsymbol{z}}\|_{2}=1$ there exists ${\boldsymbol{a}}\in\Delta^{r}$ such that

[TABLE]

Note that ${\boldsymbol{H}}^{\sf T}$ has linearly independent columns. Multiplying the previous equation by $({\boldsymbol{H}}^{\sf T})^{\dagger}$ the left inverse of ${\boldsymbol{H}}^{\sf T}$, we get

[TABLE]

Let

[TABLE]

where ${\boldsymbol{v}}$ is the right singular vector corresponding to the largest singular value of $({\boldsymbol{H}}^{\sf T})^{\dagger}{\boldsymbol{U}}$. Therefore, we have

[TABLE]

Thus, for ${\boldsymbol{a}}_{1},{\boldsymbol{a}}_{2}\in\Delta^{r}$

[TABLE]

Note that

[TABLE]

Thus,

[TABLE]

Hence,

[TABLE]

∎

The following lemma states an important property of $\widehat{\boldsymbol{H}}$ the optimal solution of problem (2.5).

Lemma B.5**.**

If $\max_{i}\|{\boldsymbol{Z}}_{i,.}\|_{2}\leq\delta$ and $\widehat{\boldsymbol{H}}$ is the optimal solution of problem (2.5), then we have

[TABLE]

Proof.

First note that since $\delta\geq\max_{i}\|{\boldsymbol{Z}}_{i,.}\|_{2}$, we have

[TABLE]

Hence, ${\boldsymbol{H}}_{0}$ is a feasible solution for the problem (2.5). Therefore, we have

[TABLE]

Letting $\tilde{\boldsymbol{\alpha}}_{i}=\arg\min_{{\boldsymbol{\alpha}}\in\Delta^{n}}\|{\widehat{\boldsymbol{H}}_{i,.}}^{\sf T}-{\boldsymbol{X}}^{\sf T}{\boldsymbol{\alpha}}\|_{2}$, we have

[TABLE]

Using the fact that (by triangle inequality) $\|{\boldsymbol{Z}}^{\sf T}\tilde{\boldsymbol{\alpha}}_{i}\|_{2}\leq\delta$, we have

[TABLE]

where $U^{2}=\sum_{i=1}^{r}\|{\widehat{\boldsymbol{H}}_{i,.}}^{\sf T}-{\boldsymbol{X}}_{0}^{\sf T}\tilde{\boldsymbol{\alpha}}_{i}\|_{2}^{2}$. Note that $\mathscrsfs{D}(\widehat{\boldsymbol{H}},{\boldsymbol{X}})\geq 0$ and for $U\geq 2\delta\sqrt{r}$, the function $U^{2}-2\delta\sqrt{r}U$ is increasing. Hence, since

[TABLE]

we have

[TABLE]

Therefore,

[TABLE]

In addition,

[TABLE]

Hence,

[TABLE]

Combining equations (B.103), (B.109), and (B.95), we get

[TABLE]

This completes the proof of lemma. ∎

Lemma B.6**.**

Let ${\boldsymbol{X}}_{0}$ be such that the uniqueness assumption holds with parameter $\alpha>0$, and ${\rm{conv}}({\boldsymbol{X}}_{0})$ has internal radius at least $\mu>0$. In particular, we have ${\boldsymbol{z}}_{0}+{\boldsymbol{U}}{\sf B}_{r-1}(\mu)\subseteq{\rm{conv}}({\boldsymbol{X}}_{0})$ for ${\boldsymbol{z}}_{0}\in\mathbb{R}^{d}$, and an orthogonal matrix ${\boldsymbol{U}}\in\mathbb{R}^{d\times(r-1)}$. Finally assume $\max_{i\leq n}\|{\boldsymbol{Z}}_{i,.}\|_{2}\leq\delta$. Then for $\widehat{\boldsymbol{H}}$ the optimal solution of problem (2.5), we have

[TABLE]

where ${\boldsymbol{P}}_{0}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the orthogonal projector onto ${\rm{aff}}({\boldsymbol{H}}_{0})$ (in particular, ${\boldsymbol{P}}_{0}$ is an affine map).

Proof.

Let $\widetilde{\boldsymbol{H}}$ be such that ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{conv}}(\widetilde{\boldsymbol{H}})$. The uniqueness assumption implies

[TABLE]

Note that Lemma B.5 implies

[TABLE]

Therefore,

[TABLE]

Hence,

[TABLE]

In addition, for a convex set $S$, by triangle inequality we have

[TABLE]

Therefore, using

[TABLE]

we have

[TABLE]

Hence,

[TABLE]

and

[TABLE]

In addition, similarly to the proof of Lemma B.5, we can write

[TABLE]

Thus,

[TABLE]

Therefore, combining (B.115), (B.119), (B.120), (B.127), we get

[TABLE]

Now, we would like to bound the terms $\|\widetilde{\boldsymbol{H}}-\widehat{\boldsymbol{H}}\|_{F}$, $\max_{i\in[r]}\|\widetilde{\boldsymbol{H}}_{i,.}-\widehat{\boldsymbol{H}}_{i,.}\|_{2}$. Note that using the fact that $\widehat{\boldsymbol{H}}$ is feasible for Problem (2.5), we have

[TABLE]

Thus,

[TABLE]

In addition, we know that $({\boldsymbol{X}}_{0})_{i,.}\in{\rm{aff}}({\boldsymbol{H}}_{0})$, where ${\rm{aff}}({\boldsymbol{H}}_{0})$ is a $r-1$ dimensional affine subspace. Therefore, ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{aff}}({\boldsymbol{H}}_{0})$ and, by convexity of ${\sf B}_{d}(2\delta,\widehat{\boldsymbol{H}})$, we get

[TABLE]

First consider the case in which $\widehat{\boldsymbol{H}}_{i,.}\in{\rm{aff}}({\boldsymbol{H}}_{0})$. for all $i\in\{1,2,\dots,r\}$. By a perturbation argument, we can assume that the rows of $\widehat{\boldsymbol{H}}$ are linearly independent, and hence ${\rm{aff}}(\widehat{\boldsymbol{H}})={\rm{aff}}({\boldsymbol{H}}_{0})$. Consider $\tilde{\boldsymbol{Q}}\in\mathbb{R}^{r\times d}$ defined by

[TABLE]

where $\xi=2r\delta_{0}$. Note that for every ${\boldsymbol{y}}\in{\sf B}_{d}(2\delta_{0};{\boldsymbol{E}}_{r,d})\cap{\rm{aff}}({\boldsymbol{E}}_{r,d})$, we have $\mathscrsfs{D}({\boldsymbol{y}},{\boldsymbol{E}}_{r,d})^{1/2}\leq 2\delta_{0}$. In addition, since ${\boldsymbol{y}}\in{\rm{aff}}({\boldsymbol{E}}_{r,d})$, $\langle{\boldsymbol{y}},{\bf 1}\rangle=1$. Hence, for ${\boldsymbol{y}}\in{\sf B}_{d}(2\delta_{0};{\boldsymbol{E}}_{r,d})\cap{\rm{aff}}({\boldsymbol{E}}_{r,d})$, we can write

[TABLE]

where ${\boldsymbol{\pi}}\in{\rm{conv}}({\boldsymbol{E}}_{r,d})$, ${\boldsymbol{x}}\in\mathbb{R}^{d}$, $\langle{\bf 1},{\boldsymbol{x}}\rangle=0$, $\|{\boldsymbol{x}}\|_{2}\leq 2\delta_{0}$. It is easy to check that for this ${\boldsymbol{y}}$ we have

[TABLE]

where ${\boldsymbol{\beta}}\in\mathbb{R}^{r}$ is such that for $i=1,2,\dots,r$,

[TABLE]

Further, note that since ${\boldsymbol{\pi}}\in{\rm{conv}}({\boldsymbol{E}}_{r,d})$, $\pi_{i}\geq 0$ and $x_{i}\geq-\|{\boldsymbol{x}}\|_{2}\geq-2\delta_{0}$, we have $\pi_{i}+x_{i}\geq-2\delta_{0}$. Hence, for $i\in\{1,2,\dots,r\}$,

[TABLE]

In addition,

[TABLE]

Therefore, every ${\boldsymbol{y}}\in{\sf B}_{d}(2\delta_{0};{\boldsymbol{E}}_{r,d})\cap{\rm{aff}}({\boldsymbol{E}}_{r,d})$ can be written as a convex combination of the rows of $\widetilde{\boldsymbol{Q}}$. Hence,

[TABLE]

Let $\widehat{\boldsymbol{H}}={\boldsymbol{E}}_{r,d}{\boldsymbol{M}}$, ${\boldsymbol{M}}\in\mathbb{R}^{d\times d}$. Since ${\rm{aff}}(\widehat{\boldsymbol{H}})={\rm{aff}}({\boldsymbol{H}}_{0})$, by taking $\widetilde{\boldsymbol{H}}=\widetilde{\boldsymbol{Q}}{\boldsymbol{M}}$, we have

[TABLE]

provided that $\delta_{0}=\delta/\sigma_{\min}({\boldsymbol{M}})=\delta/\sigma_{\min}(\widehat{\boldsymbol{H}})$. Hence, using (B.131) for this $\delta_{0}$, ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{conv}}(\widetilde{\boldsymbol{H}})$. Note that for $\widetilde{\boldsymbol{Q}}$, we have $\|\widetilde{\boldsymbol{Q}}_{i,.}-{\boldsymbol{e}}_{i}\|_{2}\leq 2r\delta_{0}$. Thus,

[TABLE]

Therefore, there exists $\widetilde{\boldsymbol{H}}\in\mathbb{R}^{r\times d}$ such that ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{conv}}(\widetilde{\boldsymbol{H}})$ and

[TABLE]

Now consider the general case in which ${\rm{aff}}(\widehat{\boldsymbol{H}})\neq{\rm{aff}}({\boldsymbol{H}}_{0})$. Let ${\boldsymbol{H}}^{\prime}\in\mathbb{R}^{r\times d}$ be such that ${\boldsymbol{H}}^{\prime}_{i,.}$ is the projection of $\widehat{\boldsymbol{H}}_{i,.}$ onto ${\rm{aff}}({\boldsymbol{H}}_{0})$. Assuming that the rows of ${\boldsymbol{H}}^{\prime}$ are linearly independent, ${\rm{aff}}({\boldsymbol{H}}^{\prime})={\rm{aff}}({\boldsymbol{H}}_{0})$. Note that since ${\rm{conv}}({\boldsymbol{X}}_{0})\in{\rm{aff}}({\boldsymbol{H}}_{0})$, for every point ${\boldsymbol{x}}\in{\rm{conv}}({\boldsymbol{X}}_{0})$, $\mathscrsfs{D}({\boldsymbol{x}},{\boldsymbol{H}}^{\prime})^{1/2}\leq\mathscrsfs{D}({\boldsymbol{x}},\widehat{\boldsymbol{H}})^{1/2}\leq 2\delta$. Thus,

[TABLE]

Therefore, using the above argument for the case where ${\rm{aff}}(\widehat{\boldsymbol{H}})={\rm{aff}}({\boldsymbol{H}}_{0})$, we can find $\widetilde{\boldsymbol{H}}$ such that ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\rm{conv}}(\widetilde{\boldsymbol{H}})$ and

[TABLE]

Hence, for every $i=1,2,\dots,r$,

[TABLE]

where ${\boldsymbol{P}}_{0}$ orthogonal projection onto ${\rm{aff}}({\boldsymbol{H}}_{0})$. We next use the assumption on the internal radius of ${\rm{conv}}({\boldsymbol{X}}_{0})$ to upper bound the term $\|{\boldsymbol{P}}_{0}(\widehat{\boldsymbol{H}}_{i,.})-\widehat{\boldsymbol{H}}_{i,.}\|_{2}$. Note that since ${\rm{conv}}({\boldsymbol{X}}_{0})\subseteq{\sf B}_{d}(2\delta,\widehat{\boldsymbol{H}})$, letting $\bar{\boldsymbol{H}}=\widehat{\boldsymbol{H}}-{\bf 1}{\boldsymbol{z}}_{0}^{\sf T}$, for some orthogonal matrix ${\boldsymbol{U}}\in\mathbb{R}^{d\times r^{\prime}}$, $r^{\prime}=r-1$, we have

[TABLE]

Now, using Cauchy-Schwarz inequality we can write

[TABLE]

Note that,

[TABLE]

Hence, using (B.151)

[TABLE]

In particular, for $\rho=0$ we get

[TABLE]

Taking $\bar{\boldsymbol{H}}=\tilde{\boldsymbol{U}}{\boldsymbol{\Sigma}}\tilde{\boldsymbol{V}}^{\sf T}$, the singular value decomposition of $\bar{\boldsymbol{H}}$, we have $\sigma_{\max}(\bar{\boldsymbol{H}})=\sigma_{\max}(\widehat{\boldsymbol{H}}-{\bf 1}{\boldsymbol{z}}_{0}^{\sf T})=\max_{i}\Sigma_{ii}$. Letting ${\boldsymbol{U}}^{\sf T}\tilde{\boldsymbol{V}}={\boldsymbol{Q}}$, we get

[TABLE]

Letting $q=\sigma_{\min}({\boldsymbol{Q}})$, this results in

[TABLE]

In addition, note that, by the internal radius assumption, for any ${\boldsymbol{z}}\in\mathbb{R}^{r^{\prime}}$, ${\boldsymbol{z}}_{0}+{\boldsymbol{U}}{\boldsymbol{z}}\in{\rm{aff}}({\boldsymbol{H}}_{0})$. Further, since ${\boldsymbol{z}}_{0}\in{\rm{aff}}({\boldsymbol{H}}_{0})$,

[TABLE]

where ${\boldsymbol{P}}_{\boldsymbol{U}}$ is the projector onto the column space of ${\boldsymbol{U}}$. Note that,

[TABLE]

where the last inequality follows from (B.157). This results in

[TABLE]

Therefore, $\|{\boldsymbol{P}}_{0}(\widehat{\boldsymbol{H}})-\widehat{\boldsymbol{H}}\|_{F}\leq 2\sigma_{\max}(\widehat{\boldsymbol{H}}-{\bf 1}{\boldsymbol{z}}_{0}^{\sf T})\delta\sqrt{r}/\mu$. Hence, using (B.148) we get

[TABLE]

Replacing this in (B.128) completes the proof. ∎

B.3.2 Proof of Theorem 2

For simplicity, let $\mathscrsfs{D}=\alpha(\mathscrsfs{D}(\widehat{\boldsymbol{H}},{\boldsymbol{H}}_{0})^{1/2}+\mathscrsfs{D}({\boldsymbol{H}}_{0},\widehat{\boldsymbol{H}})^{1/2})$. First note that under the assumption of Theorem 2 we have

[TABLE]

Therefore, using Lemma B.4 with ${\boldsymbol{H}}={\boldsymbol{H}}_{0}$ and $\delta=0$, we have

[TABLE]

In addition, since ${\boldsymbol{z}}_{0}\in{\rm{conv}}({\boldsymbol{H}}_{0})$ we have ${\boldsymbol{z}}_{0}={\boldsymbol{H}}_{0}^{\sf T}{\boldsymbol{\alpha}}_{0}$ for some ${\boldsymbol{\alpha}}_{0}\in\Delta^{r}$. Therefore,

[TABLE]

Note that

[TABLE]

Therefore, using Lemma B.6 we have

[TABLE]

In addition, Lemma B.2 implies that

[TABLE]

Further, let ${\boldsymbol{P}}_{0}$ denote the orthogonal projector on ${\rm{aff}}({\boldsymbol{H}}_{0})$. Hence, ${\boldsymbol{P}}_{0}$ is a non-expansive mapping: for ${\boldsymbol{x}},{\boldsymbol{y}}\in\mathbb{R}^{d}$, $D({\boldsymbol{P}}_{0}({\boldsymbol{x}}),{\boldsymbol{P}}_{0}({\boldsymbol{y}}))\leq D({\boldsymbol{x}},{\boldsymbol{y}})$. Therefore, since ${\rm{conv}}({\boldsymbol{H}}_{0})\subset{\rm{aff}}({\boldsymbol{H}}_{0})$, for any ${\boldsymbol{h}}\in\mathbb{R}^{d}$

[TABLE]

Therefore,

[TABLE]

First consider the case in which

[TABLE]

Note that in this case $\delta\leq\mu/2$. Hence, using Lemma B.3 to upper bound $\sigma_{\max}(\widehat{\boldsymbol{H}})$, $\sigma_{\max}({\boldsymbol{P}}_{0}(\widehat{\boldsymbol{H}}))$ and Lemma B.4 to lower bound $\sigma_{\min}({\boldsymbol{P}}_{0}(\widehat{\boldsymbol{H}}))$, by (B.177), we get

[TABLE]

Replacing these in (B.174) we have

[TABLE]

Therefore,

[TABLE]

Notice that condition (B.178) implies that $\mu-2\delta\geq\mu/2$ and

[TABLE]

Using the previous two equations, under condition (B.178) we have

[TABLE]

Combining this with (B.175), and using the fact that $1+2\alpha\leq 3$, we have under condition (B.178)

[TABLE]

Note that using (B.171), (B.172) and since $\alpha\geq 0$

[TABLE]

Therefore,

[TABLE]

Thus,

[TABLE]

where $C_{*}$ is defined in Theorem 2.

Next, consider the case in which

[TABLE]

Note that using (B.171),(B.172) and since $1+2\alpha\leq 3$, this condition on $\delta$ implies that

[TABLE]

In particular, condition (B.178) holds. Hence, using equation (B.184) we get

[TABLE]

Further, note that since ${\boldsymbol{P}}_{0}$ is a projection onto an affine subspace, for ${\boldsymbol{x}}\in\mathbb{R}^{d}$, ${\boldsymbol{P}}_{0}({\boldsymbol{x}})=\widetilde{{\boldsymbol{P}}_{0}}{\boldsymbol{x}}+{\boldsymbol{x}}_{0}$ for some $\widetilde{{\boldsymbol{P}}_{0}}\in\mathbb{R}^{d\times d},{\boldsymbol{x}}_{0}\in\mathbb{R}^{d}$. Hence, for any ${\boldsymbol{\pi}}\in\Delta^{r}$, ${\boldsymbol{h}}=\widehat{\boldsymbol{H}}^{\sf T}{\boldsymbol{\pi}}\in{\rm{conv}}(\widehat{\boldsymbol{H}})$, we have

[TABLE]

where ${\boldsymbol{e}}_{i}$ is the $i$’th standard unit vector. Hence,

[TABLE]

Thus, for ${\boldsymbol{h}}_{0}\in\mathbb{R}^{d}$ an arbitrary row of ${\boldsymbol{H}}_{0}$, we have

[TABLE]

In addition, using non-expansivity of ${\boldsymbol{P}}_{0}$, we have

[TABLE]

This implies that

[TABLE]

Therefore, using (B.177), (B.197) and (B.192) we get

[TABLE]

Hence, in this case Lemma B.3 implies that

[TABLE]

Replacing this in (B.174), we have

[TABLE]

Hence, using (B.175) under assumption (B.190), we have

[TABLE]

Hence, for $C_{*}^{\prime\prime}$ as defined in the statement of the theorem, we get

[TABLE]

This completes the proof.

Appendix C Proof of Proposition 4.2

The proof follows immediately from the following two propositions.

Proposition C.1**.**

Let ${\boldsymbol{X}}\in\mathbb{R}^{n\times d}$ and $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$. Then the gradient of the function ${\boldsymbol{u}}\mapsto\mathscrsfs{D}({\boldsymbol{u}},{\boldsymbol{X}})$ is given by

[TABLE]

Proof.

Note that $\mathscrsfs{D}({\boldsymbol{u}},{\boldsymbol{X}})$ is the solution of the following convex optimization problem.

[TABLE]

The Lagrangian for this problem is

[TABLE]

The KKT condition implies that at the minimizer $({\boldsymbol{y}}^{*},{\boldsymbol{\pi}}^{*},{\boldsymbol{\rho}}^{*},\tilde{\rho}^{*},{\boldsymbol{\lambda}}^{*})$, we have

[TABLE]

and therefore

[TABLE]

and the dual of the above optimization problem is

[TABLE]

Note that since (C.2) is strictly feasible, Slater condition holds and by strong duality the optimal value of (C.6) is equal to $f({\boldsymbol{u}})$. Hence, we have written $f({\boldsymbol{u}})$ as pointwise supremum of functions. Therefore, subgradient of $f({\boldsymbol{u}})$ can be achieved by taking the derivative of the objective function in (C.6) at the optimal solution (see Section 2.10 in [MN13]). Note that the derivative of this objective function at the optimal solution is equal to ${\boldsymbol{\rho}}^{*}=2({\boldsymbol{u}}-{\boldsymbol{y}}^{*})=2({\boldsymbol{u}}-{\boldsymbol{\Pi}}_{{\rm{conv}}({\boldsymbol{X}})}({\boldsymbol{u}}))$ (where we used Eq. (C.5)). Since the dual optimum is unique (by strong convexity in ${\boldsymbol{\rho}}$), the function ${\boldsymbol{u}}\mapsto\mathscrsfs{D}({\boldsymbol{u}},{\boldsymbol{X}})$ is differentiable with gradient given by Eq. (C.1). ∎

Proposition C.2**.**

Let ${\boldsymbol{u}}\in\mathbb{R}^{d}$ and $D({\boldsymbol{x}},{\boldsymbol{y}})=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}^{2}$, and assume that the rows of ${\boldsymbol{H}}_{0}\in\mathbb{R}^{r\times d}$ are affine independent. Then the function ${\boldsymbol{H}}\mapsto\mathscrsfs{D}({\boldsymbol{u}},{\boldsymbol{H}})$ is differentiable at ${\boldsymbol{H}}_{0}$ with gradient

[TABLE]

Proof.

We will denote by ${\boldsymbol{G}}$ the right hand side of Eq. (C.7). For ${\boldsymbol{V}}\in\mathbb{R}^{r\times d}$, we have

[TABLE]

Note that $({\boldsymbol{H}}_{0}+{\boldsymbol{V}})$ has affinely independent rows for ${\boldsymbol{V}}$ in a neighborhood of ${\boldsymbol{0}}$, and hence has a unique minimizer there, that we will denote by ${\boldsymbol{\pi}}_{{\boldsymbol{V}}}$. By optimality of ${\boldsymbol{\pi}}_{{\boldsymbol{V}}}$, we have

[TABLE]

On the other hand, by optimality of ${\boldsymbol{\pi}}_{0}$,

[TABLE]

Letting $R({\boldsymbol{V}})=|\mathscrsfs{D}({\boldsymbol{u}},{\boldsymbol{H}}_{0}+{\boldsymbol{V}})-\mathscrsfs{D}({\boldsymbol{u}},{\boldsymbol{H}}_{0})-\langle{\boldsymbol{G}},{\boldsymbol{V}}\rangle|$ denote the residual, we get

[TABLE]

Note that ${\boldsymbol{\pi}}_{{\boldsymbol{V}}}$ must converge to ${\boldsymbol{\pi}}_{0}$ as ${\boldsymbol{V}}\to 0$ because ${\boldsymbol{\pi}}_{0}$ is the unique minimizer for ${\boldsymbol{V}}={\boldsymbol{0}}$. Hence we get $R({\boldsymbol{V}})/\|{\boldsymbol{V}}\|_{F}\to 0$ as $\|{\boldsymbol{V}}\|_{F}\to 0$, which proves our claim. ∎

Appendix D Proof of Proposition 4.1

We use the results of [BST14] to prove Proposition 4.1. We refer the reader to [BST14] for the definitions of the technical terms in this section. First, consider the function

[TABLE]

Note that using the main theorem of polytope theory (Theorem 1.1 in [Zie12]), we can write

[TABLE]

for some ${\boldsymbol{a}}_{i}\in\mathbb{R}^{d}$, $b_{i}\in\mathbb{R}$ and a finite $m$. Hence, using the definition of the semi-algebraic sets (see Definition 5 in [BST14]), the set ${\rm{conv}}({\boldsymbol{X}})$ is semi-algebraic. Therefore, the function $f({\boldsymbol{H}})$ which is proportional to the sum of squared $\ell_{2}$ distances of the rows of ${\boldsymbol{H}}$ from a semi-algebraic set, is a semi-algebraic function (See Appendix in [BST14]). Further, the function

[TABLE]

is the sum of indicator functions of semi-algebraic sets (Note that using the same argument used for ${\rm{conv}}({\boldsymbol{X}})$, $\Delta^{r}$ is semi-algebraic). Therefore, the function $g$ is semi-algebraic (See Appendix in [BST14]). In addition, the function

[TABLE]

is a polynomial. Hence, it is semi-algebraic. Therefore, we deduce that the function

[TABLE]

is semi-algebraic. In addition, since $\Delta^{r}$ is closed, $\Psi$ is proper and lower semi-continuous. Therefore, $\Psi({\boldsymbol{H}},{\boldsymbol{W}})$ is a KL function (See Theorem 3 in [BST14]).

Now, we will show that the Assumptions 1,2 in [BST14] hold for our algorithm. First, note that since $\Delta^{r}$ is closed, the functions $f({\boldsymbol{H}})$ and $g({\boldsymbol{W}})$ are proper and lower semi-continuous. Further, $f({\boldsymbol{H}})\geq 0,\,g({\boldsymbol{W}})\geq 0,\,h({\boldsymbol{H}},{\boldsymbol{W}})\geq 0$ for all ${\boldsymbol{H}}\in\mathbb{R}^{r\times d},\,{\boldsymbol{W}}\in\mathbb{R}^{n\times r}$. In addition, the function $h({\boldsymbol{H}},{\boldsymbol{W}})$ is $C^{2}$. Therefore, it is Lipschitz continuous over the bounded subsets of $\mathbb{R}^{r\times d}\times\mathbb{R}^{n\times r}$. Also, the partial derivatives of $h({\boldsymbol{H}},{\boldsymbol{W}})$ are

[TABLE]

It can be seen that for any fixed ${\boldsymbol{W}}$, the function ${\boldsymbol{H}}\mapsto\nabla_{\boldsymbol{H}}h({\boldsymbol{H}},{\boldsymbol{W}})$ is Lipschitz continuous with moduli $L_{1}({\boldsymbol{W}})=2\|{\boldsymbol{W}}^{\sf T}{\boldsymbol{W}}\|_{F}$. Similarly, for any fixed ${\boldsymbol{H}}$, the function ${\boldsymbol{W}}\mapsto\nabla_{\boldsymbol{W}}h({\boldsymbol{H}},{\boldsymbol{W}})$ is Lipschitz continuous with moduli $L_{2}({\boldsymbol{H}})=2\|{\boldsymbol{H}}{\boldsymbol{H}}^{\sf T}\|_{F}$. Note that since in each iteration of the algorithm the rows of ${\boldsymbol{W}}^{k}$ are in $\Delta^{r}$. Hence,

[TABLE]

for some some positive constants $\lambda_{1}^{-},\,\lambda_{1}^{+}$. In addition, note that because the PALM algorithm is a descent algorithm, i.e., $\Psi({\boldsymbol{H}}^{k},{\boldsymbol{W}}^{k})\leq\Psi({\boldsymbol{H}}^{k-1},{\boldsymbol{W}}^{k-1})$ for $k\in\mathbb{N}$, and since $f({\boldsymbol{H}})\to\infty$ as $\|{\boldsymbol{H}}\|_{F}\to\infty$, the value of $L_{2}({\boldsymbol{H}}^{k})=\|{\boldsymbol{H}}^{k}{\boldsymbol{H}}^{k^{\sf T}}\|_{F}$ remains bounded in every iteration. Finally, note that by taking $\gamma_{2}^{k}>\max\left\{\left\|{\boldsymbol{H}}^{k+1}{\boldsymbol{H}}^{k+1^{\sf T}}\right\|_{F},{\varepsilon}\right\}$ for some constant ${\varepsilon}>0$, we make sure that the steps in the PALM algorithm remain well defined (See Remark 3(iii) in [BST14]). Hence, we have shown that the assumptions of Theorem 1 in [BST14] hold. Therefore, using this theorem, the sequence $\left\{{\boldsymbol{H}}^{k},\,{\boldsymbol{W}}^{k}\right\}_{k\in\mathbb{N}}$ generated by the iterations in (4.7) - (4.9) has a finite length and it converges to a stationary point $\left({\boldsymbol{H}}^{*},{\boldsymbol{W}}^{*}\right)$ of $\Psi$.

Appendix E Other optimization algorithms

Apart from the proximal alternating linearized minimization discussed in Section 4.2, we experimented with two other algorithms, obtaining comparable results. For the sake of completeness, we describe these algorithms here.

E.1 Stochastic gradient descent

Using any of the initializations discussed in Section 4.1 we iterate

[TABLE]

The step size $\gamma_{t}$ is selected by backtracking line search. Ideally, the direction ${\boldsymbol{G}}^{(t)}$ can be taken to be equal to $\nabla\mathscrsfs{R}_{\lambda}({\boldsymbol{H}}^{(t)})$. However, for large datasets this is computationally impractical, since it requires to compute the projection of each data point onto the set ${\rm{conv}}({\boldsymbol{H}}^{(t)})$. In order to reduce the complexity of the direction calculation, we estimate this sum by subsampling. Namely, we draw a uniformly random set $S_{t}\subseteq[n]$ of fixed size $|S_{t}|=s\leq n$, and compute

[TABLE]

E.2 Alternating minimization

This approach generalizes the original algorithm of [CB94]. We rewrite the objective as a function of ${\boldsymbol{W}}=(bw_{i})_{i\leq n}$, ${\boldsymbol{w}}_{i}\in\Delta^{r}$,${\boldsymbol{H}}=({\boldsymbol{h}}_{i})_{i\leq r}$, ${\boldsymbol{h}}_{i}\in\mathbb{R}^{d}$ and ${\boldsymbol{A}}=({\boldsymbol{\alpha}}_{\ell})_{\ell\leq r}$, ${\boldsymbol{\alpha}}_{\ell}\in\Delta^{n}$

[TABLE]

The algorithm alternates between minimizing with respect to the weights $({\boldsymbol{w}}_{i})_{i\leq n}$ (this can be done independently across $i\in\{1,\dots,n\}$) and minimizing over $({\boldsymbol{h}}_{\ell},{\boldsymbol{\alpha}}_{\ell})$, which is done sequentially by cycling over $\ell\in\{1,\dots,r\}$. Minimization over ${\boldsymbol{w}}_{i}$ can be performed by solving a non-negative least squares problem. As shown in [CB94], minimization over $({\boldsymbol{h}}_{\ell},{\boldsymbol{\alpha}}_{\ell})$ is also equivalent to non-negative least squares. Indeed, by a simple calculation

[TABLE]

where ${\boldsymbol{H}}_{\neq\ell}=({\boldsymbol{h}}_{i})_{i\neq\ell,i\leq r}$, $\widetilde{F}({\boldsymbol{H}},{\boldsymbol{W}},{\boldsymbol{A}})$ does not depend on $({\boldsymbol{h}}_{\ell},{\boldsymbol{\alpha}}_{\ell})$, and we defined

[TABLE]

It is therefore sufficient to minimize $f_{\ell}({\boldsymbol{h}}_{\ell},{\boldsymbol{\alpha}}_{\ell};{\boldsymbol{H}}_{\neq\ell},{\boldsymbol{W}},{\boldsymbol{A}})$ with respect to its first two arguments, which is equivalent to a non-negarive least squares problem. This can be seen by minimizing $f_{\ell}(\cdots)$ explicitly with respect to ${\boldsymbol{h}}_{\ell}$ and writing the resulting objective function.

The pseudocode for this algorithm is given below.

[TABLE]

Here ${\boldsymbol{H}}_{<\ell}=({\boldsymbol{h}}_{i})_{i<\ell}$, ${\boldsymbol{H}}_{>\ell}=({\boldsymbol{h}}_{i})_{\ell<i\leq r}$, and similarly for ${\boldsymbol{A}}$.