A Provably Correct and Robust Algorithm for Convolutive Nonnegative Matrix Factorization

TL;DR

This paper introduces a provably correct, polynomial-time algorithm for convolutive nonnegative matrix factorization (CNMF) under separability assumptions, with applications in audio and neural data analysis.

Contribution

It develops the first provably correct algorithm for CNMF that leverages existing separable NMF methods, ensuring accurate solutions in low noise conditions.

Findings

Algorithm successfully recovers underlying factors in synthetic data.

Effective in separating sources in a singing bird audio sequence.

Runs efficiently in polynomial time.

Abstract

In this paper, we propose a provably correct algorithm for convolutive nonnegative matrix factorization (CNMF) under separability assumptions. CNMF is a convolutive variant of nonnegative matrix factorization (NMF), which functions as an NMF with additional sequential structure. This model is useful in a number of applications, such as audio source separation and neural sequence identification. While a number of heuristic algorithms have been proposed to solve CNMF, to the best of our knowledge no provably correct algorithms have been developed. We present an algorithm that takes advantage of the NMF model underlying CNMF and exploits existing algorithms for separable NMF to provably find a solution under certain conditions. Our approach guarantees the solution in low noise settings, and runs in polynomial time. We illustrate its effectiveness on synthetic datasets, and on a singing bird audio sequence.

Figures

Equations

$$V=\left[\begin{array}{ccc}{W}_1& \dots & W_L\end{array}\right]\text{ and }G=\left[\begin{array}{c}{H}_1\\ ⋮\\ H_L\end{array}\right].$$

$$\underset{\ell =0,\dots ,L-1}{\max }\max \left(\cos \left(S_{\ell }^T{g}_i,g_j\right),\cos \left(g_i,S_{\ell }^T{g}_j\right)\right),$$

$$H[k,j]$$

$${\mu }_{\text{left}}(i,j)=\underset{\ell =0,\dots ,L-1}{\max }\cos (S_{\ell }^T\tilde{G}[i,:]-\tilde{G}[j,:]),$$

$${\mu }_{\text{right}}(i,j)=\underset{\ell =0,\dots ,L-1}{\max }\cos (\tilde{G}[i,:]-S_{\ell }^T\tilde{G}[j,:]).$$

