{'fraction_non_alphanumeric': 0.05613124512277008, 'fraction_numerical': 0.014717891947806835, 'mean_word_length': 4.073098876818266, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Memory-Efficient Convex Optimization for Self-Dictionary Separable Nonnegative Matrix Factorization: A Frank-Wolfe Approach

Nonnegative matrix factorization (NMF) often relies on the separability condition for tractable algorithm design. Separability-based NMF is mainly handled by two types of approaches, namely, greedy pursuit and convex programming. A notable convex NMF formulation is the so-called self-dictionary multiple measurement vectors (SD-MMV), which can work without knowing the matrix rank a priori, and is arguably more resilient to error propagation relative to greedy pursuit. However, convex SD-MMV renders a large memory cost that scales quadratically with the problem size. This memory challenge has been around for a decade, and a major obstacle for applying convex SD-MMV to big data analytics. This work proposes a memory-efficient algorithm for convex SD-MMV. Our algorithm capitalizes on the special update rules of a classic algorithm from the 1950s, namely, the Frank-Wolfe (FW) algorithm. It is shown that, under reasonable conditions, the FW algorithm solves the noisy SD-MMV problem with a memory cost that grows linearly with the amount of data. To handle noisier scenarios, a smoothed group sparsity regularizer is proposed to improve robustness while maintaining the low memory footprint with guarantees. The proposed approach presents the first linear memory complexity algorithmic framework for convex SD-MMV based NMF. The method is tested over a couple of unsupervised learning tasks, i.e., text mining and community detection, to showcase its effectiveness and memory efficiency.


I. INTRODUCTION
Nonnegative matrix factorization (NMF) aims at factoring a data matrix X ∈ R M×N into a product of two latent nonnegative factor matrices, i.e., where K ≤ min{M, N }. The NMF technique has been a workhorse for dimensionality reduction, representation learning, and blind source separation. It finds a plethora of applications in machine learning and signal processing; see, e.g., [1], [2]. In particular, NMF plays an essential role in many This work is supported in part by the Army Research Office (ARO) under Project ARO W911NF19-1-0407.
One major challenge of NMF lies in computation. It was shown in [13] that NMF is an NP-hard problem in the worst case, even without noise. In the last two decades, many approximate algorithms were proposed to tackle the NMF problem; see [1], [2], [14]. Notably, a line of work exploiting physically reasonable assumptions to come up with polynomialtime NMF algorithms has drawn considerable attention. To be specific, the so-called separable NMF approaches leverage the separability condition to design efficient and tractable algorithms. More importantly, separable NMF algorithms are often provably robust to noise [3], [15]- [18].
The separability condition is a special sparsity-related condition. Interestingly, it is nonetheless well-justified in many applications. For example, in topic modeling, the separability condition holds if every topic has an "anchor word" that does not appear in other topics [19]. In community detection, separability translates to the existence of a set of "pure nodes" whose membership is only associated with a specific community [7]. The "pure pixel assumption" in hyperspectral imaging is also identical to separability, which means that there exist pixels that only capture one material spectral signature [20], [21]. In crowdsourced data labeling, separability is equivalent to the existence of expert annotators who are specialized for one class [10].
Two major categories of algorithms exist for separable NMF. The first category is greedy algorithms. The representative algorithm is the successive projection algorithm (SPA), which was first proposed in the hyperspectral imaging community [21], and was re-discovered a number of times from different perspectives; see [3], [17], [22]- [26]. Many of these greedy algorithms have a Gram-Schmidt structure, and thus are very scalable. However, they also share the same challenge of Gram-Schmidt, i.e., error propagation. The second category formulates the separable NMF problem as all-at-once convex optimization criteria (see, e.g., [15], [16], [19], [27]- [31]), which are arguably more robust to noise and less prone to error propagation.
Among the all-at-once convex formulations of separable NMF, the self-dictionary multiple measurement vectors (SD-MMV) based framework in [15], [16], [27]- [30] has a series of appealing features, e.g., identifiability of the latent factors without knowing the model rank, computational tractability, and noise robustness. Nonetheless, this line of work has seri-ous challenges in terms of memory. These algorithms induce a variable that has a size of N × N , which is not possible to instantiate if N reaches the level of N = 100, 000-which leads to a 74.5GB matrix if the double precision is used. Consequently, these algorithms are often used together with a data selection pre-processing stage to reduce problem size (see, e.g., [30]). This may again create an error propagation issue-if the data selection stage missed some important data samples (e.g., those related to anchor words or pure nodes in topic modeling and community detection, respectively), the algorithms are bound to fail.
Contributions. In this work, we revisit convex optimizationbased SD-MMV for separable NMF. Our goal is to offer an NMF identifiability-guaranteed algorithm that is also provably memory-efficient. Our contributions are as follows: • A Memory-Efficient Algorithmic Framework. We first show that applying with the standard FW algorithm onto a special SD-MMV formulation for separable NMF admits identifiability of the ground-truth latent factors, if the noise is absent. More importantly, the memory cost of this algorithm is O(KN ) other than O(N 2 ), where K ≪ N often holds. Based on this simple observation, we further show that similar conclusions can be drawn even if noise is present-under more conditions. • Regularization-Based Performance Enhancement. To circumvent stringent conditions and to improve noise robustness, we propose a smoothed group-sparsity regularization that is reminiscent of the mixed norm regularization often used in SD-MMV formulations (see [15], [27], [29], [32], [33]). We show that the optimal solution of such regularized formulation corresponds to the desired ground-truth NMF factors in the noisy case-i.e., identifiability of the NMF holds. We further show that this regularization can better safeguard the memory consumption within the range of O(KN ) compared to the unregularized version, if the FW algorithm is initialized reasonably. To our best knowledge, the proposed FW algorithmic framework is the first convex SD-MMV method whose memory cost scales linearly in N , while existing methods such as those in [15], [16], [27], [28], [30] all need O(N 2 ) memory.
• Synthetic and Real data Evaluation. We test the proposed approach on various synthetic and real datasets. In particular, we evaluate our algorithm using a couple of text corpora, i.e., NIST Topic Detection and Tracking (TDT2) [34] and the Reuters-21578 data 1 on topic mining tasks. We also use social network data from [35], [36] to evaluate the performance on mixed membership community detection tasks. The proposed approach is benchmarked by competitive state-of-the-art baselines for solving separable NMF.
Notation. We follow the established conventions in signal processing. In particular, x ∈ R N denotes a real-valued Ndimensional vector; x n and [x] n both denote the nth element of x; X ∈ R N ×M denotes a real-valued matrix; rank(·) denotes matrix rank; superscript ⊤ is used for transpose; · p where p ≥ 1 denotes the vector ℓ p norm; · F denotes the Frobenius norm; X(n, :) and x ℓ denote the nth row and ℓth 1 http://www.daviddlewis.com/resources/testcollections/reuters21578 column of X, respectively; x n,ℓ , [X] n,ℓ and X(n, ℓ) all denote the (n, ℓ)th element of X; X(−n, :) and X(:, −ℓ) denote the submatrices constructing from X by removing the nth row and ℓth column, respectively; |K| and K c denote the cardinality and complement of the set K, respectively; X row−0 counts the number of nonzero rows of X; X ∞,1 = M i=1 X(i, :) ∞ is the mixed ℓ ∞ /ℓ 1 norm; supp(x) denotes a set of all indices of non-zero elements of vector x, i.e., supp(x) := {i | [x] i = 0}; conv{x 1 , x 2 , . . . , x k } denotes convex hull of set {x 1 , x 2 , . . . , x k }; 1 and 0 denote an allone matrix/vector and an all-zero matrix/vector, respectively, with proper sizes; e i denotes the ith unit vector; I N denotes an identity matrix with a size of N × N ; both notations x ≥ 0 and X ≥ 0 mean that the nonnegativity is applied element-wise; i ∈ [N ] means i ∈ {1, . . . , N }; σ max (X) and σ min (X) denote the largest and smallest singular value of X, respectively; λ max (X) denotes the largest eigenvalue of X; N denotes the set of natural numbers.

II. PROBLEM STATEMENT
Consider a noisy NMF model, i.e., where W ∈ R M×K and H ∈ R K×N are nonnegative latent factors as defined before, and V is a noise term. We further assume that i.e., the columns of H reside in the probability simplex. This assumption naturally holds in many applications, e.g., topic modeling, community detection, and hyperspectral unmixing [6], [19], [20], [37]. When H does not have sum-to-one columns, this assumption can be "enforced" through column normalization of X, under the condition that W is nonnegative, and there is no noise; see [1], [18]. We should mention that although our interest lies in NMF, the proposed method can also be applied to the so-called simplex-structured matrix factorization, where W is not required to be nonnegative; see, e.g., [8], [9], [25], [37]. Estimating the ground-truth W and H from X is in general an NP-hard problem [13]. However, if the so-called separability condition holds, the separable NMF algorithms are often tractable, even under noise.

A. Separable NMF
The separable NMF problem uses the following premise: The condition was first seen in [38] in the NMF literature. The remote sensing community discovered it even earlier [21], where the same condition is called the "pure pixel condition" [20]. The condition has many names in applications, e.g., = X C W H 0 0 0 X Fig. 1: Visualization of the idea of the row-sparsity-based SD-MMV in [17], [27].
Under separability, the NMF task boils down to finding K, since W ≈ X(:, K) if the noise level is not high. Then, H can be estimated by (constrained) least squares if rank(W ) = K.

B. Convex Separable NMF
A way to look at the K-finding problem is to cast it as a sparse atom selection problem. To be specific, when noise is absent, consider the following row-sparsity minimization formulation: where C row−0 counts the nonzero rows of C. For example, if K = {1, . . . , K}, then an optimal solution of (4) is under mild conditions (e.g., rank(W ) = K). More formally, we have: see a proof for the general case in [17] and an illustration in Fig. 1. Hence, K can be identified by inspecting the nonzero rows of C ⋆ . The formulation in (4) is reminiscent of the multiple measurement vectors (MMV) problem from compressive sensing [40], [41], but using the data itself as the dictionary-which is the reason why (4) is called selfdictionary multiple measurement vectors (SD-MMV) [17]. We should mention that SD-MMV can be regarded as a way of picking up the vertices of the convex hull of {x 1 , . . . , x N }, which is a popular perspective that many separable NMF algorithms take; see, e.g., [3], [23], [24], [42]. Unlike the classic vertex picking methods that often need the knowledge of K, SD-MMV can work without knowing the number of vertices.
The formulation in (4) is not easy to tackle, due to the combinatorial nature of C row−0 . One popular way is to use greedy pursuit, which selects the basis that represents X using its convex combinations from the self-dictionary X one by one. This leads to the successive projection algorithm (SPA) [17], [18], [21], [23]. The greedy pursuit methods are often effective, efficient, and robust to a certain extent. However, they also share the same challenge-the error could be accumulated and propagated through the greedy steps.
Another line of work employs the convex relaxation idea and use a convex surrogate for C row−0 -see [15], [27]- [30], [32], [33], [43] for different options. For example, the work in [15], [27] uses as their convex surrogate. When the noise is present, the work in [15], [27] also advocated the following formulation where λ ≥ 0 is a regularization parameter. The formulation is convex and thus is conceptually easy to handle. After using any convex optimization algorithm to obtain an optimal solution C, the K set can be identified via inspecting C(n, :) ∞ for all n and picking the K indices that are associated with K largest row ℓ ∞ -norms. This method was shown to have identifiability of K under noise [15]-also see similar results for close relatives of (8) in [16], [28], [30], [43].
In terms of noise robustness, the convex optimization-based SD-MMV methods are arguably more preferable over greedy methods, since they identify all elements of K simultaneously, instead of in a one by one manner that is prone to error propagation.

C. The Memory Challenge
Using convex programming to handle separable NMF is appealing, since many off-the-shelf algorithms are readily available. The convergence properties of convex optimization algorithms are also well understood. However, general-purpose convex optimization algorithms, e.g., the interior-point methods and gradient descent, may encounter serious challenges when handling Problem (8). The reason is that the C variable induces O(N 2 ) memory if a general-purpose solver is used. This often makes running such algorithms costly and slow, even when N is only moderate.
Many attempts were made towards accelerating convex SD-MMV. The early work in [28] uses a fast projection to expedite the algorithm for a variant of (8). The recent work in [30] employed the fast gradient paradigm for accelerated convergence. However, the O(N 2 ) memory issue could not be circumvented in both approaches. In this work, our objective is a convex optimization algorithm that share similar identification properties of those in [15] but has provable memory-efficiency-i.e., only O(N K) memory is needed for instantiating the optimization variable C.

A. Preliminaries on Frank-Wolfe Algorithm
Our development is based on the idea of the Frank-Wolfe (FW) algorithm that was developed in the 1950s [44]. The FW algorithm deals with problems of the following form: where C is compact and convex and f (·) : R d → R is differentiable and convex. The FW algorithm uses the following updates: Here, α t is a pre-defined sequence, i.e., The above two steps constitute a standard procedure of the FW algorithm [45], [44]. The idea of FW is intuitive: In each iteration, FW linearizes the cost function and solve it locally over the compact constraint set. For convex problems, the plain-vanilla FW algorithm converges to an ε-optimal solution using O(1/ε) iterations (which is also known as a sublinear convergence rate). Recent works showed that FW (and its variants) converges even faster (with a linear rate) under some more conditions; see, e.g., [45], [46]. When dealing with large-scale optimization, especially memory-wise costly optimization problems, FW may help circumvent memory explosion due to its special update rulewhich has already facilitated economical nuclear norm optimization paradigms that are important in recommender systems [45]. FW also features simple updates if the constraint set is In this case, u in (10) is always a unit vector. This is because (10a) is a linear program over the probability simplex C, and the minimum is always attained at a vertex of C [47]. The vertices of C are {e n } N n=1 . Hence, the solution of (10a) is e j where j = arg min In our algorithm design, we will take advantage of this fact to reduce the memory cost of solving SD-MMV. Note that one needs not always start the FW algorithm from t = 0. In [46], a warm-start based Frank-Wolfe (WS-FW) algorithm was proposed for accelerated convergence. There, one can use θ tinit = θ init with t init ≥ 1 and start the FW algorithm from the t init th iteration using the corresponding α tinit as defined in (11). By carefully choosing t init according to the quality of θ init and the problem structure (e.g., the curvature of the cost function), WS-FW is shown to converge to the global optimum more efficiently.

B. Warm-up: The Noiseless SD-MMV
Our algorithm design starts with the simple formulation as follows: Note that in general, (13) does not admit identifiability of Kthe optimal solution of (13) needs not to have the form in (6). A simple counter example is C = I N -which gives zero objective value, but is not the desired solution. That is, one cannot infer K from the nonzero rows of C = I N .

C. Noiseless Case and Simple Self-Dictionary Fitting
The first observation is that although the criterion in (13) does not admit identifiability of K, the FW algorithm guarantees finding a C as in (6), thereby pining down K.
To see this, note that the updates of c ℓ 's do not affect each other. Consider the following update rule for ℓ = 1, . . . , N : Since the formulation is convex, the above updates are guaranteed to solve the problem in (13). Before we examine what solution the updates will lead to, let us discuss its memory complexity, i.e., the reason why such updates are potentially memory-economical. To see how much memory that the algorithm needs beyond storing X (which is inevitable), let us analyze the steps in (14). First, evaluating gradient in (14a) is not memory-costing. One can always evaluate c t ℓ − x ℓ , followed by multiplying with X, followed by another left multiplying X ⊤ . These three operations produce vectors as their results and hence require O(N ) memory complexity in total. Note that one can evaluate the gradient w.r.t. C column by column, and no previously evaluated columns need to be stored. Hence, the total memory cost is O(N ). Second, the remaining memory requirement lies in storing the iterates c t+1 ℓ for ℓ = 1, . . . , N . Consider an ideal case where the j found in (14b) satisfies j ∈ K for all ℓ and all t. Then, if the initialization C 0 = 0, the updating rule in (14b) always maintains supp(c t ℓ ) ⊆ K for all ℓ = 1, . . . , N . This leads to an O(KN ) memory for storing C t .
Based on our discussion, the key to attain the above described memory efficiency is to ensure that j ∈ K for all ℓ and all t in (14b). We show that this is indeed the case under some conditions: Theorem 1 (Noiseless Case, Memory Efficiency) Suppose that Assumption 1 holds, that rank(W ) = K, that no repeated unit vectors exist in H 2 , and that the noise is absent . Assume that the following holds: before c t ℓ reaches an optimal solution. Then, the FW algorithm in (10) with initialization C 0 = 0 always outputs the desired solution in (6) using O(KN ) memory beyond storing X.
The proof is relegated to Appendix A. Note that Theorem 1 requires that Eq. (15) holds-i.e., no repeated minimum-value elements of q ℓ exist for any c t ℓ . This is not hard to meet under some mild conditions:

Fact 1 Assume that the columns of W follow any joint continuous distribution. Then, we have
Proof: Under the assumption, suppose that there exists any b such that we have Pr(w ⊤ i z = b, w ⊤ j z = b) > 0 for any vector z. That is, there is a positive probability that q ℓ could have two identical elements. However, this cannot happen. Indeed, note that y i = w ⊤ i z and y j = w ⊤ j z are also continuous random variables. Hence, and applying the above result complete the proof. Although Theorem 1 is concerned with the ideal noiseless case, which is hardly practical, the observation in Theorem 1 opens another door for self-dictionary convex optimizationbased NMF. That is, it shows that using FW to solve the selfdictionary problem may successfully circumvent the memory explosion issue.

D. The Noisy Case
A natural question lies in if the same procedure works in the presence of noise. The answer is affirmative. To understand this, we first notice that in the noiseless case, we have where c ⋆ ℓ is the ℓth column in C ⋆ defined in (6). This further leads to C(n, :) ∞ = 1, ∀n ∈ K, C(n, :) ∞ = 0, ∀n ∈ K c , which is the reason why one can easily infer K from the solution C = C ⋆ . When noise is present, FW may not be able to find the exact C ⋆ . However, an approximate C ≈ C ⋆ often suffices to identify K. In the following, Lemma 1 shows that instead of finding c ℓ such that Hc ℓ = h ℓ , seeking a c ℓ such that Hc ℓ ≈ h ℓ using FW does not hurt recovering K, under reasonable conditions: Lemma 1 Assume that the separability condition holds, and that rank(W ) = K, and that no repeated unit vectors appear in H. Suppose that a feasible solution C satisfies supp(c ℓ ) ⊆ K and Hc ℓ − h ℓ 2 ≤ ω for all ℓ ∈ [N ]. Then, we have i.e., K can still be identified by inspecting C(n, : On the other hands, since supp(c ℓ ) ⊆ K and C is initialized at 0, we have This completes the proof. Next, we show in the following theorem that even in the noisy case, applying FW onto (13) produces a solution that is a reasonable approximation of C ⋆ . To proceed, let us define the following quantities: Definition 1 Define the following terms: In particular, a small d(H) means that all h ℓ 's for ℓ ∈ K c are sufficiently different from the unit vectors. This is a desirable case, since small perturbation would not confuse such h ℓ 's and the unit vectors.
Theorem 2 (Noisy Case, Memory Efficiency) Suppose that Assumption 1 holds, that there is no repeated unit vector in H, and that rank( During FW's updates, assume a positive gap between the smallest and second smallest elements of q t exists, i.e., before the FW algorithm terminates. Also assume that for some η > 0. Then, if one terminates FW when x ℓ − Xc t ℓ 2 ≤ η + 2δ, the algorithm produces a solution C that satisfies In addition, during the process, supp(c t ℓ ) ⊆ K always holds for all t and ℓ, and thus only O(KN ) memory is taken for instantiating C t for all t.
The proof can be found in Appendix B. Theorem 2 shows that under certain conditions, the FW algorithm for solving (13) indeed gives a reasonable solution C using only O(N K) memory. Combining with Lemma 1, one can see that K can be correctly selected if the noise level is not high. Both Theorems 1 and 2 reveal an unconventional identifiability result. Note that the formulation in (13) per se does not have identifiability of C ⋆ . That is, the optimal solutions of (13) do not necessarily reveal K-as mentioned, I N is also an optimal solution. However, when one uses a particular algorithm (i.e., FW) to solve (13), the produced solution sequence converges to a K-revealing C, even if noise is present.
On the other hand, the identifiability and memory efficiency come with caveats. First, the condition in (16) is hard to check or guarantee. As seen in Fact 1, ν > 0 does exist under mild conditions-but the quantity of this parameter may change from instance to instance and is hard to acquire. Using the formulation in (13) is also not the most "natural", since we know that the desired C ⋆ should be very row-sparse-why not using this information for performance enhancement. Can we get rid of the condition in (16) and effectively use the prior knowledge on C ⋆ ? The answer is affirmative, with a re-design of the objective function. In the next section, we will discuss these aspects.

IV. PERFORMANCE ENHANCEMENT VIA REGULARIZATION
Under the formulation in (13), the gap specified in (16) was essential for the FW to pick j ∈ K in every step. In this section, our idea is to use a regularization term to help the FW algorithm to achieve the same goal while not relying on the condition in (16). In addition, since the desired C ⋆ in (6) has a row-sparse structure, it is natural to add a regularization to exploit this prior information.
We still hope to use a regularization term like C ∞,1 , which was shown to have nice identifiability properties in SD-MMV [15], [27]. To make this mixed-norm based nonsmooth row-sparsity regularization "FW-friendly", we use the following lemma: The above smoothing technique is from [49]. A proof is presented in Appendix K in the supplementary material for completeness. Building upon Lemma 2, a smooth approximation of C ∞,1 is readily obtained as Using Φ µ (C), our working formulation is as follows: The formulation can be understood as a smoothed version of those in [15], [27], [29], [32]. Note that the problem in (18) is still convex, but easier to handle by gradient-based algorithms relative to the nonsmooth version.

A. Identifiability of K
Our first step is to understand the optimal solution of (18)i.e., if one optimally solves Problem (18), is the solution still K-revealing as in the nonsmooth version? To this end, we will use the following definition [15], [16]: The quantity κ(W ) is defined as follows: The term κ(W ) in a way reflects the "conditioning" of W . If κ(W ) is large, it means that every w k is far away from the convex hull spanned by the remaining columns of W , which implies that conv{w 1 , . . . , w K } is well-stretched over all directions. This is a desired case, since such convex hulls are more resilient to small perturbations when estimating the vertices (i.e., w k for k = 1, . . . , K).
With this definition, we show that the optimal solution of (18) can reveal K under some conditions: Theorem 3 (Identifiability) Assume that Assumption 1 holds, that there is no repeated unit vector in H, and that rank(W ) = K. Also assume that v i 2 ≤ (ρ/N ) V 2 F for some ρ and v i 2 ≤ δ. Then, any optimal solution C opt of Problem (18) satisfies: The proof is relegated to Appendix C.
Theorem 3 reveals the "interplay" between the noise level and the hyperparameters λ, µ. In other words, it states that given a certain noise level, there may exist a pair of λ, µ such that K will be identified using the proposed FW algorithm. For example, when the noise level is not high, a natural choice is to concentrate more on the fitting error rather than the regularization term. This is reflected in conditions (20a) and (20b). Since V 2 F is small, δ would be also small. Then, with a small λ to suppress the term 2λK in the expression of β, a small β can be expected. Similarly, with a small µ to suppress the term µN log N in (20b), the right hand side (RHS) of (20b) would be close to 0. In addition, when the noise level is relatively high, one would want to increase λ to counter the effect of increasing V 2 F in (20b), but not to increase it to an overly large extent (in order to keep β close to 0, due to the presence of 2λK in the expression of β).
Theorem 3 asserts that with the proposed regularization, finding an optimal solution of (18) is useful for identifying K. Such optimal solutions can be attained by any convex optimization algorithm. Nonetheless, theorem 3 only safeguards the final solution of our formulation in (18), which is a similar result as in [15] for the nonsmooth regularization version (8). However, Theorem 3 does not speak for the memory efficiency. (18) In this subsection, we show that running FW to optimize the new objective in (18)  In the regularized case, the FW updates become the following:

B. Memory Efficiency of Frank-Wolfe for
where g ℓ is the ℓth . Similar as before, the key to establish memory efficiency is to show that j picked in (21) always belongs to K. The next theorem shows that the regularization helps achieve this goal: Theorem 4 (Regularized Case, Memory Efficiency) Consider the regularized case and the FW algorithm in (21). Assume that Assumption 1 holds and that rank(W ) = K. Also assume that (i) in iteration t, C t satisfies supp(c t ℓ ) ⊆ K for all ℓ ∈ [N ] and (ii) there exists at least an n ⋆ ∈ K such that C t (n ⋆ , :) is not a constant row vector. If the noise bound satisfies where The proof can be found in Appendix D. The noise bound in Theorem 4 is arguably more favorable relative to the unregularized case in Theorem 2. The reason is that λ can be tuned to compensate noise. This is also intuitive since a larger λ means that one has lower confidence in the data quality due to higher noise.
The key condition in Theorem 4 is the existence of C t (n ⋆ , :) that is not a constant-which is reflected in ψ. At first glance, this is hard to guarantee since it is a characterization of C t . Nonetheless, we show that if one can properly initialize the FW algorithm with an initial solution C init that satisfies a certain regularity condition, the non-constant condition is automatically satisfied-due to the "predictable" update rule of FW: for some t init ∈ N, t init ≥ 1. Suppose that for some n ⋆ ∈ K, there exists a pair of i ⋆ , j ⋆ such that the following regularity condition is satisfied: Running WS-FW with C t = C init starting with t = t init for T iterations, the produced solution sequence {C} t for t ≥ t init + 1 satisfies C t (n ⋆ , :) is not a constant row vector and min The proof is in Appendix J. Theorem 4 and Proposition 1 together mean that for a given problem instance and under a certain noise level, there exist proper parameters λ, µ that guarantee the recovery of K using O(KN ) memory to instantiate C t . Simply speaking, Proposition 1 asserts that, if there exists a row in C init (n ⋆ , :) where n ⋆ ∈ K and this row has two elements whose difference is not a natural number, then Theorem 4 holds with in finite iterations. Under Theorem 4 and Proposition 1, it is natural to run the WS-FW with a t = t init > 0 as in [46]. Many lightweight algorithms (e.g., the greedy methods) can be employed to provide the initialization.
The condition in Proposition 1 is fairly mild, since it boils down to the existence of two distinctive elements in any row of the initial C. The Proposition also suggests that using some existing algorithms to initialize the proposed FW algorithm may be appropriate, since one needs at least one n ⋆ ∈ supp(C init ) ∩ K such that the specified conditions are met. Any greedy algorithm, e.g., [17], [18], could help offer this n ⋆ using a couple of iterations. Although checking the initialization condition is easy, we should also mention that it may not be necessary to really check it in practice. The reason is that the regularity condition in Proposition 1 is only sufficient-which means that in practice one often needs not to enforce C init to satisfy it. In fact, the FW algorithm works well and maintains a low memory footprint even using C init = 0.
To summarize, we present the memory-efficient Frank-Wolfe based nonnegative matrix factorization (MERIT) algorithm in Algorithm 1. An implementation of MERIT can be downloaded from the authors' website 3 .

V. NUMERICAL RESULTS
In this section, we use synthetic and real data experiments to showcase the effectiveness of the proposed FW-based approach.

A. Synthetic Data Simulations
We create synthetic data matrices with different M, N and K, under the noisy signal model X = W H + V . The 3 https://github.com/XiaoFuLab/Frank-Wolfe-based-method-for-SD-MMV matrix W is drawn from the uniform distribution U(0, 1), the first K columns of H are assigned to be an identity matrix (which means that K = {1, . . . , K}), the N − K remaining columns of H are generated so that every column resides in the probability simplex (see details later). After adding zeromean σ 2 -variance Gaussian noise V to W H, the columns of the data matrix are then permuted randomly to obtain the final X-this means that K for each random trial is different. The signal-to-noise ratio (SNR) used in this section is defined as SNR= 10 log 10 ( Baselines and Metric. We use the SPA algorithm [17], [18], [21] that is the prominent greedy algorithm for Kidentification in separable NMF as our baseline. We also employ the FastGradient algorithm [30] that is designed to solve a convex self-dictionary formulation of SD-MMV. The algorithm uses accelerated gradient for fast convergence, and is considered state-of-the-art. We apply the algorithms to the problem instances and select K from their outputs as follows. For SPA, we use the first K indices output by the greedy steps to serve as K. For FastGradient, we use the authors' implementation and its default methods to pick up K. Following Theorem 3, we select indices of K rows that have K largest C(n, :) ∞ values as an estimation of K.
To select the hyperparameter λ of MERIT, we use an idea similar to that in [30], which is a heuristic that selects λ to balance the data fitting residue and the regularization term. In our case, the suggestion is to set λ = X − XC 0 F /Φ µ (C 0 ) (or simply λ = X − XC 0 F /K when K is known), where C 0 is an initial solution that can be constructed by some fast separable NMF algorithms, e.g., SPA. The parameter µ is set to be 10 −5 throughout this section unless otherwise specified, as it is fairly inconsequential. The hyperparameter of FastGradient is chosen by its default heuristic.
We use a number of metrics to evaluate the performance. In particular, we primarily use the success rate of K identification, which is defined as success rate = Pr(K = K).
In our simulations, the success rate is estimated using 50 random trials. We also adopt two complementary metrics from [30], namely, the mean-removed spectral angle (MRSA) and the relative approximation error (RAE). In a nutshell, the MRSA measures how well W is estimated via W = X(:, K) and the RAE measures how well the estimated W (together with an estimated H) can reconstruct the data X; see details in [30]. Following [30], the MRSA values are normalized to an interval of [0, 100]. In addition, the RAE values are in between 0 and 1. Lower MRSAs and higher RAEs correspond to better performance of separable NMF. Results. We first evaluate the algorithms under a setting from [30,Sec. 4.1], where H(:, K c )'s columns are middle points between the extreme points of the unit simplex. This way, noise could easily confuse the x ℓ 's associated with the middle points with the true extreme points of conv{x 1 , . . . , x N }thereby presenting a challenging case for separable NMF algorithms.   [30], where (M, N, K) = (50, 55, 10). One can see that FastGradient and MERIT both exhibit more satisfactory performance relative to SPA, which echos our comment that the all-at-once convex approaches often have better noise robustness relative to greedy pursuit. In particular, FastGradient and MERIT reach 100% success rates at SNR= 16dB and SNR= 10dB, respectively, while SPA does not reach this accuracy even when SNR= 20dB. Fig. 3 shows the MRSAs and RAEs of the algorithms under the same setting, where similar observations are made. Fig. 4 (a) shows the success rates of the algorithms under different SNRs using (M, K, N ) = (50,40,200) and H(:, K c )'s that are less special than that in the previous case; i.e., H(:, K c )'s columns are generated following the uniform Dirichlet distribution with its parameter being 1 ∈ R K . One can see that the algorithms perform similarly as in the case of Fig. 2, except that the gap between SPA and the convex approaches MERIT (with regularization) and FastGradient is larger than that in the previous case. Fig. 4 (b) shows the memory costs of the two convex optimization-based algorithms under different N 's. Here, we set (M, K) = (50, 40) and SNR= 10dB. The memory is measured in terms of maximum resident set size (RSS), which is the amount of allocated memory in RAM for a running pro-     gram. The RSS is measured using a Linux built-in command named time 4 . One can see that MERIT's memory growth along with N is very graceful, but FastGradient quickly reaches the level that is close to 5GB when N = 10, 000while MERIT uses less than 0.1GB memory under the same problem size. Table I shows the performance of the algorithms under various K's. As expected, SPA works better when K is relatively small. The performance deteriorates when K increases, showing the effect of error accumulation. The convex approaches work similarly and exhibit consistently good performance across all the K's under test. Fig. 5 shows the impact of the hyperparameters λ and µ on the MERIT algorithm. One can see that for lower SNRs, a larger λ often works better-which is consistent with our analysis and intuition. The parameter µ is less consequential. That is, the wide range of µ's tested in our simulations give almost the same success rate curves.

B. Real Data Experiment: Topic Modeling
Data. We use two popular datasets, namely, the NIST Topic Detection and Tracking (TDT2) and the Reuters-21578 cor-pora, for the evaluation. Following the settings in topic modeling papers, e.g., [4], [50], we use single-topic documents so that classic evaluation metrics (e.g., clustering accuracy [50], [51]) can be easily used. TDT2 contains D = 8, 384 singletopic documents with N = 36, 771 words as its vocabulary, while Reuters-21578 contains D = 8, 293 single-topic document with N = 18, 933 words in its vocabulary. All stop words are removed before running each trial.
Note that we follow the formulation in [3]- [5] that applies NMF in the correlation domain. That is, the wordword correlation (or co-occurrence) data matrix X has a size of N × N . The matrix X is estimated by the Gram matrix of the word-document term-frequency-inverse-documentfrequency (TF-IDF) representation of the data; see more details in [4], [5]. Following the work in [17], a pre-processing step for noise reduction is used in this subsection. Specifically, before running the separable NMF algorithms, the principal component analysis (PCA)-based dimensionality reduction (DR) is used to reduce the row dimension of the co-occurrence matrix to M ′ = 2K, which serves as an over estimate for K.
In practice, such DR method can be easily implemented with any M ′ > K if K is roughly known. After the DR process, one factor of the dimension-reduced factorization model may not be nonnegative. Thus, technically, the model is not "NMF". Nonetheless, all the SD-MMV methods can still be applied since the nonnegativity of the left factor (W ) is never used in attaining K. Baseline and Algorithm Settings. We use a number of separable NMF based topic modeling algorithms as our baselines in this experiment. In addition to SPA and FastGradient, we also use XRAY [52] and FastAnchor [3] that are both greedy algorithms and variants of SPA developed under the context of topic modeling. The LDA [53] method using Gibbs sampling [54] is also employed, as a standard baseline.
We use warm start for both MERIT and FastGradient. Particularly, we use SPA to extract an estimate of K, denoted as K SPA . Then, we compute C init by C init (K SPA , :) = arg min and let C init ( K c SPA , :) = 0. For MERIT, we set λ = 10 −6 and µ = 10 −5 in all cases. We try multiple λ's for FastGradient and present the best performance it attains. For MERIT, since we use WS-FW (cf. Sec. III-A), the t init needs to be determined. In the WS-FW paper [46], t init in theory is computed using the curvature constant associated with optimization problem and initial duality gap. However, as admitted in [46], estimating the curvature constant is still an open challenge. Hence, we use a heuristic that t init = round(1/RMSE init ), where RMSE init = 1 /N X − XC init 2 F . This reflects the idea that a better initialization should use a larger t init . Metrics. We use the three metrics from [4], [5], namely, coherence (Coh), word similarity count (SimCount), and clustering accuracy (ClustAcc). The Coh metric evaluates if a single topic's quality by measuring the coherence of the high-frequency words contained in this topic. The SimCount metric measures how diverse are the mined topics. The ClustAcc compares the estimated H with the ground-truth labels after automatic permutation removal using the Kuhn-Munkres algorithm; see more details of the evaluation process in [1], [4] and the definition of ClustAcc in [51, Sec. 5.2].. A good topic mining algorithm is expected to attain high Coh values, low SimCount and high ClustAcc values. Among the three, ClustAcc is arguably the most indicative for the quality of the mined topics if the objective is to use the topics for downstream tasks. For each trial, we randomly draw documents associated with K topics from the datasets and apply the algorithms. The results for each K are averaged from 50 trials. Results. Tables II shows the performance of the algorithms on  TDT2 and Reuters-21578, respectively. One can see that both the regularized and unregularized versions of the proposed method, i.e., MERIT and MERIT(0) (i.e., the version of MERIT with λ = 0), exhibit competitive performance. The proposed methods outperform all the baselines on TDT2 in terms of Coh and ClustAcc. The SimCount performance of the proposed algorithms is also reasonable. In particular, when K = 10, the ClustAcc of MERIT exhibits a 5% improvement compared to the best baseline, which is a remarkable margin. On Reuters-21578, the MERIT method also consistently offers the best and second best performance in terms of ClustAcc in most cases.
From these results, one can see that the convex optimization based separable NMF algorithms, i.e., FastGradient, MERIT and MERIT(0), in general work better than the greedy methods, namely, SPA, XRAY, and FastAnchor. This is consistent with our observation in the synthetic data experiments. This advocates using such all-at-once algorithms for real applications. Another observation is that MERIT slightly outperforms MERIT(0), which shows that the designed regularization is still effective on topic modeling problems. Fig. 6 shows the memory consumption of the algorithms on TDT2 and Reuters-21578, respectively. Since the data matrix's size changes in each trial due to the varying stop words, we only plot the first trial for each K. One can see that when N reaches 20, 000, FastGradient uses more than 20GB memory, while MERIT and MERIT(0) use less than 2GB. We observe that MERIT(0) works well in the topic modeling experiment, perhaps because the data model is reasonably well aligned with that in (2) without much modeling error.

C. Real Data Experiment: Community Detection
The link between mixed membership stochastic model (MMSB)-based overlapped community detection and simplexstructured matrix factorization has often been used in the literature; see [6]- [9]. By applying eigendecomposition to a nodenode undirected adjacency matrix, the membership estimation problem boils down to a noisy matrix factorization under probability simplex constraints. This leads to X ≈ W H, where X ∈ R K×N is extracted by the eigendecomposition of the N × N adjacency matrix (the rows of X are the first K eigenvectors), H's columns are the mixed membership vectors correspond to the N node. More precisely, h k,ℓ is the probability that node ℓ is associated with community k. This physical interpretation also means that 1 ⊤ H = 1 ⊤ and H ≥ 0. Hence, both SPA and convex optimization based separable NMF algorithms can be applied to tackle this problem, if the separability condition holds. In the context of community detection, separability is equivalent to the existence of "pure nodes" for each community, i.e., nodes who are only associated with a single community. Data. We use two co-authorship networks, namely, the Data Base systems and Logic Programming (DBLP) data and the Microsoft Academic Graph (MAG) data. A community in DBLP is defined as a group of conferences. The "field of study" is considered as a community in MAG. The groundtruth memberships of the nodes in these two datasets are known; see the Matlab format of this data from [7].
In our experiments, we consider the nodes who contribute to 99% of the energy (in terms of squared Euclidean norm). The remaining nodes are not used in the algorithms due to their rare collaboration with others. The detailed statistics of the used DBLP and MAG data are given in Table III. Baselines. We compare MERIT and MERIT(0) with three baselines. The first two algorithms are GeoNMF [7] and SPOC [8] as they are reportedly popular within the class of greedy method in context of community detection. And FastGradient as a candidate for the convex optimiza-  tion based approach is included. MERIT uses the same hyperparameters setting as in the topic modeling problem, i.e., λ = 10 −6 , µ = 10 −5 . Again, we try multiple λ's for FastGradient and report its best performance. Metric. Following [7]- [9], we evaluate the performance based on the averaged Spearman's rank correlation (SRC) coefficient between the learned community membership matrix H and the ground-truth H. By the definition, the SRC measures the ranking similarity between the estimated and ground-truth mixed membership vectors. The estimated H is obtained by probability simplexconstrained least squares using data X and the basis W = X(:, K). By definition, SRC can take value in the interval from −1 to 1. A higher value indicates a better alignment between H and H. Results. From Table IV, one can see that the proposed MERIT method offers competitive SRC results over all the datasets under test. Similar as before, the convex all-at-once algorithms, especially, MERIT, and MERIT(0), often output competitive SRC results. In particular, the MERIT class' SRC values are in top 2 over five out of seven datasets. The greedy methods (i.e., GeoNMF, and SPOC) also work reasonably well, but less competitive in some datasets. For example, in MAG2 and DBLP2, the performance gap between the greedy algorithms and the convex methods are particularly articulate. Fig. 7 compares the memory consumption of FastGradient and MERIT. As N reaches 6500, FastGradient consumes more than 2GB while MERIT and MERIT(0) use significant less memory, i.e., under 0.3GB. This is consistent with our observations in the topic modeling examples.
A final remark is that the theorems in this work are based on worst-case analyses, and thus the noise bounds and identification conditions are naturally pessimistic. However, note that our conditions are only sufficient (instead of sufficient and necessary). This may explain the reason why the proposed method works under many settings of the experiments where the noise level is quite high.

VI. CONCLUSION
In this work, we revisited convex optimization-based selfdictionary SD-MMV for separable NMF. This line of work emerged about a decade ago as a tractable and robust means for solving the challenging NMF problem. The method is recognized as an important category of NMF algorithms, but has serious challenges when dealing with big data. In particular, existing convex SD-MMV approaches' memory complexity scales quadratically with the number of samples, which is hardly feasible for datasets that have more than 1000 samples. We proposed a new algorithm based on the Frank-Wolfe method, or, conditional gradient. Unlike existing algorithms, our method is shown to have linear memory complexity under mild conditions-even in the presence of noise. For performance enhancement, we also offered a smoothed rowsparsity-promoting regularizer, and showed that it can provide stronger guarantees for memory efficiency against noise under the FW framework. We tested the algorithm using synthetic data and real-world topic modeling and community detection datasets. The results corroborate our theoretical analyses. PROOF OF THEOREM 1 i.e., the cost function is decomposable over the columns of X, one can consider each f (c ℓ ) = 1 /2 Xc ℓ − x ℓ 2 2 individually. The FW algorithm first finds u such that ∇ c ℓ f (c ℓ ) ⊤ u is minimized over the probability simplex. This is a linear program, and its solution is always attained at a vertex of the probability simplex [47]. Hence, the solution is an unit vector, i.e., u = e n ⋆ , n ⋆ = arg min [∇f (c ℓ )] n .
Next, we show that n ⋆ ∈ K always holds. To see this, we have If q ℓ = 0, then the lower bound can be attained if h n = e j ⋆ where j ⋆ = arg min j∈[K] q j,ℓ . Note that such h n always exists because of the separability assumption. Let us denote it as h n ⋆ , we have n ⋆ = arg min n [∇f (c ℓ )] n . Since h n ⋆ is a unit vector, and because of the assumption that non-repeated unit vectors appears in H, one can conclude n ⋆ ∈ K.
If there are more than one smallest element in q ℓ , say, there are 2 smallest elements q 1,ℓ = q 2,ℓ = min j q j,ℓ , then, any h n ⋆ in a form of h n ⋆ = z, 1 − z, 0, . . . , 0 ⊤ for any 0 ≤ z ≤ 1 also makes the lower bound in (24) attained. Such h n ⋆ might be not a unit vector and hence n ⋆ / ∈ K. Nonetheless, the assumption that no duplicate minimal values appear in the gradient [cf. Eq. (15)] assures that this case never happens.
If q ℓ = 0, which means that the objective function already reaches optimal value of 0 since Therefore, one should stop FW.

APPENDIX B PROOF OF THEOREM 2
Under the noisy model, for any c ℓ , we have Hence, Similar to the proof of Theorem 1, the key is to show that for Since the proof in the sequel holds for all ℓ, we omit the subscript and use n ⋆ instead. From the proof of Theorem 1, when ǫ n = 0, we know that h n ⋆ = e j where e j is a certain unit vector. Our goal is to show that with such n ⋆ , for any n = n ⋆ , Equivalently, we hope to show the following which will leads to that the FW algorithm selects n ⋆ at the current iteration with c ℓ . We will prove (26) using the following two Lemmas. Their proofs are provided in Appendix E and Appendix F in the supplementary material.
Lemma 4 For any ℓ ∈ [N ], n ∈ [N ], n = n ⋆ , where Using upper bound of noise δ, we have Using (27), (28), and (29), one can see that which is exactly (26). Until now, we have established that the FW algorithm will not update any C(n, :) such that n ∈ K c before it reaches the stopping criterion. Whenever FW terminates at C, we have This completes the proof.

APPENDIX C PROOF OF THEOREM 3
Without loss of generality, assume that K = {1, . . . , K} and H(:, 1:K) = I. Therefore, we have that is the desired solution. Note that C ⋆ is an optimal solution of Problem (4). Our goal is to show that Problem (18)'s optimal solutions are close to C ⋆ . We will show this step by step. a) Step 1: We first find an upper bound of the objective function associated with C ⋆ . For 1 ≤ ℓ ≤ N , it can be seen that which leads to where the first equality holds because we have x ℓ = Xc ⋆ ℓ for 1 ≤ ℓ ≤ N .
By Lemma 2, we have Combining (30) and (31), we have Consider any solution C such that Therefore, we have for any k ∈ K. b) Step 2: Following the idea in [16,Lemma 17], define τ = H(k, :)c k . It is readily seen that 0 ≤ τ ≤ 1. Suppose τ < 1, one can see that for where W (:, −k), H(−k, :) refers to a sub-matrix constructed by removing the kth column and kth row from W and H, respectively.
Denote θ = 1 1 − τ H(−k, :)c k , it is easy to verify that by the definition of κ(W ) for k ∈ K.

This leads to
where the last inequality was established in Step 1. Finally, for n / ∈ K C(n, : This completes the proof.

Lemma 6
Given vector x ∈ R K . Suppose x satisfies the followings: x ≥ 0, 1 ⊤ x = 1, x ∞ ≤ a for some a ≤ 1, then Lemma 7 For n ∈ K, m ∈ K c , the following holds: where Proofs of Lemmas 5, 6, and 7 are relegated to the supplementary material in Appendices G, H, and I, respectively. Let us assume C t = C. The gradient of (18) is Denote P = X ⊤ (XC − X) and Y = ∇Φ µ (C). One can see that Our goal is to show that there always exists n ∈ K such that c n,ℓ satisfies condition (38) in Lemma 5 for any ℓ. If this holds, then the corresponding gradient value [∇ f (C)] n,ℓ is expected to be small, since the corresponding y n,ℓ is small. To this end, denote g ℓ as the ℓth column of ∇ f (C), i.e., where p ℓ and y ℓ are the ℓth columns of P and Y , respectively. Our objective then amounts to showing that j ∈ K where j = arg min p n,ℓ + λy n,ℓ .
Again, w.o.l.g., we assume that K = [K] and h n = e n for n ∈ K. We use a contradiction to show our conclusion. Suppose that for every n ∈ K, c n,ℓ is the largest element in row C(n, :)-i.e., c n,ℓ ≥ c n,ℓ ′ for all ℓ ′ = ℓ for every n ∈ K.
Since c n,ℓ is the largest element, and since C(n, :) is not a constant by our assumption, one can always find an ℓ ′ such that c n,ℓ > c n,ℓ ′ . Then, we have The third equality holds because of the assumption that supp(c ℓ ) ⊆ K. The above is a contradiction, which means that for any given ℓ, there must be at least an n ∈ K such that c n,ℓ < c n,ℓ ′ for a certain ℓ ′ . Therefore, by Lemma 5, we have the following inequality: In the meantime, for m ∈ K c , we have because C(m, :) = 0 ⊤ . For an n ∈ K that satisfies (40) and an m ∈ K c , we have g m,ℓ − g n,ℓ = p m,ℓ + λy m,ℓ − p n,ℓ − λy n,ℓ = p m,ℓ − p n,ℓ + λ 1 N − y n,ℓ .
Using Lemma 7, we can establish an lower bound of p m,ℓ − p n,ℓ , i.e., where the first equality is by (25), the first and second inequalities are by (28) and Lemma 7, respectively. Therefore, we have (41) where the last inequality can be derived from noise bound given in (22).
Hence, one can see that arg min n p n,ℓ + λy n,ℓ ∈ K As a result, the update rule of FW will make supp(c ℓ ) ⊆ K for the next iteration. This completes the proof. By the stopping criterion, before FW terminates, η + 2δ ≤ Xc ℓ − x ℓ 2 holds. Hence, the following chain of inequalities holds: .
The last inequality holds since V c ℓ 2 ≤ N i=1 c i,ℓ v i 2 ≤ max j v j 2 = δ.
Let q ℓ = W ⊤ W (Hc ℓ − h ℓ ). In addition, w.o.l.g., let q j,ℓ , q s,ℓ be the smallest and the second smallest elements in q ℓ , respectively. By the definition of d(H), we have h ⊤ n q ℓ = K k=1 h k,n q k,ℓ = h j,n q jℓ + h s,n q s,ℓ + k =j,k =s h k,n q k,ℓ ≥ h j,n q j,ℓ + h s,n q s,ℓ + q s,ℓ k =j,k =s h k,n = h j,n q j,ℓ + h s,n q s,ℓ + q s,ℓ (1 − h j,n − h s,n ) = h j,n q j,ℓ + q s,ℓ (1 − h j,n ) ≥ d(H)q j,ℓ + (1 − d(H))q s,ℓ , where the last inequality holds because By the definition of ν in (16), we have Consequently, we have e ⊤ j W ⊤ W (Hc ℓ − h ℓ ) = e ⊤ j q ℓ = q j,ℓ .
APPENDIX G PROOF OF LEMMA 5 Let L n be a set of indices of the largest elements in row C(n, :): L n := {arg max i c n,i }.
APPENDIX H PROOF OF LEMMA 6 Suppose that x 1 = max i x i , which is without loss of generality. It is easy to verify a ≥ x 1 ≥ 1/K. Therefore, Note that Since characteristic polynomial of Z is hence λ is eigenvalue of W ⊤ W leads to −λ, λ are eigenvalues of Z.  where I t n,i , J t n,j can only be either 0 or 1. Let S t n,(i,j) = c t n,i − c t n,j , and U t n,(i,j) = I t n,i − J t n,j . The updating rule in terms of S t i,j is S t+1 n,(i,j) = t t + 2 S t n,(i,j) + 2 t + 2 U t n,(i,j) .
(47) Indeed, (47) can be shown by induction. Since (47) involves the sample tuple n, i, j on both sides, and for the sake of simplicity, we omit these subscript temporarily in the following induction proof.
In addition, This completes the proof.