Nonnegative Least Squares Algorithm (NNLS) [closed]

Question

Can anyone optimize the code below which is developed in an old version of Mathematica (2003) both in terms of efficiency and adaption to the latest versions of Mathematica?

Description of the code:

Nonnegative least squares (NNLS) is a well-known regression algorithm to fit data as in least-squares, subject to the constraint that all the coefficients are nonnegative. See the full detail at Wikipedia.

Upon my earlier question and suggestions to use LeastSquares and the need for a nonnegative version of it, I found out that Mathematica does not have a built-in command for NNLS.

Luckily, the algorithm due to Lawson and Hanson ["Solving Least Squares Problems" (Prentice-Hall, 1974, republished SIAM, 1995)] was implemented by Michael Woodhams in 2003 and its posted here.

For simplicity, I repost Michael's code below. While the code is free to the public, Michael would

appreciate it if you acknowledge [his] authorship if you use it.

(* Coded by Michael Woodhams, from algorithm by Lawson and Hanson, *)
(* "Solving Least Squares Problems", 1974 and 1995. *)
bitsToIndices[v_] := Select[Table[i, {i, Length[v]}], v[[#]] == 1 &];
NNLS[A_, f_] := Module[
      {x, zeroed, w, t, Ap, z, q, \[Alpha], i, zeroedSet, positiveSet,
        toBeZeroed, compressedZ, Q, R},
      (* Use delayed evaluation so that these are recalculated on the
fly as \
needed : *)
      zeroedSet := bitsToIndices[zeroed];
      positiveSet := bitsToIndices[1 - zeroed];
      (* Init x to vector of zeros, same length as a row of A *)
      debug["A=", MatrixForm[A]];
      x = 0 A\[LeftDoubleBracket]1\[RightDoubleBracket];
      debug["x=", x];
      (* Init zeroed to vector of ones,
        same length as x *)
      zeroed = 1 - x;
      debug["zeroed=", zeroed];
      w = Transpose[A].(f - A.x);
      debug["w=", w];
      While[zeroedSet != {} 
          && Max[w\[LeftDoubleBracket]zeroedSet\[RightDoubleBracket]]
> 0,
        debug["Outer loop starts."];
        (* The index t of the largest element of w, *)
        (* subject to the constraint t is zeroed *)
        t = 
          Position[w zeroed, Max[w zeroed], 1, 
                1]\[LeftDoubleBracket]1\[RightDoubleBracket]\
\[LeftDoubleBracket]1\[RightDoubleBracket];
        debug["t=", t];
        zeroed\[LeftDoubleBracket]t\[RightDoubleBracket] = 0;
        debug["zeroed=", zeroed];
        (* Ap = the columns of A indexed by positiveSet *)
        Ap = 
          Transpose[
            Transpose[A]\[LeftDoubleBracket]
              positiveSet\[RightDoubleBracket]];
        debug["Ap=", MatrixForm[Ap]];
        (* Minimize (Ap . compressedZ - f) by QR decomp *)
        {Q, R} = QRDecomposition[Ap];
        compressedZ = Inverse[R].Q.f;
        (* 
          Create vector z with 0 in zeroed indices and compressedZ
entries \
elsewhere *)
        z = 0 x;
        z\[LeftDoubleBracket]positiveSet\[RightDoubleBracket] =
compressedZ;
        debug["z=", z];
        While[Min[z] < 0,
          (* There is a wart here : x can have zeros, 
            giving infinities or indeterminates. They don't matter, 
            
            as we ignore those elements (not in postitiveSet) but it
will \
produce warnings. *)
          debug["Inner loop start"];
          (* 
            find smallest x\[LeftDoubleBracket]
                  q\[RightDoubleBracket]/(x\[LeftDoubleBracket]q\
\[RightDoubleBracket] - z\[LeftDoubleBracket]q\[RightDoubleBracket])
*)
          (* such that : q is not zeroed, 
               z\[LeftDoubleBracket]q\[RightDoubleBracket] < 0 *)
          \[Alpha] = Infinity;
          For[q = 1, q <= Length[x], q++,
            
            If[zeroed\[LeftDoubleBracket]q\[RightDoubleBracket] == 0
&&
                z\[LeftDoubleBracket]q\[RightDoubleBracket] < 0,
              \[Alpha] = 
                Min[\[Alpha], 
                  x\[LeftDoubleBracket]q\[RightDoubleBracket]/(x\
\[LeftDoubleBracket]q\[RightDoubleBracket] - 
                        z\[LeftDoubleBracket]q\[RightDoubleBracket])];
              debug["After trying index q=", q, " \[Alpha]=",
\[Alpha]];
              ]; (* if *)
            ]; (* for *)
          debug["\[Alpha]=", \[Alpha]];
          x = x + \[Alpha](z - x);
          debug["x=", x];
          
          toBeZeroed = 
            Select[positiveSet, 
              Abs[x\[LeftDoubleBracket]#\[RightDoubleBracket]] <
10^-13 &];
          debug["toBeZeroed=", toBeZeroed];
          zeroed\[LeftDoubleBracket]toBeZeroed\[RightDoubleBracket] =
1;
          x\[LeftDoubleBracket]toBeZeroed\[RightDoubleBracket] = 0;
          
          (* Duplicated from above *)
          (* Ap = the columns of A indexed by positiveSet *)
          
          Ap = Transpose[
              Transpose[
                  A]\[LeftDoubleBracket]positiveSet\[RightDoubleBracket]];
          debug["Ap=", MatrixForm[Ap]];
          (* Minimize (Ap . compressedZ - f) by QR decomp *)
          {Q, R} = QRDecomposition[Ap];
          compressedZ = Inverse[R].Q.f;
          (* 
            Create vector z with 0 in zeroed indices and compressedZ
entries \
elsewhere *)
          z = 0 x;
          
          z\[LeftDoubleBracket]positiveSet\[RightDoubleBracket] = 
            compressedZ;
          debug["z=", z];
          ]; (* end inner while loop *)
        x = z;
        debug["x=", x];
        w = Transpose[A].(f - A.x);
        debug["w=", w];
        ]; (* end outer while loop *)
      Return[x];
      ]; (* end module *)

Answer 1

We can compare 6 algorithms based on FindMinimum,NMinimize, ConvexOptimization, QuadraticOptimization with algorithm by Lawson and Hanson implemented by Michael Woodhams. First code (Michael)

bitsToIndices[v_] := Select[Table[i, {i, Length[v]}], v[[#]] == 1 &];
NNLS[A_, f_] := 
  Module[{x, zeroed, w, t, Ap, z, q, \[Alpha], i, zeroedSet, 
    positiveSet, toBeZeroed, compressedZ, Q, R}, 
   zeroedSet := bitsToIndices[zeroed];
   positiveSet := bitsToIndices[1 - zeroed];
   x = 0 A[[1]];
   zeroed = 1 - x;
   w = Transpose[A] . (f - A . x);
   While[zeroedSet != {} && Max[w[[zeroedSet]]] > 0, 
    debug["Outer loop starts."];
    t = Position[w zeroed, Max[w zeroed], 1, 1][[1]][[1]];
    zeroed[[t]] = 0;
    Ap = Transpose[Transpose[A][[positiveSet]]];
    {Q, R} = QRDecomposition[Ap];
    compressedZ = Inverse[R] . Q . f;
    z = 0 x;
    z[[positiveSet]] = compressedZ;
    While[Min[z] < 0, debug["Inner loop start"];
     \[Alpha] = Infinity;
     For[q = 1, q <= Length[x], q++, 
      If[zeroed[[q]] == 0 && 
        z[[q]] < 0, \[Alpha] = Min[\[Alpha], x[[q]]/(x[[q]] - z[[q]])];
       ]];
     x = x + \[Alpha] (z - x);
     toBeZeroed = Select[positiveSet, Abs[x[[#]]] < 10^-13 &];
     zeroed[[toBeZeroed]] = 1;
     x[[toBeZeroed]] = 0;
     Ap = Transpose[Transpose[A][[positiveSet]]];
     {Q, R} = QRDecomposition[Ap];
     compressedZ = Inverse[R] . Q . f;
     z = 0 x;
     z[[positiveSet]] = compressedZ;
     ]; x = z;
    w = Transpose[A] . (f - A . x);
    ]; Return[x];];

Second code by Jean - Claude Poujade

NNLSFindMinimum[A_, f_] := 
  Module[{nbx = Length[First[A]], xi, x, axf, xinit}, 
   xi = Array[x, nbx];
   axf = A . xi^2 - f;
   xinit = PseudoInverse[A] . f;
   If[And @@ (# >= 0 & /@ xinit), xinit, 
    fm = FindMinimum[Evaluate[axf . axf], 
      Evaluate[Sequence @@ Transpose[{xi, xinit}]], 
      MaxIterations -> 1000];
    xi^2 /. fm[[2]]]];

Code 3

NNLSNMin[A_, f_] := 
  Module[{nbx = Length[First[A]], xi, x, axf}, xi = Array[x, nbx];
   axf = A . xi^2 - f;
   fm = NMinimize[Evaluate[axf . axf], xi];
   xi^2 /. fm[[2]]];

Code 4

NNLSNMinCon[A_, f_] := 
  Module[{nbx = Length[First[A]], xi, x, axf}, xi = Array[x, nbx];
   axf = A . xi - f;
   fm = NMinimize[{axf . axf, Table[xi[[i]] >= 0, {i, Length[xi]}]}, 
     xi];
   xi /. fm[[2]]];

As it well known NNLS problem is equivalent to a quadratic programming problem. With this method we can test

Code 5

NNLSN2Min[A_, f_] := 
  Module[{nbx = Length[First[A]], X, x, axf}, X = Array[x, nbx];
   axf = 1/2 X . (Transpose[A] . A) . X - (Transpose[A] . f) . X;
   sol = NMinimize[{axf, Table[X[[i]] >= 0, {i, nbx}]}, X];
   X /. sol[[2]]];

Code 6

 NNLSNquadratic[A_, f_] := 
  Module[{nbx = Length[First[A]], X, x, axf}, X = Array[x, nbx];
   axf = 1/2 X . (Transpose[A] . A) . X - (Transpose[A] . f) . X;
   sol = ConvexOptimization[axf, Table[X[[i]] >= 0, {i, nbx}], X];
   X /. sol]; 

Code 7

NNLSNquadraticOp[A_, f_] := 
  Module[{nbx = Length[First[A]], X, x, axf}, X = Array[x, nbx];
   axf = 1/2 X . (Transpose[A] . A) . X - (Transpose[A] . f) . X;
   res = QuadraticOptimization[axf, Table[X[[i]] >= 0, {i, nbx}], X];
   X /. res];

Test

SeedRandom[12345];

a = Table[Random[], {i, 100}, {j, 200}];
f = Table[Random[], {i, 100}];

OUT = Array[out, {7}];

With these data we test 7 code

out[1] = NNLS[a, f] // AbsoluteTiming
out[2] = NNLSFindMinimum[a, f] // AbsoluteTiming;
out[3] = NNLSNMin[a, f] // AbsoluteTiming;
out[4] = NNLSNMinCon[a, f] // AbsoluteTiming;
out[5] = NNLSN2Min[a, f] // AbsoluteTiming;
out[6] = NNLSNquadratic[a, f] // AbsoluteTiming;
out[7] = NNLSNquadraticOp[a, f] // AbsoluteTiming;

Compare time and absolute error

ListPlot[Table[out[i][[1]]/out[1][[1]], {i, 1, 7}], 
 ScalingFunctions -> "Log", Filling -> Axis]

ListPlot[Table[Abs[out[i][[2]] - out[1][[2]]], {i, 2, 7}], 
 PlotRange -> {10^-13, 10^-7}, PlotStyle -> PointSize[Medium], 
 Frame -> True, PlotLegends -> Table[i, {i, 2, 7}], 
 PlotMarkers -> Automatic, ScalingFunctions -> "Log"]

As we can see from this picture the fastest code is Code 1 developed by Michael Woodhams in 2003. We can improve Code 1 by replacing While and For loop with Do loop. But improvement is very small about 10% only and eliminates on a large matrix $10^4\times10^4$. Nevertheless, let consider

Code 8

bitsToIndices[v_] := Select[Table[i, {i, Length[v]}], v[[#]] == 1 &];
NNLS1[A_, f_] := 
  Module[{x, zeroed, w, t, Ap, z, q, \[Alpha], i, zeroedSet, 
    positiveSet, toBeZeroed, compressedZ, Q, R}, 
   zeroedSet := bitsToIndices[zeroed];
   positiveSet := bitsToIndices[1 - zeroed];
   x = 0 A[[1]];
   zeroed = 1 - x;
   w = Transpose[A] . (f - A . x);
   Do[If[zeroedSet != {} && Max[w[[zeroedSet]]] > 0,
      t = Position[w zeroed, Max[w zeroed], 1, 1][[1]][[1]];
      zeroed[[t]] = 0;
      Ap = Transpose[Transpose[A][[positiveSet]]];
      {Q, R} = QRDecomposition[Ap];
      compressedZ = Inverse[R] . Q . f;
      z = 0 x;
      z[[positiveSet]] = compressedZ;
      Do[If[Min[z] < 0,
         \[Alpha] = Infinity;
         Do[
          If[zeroed[[q]] == 0 && 
             z[[q]] < 0, \[Alpha] = 
              Min[\[Alpha], x[[q]]/(x[[q]] - z[[q]])];
            ];, {q, 1, Length[x]}];
         x = x + \[Alpha] (z - x);
         toBeZeroed = Select[positiveSet, Abs[x[[#]]] < 10^-13 &];
         zeroed[[toBeZeroed]] = 1;
         x[[toBeZeroed]] = 0;
         Ap = Transpose[Transpose[A][[positiveSet]]];
         {Q, R} = QRDecomposition[Ap];
         compressedZ = Inverse[R] . Q . f;
         z = 0 x;
         z[[positiveSet]] = compressedZ;, Break[]];, {Infinity}
       ]; x = z;
      w = Transpose[A] . (f - A . x);, Break[]];, {Infinity}
    ]; Return[x];]; 

Test 2

SeedRandom[1];
a = Table[Random[], {i, 10000}, {j, 12000}];
f = Table[Random[], {i, 10000}];
OUT = Array[out, {10}]; out[2] = NNLS1[a, f] // AbsoluteTiming; 
out[1] = NNLS[a, f] // AbsoluteTiming;

Table[out[i][[1]], {i, 2}]

Out[]= {180.171, 174.289}