$$(i{\le }_{L,\tilde{G}}j)=\{\begin{array}{ll}\text{True}& \text{if }{\mu }_{\text{left}}(i,j)\ge {\mu }_{\text{right}}(i,j)\\ \text{False}& \text{otherwise}\end{array}.$$

$$ϵ+t$$

$$\frac{ϵ}{t}$$

$$ϵ$$

$$\underset{1\le j\le K}{\min }\cos \left(h_j,{\tilde{h}}_{P(j)}\right)\ge 1-\frac{\delta }{2},$$

$$M_{ij}=\{\begin{array}{ll}1& \text{if }\exists k\in \{1,\dots ,K\}\text{ such that }i,j\in {\mathcal{C}}_k,\\ 0& \text{otherwise.}\end{array}$$

$$\text{score}(H,\tilde{H})=\underset{P\in \mathcal{P}}{\min }\frac{1}{K}\sum _{k=1}^K\cos (H[k,:],\tilde{H}[P(k),:]),$$

$${\text{MSE}}_{\text{rel}}=\frac{\Vert \tilde{X}-{\sum }_{\ell =1}^L{W}_{\ell }H{S}_{\ell -1}{\Vert }_F}{\Vert \tilde{X}{\Vert }_F},$$

$$\underset{1\le j\le R}{\max }\Vert \tilde{X}[:,J(j)]-V[:,P(j)]\Vert =C_2ϵ\kappa (V{)}^2,$$

$$\frac{ϵ}{t}<\frac{C_1{\sigma }_{\min }(V{D}_V)}{2\sqrt{R}\kappa (V{D}_V{)}^2}.$$

$$\underset{1\le j\le r}{\max }\Vert (\tilde{X}{D}_{\tilde{X}})[:,J(j)]-(V{D}_V)[:,P(j)]{\Vert }_2$$

$$\le 2C_2(ϵt^{-1})\kappa (V{D}_V{)}^2,\text{ and }$$

$$\underset{1\le j\le r}{\max }\Vert \tilde{X}[:,J(j)]-V[:,P(j)]{\Vert }_2$$

$$\le (\Vert V{\Vert }_{1,col}+ϵ)2C_2(ϵt^{-1})\kappa (V{D}_V{)}^2.$$

$$Y_{jj}=\{\begin{array}{ll}0& \text{ if }\Vert {\tilde{x}}_j{\Vert }_1<t\\ \Vert {\tilde{x}}_j{\Vert }_1^{-1}& \text{ otherwise}\end{array}.$$

$${ϵ}^{\prime }$$

$$\le \frac{\Vert e_j{\Vert }_2}{\Vert x_j+e_j{\Vert }_1}+\left(1-\frac{\Vert x_j{\Vert }_1}{\Vert x_j+e_j{\Vert }_1}\right)\frac{\Vert x_j{\Vert }_2}{\Vert x_j{\Vert }_1}$$

$$\le \frac{\Vert e_j{\Vert }_2}{\Vert x_j+e_j{\Vert }_1}+\left(\frac{\Vert x_j{\Vert }_1+\Vert e_j{\Vert }_1-\Vert x_j{\Vert }_1}{\Vert x_j+e_j{\Vert }_1}\right)\frac{\Vert x_j{\Vert }_2}{\Vert x_j{\Vert }_1}$$

$$\le \frac{\Vert e_j{\Vert }_2}{\Vert x_j+e_j{\Vert }_1}+\frac{\Vert e_j{\Vert }_1}{\Vert x_j+e_j{\Vert }_1}\le \frac{2ϵ}{t}.$$

$$\Vert y-\tilde{y}{\Vert }_2\le \frac{{\sigma }_{\max }(F)\Vert y{\Vert }_2+\Vert e{\Vert }_2}{{\sigma }_{\min }(A)-{\sigma }_{\max }(F)}.$$

$$\tilde{G}$$

$$\Vert g_j-{\tilde{g}}_j{\Vert }_2^2\le {\left(\frac{{\sigma }_{\max }(F)\Vert g_j{\Vert }_2+ϵ}{{\sigma }_{\min }(V)-{\sigma }_{\max }(F)}\right)}^2.$$

$$\Vert G[i,:]-\tilde{G}[i,:]{\Vert }_2^2\le T{\left(\frac{{\sigma }_{\max }(F)\gamma +ϵ}{{\sigma }_{\min }(V)-{\sigma }_{\max }(F)}\right)}^2.$$

$$\frac{\theta -3ϵ}{(1+ϵ{)}^2}\le \cos (y_1+e_1,y_2+e_2)\le \frac{\theta +3ϵ}{(1-ϵ{)}^2}.$$

$$\theta -35ϵ\le \cos (y_1+e_1,y_2+e_2)\le \theta +35ϵ.$$

$$\frac{\theta -3ϵ}{(1+ϵ{)}^2}$$

$$\ge (\theta -3ϵ)(1-ϵ{)}^2\text{ (since ϵ2n−ϵ2n+1>0)}$$

$$\ge (\theta -3ϵ)(1-2ϵ)\ge \theta -5ϵ.$$

$$\max (\cos (S_{\ell }^T{\tilde{g}}_j,{\tilde{g}}_i),\cos ({\tilde{g}}_j,S_{\ell }^T{\tilde{g}}_i))<1-\frac{\delta }{2},$$

$$\underset{1\le j\le KL}{\max }\Vert g_j-{\tilde{g}}_{P(j)}{\Vert }_2=ϵ<\frac{\delta }{70}\Vert G{\Vert }_{-2,row},$$

Full text

1 Introduction

Nonnegative matrix factorization (NMF) is a standard unsupervised learning technique for analyzing large datasets. Given an $N\times T$ matrix $X$, NMF seeks a $N\times K$ matrix $W\geq 0$ (where the inequality is to be interpreted elementwise) and a $K\times T$ matrix $H\geq 0$ such that $X\approx WH$ and $K\ll\min(N,T)$. NMF has been successfully applied to a number of practical problems; these include hyperspectral unmixing, text mining, audio source separation, and image processing; see [5, 13, 11] and the references therein.

Let us discuss three important limitations of NMF, each of which we will address in this paper.

A first limitation of NMF is that it fails to capture local correlations in the data. For example, in imaging applications a column of $X$ may represent a pixel, and adjacent columns will often correspond to pixels adjacent to one another in the image. Neighboring pixels tend to be quite similar, especially in low contrast images. In audio or neuroscience datasets, each column is a certain instant in time, and therefore neighboring columns are often highly correlated. To capture these local correlations, Smaragdis [32] proposed a convolutive variant of NMF, known as convolutive NMF (CNMF). CNMF attempts to find $L$ matrices $W_{1},\ldots,W_{L}$ of size $N\times K$ and a matrix $H$ of size $K\times T$ such that $X=\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}$, where $S_{\tau}$ is a square matrix that shifts the columns of $H$ by $\tau$ places to the right, zero-padding the leftmost $\tau$ columns. Explicitly, $S_{\tau}$ has ones on its $\tau$th upper diagonal and zeros elsewhere. It is often convenient to instead define an $L\times N\times K$ tensor $W$ and denote $W_{\ell::}=W_{\ell}$ the resulting $N\times K$ matrix when the first index is fixed to $\ell$. Then another, perhaps more intuitive, definition of CNMF is to equivalently write $X=\sum_{k=1}^{K}W_{::k}^{T}*h_{k}$, where $*$ is the 2D-convolution operator defined as

$(A*b)_{ij}=\sum_{\tau=1}^{L}A_{i\tau}b_{j-\tau}\in\mathbf{R}^{N\times T}$, with $A\in\mathbf{R}^{N\times L}$ and $b\in\mathbf{R}^{1\times T}$, $W_{::k}\in\mathbf{R}^{L\times N}$ is the resulting matrix when the third index is fixed to $k$, and $h_{k}\in\mathbf{R}^{1\times T}$ is the $k$th row of $H$. Thus, CNMF is a factorization of $X$ into a sum of 2D-convolutions. This characterization is visually apparent in Figure 1.

Since its conception, CNMF has been found particularly useful for audio source separation [32, 31, 35] and neural sequence identification [27]. In general, CNMF will provide a more concise and interpretable factorization than NMF whenever the columns of $X$ exhibit local correlations, which is often the case when the columns represent points in time or space.

A second limitation of NMF is that solving NMF problems is challenging in practice. In fact, it is in general NP-hard [33], leading researchers to rely on a number of heuristic algorithms; see, e.g., [13] and the references therein. As CNMF is a generalization of NMF (the model is exactly NMF when $L=1$), it inherits problems related to computational intractability. The first algorithm proposed was an averaged multiplicative update rule [32], which treats each pair $(W_{\ell},HS_{\ell-1})$ as an NMF and updates them using an NMF multiplicative update [23], averaging the updates for $H$. Since then, numerous other algorithms have been proposed to fit CNMF, including multiplicative updates [31, 34, 18, 9] (which are derived using majorization minimization), projected alternating least squares [30], alternating nonnegative least squares (ANLS; [7]), and hierarchical alternating least squares (HALS; [7]). Many of these algorithms are generalizations of algorithms for NMF and perform well on both synthetic and experimental data. Several of them also have convergence guarantees to stationary points. However, all the aforementioned algorithms are heuristics, in the sense that there is no guarantee that they reach a global minimum.

A third limitation of NMF is that, in general, the NMF $WH$ of a matrix $X$ is non-unique, that is, there are other $W^{\prime},H^{\prime}$ such that $WH=X=W^{\prime}H^{\prime}$, but $W\neq W^{\prime}$ and $H\neq H^{\prime}$. We therefore cannot be certain that the recovered factors $(W,H)$ correspond to the true sources that generated the data. The same problem arises with CNMF; see Section 5.1 for some examples. In practice, most researchers use additional regularization terms to promote some structure in the sought solution such as sparsity [20]. Unfortunately, identifiability results for NMF are rather scarce; see the recent survey [11]. As far as we know, the only NMF model that is both tractable and identifiable is separable NMF [3]. When trying to identify an NMF $WH=X$, separable NMF makes the additional assumption that the columns of $W$ are contained somewhere in $X$, that is, $W=X[:,\mathcal{C}]$ for some index set $\mathcal{C}$. This assumption is called the separability assumption. This terminology was introduced in the paper of Donoho and Stodden [8] where it was shown that NMF has a unique solution under the condition that $W=X[:,\mathcal{C}]$ and additional sparsity conditions on $H$. However, this NMF model dates back from the 1990’s in the hyperspectral imaging literature where the separability assumption is referred to as the pure-pixel assumption; see [26] and the references therein. Since the paper of Arora et al. [3], which proved that separable NMF model is unique, robust to noise, and computable in polynomial time, numerous algorithms for separable NMF have been proposed [14]. Moreover, this model has been used successfully in applications such as topic modeling [2], community detection [29], and the analysis of time-resolved Raman spectra [25], to cite a few. Of course, the separability assumption is rather strong and does not hold in all applications. Nevertheless, it has also be shown to be a powerful way to initialize more sophisticated NMF models, such as minimum-volume NMF models, which are identifiable but suffer from intractability [11]. One of the most popular and powerful separable NMF algorithm is the successive projection algorithm (SPA) [1] which is guaranteed to recover the columns of $W$, even in the presence of noise [16]. SPA is a greedy algorithm that sequentially identifies the columns of $W$ using successive orthogonal projections (see Algorithm 3 and Appendix A for more details.). SPA is closely related to the modified Gram-Schmidt algorithm with column pivoting, and has rediscovered many times; see the discussion in [13].

1.1 Contribution and Outline of the Paper

In this paper, we consider conditions under which the CNMF problem can be provably solved in polynomial time, even in the presence of noise, and has a unique solution (up to permutation and scaling of the rank-one a factors). These conditions generalize the separability assumption discussed in the previous paragraph. We propose an algorithm that reduces CNMF to NMF with linear constraints and takes advantage of existing separability-based methods for NMF. This algorithm provably finds the solution to CNMF in polynomial time, both for exact problems and problems with bounded noise. In particular, we utilize SPA to estimate the columns of $W_{1},\ldots,W_{\ell}$, then apply estimation, clustering, and sorting techniques to uncover the convolutive structure. We later generalize our approach to use any separable NMF algorithm and show how the choice of the separable NMF algorithm affects our assumptions, run-time, and noise tolerance. To the best of our knowledge, this is the first algorithm to provably solve the CNMF problem in either the absence or presence of noise.

Hence, under separability-like assumptions, our approach resolved the three important limitations of NMF mentioned in the introduction.

The paper is organized as follows. In Section 2, we state the CNMF problem and formally define the convolutive separability assumptions. In Section 3, we propose the Locate-Estimate-Cluster-Sort Algorithm (LECS) for recovery and show that LECS provably finds the solution to the CNMF problem under the convolutive separability assumptions. In particular, Theorem 2 guarantees that LECS will find the unique optimal solution (up to permutation and scaling) in the absence of noise. Theorem 3 generalizes this result to problems with bounded noise. The proof of Theorem 3 is deferred to the Appendix.

Notation

The $j$th entry of a vector $a$ is denoted by $a_{j}$. For matrices $A\in\mathbf{R}^{m\times n}$, we use $A[i,:]=A_{i:}$ and $A[:,j]=A_{:j}$ to denote the $i$th row of $A$ and $j$th column of $A$, respectively. We will also use lower case letters for columns and rows, but define them explicitly first. The entry in the $i$th row and the $j$th column of $A$ will be denoted $A[i,j]=A_{ij}$. A nonnegative vector $a$ or a nonnegative matrix $A$ are denoted using $a\in\mathbf{R}^{n}_{+}$ and $A\in\mathbf{R}^{m\times n}_{+}$. Define $S_{\tau}$ as a sqaure $T$-by-$T$ matrix with 1 on its $\tau$th upper diagonal and zeros elsewhere. For vectors $a\in\mathbf{R}^{n}$, the $p$-norm for $1\leq p<\infty$ is defined as $\|a\|_{p}=\left(\sum_{j=1}^{n}|a_{j}|^{p}\right)^{1/p}$. The maximum and minimum $p$-norms of a column of $A$ are denoted by $\|A\|_{p,col}=\max_{1\leq j\leq n}\|A_{:j}\|_{p}$ and $\|A\|_{-p,col}=\min_{1\leq j\leq n}\|A_{:j}\|_{p}$, respectively. We also define $\|A\|_{p,row}=\|A^{T}\|_{p,col}$ and $\|A\|_{-p,row}=\|A^{T}\|_{-p,col}$. The Frobenius norm is denoted as $\|A\|_{F}=\left(\sum_{i,j=1}^{m,n}A_{ij}^{2}\right)^{1/2}$. We denote $\sigma_{\max}(A)$ to be the largest singular value of $A$, and $\sigma_{\min}(A)$ to be the smallest singular value of $A$. The condition number of a matrix $A$ induced by the 2-norm is denoted as $\mathbf{\kappa}(A)=\sigma_{\max}(A)/\sigma_{\min}(A)$. Given a matrix $A\in\mathbf{R}^{m\times n}$, the diagonal matrix $D_{A}\in\mathbf{R}^{n\times n}$ is defined by $(D_{A})_{ii}=\begin{cases}0&\text{ if }\|A_{:i}\|_{1}=0\\ \|A_{:i}\|_{1}^{-1}&\text{ otherwise}\end{cases}$.

2 Problem Setup and the Separability Conditions

Suppose there is a matrix $X\in\mathbf{R}^{N\times T}_{+}$ generated as $X=\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}$, where $W_{\ell}\in\mathbf{R}^{N\times K}_{+}$ for $\ell=1,2,\dots,L$, and $H\in\mathbf{R}^{K\times T}_{+}$. Consider $\tilde{X}=X+E$, where $E$ is some matrix of noise. Given $\tilde{X},K,L$, the convolutive NMF problem (or CNMF problem) is to approximately recover $W_{\ell},H$ (up to scaling or permutation). When $E=0$ and hence $\tilde{X}=X$, this problem amounts to finding an exact convolutive NMF.

2.1 Reformulation as Constrained NMF

Our approach will take advantage of existing literature on NMF. In particular, we look to separable NMF algorithms and attempt to extend them to the convolutive case. In this vein, it is useful to reformulate the convolutive NMF model as an NMF model with linear constraints of $H$.

Consider the sum $\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}$. Define $H_{\ell}=HS_{\ell-1}$ for $\ell=1,2,\dots,L$. Now define

[TABLE]

It follows that $VG=\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}=X$. Now we may think of the exact CNMF problem as a restriction of the exact NMF problem; Given $X,K,L$, our goal is to find $V$ and $G$ such that $VG=X$ and $G$ is defined as in (1), that is, each block $H_{\ell}$ is a shifted version of the first block $H_{1}$. The above formulation further elucidates the relationship between NMF and CNMF. When $L=1$, CNMF reduces to NMF. On the other hand, any CNMF (given by $X,W_{1},\ldots,W_{\ell},H$) is also an NMF (given by $X,V,G$). Therefore, when the $V,G$ found from NMF have the form in (1), a CNMF has also been identified. Given the reformulation of (1), our approach is to utilize separable NMF algorithms to locate the columns of $V$ within $\tilde{X}$ and then estimate $G$. Once these matrices have been identified up to permutation and scaling, we use clustering and sorting methods to identify $H$ and $W_{1},\ldots,W_{L}$. We name this approach the Locate-Estimate-Cluster-Sort Algorithm (LECS); see Algorithm 1 which is described in details in Section 3.

