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)?

Answer 1

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.

Answer 2

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.

Answer 3

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} *)

Answer 4

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