Matrix left division with constraints?

I have an underdetermined system of equations A*x=b, and I am looking for any non-negative solution. I.e. A is an n-by-m matrix with n < m, and I need to solve the system for x subject to the following two constraints:

1. sum(x) == 1
2. all(x >= 0)

I can solve this using linprog, but it is fairly slow (I have to solve 1e6 equations of this form, and A is relatively large, on my machine this was going to take about 6 days).

On the other hand I just discovered that matrix left division will solve this problem orders of magnitude faster, but I don't know how to enforce the non-negativity constraint (the summing to 1 constraint is easy as an additional row of 1's in A and b). On my machine it took about 2 min to solve all 1e6 systems of equations (but obviously the solutions were not non-negative).

I suppose I can use lsqlin or lsqnonneg, but these seem to be unnecessary for my needs. I know that there are an infinite number of solutionns for every system of equations I want to solve, and I don't care which solution I get (i.e. doesn't have to have the smallest norm, or anything), so long as the non-negativity and summing to 1 constraints are satisfied. Is there a way to take advantage of the speed of matrix left division and also enforce non-negativity?

Answer 1

[Shashank Prasanna]

Your best bet is LSQLIN if you have constraints for a linear system.

http://www.mathworks.com/help/optim/ug/lsqlin.html

Since your have a large system, do you know if it is sparse? If it is then you can make the matrix sparse and LSQLIN will handle it.

Why do you say that they are "unnecessary for your needs" ? Did LSQLIN not work or not provide a solution?

Backslash or MLDIVIDE is optimized for unconstrained linear system solution. I am not aware of a way to modify \ to honor these constraints.

[John D'Errico]

And of course, you cannot simply modify backslash(mldivide) to honor constraints. Except that the tool which does solve that problem is lsqlin OR linprog. For the huge problem size posed by the OP, this is not a surprise. The constraints mean that any speed gains found in backslash are no longer possible.

Big problems take big time.

Answer 2

[Mohamed Abdelhameed]

I am trying to solve an overdetermenind system of equations (45 equations with 18 unknowns) in the form of A*X = B; where:

A = [8.19 0.3; 12.39 0.86; 16.15 1.68; 17.7 2.16; 24.6 8.33];

A is 5 rows, 2 columns

B = [0.85 3.8 0.225 0.175 0.05 0.05 0 1.919 0.165; 1.5 6.825 0.45 0.4 0.175 0.125 0.025 2.907 0.315; 3 7.625 0.9 1.05 0.475 0.225 0.1 3.823 0.544; 3.05 7.9 1.1 1.8 1.225 0.25 0.45 4.153 0.881; 5.1 12.95 1.975 3.65 3.075 0.25 1.075 4.045 1.217];

B is 5 rows, 9 columns

X = [a1 b1 c1 d1 e1 f1 g1 h1 i1;a2 b2 c2 d2 e2 f2 g2 h2 i2];

X is 2 rows, 9 columns

I need to solve X with the following constraints: sum (a1 to i1) = 1 & sum (a2 to i2) = 1 & all elements in X shall be higher than or equal to zero (non negative).

[Benjamin Kelm]

Same question here, haven't found an answer yet.

[Walter Roberson]

Have you experimented with using https://www.mathworks.com/help/optim/problem-based-approach.html ?

[John D'Errico]

@Benjamin Kelm - Why have you not found an answer? This is a trivial probem to solve using lsqlin, at least for a small problem. And sereral people have suggested lsqlin. So saying you have not found an answer just means you have not looked.

It is even easier to formulate using the problem based tools in the optimization toolbox, where you don't even need to understand how to call lsqlin.

[Torsten]

X = sym('X',[2 9]);
A = [8.19 0.3; 12.39 0.86; 16.15 1.68; 17.7 2.16; 24.6 8.33];
B = [0.85 3.8 0.225 0.175 0.05 0.05 0 1.919 0.165; 1.5 6.825 0.45 0.4 0.175 0.125 0.025 2.907 0.315; 3 7.625 0.9 1.05 0.475 0.225 0.1 3.823 0.544; 3.05 7.9 1.1 1.8 1.225 0.25 0.45 4.153 0.881; 5.1 12.95 1.975 3.65 3.075 0.25 1.075 4.045 1.217];
eqn = A*X==B;
[As,Bs] = equationsToMatrix(eqn);
X_without_constraints = double(As)\double(Bs);
norm(double(As)*X_without_constraints-double(Bs))

ans = 1.6303

Aeq = [ones(1,9) zeros(1,9);zeros(1,9) ones(1,9)];
beq = [1;1];
lb = zeros(18,1);
ub = Inf(18,1);
X_with_constraints = lsqlin(double(As),double(Bs),[],[],Aeq,beq,lb,ub);

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

norm(double(As)*X_with_constraints-double(Bs))

ans = 2.2074

X = reshape(X_with_constraints,[9,2]).'

X = 2×9

0.1495 0.4739 0.0504 0.0447 0.0145 0.0109 0.0102 0.2040 0.0419 0.1524 0.1159 0.0651 0.2877 0.3047 0.0000 0.0743 0.0000 0.0000

sum(X,2)

ans = 2×1

1 1