2.2 Convolutive Separability

To introduce a notion of separability for the convolutive NMF model, we need to address two relevant details that are not present in the standard NMF model. First, to handle time dependencies, we need some notion of distance up to a shift. The following definition formalizes this notion mathematically.

Definition 1** ($L$-shift similarity).**

Consider two vectors $g_{i},g_{j}\in\mathbf{R}^{T}$. The $L$-shift cosine similarity $\cos_{L}$ between $g_{i}$ and $g_{j}$ is denoted $\cos_{L}(g_{i},g_{j})$ and given by

[TABLE]

where $\cos(x,y)=(x^{T}y)/(\|x\|_{2}\|y\|_{2})$. This is exactly the cosine of the minimum angle between $g_{i},g_{j}$ over all shifts of length $\ell<L$.

The utility of the $L$-shift similarity comes from the fact that $\cos_{L}(x,y)=1$ is equivalent to the statement that either $x=\alpha S_{\ell}^{T}y$ or $y=\alpha S_{\ell}^{T}x$ for some $\ell\in\{0,\ldots,L-1\}$ and for some scalar $\alpha>0$.

The second detail is minor but essential. Separable NMF algorithms like SPA require that each column of $X$ lie in the convex hull of $V$, that is, $\sum_{i}G[i,j]\leq 1$ for all $j$ (see Appendix A for more details). Nevertheless, when this does not hold we can still scale the columns of $V$ and $X$ to sum to one, and thus apply SPA. This is not the case in the convolutive NMF model, since the magnitude of the sequence $W_{::k}$ may change over time (for example, a crescendo in an audio sequence). Thus, we will have to allow arbitrary positive scaling of $V$ by some positive diagonal matrix $A$, and apply SPA to $VA$.

With these details in mind, we can now formulate a concept of separability for the CNMF problem. We first provide a definition, then explain each condition in more details.

Definition 2** (Convolutive separable).**

Given a CNMF problem with inputs $\tilde{X}=X+E$ where $X=\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}$, $L$ and $K$, we say the problem is convolutive separable with respect to $\delta,\epsilon,A$, where $\delta>0$, $\epsilon\geq 0$, and $A\in\mathbf{R}^{KL\times KL}$ is a diagonal matrix with strictly positive diagonal entries, if the following conditions are satisfied:

(A) [Separable] For each $k,\ell$, the vector $W_{k}[:,\ell]$ appears as a scaled column of $X$. Explicitly, there are several equivalent ways to express this condition.

$\bullet$ The matrix $V$ satisfies $VA=X[:,\mathcal{C}]$ for some set $\mathcal{C}$.

$\bullet$ The matrix $G$ satisfies $G=[A\ M]\Pi$ for some $M\in\mathbf{R}^{KL\times T-KL}_{+}$ and permutation matrix $\Pi\in\{0,1\}^{T\times T}$.

(B) [Sequentially unique]

$\bullet$ For any two rows $h_{i},h_{j}$ of $H$, the $2L$-shift similarity between them is $\cos_{2L}(h_{i},h_{j})\leq 1-\delta$.

$\bullet$ For any row $h_{i}$ of $H$ and any $\ell,\tau=1,\ldots,L-1$, we have $\cos(S_{\tau}^{T}h_{i},S_{\ell}^{T}h_{i})\leq 1-\delta$ when $\tau\neq\ell$.

(C) [Full Rank] The matrix $V$ is full rank, that is, $V$ has rank $KL$.

(D) [Bounded noise] The noise matrix $E$ satisfies $\|E\|_{1,col}\leq\epsilon$.

Definition 2 guarantees a unique solution in the following sense.

Theorem 1** (Unique Solution).**

Suppose $X=\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}$ satisfies (A), (B), and (C). Then any other CNMF given by $W_{1}^{\prime},\ldots,W_{L}^{\prime},H^{\prime}$ that satisfies (A), (B), (C), and $X=\sum_{\ell=1}^{L}W_{\ell}^{\prime}H^{\prime}S_{\ell-1}$ differs at most by permutation and scaling, in the sense that $V=V^{\prime}\Lambda\Pi$ for a positive diagonal matrix $\Lambda\in\mathbf{R}^{T\times T}$ and a permutation matrix $\Pi\in\mathbf{R}^{T\times T}$.

The proof of this result follows directly from Theorem 2, since our proposed algorithm (namely, LECS described in Algorithm 1) is deterministic and guaranteed to recover the underlying CNMF up to permutation and scaling. The same approach entails the uniqueness of a solution to a noisy CNMF problem $\tilde{X}=\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}+E$ among solutions satisfying the conditions of Theorem 3.

The first condition [Separable] is inherited directly from the separable NMF literature and used only to guarantee that a separable NMF algorithm can identify the columns of $V$. One notable exception is that in separable NMF problems, one usually requires that $V=X[:,\mathcal{C}]$. In this paper, we relax this condition to $VA=X[:,\mathcal{C}]$ for some diagonal matrix $A$ with strictly positive diagonal entries. This weaker constraint allows the columns of $V$ to be arbitrarily scaled in the matrix $X$. However, this more general condition makes the problem no harder, since we can use a separable NMF algorithm to identify the columns of $VA$ and use this scaled matrix throughout the rest of our procedure. In Section 4, we show that [Separable] can actually be weakened to the more general sufficiently scattered condition.

The second condition, [Sequentially unique]  tells us that each row of $H$ must differ from all other rows and their shifted variants. Moreover, it tells us the rows must be distinct from its own shifted variants. Since we are generally unconcerned with the scale of a row of $H$, the cosine of the angle between two vectors serves as a useful scale-invariant measurement of similarity (if the vectors have mean zero, it is exactly the correlation between the two vectors). [Sequentially unique] will be an important condition after we have a permuted estimate of $G$, since we will need to rearrange the rows to obtain a matrix that satisfies the constraints of (1).

In practice, this is a very reasonable assumption; all we require is (a) that any two factors $h_{i},h_{j}$ must not be highly correlated (in which case we should just have one factor), and (b) that the factor $h_{j}$ is not periodic with a period smaller than $L$ (in which case we can choose a smaller $L$).

In our approach, [Full Rank] guarantees that SPA terminates successfully and locates the columns of $V$. We also use it to bound the error when estimating $G$ via nonnegative least squares. However, other algorithms like the Successive Nonnegative Projection Algorithm (SNPA) from [12] do not require [Full Rank]. Similarly, [24] defines least squares robustness results in the rank deficient case. Although we assume [Full Rank] in our results to simplify our analysis, it can be weakened to the condition that no column of $V$ is contained in the convex cone generated by the other columns of $V$. We elaborate on this in Section 4.

[Bounded noise] bounds the 1-norm of any column of the noise matrix $E$. When $E=0$, this is clearly satisfied. Otherwise, we will later see that if $\epsilon$ is sufficiently small, we can still recover noisy estimates of $W_{1},\ldots,W_{\ell},H$. Note that this bound differs slightly from the noise bounds in many separable NMF algorithms. In particular, most separable NMF algorithms require bounds on the 2-norm of any column of the noise matrix $E$, whereas we bound the 1-norm. Since $\|v\|_{2}\leq\|v\|_{1}$ for any vector $v\in\mathbf{R}^{n}$, we make a stronger assumption on the noise. This stronger condition is intimately related to the following. In separable NMF, we require $\|H\|_{1,col}\leq 1$. In practical applications where this might not hold, we scale the columns of $X$ so that this assumption is satisfied and show that the new problem is equivalent; see for example [16]. This scaling trick does not apply to CNMF however, since the matrix $H$ is difficult to appropriately scale in the expression $\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}$. In Appendix A, we explain this scaling trick and its limitations in more details.

Is convolutive separability reasonable in practice? In challenging real-world scenarios, the convolutive separability condition will most likely be violated, because several sources will overlap in all time windows. Also, the noise level will be typically higher than what our bounds allow (see Theorem 3). These two issues are analogous to those encountered by separable NMF. Nevertheless, as explained in the introduction, separable NMF has been used successful in many challenging real-world problems, either directly or as an initialization strategy. Accordingly, we believe many real-world signals can be approximated well under the convolutive separability assumption because, although several sources are active in all time windows, some time windows will contain a single source that is dominant. Moreover, as it will be shown in the experiments, our approach leads to efficient initialization strategies in challenging scenarios. Additionally, our approach can be generalized in a number of different ways (see Section 4). For example, the separability assumption can be relaxed to the sufficiently scattered condition by using minimum-volume NMF instead of separable NMF.

3 Algorithm Description and Recovery Guarantee

In this section, we first describe in details our proposed algorithm (LECS, Algorithm 1) in Section 3.1, and then prove its correctness under the convolutive separability assumptions in the noiseless case (Section 3.2) and in the presence of noise in (Section 3.3).

3.1 Description of the LECS Algorithm

At a high level, our approach will be as follows: we will use the successive projection algorithm (SPA, [16]) to identify the column indices $\mathcal{C}$ of $\tilde{X}$ that correspond to scaled columns of $V$, then leverage this noisy version of $V$ to obtain some estimate of $G$. Finally, we will cluster and sort the rows of $G$ to find the rows of $H$ up to some error.

The proposed Locate-Estimate-Cluster-Sort (LECS) algorithm is broken down into four main steps.

$\bullet$ Locate. First, LECS ‘locates’ the column indices $\mathcal{C}$ of $\tilde{X}$ that correspond to $V^{\prime}=VAD_{VA}$; the columns of $VA$ are present in $X$, and accordingly their noisy variants are present in $\tilde{X}$, so locating the column indices $\mathcal{C}$ such that $VA=X[:,\mathcal{C}]$ will identify a scaled, permuted variant of $V$. The algorithm used a modified version of SPA [16], which we call Oracle Conic SPA (OrConSPA); see Algorithm 3 and Appendix A for more details.

