5
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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.

non-negative least-squares solution

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Kagaratsch Aug 27 at 20:14
  • $\begingroup$ 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! $\endgroup$ – Roman Aug 27 at 20:29
  • 1
    $\begingroup$ @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. $\endgroup$ – Carl Woll Aug 27 at 21:46
  • $\begingroup$ @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... $\endgroup$ – Kagaratsch Aug 27 at 22:47
1
$\begingroup$
Minimize[{expr.expr, Thread[X > 0]}, X]
$\endgroup$
  • $\begingroup$ This runs for a looong time and never seems to finish... $\endgroup$ – Kagaratsch Aug 27 at 11:31
  • $\begingroup$ Yes it's a tough problem. Runs fine for nn=100. $\endgroup$ – Roman Aug 27 at 15:50

Your Answer

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.