{'fraction_non_alphanumeric': 0.05677632414581481, 'fraction_numerical': 0.01511494982852788, 'mean_word_length': 3.9341967829538333, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerating Non-Negative and Bounded-Variable Linear Regression Algorithms with Safe Screening

Non-negative and bounded-variable linear regression problems arise in a variety of applications in machine learning and signal processing. In this paper, we propose a technique to accelerate existing solvers for these problems by identifying saturated coordinates in the course of iterations. This is akin to safe screening techniques previously proposed for sparsity-regularized regression problems. The proposed strategy is provably safe as it provides theoretical guarantees that the identified coordinates are indeed saturated in the optimal solution. Experimental results on synthetic and real data show compelling accelerations for both non-negative and bounded-variable problems.


Introduction
Due to their fundamental importance in many fields, boxconstrained linear regression problems-including in the large sense the popular non-negative least squares (NNLS) problem-have received considerable attention for several decades (Lawson & Hanson, 1974;Stark & Parker, 1995;Chen & Plemmons, 2010). They are in particular relevant for under-determined linear regression, when the number of variables is larger than the number of measurements. For example, it has been shown that for a variety of problems, the sole non-negativity constraint can be as efficient as sparsitybased regularization (Slawski & Hein, 2013), (Boyer et al., 2019, Proposition 4.1).
Contributions and outline In this work, we propose a generic strategy to accelerate existing solvers for boxconstrained linear regression. It relies on the early identification of saturated coordinates (those attaining the box limits in the solution vector) during the course of iterations, akin to safe screening techniques for sparse-regularized regression problems (Ndiaye et al., 2017;Dantas et al., 2021). Our general framework takes as a special case the safe squeezing approach of Elvira & Herzet (2020) for ∞ -regularized problems (see the details in Appendix A). The general optimization problem addressed in the paper is formulated in Section 2. Then, we show in Section 3 how some saturated coordinates in the solution vector can be safely identified from a primal-dual feasible pair of vectors. This allows us to derive a generic dynamic safe screening procedure for box-constrained linear regression problems. In Section 4, we focus on the dual update step, which requires specific developments in the case of non-negativity constraints. Finally, numerical experiments are reported in Section 5. They demonstrate how our proposed approach can significantly accelerates various solvers of the literature.
The ith entry of a vector z ∈ R n is denoted z i (or sometimes [z] i to avoid ambiguities). Given a subset of indices S ⊆ [n] with cardinality |S| = s, z S ∈ R s denotes the restriction of z to its entries indexed by the elements of S. For a matrix A ∈ R m×n , a j ∈ R m stands for its jth column and A S ∈ R m×s is the matrix formed out of the columns of A indexed by the set S. Vector inequalities are taken coordinate-wisely, i.e., a ≤ b means that a i ≤ b i , ∀i. Given two vectors l ∈ R n and u ∈ R n , we define the set 1 [l,u] = [l 1 , u 1 ] × · · · × [l n , u n ]. We use the notation [z] + (resp., [z] − ) to refer to the positive (resp., negative) part operation defined as max(0, z i ) (resp., min(0, z i )) for all i ∈ [n]. For a convex function f : R n → (−∞, +∞], f * (u) = sup z∈R n z, u − f (z) denotes its Fenchel-Legendre transform (or conjugate function).

Box-constrained linear regression
Given a matrix A ∈ R m×n and a data vector y ∈ R m , we consider the generic box-constrained linear regression problem where l ∈ R n and u ∈R n withR = R ∪ {∞} the Alexandroff extension of the real number line. Given y, the function f (z; y) is assumed proper, lower semi-continuous, convex, and differentiable with respect to (w.r.t.) z. We also assume it has Lipschitz gradient with constant 1/α. We define F (z; y) = i f (z i ; y i ), referred to as loss function. In the rest of the paper, it will be treated as a function of z only. Problem (1) encompasses both bounded-variable linear regression (BVLR) when u ∈ R n and non-negative linear regression (NNLR) when l i = 0 and u i = ∞ ∀i ∈ [n]. Note that our framework can also account for mixed constraints where some entries of u are finite and the others are infinite.

Early identification of saturated variables
Building upon the Gap safe screening technique (Ndiaye et al., 2017) for sparse linear regression, we propose a generic approach (Algorithm 1) to accelerate solvers for (1) through the early identification of saturated coordinates (i.e., identification of j ∈ [n] such that x j = l j or x j = u j when u j < ∞).

Dual problem
Let J u ∞ := {j ∈ [n] : u j = ∞} denotes the set of indices for which the upper bound constraint in (1) is infinity. Then, the dual problem of (1) reads where the dual objective function and the dual feasible set are given by (see Appendix B) Note that the dual solution θ is unique (thanks to the differentiability of f ).
From (3) and (4), we see that when u ∈ R n (BVLR), the dual problem is unconstrained, that is F D = R m .
In contrast, when l = 0 and all the entries of u are equal to ∞ (NNLR), the dual cost function simplifies as Remark 1. Note that if u = −l = δ1 m , the second and third terms of (3) sum up to δ A T θ 1 and we obtain a form of Lasso problem. More generally, the dual problem of BVLR can be seen as a further generalization of the Lasso, where the penalization on the vector A T θ ∈ R n is weighted by l (resp., u) for its negative (resp., positive) entries.

Safe identification of saturated variables
From the primal (1) and dual (2) problems, the first-order primal-dual optimality conditions for a primal-dual solution pair (x , θ ) ∈ 1 [l,u] × F D are given by (see Appendix C): . (6) Let us emphasize that, for NNLR, the second case of (6) never occurs (as x j < ∞ = u j ).
Upon knowledge of the dual solution θ , the optimality condition (6) (known as sub-differential inclusion) constitutes a natural criterion to identify saturated coordinates of the primal solution x . More precisely, we have Although this criterion cannot be used in practice as the dual solution θ is not known in advance, a relaxed version of it can be obtained with only a partial knowledge of the location of θ . More precisely, upon knowledge of a region R ⊂ R m such that θ ∈ R-referred to as safe region-we can define a weaker version of (7) as Clearly, the smaller the region R θ , the larger the number of saturated variables that can be identified. A convenient and efficient choice for R is presented in the next section.

The Gap safe sphere
The Gap safe sphere-initially proposed by Ndiaye et al. (2017) in the context of sparse linear regression-is defined for any primal-dual pair (x, θ) ∈ 1 [l,u] where the duality gap is given by Let us recall that α in (9) is the inverse of the Lipschitz constant of the gradient of f . Equivalently, this means that −f (and hence D) is α-strongly concave. This safe region leads to state-of-the-art screening performances due, in particular, to two key properties.
• A simple geometry which allows to simplify (8) as This significantly limits the computational overhead that is introduced when testing the validity of this criterion (lines 11 and 12 in Algorithm 1).
• It vanishes when the duality gap tends to zero. In other words, if strong duality holds (i.e., P (x ) − D(θ ) = 0), then the radius of the Gap safe sphere vanishes as the iterates {x k , θ k } converge to {x , θ }.
A proof that B(θ, r) is indeed a safe region (i.e., that θ ∈ B(θ, r)) was provided by Ndiaye et al. (2017, Theorem 6) in the context of sparse linear regression. Actually, this proof can be directly applied to the linear regression problems with box constraints considered in the present paper.
Remark 2. Note that the Gap safe sphere can also be defined using local strong concavity bounds of D computed on well-chosen subsets of the domain (Dantas et al., 2021, Theorem 5). This would allow to extend the applicability of Algorithm 1 (following Dantas et al. (2021)) to a more general class of functions f such as the β-divergences with β ∈ [1, 2), that includes in particular the popular Kullback-Leibler divergence (Cichocki & Amari, 2010).

Resulting screening algorithm
The proposed safe screening approach for Problem (1) is presented in Algorithm 1. It can be deployed with any iterative solver for (1), as indicated by the generic notation This has to be understood as performing few iterations of a given primal solver on F (A · +z; y) from the initial point x.
The vector η contains the hyperparameters of the solver (e.g., step sizes) and we will make the role of z ∈ R m explicit hereafter.
In Algorithm 1, Θ : R n → F D denotes a function that computes a dual feasible point from the current primal point x.
Details on how to define this function will be provided in Section 4. The quantity A ⊆ [n] refers the preserved set, which is the complement of the set of screened coordinates. Starting from A = [n], it is dynamically reduced (line 15) by removing the components that are surely identified as being saturated in the solution vector x * (lines 11 and 12), according to the safe rule (11). At the same time, these identified saturated components in x are permanently set to their optimal value (line 13), and their contribution to the vector of measurements is stored in z ∈ R m (line 14).
It follows that the computation of can be reduced from O(mn) to O(m(|A| + 1)) in the next calls to the primal solver. As such, the more saturated components are identified, the larger is the speed improvement in the next calls to the primal solver. Remark 3. The introduction of the variable z in Algorithm 1 is convenient to present a generic algorithm that encompasses the complete class of functions f considered in this work. Yet, this additional variable can be discarded for some loss functions. For example, when F (Ax; y) = Ax − y 2 2 , line 14 can be replaced by y ← y − A S l ∪Su x S l ∪Su , thus avoiding the use of z as well as the addition operation in (12).
It is worth mentioning that the acceleration of the primal iterates provided by the screening procedure has to be balanced with the computational overhead of the screening step itself. This mainly concerns the cost of computing the inner products a T j θ at lines 11 and 12 which, all together, have a complexity of O(|A|m). Fortunately, provided a suitable choice of the dual update function Θ(x) (see Section 4), most standard primal solvers already require the computation of these inner products. They can thus be reused for free in the screening step.
Finally, let us emphasize that Algorithm 1 can be simplified in the NNLR case as shown in Algorithm 2. Indeed, as J u ∞ = [n], we always have S u = ∅. Moreover, because l = 0, the vector z remains always zero.

Computing a dual feasible point
A crucial step of the proposed screening procedure lies in the computation of a dual feasible point (function Θ : R n → F D at line 9 of Algorithm 1). One can see from the definition of the Gap safe sphere in (9) that the closer Θ(x) is to θ , the smaller the safe region is likely to be. Moreover, the computation of Θ(x) should be cheap in order to minimize the computational load of the screening step. In this section, we present ways of defining Θ so as to meet these two desirable properties. Without loss of generality, we focus only on the BVLR and NNLR cases, the case with mixed constraints being easily deduced from the latter two. r ← 2 Gap(x, θ)/α Safe radius 11: 13: x S l ← l S l , x Su ← u Su Set saturated entries 14: z ← z + A S l ∪Su x S l ∪Su Update saturated part

12:
A ← A\S l Update preserved set 13: until Gap(x, θ) < ε gap 4.1. Dual update for BVLR As pointed out in Section 3.1, the dual problem of BVLR is unconstrained, i.e., F D = R m . As such, any point θ ∈ R m is admissible to be a center for the Gap safe sphere. Following the so-called dual scaling idea (Ndiaye et al., 2017), we propose to define Θ as The notation ∇F (Ax; y) refers to the gradient of the function F (·; y) evaluated in Ax. Note that here no scaling of this gradient is required as F D = R m . The rationale behind this choice is twofold. First, we get from the primal-dual link (5) that Θ(x) → θ as x → x . Second, for any first-order primal solver, is computed during the primal update step. This vector contains nothing else than the inner products needed for the screening test (lines 11 and 12 of Algorithm 1) and can thus be reused for free.

Dual update for NNLR
Computing a dual feasible point for the NNLR problem is more involved than for BVLR. Here, F D = {θ ∈ R m : A T θ ≤ 0} and dual scaling is no longer possible as, In other words, if −∇F (Ax; y) is not a feasible dual point (i.e., ∃j s.t. −a T j ∇F (Ax; y) > 0), then any scaled version of it remains not feasible.
Instead, assuming that Int(F D ) = ∅ (i.e., that the interior of F D is nonempty, see Remark 4), we propose a dual translation strategy (by analogy with dual scaling) defined, for any vector t ∈ Int(F D ), as where Ξ t is the translation in the direction of vector t Proposition 1 shows that this Θ indeed maps x onto F D and leads to the desired convergence property. Moreover, similarly to (14), this dual translation allows one to reuse in the screening test some quantities computed during the primal update. Note that the additional inner products A T t can be pre-computed, which keeps A T Θ(x) as cheap as with the standard dual scaling procedure.
Proposition 1 (Validity of the dual translation). Let Θ be defined as in (16). Then, for any primal point where we used the fact that t ∈ Int(F D ), i.e., a T j t < 0 ∀j ∈ [n]. Finally, by continuity of ∇F (A·; y) we get that −∇F (Ax; y) → −∇F (Ax ; y) = θ as x → x .
Remark 4 (Comment on Int(F D ) = ∅). One may wonder to what extend such a condition is restrictive in practice. From the expression of F D in (4), we get that Int(F D ) = ∅ is equivalent to In other words, the columns of A must belong to the interior of a half space of R m containing the origin on its boundary. Hence, if Int(F D ) = ∅ (i.e., (18) fails), we have 0 ∈ conv{a j } n j=1 and the NNLS problem is ill-posed as it admits infinitely many solutions. If furthermore 0 ∈ Int conv{a j } n j=1 , we have cone{a j } n j=1 = R m meaning that the non-negativity constraint is useless (Slawski & Hein, 2013). Note that the latter case corresponds to F D = {0}. To conclude, relevant NNLS problems satisfy Int(F D ) = ∅.
It remains to discuss how one can determine a point t ∈ Int(F D ). Although there is no systematic approach for a general A ∈ R m×n , this can be achieved on a case-by-case basis for many relevant matrices A.
Proposition 2. The following types of matrix A ensure that Int(F D ) = ∅. Moreover, a vector t ∈ F D can be easily computed.

1.
A ∈ R m×n with rank(A) = n ≤ m. Then, for any b ∈ R n <0 , all the solutions of A T t = b (there exists at least one) satisfy t ∈ Int(F D ).

2.
A ∈ R m×m orthogonal (i.e., A T A = I). Then, any negative linear combination of the columns of A, i.e., t = j∈[n] β j a j for β ∈ R n <0 , satisfies t ∈ Int(F D ).

≥0
with non-negative entries and no column with only zeros. Then, any negative vector t ∈ R m <0 satisfies t ∈ Int(F D ).

4.
A ∈ R m×n such that A T A contains a column (say the jth) with all entries positive. Then, t = −a j satisfies t ∈ Int(F D ).
Proof. We prove each case independently.
1. The existence of solutions for A T t = b comes from the fact that the columns of A T span R n (as rank(A) = n). Then, as b ∈ R n <0 , we have A T t = b < 0 and thus t ∈ Int(F D ).
2. Due to the orthogonality of A, we have that ∀i = j, a T i a j = 0. Hence, for any β ∈ R n <0 , we have which shows that t = j∈[n] β j a j ∈ Int(F D ).

Given that ∀i ∈ [n]
, a i ∈ R m ≥0 \{0}, we have, for any negative vector t ∈ R m <0 , that ∀i ∈ [n], a T i t < 0 which shows that t ∈ Int(F D ).
We use projected gradient descent (Polyak, 2015) to solve BVLS. For NNLS, we consider both the coordinate descent method of Franc et al. (2005) and the lsqnonneg routine of MATLAB (a variant of the original algorithm of Lawson & Hanson (1974)). The algorithms are stopped when the duality gap falls below 10 −6 . For all baselines without screening, the duality gap has been computed offline so as not to impact the measured execution times.
Section 5.1 reports results that illustrate how screening performance varies with specific experimental parameters. Then, Section 5.2 is devoted to the evaluation of the proposed screening with real datasets. Figure 1 presents the typical evolution of the speedup factor achieved with screening as a function of the saturation ratio (i.e., s/n where s ∈ N is the number of saturated components in x ). As expected, the higher the number of saturated components, the best the screening performance. Yet, there is a critical value of the saturation ratio under which the computational overhead of screening dominates the acceleration of the primal update, leading to an overall "speedup" below 1.

INFLUENCE OF THE SIZE OF THE PROBLEM
We report in Table 1 execution times for a NNLS problem with increasing size and fixed saturation ratio. For both active-set and coordinate descent algorithms, we observe that the screening performance (i.e., the speed improvement) increases with the size and the level of indeterminacy of the problem. The matrix A ∈ R m×n and vector y ∈ R m have been generated according to aij ∼ N (0, 1) and yi ∼ N (0, 1). The saturation ratio is controlled by varying the size of the box b[−1, 1]. Table 1. Execution times and speedup for a NNLS problem with fixed number of rows m = 2000 and different number of columns n. For each instance of the problem, the matrix A ∈ R m×n ≥0 has been generated using aij = |η| with η ∼ N (0, 1). The data vector y ∈ R m ≥0 has been obtained as y = Ax + wherē x ∈ R n ≥0 is such that x 0/n = 0.05 with non-zero entriesxj distributed like aij and ∈ R m is such that i ∼ N (0, 1). In Figure 2, we report for a NNLS problem the screening ratio (i.e., the number of identified saturated components relatively to the size n) as a function of the iteration number, for different choices of the dual translation vector t ∈ Int(F D ). Clearly, the choice of t affects the screening performance. Although the existence of an optimal choice remains an open question, the reported results allow for some intuition. Indeed, denoting a + (resp., a − ) the column of A ∈ R m×n ≥0 that correlates the most (resp., the least) with all other columns, we can see that setting t = −a + leads to significantly better screening performance than t = −a − . Other relevant choices for this example are t = −1 (used in all remaining NNLS examples) and t = − 1 n j a j . This suggests the conjecture that a relevant t should be close to the "central axis" of cone{a j } n j=1 .

BVLS FOR HYPERSPECTRAL UNMIXING
In this experiment, the data vector y is the observed light reflectance spectrum vector of a random pixel from the Cupitre hyperspectral image (Jia & Qian, 2007). The columns of A ∈ R 188×342 ≥0 are set as the reflectance spectra of pure materials from the USGS High Resolution Spectral Library (Clark et al., 2007), using the same pre-processing as in (Bioucas-Dias et al., 2012). The goal of the regression problem is thus to identify and determine the proportions (so-called abundances) of the materials that compose y. Physical observation constraints dictate that the abundances should belong to the interval [0, 1], leading to a BVLS problem (Bioucas-Dias et al., 2012).
The convergence and screening ratio curves of the projected gradient algorithm with and without screening are presented in Figure 3. We observe that, as the screening ratio progressively grows, the iterations become faster and the convergence curve eventually detaches from the baseline.

NNLS FOR ARCHETYPAL ANALYSIS
The NIPS papers dataset contains word counts from 2484 papers published at the NIPS conference between 1988 to 2003 (Globerson et al., 2007). We removed any allzero columns or rows from the original data matrix and normalized its columns. The input data y ∈ R 2483 ≥0 is taken as a random sample of the dataset and all remaining samples form the columns of A ∈ R 2483×14035
The convergence and screening ratio curves of the active set and coordinate descent algorithms are depicted in Figure 4. While we observe a substantial acceleration for the coordinate descent, its counterpart for the active set method is more subdued. Yet, in both cases, the proposed screening strategy allows to accelerate the considered solvers.
Oracle (%) Screened (%) Figure 3. Convergence and screening ratio for the BVLS hyperspectral experiment. A speedup of 1.68 is achieved with the screening procedure.

Conclusion
In the paper we extended the recently fruitful safe screening framework to the intrinsically distinct family of boxconstrained problems. Instead of identifying zero coordinates induced by a sparse regularization term, we manage to safely identify saturated coordinates induced by the constraints. The main technical challenge in the proposed approach lies in the choice of a dual feasible point which is non-trivial when the box limit is allowed to be unbounded (an example being the widespread NNLS problem). We proposed the simple and efficient dual translation procedure to tackle this problem and suggested some practical choices for the translation direction. Determining an optimal translation direction is actually a challenging problem that deserves further studies-inasmuch as the quality of the dual point can decisively influence a gap-based screening approach (Massias et al., 2020).
A. Relation to ∞ -regularization A particular case of (1), obtained by setting u = −l = c1, c > 0, is the ∞ -constrained optimization problem The penalized counterpart of the above problem leads to the ∞ -regularized linear regression problem considered by Elvira & Herzet (2020) x ∈ argmin where the ∞ -norm penalization is controlled by the parameter λ > 0, to which corresponds a certain value of parameter c in the previous constrained formulation. Note that, in (Elvira & Herzet, 2020), only the least-squares case is addressed.

B. Dual problems derivation
As a direct application of the Fenchel duality formalism (Rockafellar, 1970, Theorem 31.3) and the coordinate-wise separability property of the convex conjugate (Hiriart-Urruty & Lemaréchal, 1993, Ch. X, Prop. 1.3.1) we have the following pair of generic primal-dual problems where Ω(x) = 1 x∈ [l,u] . To complete the demonstration, we need to calculate the conjugate function Ω * of Ω.
From the definition of the conjugate function and the sepa- x j w j + j∈J u ∞ sup xj >lj x j w j .
The first term is a sum of sup of linear functions over compact sets. As such, each sup is attained on the boundary. If w j < 0, the sup is attained at x j = l j whereas if w j > 0, it is attained at x j = u j . Overall, the value of each sup is given by l j (w j ) − + u j (w j ) + .
The second term is a sum of sup of linear functions over umbounded sets. Each sup is either attained at x j = l j if w j < 0 or equals +∞ if w j > 0. Each sup is thus equal to l j (w j ) − + 1 wj ≤0 .
Combining all these expressions, we obtain which, with (22), completes the proof.
Finally, the proof is completed with the expression of the sub-differential of Ω j (x j ) = 1 xj ∈ [lj ,uj ] , where the second case only occurs when u j < ∞.