This method first removes the columns of $\tilde{X}$ that have a low signal to noise ratio using a hyperparameter $t$. Then, the algorithm rescales the columns of $\tilde{X}$ in order to approximately project the columns of $\tilde{X}$ from the conic hull of the columns of $VA$ onto the convex hull of the columns of $VA$ (the projection is approximate because of the noise). Finally, the algorithm applies SPA to recover the index set $\mathcal{C}$. The name Oracle Conic SPA comes from the ‘oracle’ used to select the hyperparamter $t$ and from the fact that the algorithm allows $X$ to be in the conic hull of $VA$, in contrast to standard SPA which requires $X$ to be in the convex hull of $VA$. This method requires $\mathcal{O}(NTKL)$ operations [16], applying SPA being the most expensive step.

$\bullet$ Estimate. In this step, LECS estimates the rows of $G^{\prime}=D_{VA}^{-1}A^{-1}G$ using nonnegative least squares (NNLS). Estimating $G^{\prime}$ is essential because the convolutive structure is contained in its rows, which are shifted variants of the rows of $H$. Solving a convex NNLS problem up to any given precision can be performed in polynomial time using an interior point method (IPM). However, IPMs are second-order methods and hence are computationally demanding. We therefore instead use the block pivot method from [22] which requires one least squares solve per iteration. Each least squares solve takes $\mathcal{O}(NK^{2}L^{2}+NTKL)$ operations, and the algorithm almost always converges after a few iterations.

$\bullet$ Cluster. LECS then clusters the rows of $G^{\prime}$ into $K$ groups $\mathcal{C}_{1},\ldots,\mathcal{C}_{K}$ according to which row of $H$ they are shifted variants of. One cluster should contain the $L$ shifted variants of one row of $H$. The algorithm achieves this by computing the $L$-shift similarity between every pair of rows in $G^{\prime}$, then greedily constructing clusters by adding the available row with the highest average similarity to the rows in the cluster. This simple greedy procedure requires $\mathcal{O}(TK^{2}L^{3})$ operations; see Algorithm 4 and Lemma A.8.

$\bullet$ Sort. Finally, within each cluster $C_{k}$, LECS sorts the rows of $G^{\prime}$ based on their shifted similarity to the other rows in the cluster. The algorithm uses a comparison-based sorting algorithm with a comparison operator $\leq$ defined by shifted angle scores. In particular, for two indices $i,j$, we have $i\leq j$ if $G[j,:]$ can be better expressed as a shifted copy of $G[i,:]$ than the other way around, measured via the cosine between the two vectors. This comparison operator is guaranteed to produce a strict, consistent ordering when the criteria of Theorem 3 are satisfied. However, note that, when these conditions fail to hold there is no guarantee that the comparison operator will produce a coherent ordering. Once this comparison operator has been defined, we can use any comparison-based sorting algorithm (such as merge sort or quick sort) to order the cluster $C_{k}$; in our implementation, we use a simple selection sort. When the operator is not consistent (that is, $i\leq j$, $j\leq k$ but $i\nleq k$ for some indices $i,j,k$), we select at each step the index which is less than or equal to the most other vectors (breaking ties arbitrarily). This sorting procedure requires $\mathcal{O}(TL^{3})$ operations; see Algorithm 5 and Lemma A.9.

Finally, using the clustering and sorting, the algorithm reconstructs $W_{1},\ldots,W_{\ell}$ rearranging the columns of $V$. Each row of $H$ is constructed by taking the de-shifted average of the rows of $G$ in a particular cluster. The total computational cost of Algorithm 1 is $\mathcal{O}(NTKL+TK^{2}L^{3})$ operations plus the time to solve the nonnegative least squares problem.

3.2 Guarantee for the Exact Problem

The exact problem is when $E=0$. Then [Bounded noise] is satisfied with $\epsilon=0$. In this case, Algorithm 1 provably recovers $H$ and each $W_{\ell}$ for any $\delta>0$ and any $A$ such that $\min_{i}A_{ii}>0$.

Theorem 2** (Exact Recovery).**

Given a CNMF problem that is convolutive separable with respect to $\epsilon=0$ and some $\delta,A$, Algorithm 1 with inputs $\tilde{X}=X,K,L$ recovers the unique factorization $H,W_{1},\ldots,W_{\ell}$ satisfying (A), (B), and (C), up to a permutation and scaling, in polynomial time; more precisely $\mathcal{O}(NTKL+TK^{2}L^{3})$ flops plus one NNLS solve.

Proof of Correctness.

Because of conditions [Full Rank] and [Separable], SPA is guaranteed to find the columns of $V$. [Full Rank] also guarantees nonnegative least squares will return the rows of $G$. All that remains is to show the grouping correctly recovers $W$ and $H$.

Given any two rows $g_{i},g_{j}$ that are both shifted versions of some row of $H$, their similarity is $\cos_{L}(g_{i},g_{j})=1$. Rows that are not shifted variants of one another cannot have similarity greater than or equal to $1-\delta$ by condition [Sequentially unique]. This means the grouping given by the clustering step is necessarily correct.

[Sequentially unique] also entails that the ordering from the sort step is also correct. Since the ordering is correct, the construction of $H$ and each $W_{\ell}$ must also be correct up to a permutation and scaling of the rows of $H$.

Uniqueness follows from the success of the algorithm. Suppose we have two factorizations $W_{1},\ldots,W_{L},H$ and $W_{1}^{\prime},\ldots,W_{L}^{\prime},H^{\prime}$ such that $X=\sum_{\ell=1}^{L}W_{\ell}HS_{\ell-1}=\sum_{\ell=1}^{L}W_{\ell}^{\prime}H^{\prime}S_{\ell-1}$. Then given $X$, Algorithm 1 recovers the first factorization up to permutation and scaling. However, Algorithm 1 also recovers the latter factorization up to permutation and scaling. Hence it must be that one factorization is a permute, scaled version of the other, that is, $H=\Pi\Lambda H^{\prime}$ for some permutation matrix $\Pi\in\{0,1\}^{K\times K}$ and some diagonal matrix $\Lambda\in\mathbf{R}^{K\times K}$ with positive diagonal elements. ∎

3.3 Guarantee for the Noisy Problem

Even when $E\neq 0$, we can still recover noisy estimates of $H,W_{1},\ldots,W_{\ell}$ given sufficiently low noise levels.

Theorem 3** (Noisy Recovery).**

Suppose a CNMF problem with inputs $\tilde{X},L,K$ is convolutive separable with respect to $\delta,\epsilon,A$ and consider some parameter $t>0$. Let $V^{\prime}=VAD_{VA}$ and $G^{\prime}=D_{VA}^{-1}A^{-1}G$ and suppose

[TABLE]

where $\rho=C_{a}\left(\frac{\sigma_{\min}(V^{\prime})}{\sqrt{KL}\mathbf{\kappa}(V^{\prime})^{2}}\right)-\frac{\epsilon}{t}$ and $C_{a},C_{b}$ are universal constants independent of all other terms. Furthermore, without loss of generality assume that $\|VA\|_{2,col}>1$ and $t>1$ (see below). Then in polynomial time, more precisely $\mathcal{O}(NTKL+TK^{2}L^{3})$ flops plus one NNLS solve, Algorithm 1 finds $\tilde{H}\in\mathbf{R}^{K\times T}$ with bounded error, in the sense that there exists some permutation $P$ such that

[TABLE]

where $h_{j},\tilde{h}_{P(j)}$ are the $j$th and ${P(j)}$th rows of $H$ and $\tilde{H}$, respectively.

The proof of this theorem is deferred to Appendix A.3, where we prove a more general statement in Lemma A.10. The assumptions $\|VA\|_{2,col}>1$ and $t>1$ are only used to simplify (4) and make the bound easier to read. These assumptions are also made without any loss of generality—we can always consider the equivalent problem with $\alpha V$ and $\alpha\tilde{X}=\alpha X+\alpha E$ for an arbitrarily large $\alpha>0$; in particular for $\alpha=\max(t^{-1},\|VA\|_{2,col}^{-1})$. This scaling does not impact the output of the algorithm (beyond scaling), and cancelling terms shows it also does not affect inequality (2) and inequality (3). However, one should note this scaling does impact $\|G^{\prime}\|_{2,col}$ since $G^{\prime}$ depends on $D_{VA}^{-1}$, which changes inequality (4).

Algorithm 1 depends on selecting a good value for the parameter $t$; the best error bounds is obtained as $t\rightarrow\|VA\|_{-1,col}-\epsilon$. This parameter is used to threshold columns of $\tilde{X}$ with a small norm, since these columns will have relatively low signal-to-noise ratios and cause SPA to fail after scaling $\tilde{X}$. In practice, one can run the algorithm several times using different values of $t$ and keep the best fit according to some heuristic (for example, mean square reconstruction error).

In fact, since the conditions of the Theorem 3 imply that noisy versions of the columns of $VA$ are present in $\tilde{X}$, we can be certain that $\|\tilde{X}[:,p]\|_{-1,col}$ is within $\epsilon$ of $\|VA\|_{-1,col}$ for some $p\in\{1,\ldots,T\}$. So we could run the algorithm $T$ times with parameter $t=\|X[:,p]\|-2\epsilon$ for $p=1,2,\dots,T$ and guarantee that one of the runs yields the desired solution.

4 Generalizations

One strength of our approach is that any of the subroutines used for the four steps (location, estimation, clustering, and sorting) can be substituted in exchange for different subroutines. For example, we use SPA to locate the columns of $V$, but could substitute any separable NMF algorithm in its place. In this section, we provide several examples of improvements to LECS obtained by appropriate substitutions and highlight the generality of the high level algorithm.

Rank-deficient identification of $V$

