least-squares method with a constraint

Question

I have 37 linear equations with 36 variables in the form of matrix: A x = b. (A has 37 rows and 36 columns.) The equations don't have an exact solution so I have used Matlab to find the closest answer using x = A \ b.

The problem is that I also have a condition, all elements of x should be positive: x_i > 0 for all i. x = A \ b gives negative values for some elements. How can I apply this constraint ?

Here are the concrete values of A and b that I'm working with:

A = [0.83   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0.02   0.63    0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0.02    -0.2    0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0.15   0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0.15    0   0   0   0   0.02    -0.33   0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.18    0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.018   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.17
     1  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1];
b = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]';

If you have 37 equations for 36 unknown, that is an overdetermined system, meaning there is unlikely to have multiple solution (and, actually, highly likely to have no solutions). — , May 20 '15 at 16:02
I know, but can't I at least look for an approximate solution with least error? — nasim, May 20 '15 at 16:05
If you want to find the `X0` that minimizes some error, first you need to define that error. :-) What type of norm would you use for the error `E = B - A*X0`? Euclidean, taxicab, maximum? — , May 20 '15 at 16:12
I don't even know those terms ! (I am reading about them right now) My variables are transition probabilities in a queuing network and usually range between 0.0001 and 0.1. Which one do you recommend? — nasim, May 20 '15 at 16:22
@CST-Link, the title says least squares, this specifies the cost function — A. Donda, May 20 '15 at 19:53
@nasim, can you post the exact values of A and b? This would make it easier to debug a possible answer. In the mean time, check out my solution. — A. Donda, May 20 '15 at 21:01
@nasim, no, not the comments. Edit the question and include it as code. — A. Donda, May 20 '15 at 21:06
@nasim, thanks, but this doesn't make sense. This A is 36x36, and this b is 1x38. Please fix this. — A. Donda, May 20 '15 at 21:11

A. Donda · Accepted Answer · 2015-05-20T21:38:07.043

You want to find an approximate solution x to A x = b in the least-squares sense, i.e. you want to minimize

||A x - b||² = x^T A^T A x + b^T b - 2 x^T A^T b.

Disregarding the constant term b^Tb and dividing by a factor 2, this fits the form of a quadratic programming problem, which is to minimize

1/2 x^T H x + f^T x,

if we choose H = A^T A and f = - A^T b.

The corresponding use of quadprog is:

H = A' * A;
f = - A' * b;
x = quadprog(H, f)

You also want the elements of x to be positive. A non-negativity constraint can be introduced using the additional parameters to quadprog, A and b (not to be confused with your matrices!):

n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))

A positivity constraint does not make sense, because if the optimal solution involves one or more elements of x to be exactly 0, then a strictly positive solution will be the better the smaller the corresponding elements are: 0.01 will be better than 0.1, 0.001 will be better than 0.01, etc. etc. – there is no natural bound. If you want to make sure that the solution is all-positive, you have to set a finite bound yourself:

x = quadprog(H, f, -eye(n), zeros(n, 1) + 0.001)

Now the smallest possible value of an element of x is 0.001.

Update after the question was supplemented with the actual data of A and b: Using the code

H = A' * A;
f = - A' * b;
n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))

I get the result:

Minimum found that satisfies the constraints.

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

<stopping criteria details>

x =
      0.000380906335150292
      3.90638261088393e-05
        0.0111196970167585
        0.0227055107206744
        0.0318402514628274
        0.0371743514880516
      0.000800900221354844
       0.00746652476710186
        0.0180511534370576
        0.0282423767946842
        0.0362606972021829
        0.0417582260990786
       0.00860220929402253
        0.0174105435824309
        0.0265771677458008
        0.0343071472371469
        0.0395176470725881
        0.0419494410289298
        0.0187719294637544
        0.0268976053211278
        0.0336818044612046
        0.0382365751296441
        0.0398823076542831
        0.0391016682549663
        0.0279383031707377
        0.0339393563379992
        0.0377917413001034
        0.0382731422972829
        0.0338557405807941
        0.0217568643500703
        0.0343698083354502
        0.0381554349806972
        0.0392353941260779
        0.0368010570888738
         0.031271868401718
        0.0258232230013864

Warning: Trust-region-reflective algorithm does not solve this type of problem, using active-set algorithm. For more help, see Choosing the Algorithm in the documentation. — nasim, May 20 '15 at 21:03
@nasim, I checked this code with some randomly generated data, so I cannot check this warning myself. As I commented on the question, better if you post your matrices. — A. Donda, May 20 '15 at 21:05
I will be happy to, but the A matrix does not fit in the comments. Where can I post ? — nasim, May 20 '15 at 21:06
@nasim, click the "edit" link below your question, and copy the data into the question. — A. Donda, May 20 '15 at 21:07
@nasim, thanks for the fixed data. But with these values I get a proper solution without warning. The solution is all-positive even without specifying a finite bound; the smallest x is 3.9e-5. — A. Donda, May 20 '15 at 21:23
Really?! Can I see the 36 values please? And can I please have the exact piece of code you used as an answer? — nasim, May 20 '15 at 21:25
Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/78358/discussion-between-nasim-and-a-donda). — nasim, May 20 '15 at 21:27

score 1 · Answer 2 · answered May 20 '15 at 18:42

This would be a linear optimization problem when you have x>0 condition. The best algorithm to solve that is simplex algorithm. The idea is that each linear equation aix=bi provide a line and the combination of these lines provide a polygon/polyhedron and the answer is one of the vertices of this polyhedron/polygon. Simplex algorithm is pretty standard and there are many available functions and libraries tha can calculate that.

BayesianMonk · Answer 3 · 2018-08-30T15:13:07.070

0

The Matlab functions lsqlin or lsqnonneg can be used to solve your issue. e.g:

x=lsqnonneg(A,b)

will give you what you are looking for.

edited Aug 30 '18 at 15:13

answered Aug 29 '18 at 14:29

BayesianMonk

540
7
23

least-squares method with a constraint

3 Answers3