Linear problem, linprog does not give realistic answer

I need to solve a simple Ax=b problem, where A=1x5 and b=1x1. All elements in A and b are positive and so should the elements in x be. Hence x has a lower and upper bound. I use Matlab to solve this.

When I solve this with A\b or linprog, it gives me an answer which is not realistic in the context x is in(x=[0,0,0,x4,0]^T where x4 is some number). A more realistic answer, coming from what these numbers represent, would be that all elements of x are unequal to zero and are of roughly the same order. I don't know exactly what this order is and putting in artificial inequalities for the elements in x is not an option.

I noticed that when I use fmincon for this problem with the active set algorithm, gives me what I am looking for. Can someone explain me why this is, and how the active set algorithm works? Thanks in advance.

Answer 1

[Matt J]

You should really be using lsqlin

x = lsqlin(A,b,[],[],[],[],zeros(size(b)));

A\b does not know about your upper and lower bounds on x. The fact that A and b are positive does not imply that x will be.

It's not clear why you think linprog would solve A*x=b.

I'm assuming fmincon works because you've implemented the same least squares solution that lsqlin does. But, it is overkill. lsqlin should be sufficient.

[Matt J]

You should also realize that your problem is ill-posed. Consider

 A=[1 1 1 1 1]
 b=5;

There are infinite positive x that will solve this problem. A subset of such solutions is x=[c, 2-c, 1,1,1] for any 0<=c<=2.

[MQMST66]

Thank you for your answer. I know there are infinitely many solutions, but I want to know why the fmincon with active-set algorithm gives this specific one. I also implemented the lsqlin algorithm, however this gives again another solution which is not realistic.

Do you happen to know how this specific algorithm in fmincon works?

[Matt J]

but I want to know why the fmincon with active-set algorithm gives this specific one.

There's nothing special about the active-set algorithm that you're benefiting from here.

You're getting a solution you like better as a result of luck and the initial guess x0 that you supplied. If you take the result that lsqlin gives you and supply that to fmincon as the initial x0, then fmincon will give the same result as lsqlin.

Answer 2

[John D'Errico]

What matters is NOT the algorithm in fmincon, since the solution is strongly defined by your starting values. As well, you have not even said what objective function you set for fmincon, surely an important factor! Were we to spend a great deal of time trying to explain the algorithm, it would not help you in the least.

As well, using linprog to solve the problem is also a poor choice of tool.

If these tools result in something that is unrealistic for your context, then you have failed to bring the additional information that you hold to the problem. How can mathematics (or MATLAB) read your mind?

Really, your problems seems to be:

Solve A*x = b, where A is 1x5, x 5x1, and b a scalar, subject to the goal that the elements of x are as close as possible to each other as possible, AND that all are positive.

Formulating the proper mathematical problem is half way there. Next, how does one solve the problem? Best is to use LSQLIN. It requires no starting values, is not overtly iterative, is fast, and can solve your problem. What else can you ask for?

First, I'll pick some random vector to serve as A, and an arbitrary value for b.

A = rand(1,5)
A =
      0.16261        0.119      0.49836      0.95974      0.34039

b = 1;

We will use A*x = b as an equality constraint.

lb = zeros(1,5);
ub = [];

The design matrix here will NOT be A*x = b. Instead, it will be a mostly zero matrix that will attempt to force all of the elements of A to be as close as possible to each other.

ind = nchoosek(1:5,2);
n = nchoosek(5,2);
C = full(sparse(repmat((1:n)',1,2),ind,repmat([-1 1],n,1),n,5));
d = zeros(n,1);

Spy is a good tool to visualize what C looks like here. I really did not need to make it a sparse matrix, as it is not very sparse. However, if your problem were larger, then it might be useful to have C sparse. Of course, that might be a problem for LSQLIN. (I can't recall if it handles sparse problems like this, and I'm feeling too lazy to look.)

C
C =
    -1     1     0     0     0
    -1     0     1     0     0
    -1     0     0     1     0
    -1     0     0     0     1
     0    -1     1     0     0
     0    -1     0     1     0
     0    -1     0     0     1
     0     0    -1     1     0
     0     0    -1     0     1
     0     0     0    -1     1

And now solve.

x = lsqlin(C,d,[],[],A,b,lb,ub)
Warning: Large-scale algorithm can handle bound constraints only;
    using medium-scale algorithm instead. 
> In lsqlin at 269 
Optimization terminated.
x =
      0.48075
      0.48075
      0.48075
      0.48075
      0.48075

I would point out that the solution has all its values equal, something we could probably have predicted earlier, had we looked at the problem more carefully.

I should also have specified the medium-scale algorithm in the options to avoid the warning.