We previously mentioned that [Full Rank] is required for two reasons: to guarantee the success of the SPA algorithm and to bound the error when estimating $G$ via nonnegative least squares. However, if we replace SPA with another separable NMF algorithm that does not require $V$ to have full column rank, then the first reason is no longer necessary. SNPA in particular essentially requires the weaker condition that no column of $V$ can be written as a conic combination of the other columns of $V$. By using this method in place of SPA, we can complete the location step with weaker condition than [Full Rank]. This condition is also not necessary for estimating $G$, which we discuss later in the section.

Sufficiently Scattered Condition

Since [Separable] is only used in identifying the columns of $V$, we could also use the sufficiently scattered [21, 11], which is more general and allows one to determine $V$ by identifying the minimum-volume NMF of $X$. Unfortunately, this more general condition suffers from two drawbacks: first, it has yet to be proven robust to noise, and, second, there are currently no algorithms that are guaranteed to find the minimum-volume NMF in polynomial time [11] (this problem is NP-hard in general, but it remains unknown whether this is still true under the sufficiently scattered condition).

Rank-deficient Estimation of $G$

Under some assumptions, the matrix $V$ need not be full rank to properly estimate $G$. In particular, [24] gives perturbation bounds in the rank-deficient case in a more general setting; these bounds apply to our problem when the rank of $\tilde{V}$ is the same as the rank of $V$.

Robust Least Squares

When we estimate $G$, we solve the problem $\min_{G\geq 0}\|\tilde{V}G-\tilde{X}\|_{F}^{2}$ in an attempt to find some $\tilde{G}$ with rows similar to those of $G$. If we know the noise level $\epsilon$ in advance or have some estimate of it, the optimal estimate is given by the robust optimization problem: $\min_{\tilde{G}\geq 0}\max_{(V,G,X)\in\mathcal{C}}\|\tilde{G}-G\|_{2,row}$, where $\mathcal{C}=\{(V,G,X):\|V-\tilde{V}\|_{2,col}\leq\epsilon_{V},\|X-\tilde{X}\|_{1,col}\leq\epsilon_{X},V,G,X\text{ satisfy (\ref{eq: block convsum}) and (A)-(D)}\}$, where $\epsilon_{V},\epsilon_{X}$ depend on $\epsilon$. In general, a number of least squares-like procedures may be used in place of standard nonnegative least squares in order to minimize the worst case error $\|\tilde{G}-G\|_{2,row}$.

Generalized Similarity Metric

The $L$-shift similarity measure utilizes the cosine of the angle between two vectors because it is scale invariant, allowing us to consider $a$ and $\alpha a$ equivalent, where $a\in\mathbf{R}^{n}$ and $\alpha\in\mathbf{R},\alpha>0$. The choice of similarity measure heavily influences the clustering and sorting steps, as seen in the proofs of Lemma A.8 and Lemma A.9. Therefore, a wise adjustment of this measure could lead to significant improvements in the final two steps of LECS. This might include a simple post-processing step, such as exponentiation, or an entirely different similarity measure. See [28, 10] for a more comprehensive review.

Spectral Clustering

Our greedy clustering algorithm is simple but suboptimal, in the sense that it does not take advantage of the global cluster structure. One more advanced approach is to use a spectral clustering method [28]. Let the matrix $M$ be defined as

[TABLE]

Let $v_{1},\ldots,v_{k}$ be the eigenvectors corresponding to the $K$ largest eigenvalues of $M$. Orient each eigenvector so that $|\max_{j}(v_{i})_{j}|\geq|\min_{j}(v_{i})_{j}|$ (that is, ensure the entry with the largest magnitude is positive). Then the indices of the largest $L$ entries in $v_{i}$ form the representatives for cluster $i$. The challenge is to construct a similarity matrix $\tilde{M}$ close to $M$ using only the observed data. One option is to construct $\tilde{M}$ using the $L$-shift similarity, then set the largest $KL^{2}$ entries to 1 and the remaining entries to 0. Under the conditions of Theorem 3, $\tilde{M}=M$ and therefore recovery is guaranteed. Results from matrix perturbation theory, such as the Davis-Kahan bound [6], suggest there is a strong theoretical justification for using a spectral clustering method, given the right construction of $\tilde{M}$.

5 Numerical Experiments

We first test the LECS algorithm on synthetic data and verify it correctly finds the ground truth in low noise settings. Then, we demonstrate that LECS finds reasonable results on the spectrogram of a songbird and can be used as an effective initialization for other algorithms like multiplicative updates [31] and alternating nonnegative least squares (ANLS) [7].

All code is written in Julia [4] and available on GitHub at github.com/degleris1/CMF.jl. In particular, the code used at the time of publication is available in release v0.1. We use the code from github.com/ahwillia/NonNegLeastSquares.jl to solve NNLS problems via the pivot method. The figures for Sections 5.1 and 5.2 are produced by the Jupyter notebooks figures/sep_synth.ipynb and figures/sep_song.ipynb, and the LECS algorithm itself is available file src/algs/separable.jl. Note that this repository is available to other researchers as a tool for fitting CNMF models, as well as rapidly developing and testing new CNMF algorithms.

5.1 Unstructured Synthetic Data

In this experiment, we consider synthetic data that is separable by construction, but lacks structure in the sense that matrices are generated randomly without constraints on important parameters like condition number. Since angle determines the uniqueness of a vector up to scaling and plays an important role in our theoretical analysis, we use it to measure performance. Specifically, after the algorithm estimates $\tilde{H}$, the score is computed as the average cosine of the angles between the rows of $\tilde{H}$ and the rows of $H$ (minimized over all row permutations). This is formally given by

[TABLE]

where $\mathcal{P}$ is the set of all permutations. We test four different algorithms:

$\bullet$ LECS. The LECS algorithm as presented in Algorithm 1.

$\bullet$ LECS-Pre. The LECS algorithm with SPA replaced by heuristically preconditioned SPA from [17]. Heuristically preconditioned SPA computes the SVD $U\Sigma V^{T}$ of $\tilde{X}$, then runs SPA on $\Sigma^{-1}U^{T}\tilde{X}=V^{T}$. Under the assumption that the data points are evenly distributed within the convex cone generated by $V$, it is more robust than SPA [15]. This illustrates the flexibility of LECS with respect to the choice of the different building blocks; see Section 4.

$\bullet$ Mult. The multiplicative updates algorithm for minimizing the Frobenius norm used in [31, 27]. This also corresponds to the $\beta$-divergence update rule from [18, 9] when $\beta=2$. The algorithm is initialized randomly and scaled using the method from [19]. The alternating update rules runs for 5 seconds.

$\bullet$ LECS-Mult. The same multiplicative updates algorithm as before, except using LECS as its initialization.

For each $\ell\in\{1,\ldots,L\}$, we let $(W_{\ell})_{ij}\sim\mathbf{Unif}(0.5,1.5)$ for all $i,j$, where $\mathbf{Unif}(a,b)$ is the uniform distribution in the interval $[a,b]$. We let $H_{ij}\sim\mathbf{Unif}(0,1)\mathbf{Bern}(1-p)$ for all $i,j$, where $p$ is our sparsity parameter and $\mathbf{Bern}(1-p)$ is the Bernoulli distribution with parameter $1-p$. To enforce the separability condition, we then set some entries of $H$ to zero as follows. For each row $H[k,:]$ of $H$, choose two random $t_{k},s_{k}\in\{0,\ldots T-L\}$. For all $t$ satisfying either $t_{k}-L\leq t\leq t_{k}+\lfloor L/2\rfloor$ or $s_{k}-\lfloor L/2\rfloor\leq t\leq s_{k}+L$, we set $H[:,t]=0$. Then we let $H[k,t_{k}],H[k,s_{k}]\sim\mathbf{Unif}(0.5,1.5)$. Finally, we ensure each $t_{1},s_{1},\ldots,t_{K},s_{K}$ are sufficiently far apart from one another so that the same entries are never set to zero twice. This construction ensures the first half of each sequence ($W_{::k}$) is separable at one location, and the second half at another location. We choose “half-sequences” to demonstrate the full sequence need not be separable at the same location. To generate the noise, we vary a noise parameter $\beta$ and use one of three procedures: for all $i,j$

$\bullet$ Uniform. Set $E_{ij}\sim\mathbf{Unif}(0,\beta)$.

$\bullet$ Gaussian. Set $E_{ij}=\max(-X_{ij},M_{ij})$ where $M_{ij}\sim\mathbf{Normal}(0,\beta)$.

$\bullet$ Exponential. Set $E_{ij}\sim\mathbf{Exponential}(\beta)$.

We generate $n_{\text{trials}}=10$ random matrices $X$ according the above construction with $N=100$, $T=250$, $K=3$, $L=5$, $p=0.75$. Then, for each matrix, the noise level $\beta$ is varied across $n_{\text{noise}}=13$ different values spaced logarithmically between $\beta=10^{-3}$ and $\beta=10^{3}$.

The results of these experiments are displayed in Figure 2. We observe that for all three noise types, LECS finds the true $H$ for small noise levels, as expected by Theorem 3. Additionally, as expected, the preconditioning improves its tolerance to noise. In comparison, Mult fail to find the true $H$ even in the lowest noise settings. In fact, although Mult finds a solution with low reconstruction error, it is not able to recover the ground truth factor $H$: this is due to the non-uniqueness problem with CNMF (see Section 1). As a matter of fact, Mult with LECS as an initialization keeps the true solution and outperforms random initializations. At higher noise levels, the effects are varied. In particular, LECS-Mult exhibits different performance depending on the type of noise, and notably performs the best (relative to random initialization) when the noise is Gaussian. This is likely due to the fact this update rule attempts to minimize the reconstruction error using a Frobenius loss function, which is a natural objective for Gaussian noise but not uniform nor exponential noise. A different loss function is likely necessary to identify the ground truth given these noise models; for example, entrywise $\infty$-norm and $1$-norm losses may be more appropriate for uniform and exponential noise, respectively.

