A Fast Scale-Invariant Algorithm for Non-negative Least Squares with Non-negative Data
Authors: Jelena Diakonikolas
Chenghui Li
Swati Padmanabhan
Chaobing Song
Venue: NeurIPS 2022 Accept (2022)
Keywords: acceleration
non-negative least squares
scale invariance

Abstract:  Nonnegative (linear) least square problems are a fundamental class of problems that is well-studied in statistical learning and for which solvers have been implemented in many of the standard programming languages used within the machine learning community. The existing off-the-shelf solvers view the non-negativity constraint in these problems as an obstacle and, compared to unconstrained least squares, perform additional effort to address it. However, in many of the typical applications, the data itself is nonnegative as well, and we show that the nonnegativity in this case makes the problem easier. In particular, while the worst-case dimension-independent oracle complexity of unconstrained least squares problems necessarily scales with one of the data matrix constants (typically the spectral norm) and these problems are solved to additive error, we show that nonnegative least squares problems with nonnegative data are solvable to  multiplicative error and with complexity that is independent of any matrix constants. The algorithm we introduce is accelerated and based on a primal-dual perspective. We further show how to provably obtain linear convergence using adaptive restart coupled with our method and demonstrate its effectiveness on large-scale data via numerical experiments. 

A Fast Scale-Invariant Algorithm for Non-negative Least Squares with Non-negative Data
Jelena Diakonikolas University of Musconsin-Madissee jelena@cs. wisc.
Chenghui Li
University of Wisconsin-Madison cli539cwisc. edu
Swati Padmanabhan University of Washington, Seattle pswati@ur.c edu
Chaobing Song
University of Wisconsin-Madison chaobing. songgwisn.....cc.ooo edu
Abstract
Nonnegative (linear) least square problems are a fundamental class of problems that is well-studied in statistical learning and for which solvers have been imple- mented in many of the standard programming languages used within the machine learning community. The existing off-the- shetlf solvers view the non-negativity constraint in these problems as an obstacle and, compared to unconstrained least squares, perform additional effort to address it. However, in many of the typical applications, the data itself is nonnegative as well, and we show that the nonneg- ativity this case makes the problem easier. In particular, while the worst-case dumension-indeeeddee oracle complexity for unconstrained least squares problems necessarily scales with one of the data matrix constants (typically the spectral norm) and these problems are solved to additive error, we show that nonnegative least squares problems with nonnegative data are solvable multiplicative error and with complexity independent of any matrix constants. The algorithm we introduce is accelerated and based on a primal-dual perspective. We further show how to provably obtain linear convergence using adaptive restart coupled with our method and demonstrate effectiveness on large-scale data via numerical experiments.
Introduction
Nonnegative least squares (NNLS) problems, defined by minx>0 ! | - blž where E mxx and b Th are fundamental problems and have been studied for decades in optimization and statis- tical learning [26, 7, 23), with various off-the-shelf solvers available standard packages of Python (as optimize.... nnls in the SeiPy package), Julia (as anls.jl) and MATLAB (as lsqnonneg). Within machine learning. NNLS problems arise whenever having negative labels not meaningful, for example, when representing prices, age, pixel intensities, chemical concentrations, frequency counts. NNLS is also widely used as a subroutine in nonnegative matrix factorization [10, 18, 25] to extract sparse features in applications like clustering, collaborative filtering, and community detection.
From a statistical perspective, NNLS problems can be shown to possess a regularization property that enforces sparsity similar to LASSO [46], while being comparatively simpler, without the need to tune regularization parameter or perform cross-validation [44, 8, 17, 24, 49, 43].
From algorithmic standpoint, the nonnegativity constraint in NNLS problems is typically viewed as an obstacle: most NNL.S algorithms perform additional work to handle it, and the problem
Authors ordered alphabetically.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

considered harder than unconstrained least squares. However, in many important applications of NNLS, such as text mining [6], functional MRI [3, 20]. EEG data analysis [33], pulse oximetry [22. 50], statistical procedures in observational astronomy [19], and those traditionally addressed using nonnegative matrix factorization [10], the data is also nonnegative. We argue in this paper that when the data for NNLS is nonnegative, it is in fact possible to obtain stronger guarantees than for traditional least squares.
Our Contributions. We study NNL.S problems with (element- nonnegative data matrix A, to which we refer as the NNLS+ problems, through the lens of the (equivalent) quadratic problems:
min { f(x) del' 1 A x20
- c x
where c A Tb may be assumed element-wise positive. This assumption is without loss of generality since if Cg < ( for some j. then V, (x) 2 0, implying that the jth coordinate of the optimal solution is zero?. Hence, we could fix 2 = 0 and optimize over only the remaining coordinates.
We further assume that the matrix A is non-degenerate: none of its rows columns has all elements equal to zero. This assumption is without loss of generality because (1) if such a row exists, we could remove it without affecting the objective, and (2) if the jåh column had all elements equal to zero, the optimal value of (P) would be -oo, obtained for x with x + oc. Having established our assumptions and setup, we now proceed to state our contributions, which are three-fold.
(1) A scale-invariant, E-multiplicative algorithm. We design an algorithm based on coordinate descent that, in total cost O( (A) ), constructs an E-multiplicative approximate solution to (P). Our
algorithm capitalizes on structural properties of (P) that arise as a result of the nonnegativity of A.
Theorem 1.1 (Informal; see Theorem 3.5). Given a matrix A E R and > O, define f(x) = || - x and E argminx20 fxx). Then, there exists an algorithm that in K = OO log n- ve) iterations and O ((nz(A))))))))))))))))))())))( log n + (e)) arithmetic operations returns XK E 12T such that ( - x")] S =l (xx)).
The application of our structural observations on (P) to Theorem 4.6 of [13] enables the recovery of our guarantee on the optimality gap; however, we provide guarantee the primal-dual gap, and this is stronger than the one on the optimality gap stated Theorem 1.1. What is significant about Theorem 1.1 is the invariance of the computational complexity to the scale of A-it does not depend on any matrix constants. This cost stands in stark contrast to that of traditional least squares, where the dependence of (oracle) complexity on matrix constants (specifically, the spectral norm of A in the Euclidean case) is unavoidable [35, 36], and multiplicative approximation is not possible in general. 3 In general, ocale-invariance is crucial feature in problems with data matrices, since dependence on the width implies that the algorithm is technically not polynomial-time time. This feature has, in fact, been an object of extensive study in the long line of works on packing and covering linear programs [48, 1] and its variants such as a fair packing [14]. Conceptually, our algorithm is a new acceleration technique inspired by VRPDA [45].
(2) Linear convergence with restart. By incorporating adaptive restart (P), we improve the guarantee of Theorem 1.1 to one with linear convergence (with log(1/s) complexity). Thus, we establish the first theoretical guarantee for NNLS$ that simultaneously satisfies the properties of being ecale-invarian, accelerated, and tinearly-convergent.
Theorem 1.2 (Informal; see Theorem 4.1). Consider the setup of Theorem 1.1. Then, there is an algorithm that in expected ((nnz(A)(log) + vi log(()) arithmetic operations returns 3g E RT with f(xx) - f(x*) & =lf(x*)), ), where μ is the constant in a local error bound for (P).
Proving this bound requires bounding the expected number of iterations between restarts in conjunc- tion with careful technical work in identifying an appropriate local error bound for NNLS+.
TTo see this, note that by first-order optimality condition vf(x*) x*) 2 0 for all x 2 0. Choosing with x = x" for f j and X, = in the first-order optimality condition gives x] = 0. "To see why, consider a case in which the optimal objective value equals zero. Then any problem with a multiplicative guarantee of the form in Theorem 1.1 would necessarily return an optimal solution.

