different results from fmincon with nonlcon

6 views (last 30 days)
Hi, my problem is that an optimization with fmincon (see below) is not stable, i.e. various runs deliver various results. The optimization problem is described as below:
Many thanks in advance for your help!
Best Wishes, Alex
The function to be optimized:
fun=@(x)x*meanOrRegressed
The non-linear constraint (Sum of all positive weights not above 1.3 and sum of all negative weights not below -0.3:
function [c,ceq,gradc,gradceq]=nonlcon_g(x,1.3,0.3)
nrAssets=size(x,1);
weightSumMax = sum(max(0,x(1:nrAssets-1)));
c(1) = weightSumMax-1.3;
weightSumMin = sum(min(0, x(1:nrAssets-1)));
c(2) = -weightSumMin-0.3;
gradc1=1/2*(sign(x(1:nrAssets-1))+ones(nrAssets-1,1));
gradc1=vertcat(gradc1,0);
gradc2=1/2*(sign(x(1:nrAssets-1))-ones(nrAssets-1,1));
gradc2=vertcat(gradc2,0);
gradc=horzcat(gradc1,gradc2);
ceq =[];
gradceq=[];
end
[weights,optimizedValue,exitflag] = fmincon(@fun,initialWeights',A,b,Aeq,beq,lowerBound,upperBound,@nonlcon_g,options)
A = -1 -1 -1 -1 -1 0
1 1 1 1 1 0
0 0 0 0 0 -1
0 0 0 0 0 1
b = 0.3 1.3 0.3 1.3
Aeq = 1 1 1 1 1 1
beq = 1
lowerBound = -0.3 -0.3 -0.3 -0.3 -0.3 -0.3
upperBound = 1.3 1.3 1.3 1.3 1.3 1.3 1 3
meanOrRegressed =
2.349096891796729e-004
-7.582259013820250e-005
1.190461785891006e-003
2.529756213317396e-003
1.066862350689632e-003
5.133561643835617e-005
options =
Display: []
MaxFunEvals: 60000
MaxIter: 40000
TolFun: 1.000000000000000e-010
TolX: []
FunValCheck: []
OutputFcn: []
PlotFcns: []
ActiveConstrTol: []
Algorithm: 'active-set'
AlwaysHonorConstraints: []
BranchStrategy: []
DerivativeCheck: []
Diagnostics: []
DiffMaxChange: []
DiffMinChange: []
FinDiffRelStep: []
FinDiffType: 'forward'
GoalsExactAchieve: []
GradConstr: 'on'
GradObj: []
HessFcn: []
Hessian: []
HessMult: []
HessPattern: []
HessUpdate: []
InitialHessType: []
InitialHessMatrix: []
InitBarrierParam: []
InitTrustRegionRadius: []
Jacobian: []
JacobMult: []
JacobPattern: []
LargeScale: []
LineSearchType: []
MaxNodes: []
MaxPCGIter: []
MaxProjCGIter: []
MaxRLPIter: []
MaxSQPIter: 200000
MaxTime: []
MeritFunction: []
MinAbsMax: []
NodeDisplayInterval: []
NodeSearchStrategy: []
NoStopIfFlatInfeas: []
ObjectiveLimit: []
PhaseOneTotalScaling: []
Preconditioner: []
PrecondBandWidth: []
RelLineSrchBnd: []
RelLineSrchBndDuration: []
ScaleProblem: []
Simplex: []
SubproblemAlgorithm: []
TolCon: 1.000000000000000e-008
TolConSQP: []
TolGradCon: []
TolPCG: []
TolProjCG: []
TolProjCGAbs: []
TolRLPFun: []
TolXInteger: []
TypicalX: []
UseParallel: 'always'

Accepted Answer

Matt J
Matt J on 14 Apr 2014
Edited: Matt J on 14 Apr 2014
Another solution, one which avoids the need for epsilon-approximation, is to recognize that your original "nonlinear" constraints are really linear ones. In particular, one can show that
sum(max(x,0))<=K
is equivalent to the set of linear inequalities
sum(x(J))<=K
for every index subset, J. This can be expressed by adding more rows to your A,b data. I do so as follows and get a highly accurate result.
load AlexData
N=length(initialWeights);
AA=dec2bin(0:2^N-1,N)-'0';
AA(sum(AA,2)<=1,:)=[];
M=size(AA,1);
bb(1:M,1)=maxLong;
bb(M+1:2*M)=maxShort;
A=[A;AA;-AA]; %Append new constraints
b=[b(:);bb];
lowerBound=max(lowerBound,-maxShort);
upperBound=min(upperBound, maxLong);
options.algorithm='interior-point';
options.TolX=1e-20;
options.TolFun=1e-20;
options.UseParallel=[];
[weights,optimizedValue,exitflag,output,lambda,grad,hessian] = ...
fmincon(@(x)dot(-x,ret), initialWeights', A, b, Aeq, beq,...
lowerBound, upperBound,[],options);
weightError=weights(:).' - [0, 1.3, -.3 0 0 0],
Unfortunately, this means that your constraint data A,b consume O(N*2^N) memory where N is the number of unknown weights. It should be fine, though, for small problems (N<=15 or so).
  5 Comments
alexander
alexander on 22 Apr 2014
Dear Matt, thanks again for all your help. The least I can do is to click this accept answer button, but I'm looking the third time and still cannot find it! Hope you can help! Thanks & Kind Regards, Alex
Matt J
Matt J on 23 Apr 2014
Maybe because you're not logged in to your mathworks account?

Sign in to comment.

More Answers (2)

Matt J
Matt J on 28 Mar 2014
The nonlinear constraints are not differentiable, so that could well be the cause of the instability.
  14 Comments
Alexander
Alexander on 8 Apr 2014
Hi Matt, is there anything else which might be useful for you to tackle the problem? I'd be glad to be of any assistance. Thanks & Kind Regards, Alex
Matt J
Matt J on 11 Apr 2014
Edited: Matt J on 11 Apr 2014
Hi Alex,
I've run your code, and also made some modifications of my own, which I've attached.
I do not see the instability that you describe when I run your version. I ran 10 consecutive trials and verified that the answers were identical across all of them. However, I vaguely wonder if you might get slightly different results because of UseParallel being turned on. Parallel processing routines can cause things to be executed in a somewhat random order. In theory the order of parallel tasks shouldn't matter, but in practice, for parallel reduction operations , there can be slight order-dependencies in the results due to floating point precision issues. That's just a guess, though -- I'm not even sure there are any reduction operations when fmincon runs in parallel mode. In any case, I don't really see that UseParallel would help you in terms of speed for such a small computation as this. I think it might even hurt you.
I do see limitations on the accuracy of the result. However, because you are not solving the original problem, but rather an epsilon-approximation of it, you can only expect to get approximate results. There is also an influence on accuracy of doing finite difference derivative computations.
In my version, I turned on the 'GradConstr' option and also switched back to the 'active-set' method. I also introduced a refinement step. It basically takes the solution of the epsilon-problem and uses that to decide which weights should be positive and which weights should be negative. It then constrains the search to that specific pattern of positives/negatives (similar to what I suggested earlier that you do combinatorically). This allows you to get rid of the nonlinear constraints, whereupon the solution becomes highly well-conditioned and accurate.
Finally, I corrected a bug in your implementations of max_smooth/min_smooth. You set nrAssets=length(x)-1 for some reason, when it really needs to be nrAssets=length(x).

Sign in to comment.


Alexander
Alexander on 23 Apr 2014
Hi Matt, it finally worked!! Thanks again, and all the best! Alex

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!