5.2 Songbird Spectrogram

In practice, noise levels are often too high for the LECS algorithm to identify the ground truth, even when the separability assumptions are satisfied. However, LECS can be used as an initialization to other algorithms for CNMF to achieve faster convergence than random initialization, and to obtain better solutions.

In this experiment, we fit a spectrogram of a singing bird from [27]; the authors have kindly made this dataset publicly available at github.com/FeeLab/seqNMF. The spectrogram matrix $\tilde{X}$ has 141 rows (DFT bins) and 4440 columns (timebins). We fit the spectrogram using the LECS algorithm with $K=3,\ L=20$, and a threshold value $t=10$. Then, we use this initialization to run Mult [31] and ANLS [7] algorithms for 15 and 60 iterations, respectively. Note that ANLS, as for NMF [22], solves the subproblem exactly for $W$ and $H$ using an active set method, as opposed to Mult that is a gradient-like descent method. The threshold value $t$ was chosen after sweeping over all values from 0 to 50 (approximately the maximum 1-norm of a column of $\tilde{X}$).

Figure 3 displays the loss curves for Mult and ANLS using both random initialization and the output of LECS. The loss is measured as the relative mean square error, defined as

[TABLE]

where $W_{1},\ldots,W_{L},H$ are the estimated matrices.

We observe that when either Mult or ANLS are initialized using the output of LECS, they converge to a local minimum in fewer iterations. This suggests that LECS finds a solution near some local minimum and could be used to accelerate algorithms on very large datasets. In fact, when used as an initialization for ANLS, LECS allows to obtain a high quality solution faster than random initializations. For Mult, using LECS as an initialization does not perform as well because, on average and after sufficiently many iterations, the solution generated by random initializations have lower reconstruction error. We suspect that the reason is that Mult does not deal well with input matrices with many small entries (which LECS generates) because of the zero locking phenomenon–zero entries cannot be modified by Mult while it takes many iterations for small entries to get large; see for example the discussion in [13] and the references therein. In fact, with the LECS initialization, ANLS reduces the relative error to 56.6% after 15 iterations, while MU reduces it to only 58.4% after 60 iterations. Note also that ANLS performs on average better than MU with an average relative error of 56.6% for ANLS vs. 57.7% for MU; this was already observed [7].

Moreover, despite a high relative mean square error, Figure LABEL:fig:_song-qual demonstrates that LECS outputs factors and a reconstruction with qualitative similarities to the ground truth. In fact, after just a single iteration of the ANLS algorithm, the factors visibly mimic the sequential structure in the dataset. Note that, as opposed to the experiments on the synthetic data set performed in the previous section, we cannot assess the quality of the generated factors as the ground truth sources are not available.

6 Conclusion

In this paper we have presented a provably correct and robust algorithm for CNMF under separability conditions. By addressing the questions of identifiability and provable recovery of the solution, this paper offers a novel theoretical perspective on the CNMF problem. Our algorithm draws directly from methods for separable NMF and is easily generalizable—the high level algorithm does not require any particular method for each of its four steps.

Future work includes addressing two significant weakness of our approach. First, the location step does not leverage the convolutive structure of CNMF (and could therefore possibly be made more tolerant to noise). Second, the method cannot handle the case when some $W_{::k}$ has repeated columns, which is relevant to processing music datasets that often contain pure harmonic tones.

Appendix A Appendix

A.1 Conic SPA Procedure

As mentioned in the introduction, SPA is a sequential algorithm for separable NMF that identifies the columns of $V$ among the columns of $X=VG$. At each step, it first extracts the column of $X$ that has the largest $\ell_{2}$ norm, and then project all columns of $X$ onto the orthogonal complement of the extracted column; see Algorithm 2. Under the assumptions that $V=X[:,\mathcal{C}]$ is full column rank and the columns of $G$ have unit $\ell_{1}$ norm, SPA recovers the set $V$, even in the presence of noise. In this section we generalize SPA from an algorithm for finding the vertices of a polytope (because the columns of $G$ are required to have unit $\ell_{1}$ norm) to an algorithm for finding the extreme rays of a cone. In the noiseless case, it is relatively straightforward to show this generalization by an appropriate rescaling of $X$. However, when noise is present, extra steps must be taken to ensure points that are primarily noise do not grow too large and are mistaken for a vertex.

In fact, in the robustness analysis of SPA, it is assumed from the scratch that the input noisy data matrix has the form $\tilde{X}=VG+E$, where the columns of $G$ are normalized which is key for SPA to succeed. If $G$ is not normalized, then normalization of the input matrix is necessary and this might increase the noise drastically; for example for the columns of $X$ with very small norm. Therefore, to make SPA robust to noise in this scenario, we will need to first remove the columns of $\tilde{X}$ with small norm.

Lemma A.1** (SPA).**

[16*, Theorem 3]**

Let $\tilde{X}=X+E=VG+E\in\mathbf{R}^{N\times T}$, where $V\in\mathbf{R}^{N\times R}$ has rank $R$ and $G\in\mathbf{R}^{R\times T}_{+}$ satisfies $\sum_{i}G[i,j]\leq 1$ for all $1\leq j\leq T$. Additionally, suppose that $G=[I\ M]\Pi$ for some permutation matrix $\Pi$ and $M\in\mathbf{R}^{R\times T-R}_{+}$.*

Let $\|E[:,i]\|\leq\epsilon$ for all $i$ with $\epsilon<\frac{C_{1}\sigma_{\min}(V)}{\sqrt{R}\mathbf{\kappa}(V)^{2}}$, where $C_{1}$ is a global constant. Let $J$ be the index set extracted by SPA, which takes time $\mathcal{O}(NTKL)$. Then there exists a permutation $P$ of $\{1,\ldots,R\}$ such that

[TABLE]

where $C_{2}$ is a global constant.

Note that Lemma A.1 makes two strong assumptions on $G$; first, that $\sum_{i}G[i,j]\leq 1$ for all $1\leq j\leq T$, and second, that $G$ contains the identity matrix as a submatrix. Because of the additional structure of the convolutive NMF problem, it is unreasonable to have both these assumptions simultaneously. Instead, we drop the first assumption and generalize the SPA algorithm to any $G$ such that $G=[I\ M]\Pi\in\mathbf{R}^{R\times T}_{+}$, where $\Pi$ is a permutation matrix and $M\in\mathbf{R}^{R\times T-R}_{+}$. Note that this formulation is equivalent to $G=[\Lambda\ M]\Pi$, where $\Lambda$ is some full rank diagonal matrix, since we can always scale the $j$th row of $G$ by $\lambda_{j}$ and scale the $j$th column of $V$ by $\lambda_{j}^{-1}$ without changing $VG$. We first demonstrate conditions under which we can recover $V$ if some oracle gives us a threshold value $t$ which is greater than the noise level $\epsilon$ but smaller than $\|V\|_{-1}-\epsilon$.

Lemma A.2** (Oracle Conic SPA).**

Let $\tilde{X}=X+E=VG+E\in\mathbf{R}^{N\times T}_{+}$, where $V\in\mathbf{R}^{N\times R}_{+}$ has rank $R$ and $G=[I\ M]\Pi\in\mathbf{R}^{R\times T}_{+}$ for some permutation matrix $\Pi$ and $M\in\mathbf{R}^{R\times T-R}_{+}$. Denote $e_{j},v_{j}$ the $j$th column of $E$ and $V$, respectively. Let $\|e_{j}\|_{1}<\epsilon$ for all $j$ and suppose we know some $t>0$ such that $t+\epsilon<\|V\|_{-1,col}$. Assume

[TABLE]

Let $J$ be the index set extracted by OrConSPA with inputs $\tilde{X},R,t$. Then there exists a permutation $P$ of $\{1,\ldots,R\}$ such that

[TABLE]

[TABLE]

Proof.

Let $v_{j}$ be the $j$th column of $V$, and allow the same notation for all other matrices. Since the oracle has given us $t\leq\|V\|_{-1}-\epsilon_{1}$ (where $t>0$), we know that any column $j$ with $\|\tilde{x}_{j}\|_{1}<t$ cannot satisfy $x_{j}=v_{i}$ for any $i$. Therefore, we can define $Y\in\mathbf{R}^{T\times T}$ as a diagonal matrix given by

[TABLE]

Now define $\tilde{X}^{\prime}=\tilde{X}Y$. Consider the matrices $V^{\prime}=VD_{V}$, $X^{\prime}=XD_{X}=V^{\prime}G^{\prime}$, and $G^{\prime}=D_{V}^{-1}HD_{X}$. It must be that the columns of $G^{\prime}$ sum to one, since the columns of $X^{\prime}$ and $V^{\prime}$ sum to one. Furthermore, columns $g_{j}$ of $G$ with only a single non-zero entry will only be scaled by some scalar $\alpha$, so $G^{\prime}=[I\ M^{\prime}]\Pi$ for some $M^{\prime}\in\mathbf{R}^{N\times T-R}_{+}$ and permutation matrix $\Pi$. Therefore, $X^{\prime},V^{\prime},G^{\prime}$ all satisfy the conditions of Lemma A.1; if we wish to apply SPA to identify $V^{\prime}$ from $\tilde{X}^{\prime}=X^{\prime}+E^{\prime}$, all that remains is to bound the noise $E^{\prime}=\tilde{X}^{\prime}-X^{\prime}$. To bound the noise, we note that for any non-zero column of $\tilde{X}^{\prime}$, we have $e^{\prime}_{j}=\tilde{x}^{\prime}_{j}-x^{\prime}_{j}=\frac{x_{j}+e_{j}}{\|x_{j}+e_{j}\|_{1}}-\frac{x_{j}}{\|x_{j}\|_{1}}=\frac{e_{j}}{\|x_{j}+e_{j}\|_{1}}-\left(1-\frac{\|x_{j}\|_{1}}{\|x_{j}+e_{j}\|_{1}}\right)\frac{x_{j}}{\|x_{j}\|_{1}}$, and therefore we have the following bound for $\epsilon^{\prime}=\|e^{\prime}_{j}\|_{2}$:

[TABLE]

Concerning the zero columns, they cannot be noisy versions of a column of $V$ and therefore setting them to zero is irrelevant to the performance of SPA. Therefore, if we apply SPA to $X^{\prime}$, then Lemma A.1 and (5) tells us our assumptions are satisfactory to recover a noisy estimate of $V^{\prime}$ with error $C_{2}(2\epsilon t^{-1})\mathbf{\kappa}(VD_{V})^{2}$. Scaling each column by $D_{\tilde{X}}^{-1}$ gives the error in (7). ∎

Tighter bounds are achieved as $t$ approaches $\|V\|_{-1}-\epsilon$. In real applications, this oracle is not available and $t$ is unknown. However, we do know that if such a $t$ exists, then $0\leq t\leq\|\tilde{X}\|_{1}$ and that $t=\|X_{:j}\|$ for some $j$. Since SPA is very fast (essentially $2NTR$ flops), it is reasonable to run OrConSPA for multiple values of $t$ and be certain that one of the outputs locates the columns of $V$ (that is, if such a value of $t$ exists at all). In practice, a measure like mean square error can be used to select the “best” value of $t$. The previous proof specifically uses SPA to the thresholded and scaled matrix $\tilde{X}^{\prime}=\tilde{X}Y$. However, there is nothing special about our choice of separable NMF algorithm. In fact, other algorithms such as SNPA [12] or AnchorWords [2] will yield more robust noise tolerances. For a more complete list of separable NMF algorithms, see [14]. In this vein, we can generalize Lemma A.2 to an arbitrary separable NMF algorithm.

Lemma A.3** (Oracle Conic Location).**

Let $\tilde{X}=X+E=VG+E\in\mathbf{R}^{N\times T}_{+}$, where $V\in\mathbf{R}^{N\times R}_{+}$ and $G=[I\ M]\Pi\in\mathbf{R}^{R\times T}_{+}$ for some permutation matrix $\Pi$ and $M\in\mathbf{R}^{R\times T-R}_{+}$.

Consider some arbitrary separable NMF algorithm ALG with a maximum input noise tolerance $\epsilon_{\text{in}}<m_{1}(V)$ and a maximum output noise $m_{2}(V,\epsilon_{\text{in}})$. Furthermore, suppose $VD_{V}$ satisfies any other conditions of ALG (e.g. a condition on the rank). Denote $e_{j},v_{j}$ the $j$th column of $E$ and $V$, respectively. Now suppose $\|e_{j}\|_{1}\leq\epsilon$ for all $j$ and suppose we know of some $t>0$ such that $t+\epsilon<\|V\|_{-1}$ and $\frac{\epsilon}{t}<\frac{1}{2}m_{1}(VD_{V})$. Let $J$ be the index set extracted by OCL with inputs $\tilde{X},R,t$, and ALG. Then there exists a permutation $P$ of $\{1,\ldots,R\}$ such that $\max_{j}\|\tilde{X}[:,J(j)]-V[:,P(j)]\|_{2}\leq(\|V\|_{1}+\epsilon)m_{2}(VD_{V},2\epsilon t^{-1})$.

A.2 Relevant Lemmas

A.2.1 Perturbation bounds for least squares

Lemma A.4** (Least Squares Perturbation, [24] Theorem 3, Remark 2).**

Consider full rank matrices $A,\tilde{A}\in\mathbf{R}^{T\times R}$ and vectors $b,\tilde{b}\in\mathbf{R}^{T}$. Define $F=A-\tilde{A}$ and $e=b-\tilde{b}$ Suppose there is a convex set $\mathcal{C}$ and a vector $y\in\mathcal{C}$ such that $Ay=b$ and let $\tilde{y}=\mathop{\rm argmin}_{y^{\prime}\in\mathcal{C}}\|\tilde{A}y^{\prime}-\tilde{b}\|_{2}^{2}$. If $\sigma_{\min}(A)>\sigma_{\max}(F)$, then

[TABLE]

Corollary A.5**.**

Suppose we are given $\tilde{X},X\in\mathbf{R}^{T\times N}$ such that $E=X-\tilde{X}$ satisfies $\|E[:,j]\|\leq\epsilon$ for all $j$. Additionally, suppose we also have $\tilde{V},V\in\mathbf{R}^{N\times KL}$ such that $F=V-\tilde{V}$ satisfies $\sigma_{\min}(V)>\sigma_{\max}(F)$. Suppose there is a $G\geq 0$ such that $VG=X$ and let

[TABLE]

Also let $\gamma=\|G[:,j]\|_{2,col}$ be the maximum norm of any column of $G$. Then the distance between rows of $G^{*}$ and $\hat{G}$ is bounded by $\max_{i}\ \|G[i,:]-\hat{G}[i,:]\|_{2}\leq\sqrt{T}\left(\frac{\sigma_{\max}(F)\gamma+\epsilon}{\sigma_{\min}(V)-\sigma_{\max}(F)}\right)$.

Proof.

For the $j$th column $g_{j}$ of $G$, we have the bound

[TABLE]

Then for a row $G[i,:]$ we must have

[TABLE]

Taking the square root yields the desired bound. ∎

A.2.2 Cosine Similarity

Lemma A.6**.**

Consider the vectors $y_{1},y_{2},e_{1},e_{2}\in\mathbf{R}^{T}$. Let $\theta=\cos(y_{1},y_{2})>0$. Suppose $\|e_{1}\|_{1}\leq\epsilon\|y_{1}\|_{1}$ and $\|e_{2}\|_{2}\leq\epsilon\|y_{2}\|_{2}$ for $0<\epsilon<1$. Then we have

[TABLE]

Proof.

First consider when $\|y_{1}\|_{2}=\|y_{2}\|_{2}=1$. We have $\frac{(y_{1}+e_{1})^{T}(y_{2}+b_{2})}{\|y_{1}+e_{1}\|_{2}\|y_{2}+b_{2}\|_{2}}=\frac{\theta+y_{1}^{T}e_{2}+y_{2}^{T}e_{1}+e_{1}^{T}e_{2}}{\|y_{1}+e_{1}\|_{2}\|y_{2}+e_{2}\|_{2}}\geq\frac{\theta-2\epsilon-\epsilon^{2}}{\|y_{1}+e_{1}\|_{2}\|y_{2}+e_{2}\|_{2}}\geq\frac{\theta-3\epsilon}{(1+\epsilon)^{2}}$. Now consider when $\|y_{1}\|_{2},\|y_{2}\|_{2}\neq 1$. Since cos is scale invariant, we scale the vectors to $(y_{1}+e_{1})/\|y_{1}\|_{2}$ and $(y_{2}+e_{2})/\|y_{2}\|_{2}$ and apply the above statement to obtain the lower bound in our lemma. The upper bound follows similarly. ∎

Corollary A.7**.**

Under the same assumptions as in Lemma A.6 and $\epsilon<\frac{1}{2}$, we have

[TABLE]

Proof.

For the lower bound, use the Taylor expansion of $(1+\epsilon)^{-1}$ to show that

[TABLE]

For the upper bound, consider the Taylor expansion of $(1-\epsilon)^{-1}=\left[1+\epsilon+\epsilon^{2}+\ldots\right]=\left[1+\epsilon(1+\epsilon+\ldots)\right]=\left[1+\epsilon\frac{1}{1-\epsilon}\right]\geq(1+2\epsilon),$ hence $\frac{\theta+3\epsilon}{(1-\epsilon)^{2}}\leq(\theta+3\epsilon)(1+2\epsilon)^{2}\leq(\theta+3\epsilon)(1+8\epsilon)\leq\theta+35\epsilon$. ∎

A.2.3 Clustering and Sorting

Lemma A.8** (Clustering).**

Suppose $G\in\mathbf{R}^{KL\times T}$ is produced from (1) for a CNMF problem with some $H$ that satisfies [$\delta$-sequentially unique]. Consider $\tilde{G}$ such that $\max_{1\leq j\leq KL}\ \|g_{j}-\tilde{g}_{{P(j)}}\|_{2}=\epsilon<\frac{\delta}{70}\|G\|_{-2,row}$, for some permutation $P$, where $g_{j}=G[j,:]$ and $\tilde{g}_{j}=\tilde{G}[j,:]$. Then in $O(TK^{2}L^{3})$ time we can recover $K$ sets $\mathcal{C}_{k}$ of $L$ indices such that for each row $h_{k}$ of $H$, there is some group $\mathcal{C}_{k}$ where for all $\ell=0,\ldots,L-1$, there exists some $j\in\mathcal{C}_{k}$ such that $g_{P^{-1}(j)}=S_{\ell}^{T}h_{k}$.

Proof.

Fix some $k$. Consider $i,j$ such that $g_{i}=S_{\tau_{1}}^{T}h_{k}$ and $g_{j}=S_{\tau_{2}}h_{k}$ for $\tau_{1},\tau_{2}\in\{0,\ldots L-1\}$. Then we must have one of $g_{i}=S_{\ell}g_{j}$ or $g_{j}=S_{\ell}g_{i}$ for some $\ell\in\{0,\ldots,L-1\}$. Without loss of generality, assume $g_{i}=S_{\ell}g_{j}$. Let $r_{i}=\tilde{g}_{{P(i)}}-g_{i}$ and $r_{j}=S_{\ell}^{T}\tilde{g}_{{P(j)}}-S_{\ell}^{T}g_{j}$. By assumption, $\max(\|r_{j}\|,\|r_{i}\|)<(\delta/70)\|G\|_{-2,row}$. Now we apply Corollary A.7 to $\tilde{g}_{i}=g_{i}+r_{i}$ and $S_{\ell}^{T}\tilde{g}_{j}=S_{\ell}^{T}g_{j}+r_{j}$ and consider that $\cos(g_{i},S_{\ell}^{T}g_{j})=1$ to obtain $\cos(S_{\ell}^{T}\tilde{g}_{j},\tilde{g}_{i})>1-\frac{\delta}{2}$, and hence $\cos_{L}(\tilde{g}_{j},\tilde{g}_{i})\geq 1-\frac{\delta}{2}$.

