Approximate strictly positive solution to a linear set of equations?

Question

Consider a positive matrix M and a positive vector b, e.g.

nn = 1000;
M  = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}];
b  = Table[RandomReal[{0, 100}], {i, 1, nn}];

I would like to find a positive vector X

X  = Array[x,nn];

(each x[i]>0) such that given

expr = M.X-b;

the quantity expr.expr is minimized. Is it possible to do that in Mathematica efficiently (so that it finishes within a few seconds/minutes)?

score 9 · Answer 1 · answered Aug 27 '19 at 19:44

You can use the new in M12 function QuadraticOptimization, which minimizes functions of the form:

$$\frac{1}{2} x . q. x + c . x$$

subject to linear constraints on $x$. So, the first step is to figure out what $q$ and $c$ are for your example. To do this we expand your expr, but before doing so, note that for a vector u and a matrix M we have:

$$u.M=M^T.u$$

Then, we can expand expr yielding:

$$(M.X-b).(M.X-b) = X.M^T.M.X-2 b.M.X+b.b$$

Hence we have:

q := 2 Transpose[M] . M
c := -2 b . M

Let's do a small example and compare to @Roman's answer. Setup:

nn = 10;
SeedRandom[10];
M = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}];
b = Table[RandomReal[{0, 100}], {i, 1, nn}];
X = Array[x, nn];
expr = M . X - b;

Comparison:

r1 = X /. Last @ Minimize[{expr.expr, X > 0}, X]
r2 = QuadraticOptimization[{q, c}, {IdentityMatrix[nn], ConstantArray[0, nn]}]

{0.221992, 0.188374, 0.131969, 0., 0., 0., 0., 0., 0.0849646, 0.028771}

{0.221992, 0.188374, 0.131969, 0., 0., 0., 0., 0., 0.0849646, 0.028771}

where I used constraints that ensured that each vector element is nonnegative. For your larger example:

nn = 1000;
SeedRandom[1];
M = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}];
b = Table[RandomReal[{0, 100}], {i, 1, nn}];

res = QuadraticOptimization[{q, c}, {IdentityMatrix[nn], ConstantArray[0, nn]}]; //AbsoluteTiming
res[[;;20]]

{1.48153, Null}

{0., 0., 0., 0.015694, 0., 0.00561439, 0., 0.0157487, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}

So, it finishes on the order of a couple seconds, as requested.

Wow, this is great! Unfortunately, I am on MA 11.3 right now. Any chance you could give me a hint on how to write a custom code similar to what QuadraticOptimization is doing? — Kagaratsch, Aug 27 '19 at 20:14
Carl I was naïvely assuming that internally, my solution's Minimize would be auto-mapped onto QuadraticOptimization for fastest execution. Seeing that it's not so is a bit disappointing! — Roman, Aug 27 '19 at 20:29
@Kagaratsch You could use CloudEvaluate, e.g., define qo[a__] := CloudEvaluate[System`QuadraticOptimization[a]] and then use qo. I think FindMinimum has an undocmented "QuadraticOptimization" method you could try instead if you don't want to or are unable to use the cloud. — Carl Woll, Aug 27 '19 at 21:46
@CarlWoll The cloud returns $Failed, but for FindMinimum you probably mean FindMinimum[{expr.expr, Thread[X > 0]}, X, Method -> "QuadraticProgramming"]? This works fairly quickly for nn=10 and nn=100, but for nn=1000 it still never seems to finish... — Kagaratsch, Aug 27 '19 at 22:47

score 5 · Accepted Answer · answered Aug 27 '19 at 22:49

5

Perhaps you could use the non-negative least-squares algorithm (NNLS) of Lawson and Hanson, discussed here. See the links in the comments.

NNLS was ported to Mathematica by Michael Woodhams, and used here as NNLS[M,b]. I am running Mma 11.3.0 for 64-bit Linux.

nn = 1000;
SeedRandom[10];
M = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}];
b = Table[RandomReal[{0, 100}], {i, 1, nn}];
ls = NNLS[M,b];

Timing was about 10 seconds for nn=1000.

answered Aug 27 '19 at 22:49

KennyColnago

15,209
26
62

This is brilliant! Can be found here: http://forums.wolfram.com/mathgroup/archive/2003/Oct/msg00053.html – Kagaratsch Aug 27 '19 at 23:00

score 2 · Answer 3 · answered Jun 22 '22 at 17:02

THIS IS AN EXTENDED COMMENT RATHER THAN AN ANSWER.

It is inefficient to use Table to generate random numbers.

Clear["Global`*"]
nn = 1000;
SeedRandom[10];
t1 = AbsoluteTiming[
   {M1 = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}],
     b1 = Table[RandomReal[{0, 100}], {i, 1, nn}]};][[1]]
(* 0.202688 *)
SeedRandom[10];
t2 = AbsoluteTiming[
   {M2 = RandomReal[{0, 100}, {nn, nn}],
     b2 = RandomReal[{0, 100}, nn]};][[1]]
(* 0.005653 *)

The values are identical

{M1 === M2, b1 === b2, t1/t2}
(* {True, True, 35.8549} *)

score 1 · Answer 4 · answered Aug 27 '19 at 09:22

1

Minimize[{expr.expr, Thread[X > 0]}, X]

answered Aug 27 '19 at 09:22

Roman

47,322
2
55
121

This runs for a looong time and never seems to finish... – Kagaratsch Aug 27 '19 at 11:31
Yes it's a tough problem. Runs fine for nn=100. – Roman Aug 27 '19 at 15:50

Approximate strictly positive solution to a linear set of equations?

4 Answers4

Linked