Fsolve with constraint for positivity

I have the following function:

   function y = fun(x)
   b=1:5;
  for i=1:5
  s(i)=sum(1./b(i:end).^2);
  end
   y(1)=x(1)+1/s(1)-x(2)-1/s(2);
   y(2)=x(2)+1/s(2)-x(3)-1/s(3);
   y(3)=x(3)+1/s(3)-x(4)-1/s(4);
   y(4)=x(4)+1/s(4)-x(5)-1/s(5);
   y(5)=sum(x)-5*5;
   end

And then I use fsolve:

x0=[5 5 5 5 5];
z = fsolve(@f, x0);

It works OK, but I want to add constraint for all x to be positive. Is there any possibility to do that? I tried to reformulate and run it through fmincon but error appears.

[John D'Errico]

Why in the name of god and little green apples would you use fsolve to solve a problem that is linear in the parameters? You have a linear system of 5 equations in 5 unknowns.

Answer 1

[John D'Errico]

You have a system of 5 linear equations, in 5 unknowns. As long as the system is non-singular, there is only ONE solution. Using fmincon won't give a better answer.

Write s using MATLAB...

s = flip(cumsum(1./flip(b.^2)));

Write the linear system of equations. So the right hand side is:

rhs = [diff(1./s)';5*5];
A = [1 -1 0 0 0;0 1 -1 0 0;0 0 1 -1 0;0 0 0 1 -1;1 1 1 1 1];

A
A =
   1    -1     0     0     0
   0     1    -1     0     0
   0     0     1    -1     0
   0     0     0     1    -1
   1     1     1     1     1

rank(A)
ans =
     5

b = [diff(1./s)';5*5]

b =
     1.4737
     2.5244
     5.0747
     15.244
         25

So, we need to solve for the 5x1 vector x, such that A*x=b. This is done in MATLAB using backslash most easily.

A\b
ans =
     12.772
     11.299
     8.7741
     3.6994
    -11.544

Now the problem becomes obvious. When you used fsolve, you got a negative solution in x. But the problem is this solution is UNIQUE. In fact, the matrix A is very well conditioned, not at all singular. So there is no alternative solution that exists, or can exist where all of the elements of x are nonnegative, at least none that solves the problem exactly.

At the same time, you CAN find the solution that best solves the problem, where all 5 elements of x are non-negative. It won't solve the problem exactly.

Xnonneg = lsqnonneg(A,b)
Xnonneg =
     5.6967
     6.0849
     7.2845
     7.7959
          0

LSQNONNEG is actually part of MATLAB itself. You don't need any special toolbox. How good is that solution?

[A*Xnonneg,b]
ans =
   -0.38827       1.4737
    -1.1996       2.5244
   -0.51133       5.0747
     7.7959       15.244
     26.862           25

It is not that great. You can clearly see that it does not come that close to

But it is the best possible solution that has all elements of x non-negative.

As a followup. I just noticed that you had written:

   y(1)=x(1)+1/s(1)-x(2)-1/s(2);
   y(2)=x(2)+1/s(2)-x(3)-1/s(3);
   y(3)=x(3)+1/s(3)-x(4)-1/s(5);  **********
   y(4)=x(4)+1/s(4)-x(5)-1/s(5);

I see that the 3rd equation seems to be atypical. Did you intend to subtract 1/s(5) there? Or did you really mean to write that line as

y(3)=x(3)+1/s(3)-x(4)-1/s(4);

which would seem more consistent. As I solved it above, I assumed you meant the latter form, and what you wrote in this post was a typo. It looks like a copy/paste error.

But if you REALLY wanted the solution to what you wrote, and it was NOT a typo, then we would have b as:

b = [1./s(2 3 5 5])' - 1./s(1 2 3 4])';5*5]

b =
       1.4737
       2.5244
       20.319
       15.244
           25
>> A\b
ans =
        18.87
       17.396
       14.872
      -5.4469
      -20.691
>> lsqnonneg(A,b)
ans =
       6.1884
       8.0519
       12.202
        1.895
            0

My guess is you wrote it as a typo, so the earlier solution is what you want to see.

[Miroslav Mitev]

Thank you for this answer. Yes, it was a typo. I was just wondering, how can I do the following: If lets say x(5) is negative replace it with 0 and solve the system again where sum(x) is still 25. Or if x(5) and x(4) are negative replace both with 0 but, we still have x(1)+x(2)+x(3)+0+0=25 and so on.

[Torsten]

And you don't care about the resulting errors in the first four equations after you replaced the negative x(i)'s with zeros ?

[Miroslav Mitev]

I care. Lets say if x(5) is negative we go 1 step back and solve the system for x(1), x(2), x(3), x(4) only knowing that x(5)=0, if we have negative again we go back until all we reach the case where all of them are positive. Which could be just x(1)=25 and x(2)=x(3)=x(4)=x(5)=0.

[Torsten]

You don't care. Setting x5=0 leaves you with a system of 5 equations in 4 unknowns. This means that an exact solution for x1,...,x4 that satisfies all 5 equations simultaneously will in general not exist. Thus each equation will have an error, and the sum of the 5 errors squared will in general be larger as if you had solved for all 5 variables under the constraint that all of them are nonnegative.

I think what you try to do is to minmize the sum of errors squared of the first four equations under the constraint that x1,...,x5 are nonnegative and x1+x2+x3+x4+x5=25. Use "lsqlin" to do so.

Best wishes

Torsten.

[John D'Errico]

Again, there is a UNIQUE solution. So if you choose not to accept that unique solution (presumably because one of the elements was negative) then the result will not be exact.

If A\b returns a result with one or more negative, then there is no solution that has all of them positive. This is true as long as rank(A) is 5.

The scheme you suggest is sort of what lsqnonneg uses. However it may be that the optimal non-negative solution need not always result from just setting a negative element to zero and then resolving. And that is why you use lsqnonneg. lsqlin is fine as an alternative, as long as you do have the optimization toolbox.

As Torsten points out, you can use lsqlin, if you really want to absolutely insure that the sum is 25, while allowing the other 4 equations to be solved approximately. That would have you pose the sum of the 5 variables as an equality constraint.