---

(3) Numerical experiments. We consolidate our theoretical contributions with demonstration of the empirical advantage of our restarted method over tfate-of-the-art solvers via numerical experiments on datasets from LibSVM with sizes up to 19996 x 1355191. Figure I shows that, when combined with the restart strategy, our algorithm significantly outperforms the compared algorithms.
Related work. NNLS has seen a large body of work on the empirical front. The first method that was widely adopted in practice (including in the lsquonneg implementation of MATLAB) due to the seminal work of [26] (originally published in 1974). This method, based on active sets, solves NNLS via sequence of (unconstrained) least squares problems and has been followed up by [7. 47, 34, 11] with improved empirical performance. While these variants are generally effective on small mid-scale problem instances, they are not suitable for extreme-scale problems ubiquitous in machine learning. For example, in the experiments reported in [34], Fast NNLS [7] took 6.63 days to solve a problem of size 45000 x 45000, while the TNT-NN algorithm [34] took 2.45 hours. However, the latter requires computing the Cholesky decomposition of A A at initialization, which can be prohibitively expensive both in computation and in memory. Another prominent work on the empirical front is that of [23], which performs projected gradient descent with modified Barzilai- Borwein steps [4] and step sizes carefully designed sequence of diminishing scalars.
Another, separate, line of work concerns optimization algorithms with multiplicative error guarantees. Of interest to us are standard first-order algorithms that run in time that is near-linear in the input size and are thus applicable to large-scale setting (for multiplicative-eroo algorithms applicable to broad classes of problems but that run in time that is superlinear in the input size, see [41, Chapter 7)). Most of the existing literature in this domain concerns positive (packing and covering) linear programs (e.g., [29, 51, 1, 30). Results also exist for positive semidefinite programs [2], (nonlinear) fair packing and covering problems [32, 14], and fair packing problems under Schatten norms for matrices [21]. With the exception of [14] (discussed below), none of these results are directly applicable to (P).
To the best of our knowledge, theoretical guarantees explicitly for (P) have been scarce. For in- stance, [23] mentioned in the preceding paragraph provides only asymptotic guarantees. Orthogonally, the result on 1-fair covering by [14] solves the dual of NNLS+ which also gives a multiplicative
guarantee for NNLS+ but with overall complexity Ō mnx(A) ), where (-) hides poly-log factors.
Since our algorithm is based on the coordinate descent algorithm, we highlight some results of other coordinate descent algorithms when specialized to the closely related problem of uncon- strained linear regression. The pioneering work of [37] proposed a coordinate descent method
called RCDM, which in our setting has iteration cost 0 Ei= ||A., -||xx=x'''"" where ||A [| is the Euclidean norm of the jth column of A. This was improved by [27]. in an algorithm termed ACDM, by combining Nesterov's estimation technique [39] and coordinate sampling, giving an it- (V" " [1,, ) for solving (P). The latest results in this line of 40] perform non-uniform sampling atop a framework of [37] and achieve iteration complexity of O / /اك [lA,, -EIxx with [12] dropping the dependence on maxisisn | 22 Additionally, the work of [28] develops an accelerated randomized proximal coordinate gradient (APCG) method to minimize composite convex functions.
eration complexity of O(* work by [2, 42,
As remarked earlier, [12], coupled with insights on NNLS+ problems provided in this work, can recover our guarantee for the optimality gap from Theorem 3.5. However, our work the first to bring to the fore the properties of NNLS+ required to get such guarantee, and our choice of primal-dual perspective allows for a stronger guarantee terms of an upper bound on the primal-dual gap. Further, our algorithm is novel type of acceleration, with our primal-dual perspective transparently illustrating our use of the aforrmmnnionee propertics. We believe that these technical contributions, along with our techniques to obtain vastly improved theoretical guarantees with the restart strategy applied to this problem, are valuable to the broader optimization and machine learning communities.
2 Notation and preliminaries
Throughout the paper, we use bold lowercase letters to denote vectors and bold uppercase letters for matrices. For vectors and matrices, the operator' 2 is applied element-wise, and IR., is the non-negative part of the real line. We use (a, b) to denote the inner product of vectors and b and V
3

---

for gradient. Given matrix A. we use A.j for its jth column vector, and for vector 2, denotes its jth coordinate. We use nnz(AA for the number of non-zero entries of A. We use X5, for the vector in the hth iteration and, to disambiguate indexing, use [xx]; to mean the jth coordinate of X The 7th standard basis vector is denoted by e;. For an n-dimensional vector x and A E RRxx,, we define diag(()) A.1 ..... ||.... (121) and | Σ A Finally, [n] def (1.2.....
differentiable function / R" + R is convex if for any x, & E R"; we have f(x) 2 f(x) + (V[(x) & x). A differentiable function f RR + R is said to be jestrronlly convex w.r.t the f2-norm if for any x, & E IR", we have that /(x) 2 (x) + (VI(x), X - x) + 2 x 8 We have analogous definitions for concave and strongly concave functions, which flip the inequalities noted.
Given a convex program milxxe fx), where f R" - IR differentiable and convex and X S IR" closed and convex, the first-order optimality condition of solution x* E argmin f(x) is
(Vx E x): (V/(x"), xx) 2 0.
Problem setup. As discussed in the introduction, our goal is to solve (P), with E IR Rmxx" For notational convenience, work with the problem in the following scaled form:
[(x)
der
Ax 1 1'x)
min x
This assumption w.l.o.g, since any (P) can be brought to this form by a simple change of variable i = Cjzry (see also e.g.. [1, 14] for similar scaling ideas). The scaling need not be explicit in the algorithm since the change of variable 2 = Cj is easily reversible.
Properties of the objective. To kick off our analysis, we highlight some properties inherent to the objective defined in (2). These properies, which strongly need the non-negativity of A and x, are central to obtain a acale-invariant algorithm for (P).
Proposition 2.1. Given f TR -> TR as defined in (2) and x* E Cargmin f(x), the following statements all hold.
a) Vf(x)) 0.
b) f(x')
Ax 11--1=
c) for all j E [n], we have x E [0, T
d) - 2 ΣjEnn |A I f(x") 2inngela |A
The validity of division by ||A واار in the preceding proposition is by the non-degeneracy of A discussed in the introduction. We prove this proposition in Section 6.1.
An important consequence of Proposition 2.1 (C) is that (P) can be restricted to the hyperrectangle
X' {x E R 0 < æÿ the search space. Thus, going forward, we replace the constraint x 2 0 in (P) by x E X.
||A, without affecting its optimal solution, but effectively reducing
1
Primal-dual gap perspective. As alluded earlier, our algorithm is analyzed through a primal-dual perspective. For this reason, is useful to consider the Lagrangian
L(x. y) (Ax,y) y) | 1
from which we can derive our rescaled problem (2) as the primal problem mingex P(x),
where 1 +
P(x)=
max L(x,y) - 1 max y2o
1 +(x -
x -x
y
Thus, the Lagrangian is constructed in a way that one can derive the primal problem from it while also localizing the matrix in the bilinear term and ensuring coordinate-wise separability of the terms involving either only the dual or only the primal variable since this greatly simplifies our analysis. Similar to [45], we use Eq. (3) to define the following relaxation of the primal-dual gap, for arbitrary but fixed X' and v mm
Gaptu. (u.v) (x, L(x,v) L(u,y)
(4)

---

The significance of this relaxed gap function is that for a candidate solution x and an arbitrary (x y Ilm,, a bound on 2 (x, y) translates to one the primal error, as follows. First select Axx) = f(x"), ), we have u x AR. Then, by observing that [ A = f(x) and ,
((x) [(x) L(x, AX) L(x Ax*)
For a fixed x, [(x, .) is I-strongly concave and minimized at Ax. Thus, L(x", Ax) & L(x", Ax*)
A x*)11 Hence, we have the following primal bound:
((x) f(x'
A(X x*) ((X L(x".
Gap (x' .Ax)
y).
(5)
In light of this connection, our algorithm for bounding the primal error is one that generates iterates that can be used to construct bounds Gap (u.v) (x, ý), we detail next.
3 Our algorithm and convergence analysis
Our algorithm, Scale Invariant NNLS+ (SI-NNLS+), is an iterative algorithm using the estimate sequences (x) and (y)) for k 2 1 (see Section 3.1) giving the primal and dual updates
X = argmin (x(x)
and Yk = argmax 3/ (y). YER"
(6)
xEA
(u.v) (u.
We use our algorithm's 's iterates from Eq. (6) to construct Gg. an upper estimate of Gapc where Xg, Yk are convex combinations of the iterates and u E X, v RT. where is fixed and is arbitrary. Our motivation for constructing Gx(u, v) 2 Gap (u.v) (Kx, F&) (x) is obtain a bound on the primal error, using Inequality (5). The main goal in the analysis
show that G5 Ar where As is sequence of positive numbers and Q is bounded. To obtain the claimed multiplicative approximation, we use Inequality (5) and argue that for u = x", V AXt, we have Q S (lf(xx))) ))).
Q,
Our algorithm may be interpreted as a variant of the VRPDA algorithm [45], with our analysis inspired by the approximate duality gap technique [13]; however, unlike VRPDA? which uses estimate sequences, our analysis directly bounds the primal-dual gap. Another difference from VRPDA is in our choice of primal regularizer (described shortly) and our lack of dual regularizer. variant of SI-NNLS+ suitable for analysis shown in Algorithm 1, with proofs in Section 6.2 and Section 6.3. We give equivalent implementation version ("Lazy SI-NNLS+", which updates few coordinates as is possible to improve cost per iteration) in Algorithm 2 and its analysis in Section 7.
Algorithm 1 Scale Invariant Non Regatvee Least Squares with Non- Non-negative Data (SI-NNLS+)
1: Input: Matrix A E mm with 2 4, accuracy initial point X
2: Initialize: Xo Xo. Yo Yo Axo. K n n +
Ve"
Gn
01
a2 T Ao
A1 01 (x) ! |x x
3: for k= 1, 1.2..... K do
4: Sample jk uniformly at random from 11.2....,12 ,n]
5:
X t arg minxex (x))), for (x) defined by Eq. (12) and Eq. (15)
6:
Yı - arg maXyeRm (y),,, for (y) defined by Eq. (9)
7:
x; I [Ax-1x + ax (nxx- - (n 1)xx-1)]
Te
8:
FA + Ya + au (yk - YYk-1)
a2+1
9:
A As + actl, 0k42 minf nap41
Ar+1
10: end for
11:
Return X
3.1
Gap estimate construction
The gap estimate Gy is constructed as the difference G& (u, v) Uv(v) Lا (u) Un(v) > [(x,, v) and La (u) < L(u, yx) are, respectively, upper and lower bounds we construct on the Lagrangian. It then follows by Eq. (4) that G& (u, v) is valid upper estimate of Gap
(u.v) (x,,,))
5

---

We first introduce technical component our constructions L; and U3 crucially hinge on: we define two positive sequences of numbers {a:h21 and {a sigs. with one of their properties being that both sum up to A >> 0 for k 2 1. Specifically, we define Ao 0 and {a:):21 as A; A;-1- The sequence {a²² } changes with & and for k = 1 is defined by af a1, while for l 2 2:
a1 (n 1) 1)
if 1, if - i i 1, i k.
a. = na, - (n 1)
(7)
nak:
Summing over i E [k] verifies that As- _i=1 at. For the sequence {a Jug to be non-negative, we further require that a1 - (n - 1)a2 2 0 and Vi 2 2, na, (n 1)0,+1 2 O.
The significance of these two sequences lies in defining the algorithm's 's primal-dual output pair by
&
x and
yx = Ak
Σ a.y. KE(M)
(8)
AR
RE(M)
The intricate interdependence of {a;} and {a enables expressing X in terms of only {a:}. This expression further simplifies to a cheaper recursive one, which is used Algorithm 1.
With the sequences {ai)e2 and {a hesck in tow, we are now ready to show the construction of an upper bound (v) [(x,, v) and lower bound (u) on C(u, yk)-
Upper bound.
To construct an upper bound, first observe that by Eq. (3) and Eq. (8),
[(x;, v) = (AXg, v)
2 vil 1 xx
(Ax, v) _ | 1 x;
Ak iE(A]
Consider the primal estimate sequence defined for k 0 vh and for k 2 1 by
(v) def
Σ at (Axi, v) [E[(]
1 2 MIl2 - 1 1
which ensures that L(x,, v) tt vux(v)) ). A key upshot of constructing (v) as in Eq. (9) that the quadratic term implies Ag-strong concavity of th for k 2 1, which in turn ensures that the vector yı = arg m 44 (y) from Eq. (6) is unique. This property, coupled with the first-order optimality condition in Inequality (1), gives that for any y E R vu(y)) I vy(Y))) Ax lly yu112. We are now ready to define the following upper bound by:
U&(v) def
A (v()
1 , llv y
(10)
The preceding discussion immediately implies that U2 is a valid upper bound for the Lagrangian. Lemma 3.1. For Us as defined in Eq. (10), Lagrangian defined in Eq. (3) and & E IR. in Eq. (8). we have, for all y E Rm the upper bound Ux(y) 2 L(x,,y).
Lower bound. Analogous to the preceding section, we now obtain lower bound on the Lagrangian, completing the bound on the gap estimate. However, the construction becomes more technical. We start with the same approach as for the upper bound. Since L(u, V&) is concave in Yk. by Jensen's s
inequality: L(u, FA) Z A Eielo a yA) - 1 Tu u - 1 |1: Were we to define the dual estimate sequence Φ& in the same way as we did for the primal estimate sequence thk, we would now simply define it as Aş times the right-hand side in the last inequality. However, doing so would make Φe depend on Yı. which is updated after Xt, which in Eq. (6) defined as the minimizer of the oz.
To avoid such a circular dependency, we add and subtract a linear term Σie (A (ATy Ty--1, u), where V.-1. defined later, are extrapolation points that depend only on y .i-1. We thus have
((u Z
1 Σ as (Au,
T
1ly
Σ a (A - Y1-1)
AR
2
A

---

If we now defined OR based on the first term in the above inequality, we run into another obstacle: the linearity of the resulting estimate sequence is insufficient for cancelling all the error terms in the analysis. Hence, as is common, we introduce strong convexity by adding and subtracting an appropriate strongly convex function. Our chosen strongly convex function is motivated by the box-constrained property of the optimum from Proposition 2.1 (c) and crucial in bounding the initial gap estimate. coincides with Φo: for any x E IR,, define the function
(x)
2 -xoli
(11)
This function is 1-strongly convex with respect to || a| and used in defining Φ1
(u) a1 (A yo 1, u) (u).
(12)
The definition of Φi is driven by the purpose of cancelling initial error terms. Next, we choose Φx so that for any fixed X, have
(l)) = E Σ a a 1. + do(m){
(13)
(E
where the expectation is with respect all the randomness in the algorithm. This construction is used to reduce the per-iteration complexity, for which we employ randomized coordinate update on X3 for k 2 To support such updates, we relax the lower bound to hold only in expectation.
Concretely, let ji be the coordinate sampled uniformly at random from [n] in the 5th iteration of SI-NNLS+, independent of history. Fix y; for = ....., k 1 and for k > and X, define
(x(x) Φ1 (x)
Σ na: (A
V.-1 (
1,21,
(14)
For A: 2 2, (u) can also be defined recursively via
Φx (x) ( + nan (A y 1, 1, 0ja)
(15)
The function Φ% inherits the strong convexity of Φv. This property, together with Eq. (6) and first-order optimality from Inequality (1), give
(x(x) 2 (xx(x)) + 2|1x xx
(16)
Along with strong convexity, choice of de in Eq. (15) leads to the following properties essential to our analysis: (1) Φu is separable its coordinates; (2) the primal variable x₁ is updated only at its
1 coordinate; (3) Eq. (13) is true. These formally stated in Proposition 6.2.
With the dual estimate sequence Φb defined Eq. (15), we now define the sequence L2 by
(x) (x(xx)) + 1 ||x xx A
oo(x) Eieny dr((Ax, y 1 1x:1
(17)
A
We conclude this section by justifying our choice of T.g as a valid expected lower bound on EL(( ÿx). Lemma 3.2. For T. defined in Eq. (17), for the Lagrangian in Eq. (3) and Y& in Eq. (8), have, for a fixed u E X, the lower bound (L,(u, 2 EEI, where the expectation is with respect to all the random choices of coordinates in Algorithm 1.
3.2 Bounding the gap estimate
With the gap estimate Gg constructed in the preceding section and combining Eq. (10) and Eq. (17), now achieve goal of bounding A G& (to obtain convergence rate of the order 1/Ap) by bounding the change in G2 and the initial scaled gap A
Lemma 3.3. Consider the iterates {xx} and {yk} evolving according to Algorithm 1. Let and assume that a1
and 01 2 (n 1)a2. while for k 2 3,
1.5...5
ag < min naa= V A
(18)
2n
7

---

Then, for fixed u E X. any R" and all lo 2 2, the gap estimate Gj = UR LL satisfies
B(AgGr(Xx y) 1 - yx)
S -E Ak lly vxll2 Ak-1
_ 1 BB|x 1 xx11 + _ Xk- 111
lyy - yx-118 My----118
2
(AE(A - xx). -Yk-1) ax |E(A( Xk-1), yk-2)
Al-y -vk-1l2 + Ak-2BBlly- yk-
Lemma 3.4. Given a fixed u E X, any v E Rm yo yo, and X1 and y1 from Algorithm 1, we have
A
A, v) ( yo), u} + Φo ()) do (x 2 u- x 11
|| vill
Combining the two lemmas, we now bound GK and deduce our final result on the primal error.
Theorem 3.5. [Main Result) Assume that 2 4. Given matrix A R >>> O, an arbitrary Xo E X and Yo = Yo = Axo, let Xx and Aş evolve according to SI-NNLS+ (Algorithm for k 2 1. For defined in (2), define x* E argmin [xx. Then, for all K 2 have
E (V & x*) 1 1|A((Kx - [[XXK-xx)11] I
200(x*)
X -0x11 13
Ax
Ak
When K 2 Inlog n, we have Ax 2 (K- (ng3nlogm) If ff(x) S 1f(xx)),, then for K 2 in logn
36n
Ve we have [xf(x) f(x') S elf(x*)]
).
The total cost O (nnz(A) (log + ++)).
The assumption do (x") < [[(x*) above is satisfied by Xo = 0 (c.f. Proposition 2.1 and Eq. (11)).
We reiterate that the reason All does not show up in the final bounds (thereby rendering our algorithm "scale-invariant") is because Proposition 2.1 allows bounding ||xn x 1 by [f(x*)], where we crucially used the non-negativity of A and x; this does not seem possible for general A. Remark 3.6. SI-NNLS+ (Algorithm 1) and Theorem 3.5 also generalize Increasing the batch size grows our bounds and number of data passes by a factor of at most square- root of the batch size 8, by relating the spectral norm the a columns of A corresponding to a batch to the Euclidean norms of individual columns of A from the same batch. However, due to efficient available implementations of vector operations, mini-batch variants our algorithm with small batch sizes can have lower total runtimes on some datasets (see Section 5).
mini-batch version.
4 Adaptive restart
We now describe how SI-NNLS+ can be combined with adaptive restart to obtain linear convergence rate. To apply the restart strategy, need suitable upper and lower bounds the measure of convergence rate. Our measure of optimality is the natural residual r(x) = |[[(x) la [31] for
R(x) = - UR (x - - f()) = x -VV
(19)
where Tgg is the projection operator onto R and A is as defined in Section 2. For A = I, R(x) is the natural map as defined in, e.g., [15]. Due to space constraints, we only state the main result of this section in the following theorem, while full technical details are deferred Section 6.4.
Theorem 4.1. Given an error parameter >> O and x O, consider the following algorithm A:
A: SI-NNSS+ with Restarts Initialize: le = 1. Initialize Lazy SI-NNLS+ at X:-1- Run Lazy SI-NNLS+ until the output X satisfies r(*%) S (x(xx-1). Restart Lazy SI-NNLS+ initializing at X Xk.
Increment
Repeat until r(x) E.

---

1G4 SAA TOSER SETMA
11. TOSEK TERM
1250
190
200
2000
2200
(a) real- -iim
(b) real-sim
(c) real-sim
(d) real-sim
revid 1- 19966 1250380)
EIUAGE te 10906 IY 1205380)
reNGS Io 3896E 13250L
revid 1996E n 1353340
(e) news20
ne
(g) news20
(h) news20
17009-UH TEOST, 150MME
POO 1669T, 1901MM
FIEORLAY 16097, 1909ME
#0 18507
(i) E2006train
i E2006train
(k) E2006train
E2006train
Figure NNLSS algorithms with restart on real sim, news20 and E2006train with spectral
norm 3 10-3 10-3,5 5. 10-6 and condition number 6 10 5 1020, 6 10 respectively.
Then, the expected number of arithmetic operations of A is (nnzz(A) + Ln log (~(X²) As a consequence, given E >> 0, the total expected number of arithmetic operations until a point with [(x) - [(x*) < {x((())) can be constructed by A is O (nnz(A) A) llogn + yn p
log n
p2
5 Numerical experiments and discussion
We conclude our paper by presenting the numerical performance of SI-NNLS+ and its restart versions (see the efficient implementation version in Algorithm 2) against FISTA with restart [5, 38], a general- purpose large-scale optimization algorithm, OA+DS with restart designed by [23] specifically for large-scale NNLS problems, and lazy implemented APCG [16] with restart. We use the same restart strategy for all the algorithms, proposed in Section 4. As an accelerated algorithm, FISTA has the optimal 1/k² convergence rate; OA+DS, while often efficient in practice, has only an asymptotic convergence guarantee. For FISTA, we compute the tightest Lipschitz constant (i.c., the spectral norm ||A]); for OA+DS, we follow the best practices laid out by [23]. For our SI-NNLS+ algorithm and its restart version with batch size bs Section 4.* For the restart version with batch size larger than 1, we choose the best batch size in {10, 50, 300, 500} and compute the block coordinate Lipschitz constants as the spectral norms of the corresponding block matrices. All algorithms were implemented in Julia and run on a server with 32 Intel(R) Xeon(R) Silver 41 10 32-Core Processors.
follow Algorithm 2 and the restart strategy in
We evaluated the performance of the algorithms on the large-scale sparse datasets real-sin news20, and E2006train from the LibSVM library [9]. Both real-sim and news20 datasets have non- negative data matrices, but the labels may be negative. When there exist negative labels, is possible for the elements of A b to be negative. In such cases, per the discussion from the introduction, we can simply remove the corresponding columns of A and solve an equivalent problem with smaller dimension. On the other hand, the data matrix in E2006train dataset is not nen-neeative, which means that this dataset does not satisfy the assumption required for the analysis of Algorithm 1.
'The algorithm is implemented for the non-scaled version of the problem. (P); see Section 2. An implemen- tation is in https///gihub..oo/aatuurr 1/nnlr-2222222

---

HORAA 1240 1. 201000
CORMN
17SK TERMA
1220 (a) real- - sim
25
150
2011
(b) real- -sim
(c) real-sin
(d) real-sim
reviGE 1-- 10906 1290880)
ERUNGE to 10906 n 1205380)
FAVIGE (1) 1090 1253180)
reviSS 1006 1353340
(e) news20
ne
(g) news20
(h) news20
#7009/UP TEOST, 15014R
POR 16697, 1901MM
F7E0 16597, 1909M
#0 18307
(i) E2006train
i
E2006train
(k) E2006train
() E2006train
Figure 2: NNI.S+ algorithms without restart on real-sim, news20 and E2006train with spectral
norm 3 10-3 10-$,5 10-6 and condition number 6 1018 5 1020 6 10 respectively.
However, Algorithm can still be run by keeping only the non-negativity constraints for primal updates. This example is provided solely for illustration of empirical performance.
Results. To compare all implemented algorithms, plot the natural residual/objective value gap versus number of data passes/time in Figure for all algorithms implemented with restart and in Figure 2 for all algorithms implemented without restart. As can be observed from the two figures, our proposed restart speeds up all the algorithms. Figure ((a)-(d) shows that SI-NNLS+ is better than FISTA and OA+DS in terms of number of data passes on the real-sim dataset and better than APCG in terms of time. While the proposed restart strategy speeds up all the algorithms to linear convergence, variants of SI-NNLS+ remain competitive in all the settings. In terms of the performance of different variants of SI-NNLS+, with bs = 1, we have a much better coordinate Lipschitz constant than for data passes. As FISTA and OA+DS take less time accessing the full dataset once, they have lower runtimes than SI-NNLS+ but are beaten by SI-NNLS+ with restart and bs 1. In Figure ((e)-(h), on the news20 dataset, in terms of number of data passes, restarted APCG and SI-NNLS+ with bs = 1 are dominant. However, as news20 a very sparse dataset, letting bs = 1 significantly increases the total time to access the full data once due to the overhead per iteration. As result, single-coordinate methods have the worst runtimes, while restarted OA+DS the fastest but SI-NNLS+ with = 10 remains competitive. Figure (())-l) shows the performance comparison on the E2006train dataset. On this dataset, both restarted FISTA and restarted OA+DS ran for hours without visibly reducing the function value. SI-NNLS+ outperforms FISTA and OA+DS in both number of data passes and time, and outperforms APCG in terms of time. Further experiments are left for future work.
10 and thus the case of bs
dominates bs 10 in terms of
Acknowledgements
CS was supported in part by the NSF grant 2023239. JD was supported in part by the NSF grant 2007757, by the Office of Naval Research under contract number N00014--2-1-2345 and by the Wisconsin Alumni Research Foundation. Part of this work was done while JD was visiting Simons Institute for the Theory of Computing. Part of this work was done when SP was visiting JD at the University of Misconsin-Madison
10

---

References
[1] Zeyuan Allen-Zhu and Lorenzo Orecchia. Nearly linear-time packing and covering Ip solvers. Mathematical Programming, 175(1):307-353, 2019.
[2] Zeyuan Allen-Zhu, Zheng Qu, Peter Richtárik, and Yang Yuan. Even faster accelerated coordinate descent using non-uniform sampling. In International Conference Machine Learning, pages 1 1110-1119. PMLR, 2016.
[3] Anders H Andersen and William S Rayens. Structure-secking multilinear methods for the analysis of fmri data. Neurolmage, 22(2):728-7777 2004.
[4] Jonathan Barzilai and Jonathan M Borwein. Two-point step size gradient methods. IMA journal of numerical analysis, 8(1):141-148, 1988.
[5] Amir Beck and Marc Teboulle. A fast iterative shrinkage-threeholdiin algorithm for linear inverse problems. SIAM journal on imaging sciences, 2(1):183-2022 2009.
[6] Michael W Berry, Murray Browne, Amy N Langville, V Paul Pauca, and Robert J Plemmons. Algorithms and applications for approximate nonnegative matrix factorization. Computational statistics & data analysis, 52(1):155-173, 2007.
[7] Rasmus Bro and Sijmen Jong. fast oon-negaaivityy-osstaaneee least squares algorithm. Journal of Chemometrics, 11:393-400,, 09 1997.
[8] Alfred M Bruckstein, Michael Flad, and Michael Zibulevsky. On the uniqueness of nonnegative sparse solutions to underdetermined systems of equations. IEEE Transactions Information Theory, (4((1))48133 4820, 2008.
[9] Chih-Chung Chang and Chih-Jen Lin. Libsvm: library for support vector machines. ACM transactions on intelligent systems and technology (TIST), 2(3):1277, 2011.
[10] Andrzej Cichocki, Rafal Zdunck, Anh Huy Phan, and Shun-ichi Amari. Nonnegative matrix and tensor factorizaaions: applications to exploratory multi-way data analysis and blind source separation. John Wiley & Sons, 2009.
[11] Monica Dessole, Marco Dell' Orto, and Fabio Marcuzzi. The lawson-hanson algorithm with de- viation maximization: Finite convergence and sparse recovery. arXiv preprint 1:Xiv:2108.05345,
2021.
[12] Jelena Diakonikolas and Lorenzo Orecchia. Alternating randomized block coordinate descent. In International Conference on Machine Learning, pages 1224-1232. PMLR, 2018.
[13] Jelena Diakonikolas and Lorenzo Orecchia. The approximate duality gap technique: A unified theory of first-order methods. SIAM Journal on Optimization, 29(1):660-689 2019.
[14] Jelena Diakonikolas, Maryam Fazel, and Lorenzo Orecchia. Fair packing and covering on a relative scale. SIAM Journal on Optimization, 30(4):3284-3314 2020.
[15] Francisco Facchinei and Jong-Shi Pang. Fimite-dimensional variational inequalities and com- plementarity problems. Springer Science & Business Media, 2007.
[16] Olivier Fercoq and Peter Richtárik. Accelerated, parallel, and proximal coordinate descent. SIAM Journal Optimization, 25(4):1997-2023 2015.
[17] Simon Foucart and David Koslicki. Sparse recovery by means of nonnegative least squares. IEEE Signal Processing Letters, 204 2014.
[18] Nicolas Gillis. The why and how of nonnegative matrix factorization. Connections, 12:2-2, 2014.
[19] Takashi Isobe, Eric D Feigelson, Michael G Akritas, and Gutti Jogesh Babu, Linear regression in astronomy. The astrophysical journal, 364:104-113, 1990.

---

[20] Ethan C Jackson, James Alexander Hughes, and Mark Daley. On the generalizability of linear and non-linear region of interest-based based multivariate regression models for finri data. In 2018 IFEE Conference on Computational Intelligence in Bioinformatics and Computational Biology (CIBCB), pages 1-8. IEEE, 2018.
[21] Arun Jambulapati, Jerry Li, and Kevin Tian, Robust sub- sub-gaussian principal component analysis and width- indtt-inndeeendent schatten packing. Advances in Neural Information Processing Systems, 33:15689-15701, 2020.
[22] Michael S Jennis and Joyce L Peabody. Pulse oximetry: an alternative method for the assessment of oxygenation in newborn infants. Pediatrics, 79(4):524-55) -524, 1987.
[23] Dongmin Kim, Suvrit Sra, and Inderjit Dhillon. A non-monotonic method for large-scale non-negative least squares. Optimization Methods and Software, 286)):012-10399 2013.
[24] Richard Kueng and Peter Jung. Robust nonnegative sparse recovery and the nullspace property of 0/1 measurements. IEEE Transactions on Information Theory, 64(2):689-70) 2018. doi: 10.1109/TT122117274466
[25] Abhishek Kumar, Vikas Sindhwani, and Prabhanjan Kambadur. Fast conical hull algorithms for near-separable non-negative matrix factorization. In International Conference on Machine Learning, pages 231-239. PMLR, 2013.
[26] Charles I.. Lawson and Richard J. Hanson. Solving least squares problems, volume 15 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1995. ISBN 0-89871-356-0, Revised reprint of the 1974 original.
[27] Yin Tat Lee and Aaron Sidford. Efficient accelerated coordinate descent methods and faster algorithms for solving linear systems. In 2013 ieee 54th annual symposium on foundations of computer science, pages 147-156. IEFF, 2013.
[28] Qihang Lin, Zhaosong Lu, and Lin Xiao. An accelerated proximal coordinate gradient method. Advances in Neural Information Processing Systems, 27:3059-3067 2014.
[29] Michael Luby and Noam Nisan. A parallel approximation algorithm for positive linear pro- gramming. In Proc. ACM STOC'93, 1993.
[30] Michael W Mahoney, Satish Rao, Di Wang, and Peng Zhang. Approximating the solution to mixed packing and covering LPs in parallel Ö(e-3) 9) time. In Proc. ICALP' 16, 2016.
[31] Olvi L. Mangasarian and J Ren. New improved error bounds for the linear complementarity problem. Mathematical Programming, (6(1):241-253 1994.
[32] Jelena Marasevic, Clifford Stein, and Gil Zussman. A fast distributed stateless algorithm for alpha-fair packing problems. In Proc. ICALP' 16, 2016.
[33] Eduardo Martinez-Montes, José M Banchez-Bornot and Pedro A Valdés-Sosa. Penalized parafac analysis of spontaneous eeg recordings. Statistica Sinica, pages 1464, 2008
[34] Joe M Myre, Erich Frahm, David Lilja, and Martin O Saar. TNT-NN: A fast active set method for solving large non-negative least squares problems. Procedia Computer Science, 108:755-764, 2017.
[35] Arkadi S Nemirovsky, Information- complexity of linear operator equations. Journal of Complexity, 8(2):153-175, 1992.
[36] Arkadii Semenovich Nemirovsky and David Borisovich Yudin. Problem complexity and method efficiency in optimization. 1983.
[37] Yu Nesterov. Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization, 22(2):341-362, 2012.
[38] Yu Nesterov. Gradient methods for minimizing composite functions. Mathematical program- ming, (40((1)122511))) 2013.
12

---

[39] Yurii Nesterov. A method of solving convex programming problem with convergence rate O(1/k²). In Doklady AN SSSR (translated as Soviet Mathematics Doklady), volume 269, pages 543-547, 1983.
[40] Yurii Nesterov and Sebastian Stich. Efficiency of the accelerated coordinate descent method on structured optimization problems. SIAM Journal on Optimization, (111:110-123, 2017.
[41] Yurii Nesterov et al. Lectures on convex optimization, volume 137. Springer, 2018.
[42] Zheng Qu and Peter Richtárik. Coordinate descent with arbitrary sampling i: Algorithms and complexity. Optimization Methods and Software, 31(5):829-857, 2016.
[43] Yonatan Shadmi, Peter Jung, and Giuseppe Caire. Sparse non-negative recovery from biased subgaussian measurements using nnls. arXiv preprint arXiv:/901.05727, 2019.
[44] Martin Slawski and Matthias Hein. Non-negative least squares for high-dimensional linear models: consistency and sparse recovery without regularization, 2014.
[45] Chaobing Song, Stephen J Wright, and Jelena Diakonikolas. Variance reduction via primal-dual accelerated dual averaging for nonsmooth convex finite-sums. In International Conference Machine Learning 2021.
[46] Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267-285 1996.
[47] Mark H Van Benthem and Michael R Keenan. Fast algorithm for the solution of large-scale non-negativittyconstrra least squares problems. Journal of Chemometrics: A Journal of the Chemometrics Society, 18(10):441-450 2004.
[48] Di Wang. Fast Approximation Algorithms for Positive Linear Programs. University of California, Berkeley, 2017.
[49] Meng Wang, Weiyu Xu, and Tang. unique "nonnegative" solution to an underdetermined system: From vectors to matrices. IEEE Transactions on Signal Processing, 59(3):1007-1016, 2011. doi: 10. 1199TTSP.200.228992244
[50] Michael W Wukitsch, Michael T Petterson, David R Tobler, and Jonas A Pologe. Pulse oximetry: analysis of theory, technology, and practice. Journal of clinical monitoring, 4(4):290-300 1988.
[51] Neal Young. Sequential and parallel algorithms for mixed packing and covering. In Proc. IEEE FOCS'01, 2001.
Checklist
For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] We summarized the contributions and reviewed related work Introduction. Section and 6 contain the full technical details for our theoretical results. Section 5 provides full details for numerical experiments, while Section 7 provides full technical details for the implementation version of our algorithm, including the associated discussion about maintaining the claimed theoretical guarantees.
(b) Did you describe the limitations of your work? [Yes] The main examples are: (i) that similar to ours (although weaker) results can be obtained with existing accelerated coordinate descent algorithms once the structural properties we prove are taken into account, Introduction (Section 1); (ii) the discrepancy between the theoretical and empirical guarantee depending on the batch size as a consequence of modern processor architecture, in Remark 3.6 and Section 5; and (iii) that further experiments are left for future work, in Section 5.
13

---

(c) Did you discuss any potential negative societal impacts of your work? [N/A] This is a theoretical optimization paper for generic non-negative least squares problems, and as such it is not tied to any specific application that could lead to a negative societal impact.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] We stated all the assumptions in Section 2.
(b) Did you include complete proofs of all theoretical results? [Yes] We provided a detailed technical overview in Sections and 4, and included complete proofs in Sections 6 and
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL.? [Yes] We included the link of our code in Section 5, and the data source used in this paper is also cited in Section 5.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? (Yss] We explained the data preprocessing
Section 5.
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [No] The performance of all the algorithms
consistent many runs and thus the error bars would not have been particularly informative.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] We included this information in Section 5. 4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you cite the creators? [Yes] Yes, we cited [9]; see Section 5.
(b) Did you mention the license of the assets? [N/A] The datasets which is an open source library under the 3-clause BSD license.
from LibSVM,
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] We included the link of the code.
(d) Did you discuss whether and how consent was obtained from people whose data you'l 're using/curating? [N/A] The data use is from LibSVM, which is open resource database.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] According the data's 's description, there is no such issue.
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "tvqq39bsURR", "reviewer": "Reviewer_nhFt", "summary": "This work studies Nonnegative (linear) least squares problems (NNLS) where the data itself is nonnegative. This work shows an interesting result that the nonnegativity of data makes the problem easier. In particular, they are able to present a nice and clean result of error bound (or complexity bound given an $\\epsilon$ error) that is independent of any matrix constants (for example operator norm of the matrix), in this nonnegative coefficient/data case. With a new adaptive restart strategy, a linear convergence theory is also established. As the authors already mentioned, this is the first theoretical guarantee for NNLS+ (with nonnegative coefficient matrix $A$) that simultaneously satisfies the properties of being scale-invariant, accelerated, and linearly-convergent. Numerical experiments demonstrate the superior performance of the proposed method on large-scale data. ", "questions": "1. The authors did a great job on improving the time complexity dependence from $O(1/\\sqrt{\\epsilon})$ to $O(\\log(1/\\epsilon))$ by introducing adaptive restart. I wonder if $O(\\log(1/\\epsilon))$ (linear convergence) the best one can do in theory, or, is superlinear convergence even possible in the NNLS/NNLS+ setting?\n\n2. The running time of the proposed method seems to be quite sensitive to the batch size $bs$. Ideally, different $bs$ should be used for handling different sparsity levels of matrices. Is there a way to (roughly) choose this parameter before running the algorithm? It seems to be a complicated problem, as it also depends on how the hardware works.\n", "limitations": "Yes.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. This work has a very interesting and novel idea that uses the assumption of nonegativity of data to accelerate the algorithm and establish the improved theory.\n\n2. The theory looks very nice and clean. Most importantly, the theoretical results look strong. In particular, the time complexity bounds for the proposed algorithm(s) are scale invariant. Indeed, they primarily rely on nnz($A$), and do not rely on its matrix norm. This nnz($A$) dependence also suggests that the theory and the proposed methods handle sparse matrices well (even though their dimension could be large).\n\n3. The results from the numerical experiments look impressive.\n\n4. This work is clearly motivated, well-organized, and very well-written. The proofs also look correct. I really enjoyed reading this paper.\n\nWeaknesses:\n\nGiven the current state of the paper, the only \"weakness\" I can think of is its nonegative assumption of data, which may limit its applications. However, I would NOT say that it is really a weakness, since it is this nonegative assumption that brings advantages to the proposed methods (it also motivates the entire paper). \n\nIn summary, this work has made a strong contribution to the field of NNLS. It is well-written, and its idea is novel. Thus, I highly recommend accepting this manuscript.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "fqThdIcku_I", "reviewer": "Reviewer_5CSB", "summary": "The paper considers the classical non-negative least squares problem, and makes an observation that when the \ndata itself is also non-negative, this considerably simplifies the problem.  It proposes an algorithm inspired\nby primal-dual variance reduction, and establishes scale-invariant (not dependent on matrix norms) multiplicative\nconvergence guarantees, and linear convergence using an adaptive restart technique.", "questions": "7) You report convergence in terms of natural residual.  Why is this a good measure (beyond being convenient for analysis).  It's worth \ncommenting how well it translates to other standard convergence metric -- especially for the experimental section. Does your method work especially well under the lens of natural residual? \n\n8) Algorithm 2 is not defined -- is that the algorithm in Theorem 4.1?  \n\n9) For restart versions we choose the best batch size in {10,50,300,500}.  For a solver to be used in practice -- how do you know (without \nrunning all versions) which batch size to use -- any heuristics? ", "limitations": "I would've liked to see more discussion of limitations and more thorough comparison with competing approaches. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Key strengths of the paper:\n* pointing out that NNLS is significantly easier with non-negative data (called NNLS+ in the paper). \n* Developing a new optimization algorithms for NNLS+, and its detailed convergence analysis. \n* Numerical experiments show cases with significant improvement upon existing  methods (some comments below).\n* The paper is clearly written.\n\nWeaknesses:\n\n1) I'm not fully convinced of the prevalence of NNLS+ (as opposed to plain NNLS) in applications. Taking NNLS+ and \nsimply subtracting the mean, de-trending, or removing first principal component (or any kind of residual) will in general induce \nnegative values in the data. Most practical applications of NNLS I've seen included some negative data. Even in the \nexperiments both examples included negative data (1st example is indeed trivially addressed, while 2nd is not NNLS+). \nAlso, since you applied the approach to plain NNLS in experiments -- can you comment on whether you expect the approach \nto be competitive for NNLS at large?\n\n2) As you mentioned,  an important application of NNLS+ is to NMF.  Is the proposed approach competitive w.r.t. state-of-the-art \nNMF solvers? (Are these based on alternating NNLS -- or are there other approaches as well?)\n\n3) It's worth mentioning other classical approaches to deal with non-negaitve targets, e.g. generalized linear models (e.g. Poisson \nregression). \n\n4) The paper claims that plain LS necessarily scales with matrix constants (e.g. spectral norms). In addition to iterative solution for\nLS there are also exact direct solvers, which have finite termination.  Direct solvers can be very efficient for sparse matrices\nwith good structure (small fill-in).  Also preconditioners can dramatically increase convergence of iterative methods -- especially \nfor near-chordal graphs.  So the statement seems somewhat too generic or misleading. \n\n5) There is some relevant work that would be helpful comment on / cite:  \nThere is a large body of work on multiplicative-update methods, e.g. see review by Satyan Kale, Elad Hazan and Sanjeev Arora, \n\"The Multiplicative Weights Update Method: A Meta-Algorithm and Applications\", which have a variety of applications (including NMF, \nonline convex optimization, packing and covering).\nFor NMF specifically, there's work on separable (anchored) NMF, which has nice error guarantees, (Donoho et al, and Abhishek Kumar \nand Vikas Sindhwani, \"Fast Conical Hull Algorithms for Near-separable Non-negative Matrix Factorization\".\n\n6) For numerical experiments -- you compared to plain FISTA, but it's well known that practical convergence can often be dramatically improved with adaptive restarts, e.g. \"Adaptive Restart for Accelerated Gradient Schemes\", Donoghue and Emmanuel Candes.  Since you employed adaptive restarts in your method -- it seems only fair to use them for FISTA as well.   Furthermore, it'd be helpful to improve the numerical experiments -- run on more examples, and comment why it performs well on certain problems and less well on others (e.g. in terms of runtime). Why are FISTA and OA+DS are so poor on E2006train -- is this due to ill-conditioning / poor scaling?   Also mentioning timing results w.r.t. competing code is not always informative -- as e.g. using advanced compiler techniques (SIMD/AVX512, e.t.c.) vs. say python for-loops can make orders of magnitude runtime difference in itself.  (So reporting that fast NNLS [5] took 6.63 days may mean that the algorithm is poor, or that the implementation is poor, or both). ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "f9leBAjz-sO", "reviewer": "Reviewer_mKSw", "summary": "The paper studies the important non-negative least squares problem. The authors obtain stronger guarantees than for traditional least squares by assuming that the data matrix A is also non-negative. They showed that the problem can be solved up to multiplicative error and with complexity that is independent of any matrix constants.", "questions": "Line 29: \"without the need to tune a regularization parameter or perform cross-validation [33, 6, 13]\". I would add the references https://ieeexplore.ieee.org/document/8022909, https://ieeexplore.ieee.org/document/5608519 and https://arxiv.org/abs/1901.05727 given the theoretical importance that they have for the problem.\n\nLine 32: \"This is true, for example, in problems arising in image processing, computational genomics, functional MRI\". It would be very interesting to have at least one example where this is clearly the case and/or references for each of these applications where it is crucial to have positive data since this is crucial for the paper.\n\nLina 40: Even though it seems to be correct, it is not well written. I would suggest rewriting this part.\n\nLine 48: what is nnz(A)? Is it the number of nonzero elements of A? It should be clearly defined already at line 48.\n\nLine 57: \"What is significant about Theorem 1.1 is the invariance of the computational complexity to the scale of A—it does not depend on any matrix constants\". Could you please explain the nature of this invariance? It is something that I didn't understand from the manuscript. Why does the positivity data give this invariance \"for free\" to the least squares problem? Would it be the case for other constraints as well?\n\nLine 87: \"for problems with strong convexity, a property not satisfied by underdetermined systems\". If the matrix satisfies RIP, this is a form of strong convexity, at least when applied to sparse vectors. This sentence needs to be changed and more should be commented on that issue. I think that TNT-NN still works well for problems without strong convexity, even though there are no theoretical results to support it.\n\nLine 100: \"[16] and [24] mentioned in the preceding paragraph provide only asymptotic convergence guarantees\". As far as I know, reference [24] (or its auxiliary paper https://ieeexplore.ieee.org/abstract/document/8425520) does not establish convergence of TNT for NNLS. It is a good algorithm in practice but without any theoretical guarantees.\n\nLine 104-114: It would be simpler to visualize these algorithms on a table with the different results and assumptions that each one of them needs.\n\nLine 142: I understand the change of variables. But matrix A will change as well. This should be stated in the paper. And what are the implications of this fact?\n\nI dont understand the Lagrangian defined after line 158 and its construction is a crucial part for the analysis, if I understood correctly. I understand that you recover your primal problem from it but I don't understand why it is natural or what is the connection with the \"standard\" Lagrangian. Am I missing something obvious? \n\nLine 182: Implementation version of Lazy SI-NNLS+. Is \"Lazy\" part of the name of the algorithm? If I understand correctly, in line 701, you described that the update of \\tilde(x)_k and \\tilde(y)_k would not be efficient if the matrix is sparse and, therefore, you have a lazy version of the algorithm. It would be interesting to briefly comment on this difference in the main text. I understand that on line 744 you mention the equivalence between Algorithms 1 and 2 but this should be in the main text somehow.\n\nLine 227: You employ a randomized coordinate update in order to reduce the per-iteration complexity. But can you get a deterministic result as well without the randomized coordinate update?\n\nLine 277: \"Numerical experiments and discussion\". I understand that there is a lack of space but for this kind of contribution I think it is crucial to have detailed numerical experiments. Or have more extensive experiments in the appendix. In the case of NNLS, a comparison with the Lawson-Hanson algorithm is very important given how much it is used in the literature. Maybe a comparison with TNT-NNLS as well since it was mentioned and it works really well in practice. I would also add a comparison with PGD and accelerated PGD with restart. I think that a comparison of all algorithms in situations where the matrix does not have RIP as well as matrices that are very ill-conditioned should be performed. And tests with synthetic and simple data to understand better what is going on. Also, differences between problems with sparse ground-truth versus dense ground truth. This is really important to convince a general audience of the usefulness of the algorithm in different situations.\n\nSince the matrix A has all the entries positive, the matrix obviously satisfies the M_{+} condition (see, e.g., https://ieeexplore.ieee.org/document/5608519). What happens to your algorithm when the matrix A does not satisfy the M_{+} condition but does not have all the entries positive? What is the behaviour of the algorithm for the case that the ground truth is sparse?\n", "limitations": "As mentioned before, it is not clear how important the NNLS+ problem is and how the technique developed in this paper would be useful for NNLS problems under more general assumptions.\n\nAlso, it is not clear if the results are a simple application/modification of the VRPDA algorithm or not.\n\nMoreover, the numerical section is quite limited and more extensive experiments must be performed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is clear and develops a quite fast algorithm for the NNLS+ problem. It shows convergence by using a nice combination of a relaxed gap function combined with a random coordinate update. \n\nThe problem for me is how relevant the problem NNLS+ is. I think that the argumentation should be stronger because it is not clear, in principle, how and when this particular instance of NNLS would be interesting. Very precise references linking it to the applications should be presented. For example, Proposition 2.1 is very strong and it is not so clear what would happen without it (e.g. the choice of the box-constrained property of the optimum).\nThe numerical section needs to be reformulated, in my opinion, and more extensive (and simpler) numerical results should be presented, as I discussed in the \"questions\" section. Sparse vs dense ground truths usually present different behavior for the NNLS problem. What happens with SI-NNLS+in this case?\n\nEven though something was mentioned on line 180, it would be important to stress, at least in the appendix, why this new algorithm is not a simple application/modification of the VRPDA algorithm in order to make the importance of novelty even more concrete.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "8-4aPLGx_b", "reviewer": "Reviewer_mTbb", "summary": "The paper considers NNLS (nonnegative least squares) problems with an element-wise nonnegative data matrix and argues that in this scenario it is possible to obtain stronger guarantees than for traditional least squares. To support this claim, a scale-invariant iterative algorithm is provided that constructs an epsilon-multiplicative approximate solution. To obtain their approximation guarantees, the authors use the iterates of the algorithm to construct a sequence of upper estimates of the duality gap. The gap estimate itself is computed as a difference between an upper and lower bound on the Lagrangian. The gap estimate then enables them to obtain a bound on the primal error. The authors then describe how their proposed method can be combined with adaptive restart to obtain a linear convergence rate. The paper concludes with numerical experiments where the proposed method is compared to several competing optimization algorithms. The competing algorithms are outperformed in terms of runtime on one of three datasets while being competitive on the other two.\n", "questions": "It is a bit confusing that in the plots in Figure 1, there is no consistent choice of batch sizes even for the same dataset, e.g. for the new20 dataset, the batch size 50 is used in the plot for the number of data passes (e) but batch size 10 is used for the time plot (f). How were the batch sizes chosen for each plot? The effect of the choice of batch size could be further investigated in an experiment.\n \n", "limitations": "The paper does not contain any discussion of the societal impact of the paper.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- Originality:\nThe paper provides a novel algorithm and proves its convergence. Moreover, the authors show how it can be adapted by combining it with a restart mechanism to achieve linear convergence. In the experiments it is shown to achieve competitive results.\n\n- Quality:\nThe theoretical analysis seems to be sound but I was not able to check it in detail. The experimental evaluation is a bit weak and could be extended to more datasets.\n\n- Clarity:\nThe structure of the paper could be improved, as the algorithm Lazy SI-NNLS+ is only defined in the appendix but already used in Theorem 4.1. It should at least be stated in the main paper that Lazy SI-NNLS+ is in fact equivalent to Algorithm 1, as proven in the Appendix. \n\n- Significance:\nThe paper's contribution is mostly technical and allows the authors to obtain stronger theoretical guarantees compared to earlier work. The techniques used to obtain those improved theoretical guarantees (including the use of the adaptive restart technique) may be of value for the optimization and machine learning community. On the other hand, the experimental evaluation shows that the proposed algorithm is competitive but not consistently better than competing algorithms.\n\n\n----\nEDIT: Taking into account the feedback of the authors during the discussion phase, I have increased my score.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept