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: xi > 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]';
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||2 = xT AT A x + bT b - 2 xT AT b.
Disregarding the constant term bTb and dividing by a factor 2, this fits the form of a quadratic programming problem, which is to minimize
1/2 xT H x + fT x,
if we choose H = AT A and f = - AT 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