Fsolve with constraint for positivity
2 views (last 30 days)
Show older comments
Miroslav Mitev
on 15 May 2018
Commented: John D'Errico
on 16 May 2018
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 Comment
John D'Errico
on 15 May 2018
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.
Accepted Answer
John D'Errico
on 15 May 2018
Edited: John D'Errico
on 15 May 2018
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.
5 Comments
Torsten
on 16 May 2018
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
on 16 May 2018
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.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!