Non-negative least squares

Algorithms

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. In pseudocode, this algorithm looks as follows:

• 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 j in R be the index of max(w^R) in w.
  • Add j to P.
  • Remove j from R.
  • Let A^P be A 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
    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 i in P where s_i ≤ 0.
    • Set x to x + α(s − x).
    • Move to R all indices j in P such that x_j ≤ 0.
    • Set s^P = ((A^P)^T A^P)⁻¹ (A^P)^Ty
    • Set s^R to zero.
  • Set x to s.
  • Set w to A^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)⁻¹. Variants of this algorithm are available in MATLAB as the routine lsqnonneg and in SciPy as optimize.nnls.

Many improved algorithms have been suggested since 1974. Fast NNLS (FNNLS) is an optimized version of the Lawson–Hanson algorithm. Other algorithms include variants of Landweber's gradient descent method, coordinate-wise optimization based on the quadratic programming problem above, and an active set method called TNT-NN.