Contents

Non-negative least squares

${\displaystyle \operatorname {arg\,min} \limits _{\mathbf {x} }\|\mathbf {Ax} -\mathbf {y} \|_{2}^{2}}$ subject to x ≥ 0.

Here x ≥ 0 means that each component of the vector x should be non-negative, and ‖·‖₂ denotes the Euclidean norm.

Non-negative least squares problems turn up as subproblems in matrix decomposition, e.g. in algorithms for PARAFAC^[2] and non-negative matrix/tensor factorization.^[3][4] The latter can be considered a generalization of NNLS.^[1]

Another generalization of NNLS is bounded-variable least squares (BVLS), with simultaneous upper and lower bounds α_i ≤ x_i ≤ β_i.^{[5]: 291 [6]}

Quadratic programming version[edit]

The NNLS problem is equivalent to a quadratic programming problem

${\displaystyle \operatorname {arg\,min} \limits _{\mathbf {x\geq 0} }\left({\frac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {Q} \mathbf {x} +\mathbf {c} ^{\mathsf {T}}\mathbf {x} \right),}$

where Q = A^TA and c = −A^T y. This problem is convex, as Qispositive semidefinite and the non-negativity constraints form a convex feasible set.^[7]

Algorithms[edit]

The first widely used algorithm for solving this problem is an active-set method published by Lawson and Hanson in their 1974 book Solving Least Squares Problems.^[5]: 291Inpseudocode, this algorithm looks as follows:^[1][2]

• Inputs:
  • a real-valued matrix A of dimension m × n,
  • a real-valued vector y of dimension m,
  • a real value ε, the tolerance for the stopping criterion.
• Initialize:
  • Set P = ∅.
  • Set R = {1, ..., n}.
  • Set x to an all-zero vector of dimension n.
  • Set w = A^T(y − Ax).
  • Let w^R denote the sub-vector with indexes from R
• Main loop: while R ≠ ∅ and max(w^R) > ε:
  • Let jinR be the index of max(w^R)inw.
  • Add jtoP.
  • Remove j from R.
  • Let A^PbeA restricted to the variables included in P.
  • Let s be vector of same length as x. Let s^P denote the sub-vector with indexes from P, and let s^R denote the sub-vector with indexes from R.
  • Set s^P = ((A^P)^T A^P)⁻¹ (A^P)^Ty
  • Set s^R to zero
  • While min(s^P) ≤ 0:
    • Let α = min ⁠x_i/x_i − s_i⁠ for iinP where s_i ≤ 0.
    • Set xtox + α(s − x).
    • Move to R all indices jinP such that x_j ≤ 0.
    • Set s^P = ((A^P)^T A^P)⁻¹ (A^P)^Ty
    • Set s^R to zero.
  • Set xtos.
  • Set wtoA^T(y − Ax).
• Output: x

This algorithm takes a finite number of steps to reach a solution and smoothly improves its candidate solution as it goes (so it can find good approximate solutions when cut off at a reasonable number of iterations), but is very slow in practice, owing largely to the computation of the pseudoinverse ((A^P)^T A^P)⁻¹.^[1] Variants of this algorithm are available in MATLAB as the routine lsqnonneg^[8][1] and in SciPyasoptimize.nnls.^[9]

Many improved algorithms have been suggested since 1974.^[1] Fast NNLS (FNNLS) is an optimized version of the Lawson—Hanson algorithm.^[2] Other algorithms include variants of Landweber's gradient descent method^[10] and coordinate-wise optimization based on the quadratic programming problem above.^[7]

See also[edit]

• M-matrix
• Perron–Frobenius theorem

References[edit]