# Approximate strictly positive solution to a linear set of equations?

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

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.

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? Aug 27, 2019 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! Aug 27, 2019 at 20:29
• @Kagaratsch You could use CloudEvaluate, e.g., define qo[a__] := CloudEvaluate[SystemQuadraticOptimization[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. Aug 27, 2019 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... Aug 27, 2019 at 22:47

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