Now consider $i,j$ such that $g_{i}=S_{\tau_{1}}^{T}h_{k}$ for some $\tau_{1}\in\{0,\ldots,L-1\}$, but $g_{j}\neq S_{\tau_{2}}^{T}h_{k}$ for any $\tau_{2}\in\{0,\ldots,L-1\}$. Then by [$\delta$-sequentially unique], we must have that $\max(\cos(g_{i},S_{\ell}^{T}g_{j}),\cos(S_{\ell}^{T}g_{i},g_{j}))\leq 1-\delta$ for any $\ell\in\{0,\ldots,L-1\}$. Applying Corollary A.7 with the same construction as before, we have

[TABLE]

for any $\ell\in\{0,\ldots,L-1\}$. Therefore, $\cos_{L}(\tilde{g}_{1},\tilde{g}_{2})<1-\frac{\delta}{2}$. Together, these two angle bounds gives a strict criterion for clustering the rows of $G$ into $K$ groups of $L$ vectors, according to their corresponding row of $H$. The determine the runtime, we note that we must compute the angle between $KL$ vectors $L$ times each. Computing the angle of each vector requires $\mathcal{O}(T)$ flops, so in total we require $\mathcal{O}(TK^{2}L^{3})$ operations. The grouping of the vectors can be done simply in $\mathcal{O}(K^{2}L)$ and is comparatively negligible as $K<T$. ∎

Lemma A.9** (Sorting).**

Suppose $G\in\mathbf{R}^{KL\times T}$ is produced from (1) for a CNMF problem with some $H$ that satisfies [$\delta$-sequentially unique]. Suppose there exists a permutation $P$ and a matrix $\tilde{G}$ such that

[TABLE]

where $g_{j}=G[j,:]$, $\tilde{g}_{j}=\tilde{G}[j,:]$.

Now suppose we have some set of $L$ indices $\mathcal{C}_{k}$, and that for all $\ell=0,\ldots,L-1$ there exists $j\in\mathcal{C}_{k}$ such that $g_{P^{-1}(j)}=S_{\ell}^{T}h_{k}$.

Then in $\mathcal{O}(TL^{3})$ time we can recover a bijective map $\pi:\mathcal{C}_{k}\rightarrow\{0,\ldots,L-1\}$ such that $g_{P^{-1}(j)}=S_{\pi(j)}^{T}h_{k}$.

Proof.

For any $i,j\in\mathcal{C}_{k}$, a nearly identical argument to Lemma A.8 gives us a decision criterion for finding a unique $\ell_{ij}\in\{0,\ldots,L-1\}$ such that either $g_{P^{-1}(j)}=S_{\ell_{ij}}^{T}g_{i}$ or $g_{P^{-1}(i)}=S_{\ell_{ij}}^{T}g_{j}$. We find such $\ell$ for each pair of indices, and construct a mapping $y(i,j):\mathcal{C}_{k}\times\mathcal{C}_{k}\rightarrow\{-L+1,\ldots,L-1\}$ defined by

[TABLE]

By the construction of $G$ and $\mathcal{C}_{k}$, there is only one $i$ such that $y(i,j)$ will be strictly nonnegative for all $j\in\mathcal{C}_{k}$. Then our permutation is $\pi$ is given by $\pi(j)=y(i,j)$.

Each comparison requires $L$ distance computations, and each distance computations takes $\mathcal{O}(T)$ flops. Since we must measure the distance between $L$ different vectors, this totals to $\mathcal{O}(TL^{3})$ flops. The construction of $y$ and $\pi$ requires $\mathcal{O}(L^{2})$ flops and is comparatively negligible. ∎

A.3 Proof of the Recovery Guarantee in the Presence of Noise

Lemma A.10 proves the recovery guarantee in slightly more details than in Theorem 3. It can then be used to obtain Theorem 3.

Lemma A.10**.**

Suppose a CNMF problem with inputs $\tilde{X},L,K$ is convolutive separable with respect to $\delta,\epsilon,A$. Let $V^{\prime}=VAD_{VA}$ and $G^{\prime}=D_{VA}^{-1}A^{-1}G$ and suppose we know some $t>0$ such that

[TABLE]

where $C_{1}$ and $C_{2}$ are universal constants independent of all other terms, $V$ is defined in (1), and $c=2C_{2}(\epsilon t^{-1})\mathbf{\kappa}(V^{\prime})^{2}\sqrt{KL}$. Then in $\mathcal{O}(NTKL+TK^{2}L^{3})$ time and the time for one NNLS solve, Algorithm 1 finds $\tilde{H}\in\mathbf{R}^{T\times K}$ with bounded error, in the sense that there exists some permutation $P$ such that

[TABLE]

where $h_{j},\hat{h}_{P(j)}$ are the $j$th and ${P(j)}$th rows of $H$ and $\hat{H}$, respectively.

Proof. We break the proof up into several steps.

Location

We know that $V=X[:,\mathcal{C}]A\Pi$ for some index set $\mathcal{C}$, some diagonal scaling matrix $A$, and some permutation matrix $\Pi$. Instead of identifying $V$ directly, we identify $VA=X[:,\mathcal{C}]\Pi$, which satisfies $X=VA[I\ M]\Pi^{\prime}$ for some $M\in\mathbf{R}^{N\times T-R}$ and some permutation matrix $\Pi^{\prime}$. In particular, we apply OrConSPA to inputs $\tilde{X},R,t$ to recover $V^{\prime}=VAD_{VA}$. Since $X=VA[I\ M]\Pi^{\prime}$ and the noise is bounded in (8) and (9) by $\epsilon+t\leq\|VA\|_{-1}$ and $\frac{\epsilon}{t}\leq\frac{C_{1}\sigma_{\min}(V^{\prime})}{2\sqrt{KL}\mathbf{\kappa}(V^{\prime})^{2}}$, we can apply Lemma A.2 and show that the index set $J$ output by OrConSPA satisfies

[TABLE]

Let $\tilde{V}=(\tilde{X}D_{\tilde{X}})[:,J(j)]$.

Perturbation

Given $\tilde{V}$, we estimate $G^{\prime}=D_{VA}^{-1}A^{-1}G$ using nonnegative least squares. Let $F=\tilde{V}-V^{\prime}$. We know the largest singular value is bounded above by the Frobenius norm, so we have $\sigma_{\max}(F)\leq\sqrt{KL}\left(2C_{2}(\epsilon t^{-1})\mathbf{\kappa}(V^{\prime})^{2}\right)=c$. Rearranging (9) using the latter term of the minimum gives us the bound $\sigma_{\min}(V^{\prime})>2C_{2}(\epsilon t^{-1})\mathbf{\kappa}(V^{\prime})^{2}\sqrt{KL}=c$. This allows us to apply Corollary A.5, so the recovered matrix $\tilde{G}$ must satisfy $\max_{i}\ \|G^{\prime}[i,:]-\tilde{G}[i,:]\|\leq\sqrt{T}\left(\frac{\sigma_{\max}(F)\|G^{\prime}\|_{2,col}+\epsilon}{\sigma_{\min}(V^{\prime})-\sigma_{\max}(F)}\right)$. We then have

[TABLE]

Clustering and Sorting

Since $G^{\prime}=D_{VA}^{-1}A^{-1}G$, the angle to a row of $G^{\prime}$ is the same as an angle to a row of $G$. Rearranging (10) we have

[TABLE]

Combining this with (12) lets us apply Lemma A.8 and Lemma A.9 to achieve the correct clustering and sorting. The bound (11) follows automatically. ∎

Proof of Theorem 3.

All we must do is simplify (9) and (10) to match the bounds in Theorem 3. For (9), simply set $C_{a}=\frac{\min(C_{1},C_{2}^{-1})}{2}$. For (10), we begin by recalling our assumption that $\|VA\|_{2,col}>1$ and that $t>1$. We must have that $\|A^{-1}G\|_{-2,row}>1$ since each row has an entry with a 1 in it by construction. Then since $\|VA\|_{2,col}>1$ and equivalently $(D_{VA})^{-1}_{ii}>1$, then $\|G^{\prime}\|_{-2,row}>1$. Given these assumptions, we have the sufficient condition

[TABLE]

Dividing the numerator and denominator by $2C_{2}\mathbf{\kappa}(V^{\prime})^{2}\sqrt{KL}$ leads to $\delta>\frac{70\sqrt{T}(\|G^{\prime}\|_{2,col}(\epsilon t^{-1})+\epsilon)}{\frac{\sigma_{\min}(V^{\prime})}{2C_{2}\mathbf{\kappa}(V^{\prime})^{2}\sqrt{KL}}-(\epsilon t^{-1})}$. Then we have the sufficient condition

[TABLE]

where $\rho=C_{a}\frac{\sigma_{\min}(V^{\prime})}{\mathbf{\kappa}(V^{\prime})^{2}\sqrt{KL}}-(\epsilon t^{-1})$. Since we assumed $t>1$, $\delta>\frac{140\sqrt{T}\|G^{\prime}\|_{2,col}\epsilon}{\rho}$ is also sufficient. Setting $C_{b}^{-1}=140$ gives us $\delta>C_{b}^{-1}\sqrt{T}\|G^{\prime}\|_{2,col}\epsilon/\rho$ which in turn reduces to $\epsilon<C_{b}\delta\rho/\|G^{\prime}\|_{2,col}\sqrt{T}$. ∎

Acknowledgment

The authors thank the reviewers for their feedback which helped improve the paper significantly.