4
$\begingroup$

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 *)
New contributor
Farid Shahandeh is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
14
  • 3
    $\begingroup$ What is the question? $\endgroup$
    – Syed
    Jun 21 at 9:18
  • 2
    $\begingroup$ Thanks for the effort to share useful code, but StackExchange is a questions/answers site, and does not accept arbitrary discussion. mathematica.stackexchange.com/tour I suggest you rewrite this in a Q/A format (self-answered questions are fine) to prevent the thread from getting closed. $\endgroup$
    – Szabolcs
    Jun 21 at 9:40
  • 3
    $\begingroup$ I think that we can test this code and proposed a new one used new functions introduced in time period 2003-2022. $\endgroup$ Jun 21 at 10:13
  • 4
    $\begingroup$ Adding a specific example (i.e., data and proposed model rather than just code) to the question should be essential. $\endgroup$
    – JimB
    Jun 21 at 14:27
  • 3
    $\begingroup$ Note that as an alternative to this method, there is a nice one posted by @CarlWoll in this prior MSE thread. It uses quadratic programming under the hood. $\endgroup$ Jun 21 at 16:02

1 Answer 1

10
$\begingroup$

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"]

Figure 1

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}
 
$\endgroup$
4
  • $\begingroup$ Wonderful job Alex, Thanks! $\endgroup$ 2 days ago
  • $\begingroup$ @FaridShahandeh You are welcome! $\endgroup$ 2 days ago
  • $\begingroup$ I checked all seven codes against my examples. When $f$ is a matrix, the solutions for each row form a matrix, say $S$. It turns out that only the second code returns the $S$ with minimum rank. $\endgroup$ 2 days ago
  • $\begingroup$ @FaridShahandeh Could you add this test to your post for discussion? $\endgroup$ yesterday

Your Answer

Farid Shahandeh is a new contributor. Be nice, and check out our Code of Conduct.

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.