Fsolve with constraint for positivity

Question

Miroslav Mitev 2018-5-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/400863-fsolve-with-constraint-for-positivity

评论： John D'Errico 2018-5-16

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

John D'Errico 2018-5-15

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 2018-5-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/400863-fsolve-with-constraint-for-positivity#answer_320352

编辑：John D'Errico 2018-5-15

在 MATLAB Online 中打开

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.

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.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2018-5-16

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 2018-5-16

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.

请先登录，再进行评论。