Find gradient of objective function for fmincon

Question

Arnab Samaddar-Chaudhuri 2023-1-11

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1892190-find-gradient-of-objective-function-for-fmincon

评论： Arnab Samaddar-Chaudhuri 2023-1-12

Hello,

I need to optimize 2 points. I am finding it difficult to manually calculate the gradient of the objective function The shortened version of the code is as follows:

 W1 = @(X) objfun(X,XQ,YQ,cdf1,cdf2);
 % I have a meshgrid, where XQ, YQ, cdf1, cdf2 are of mxn matrix and are known
 
 options = optimoptions('fmincon','Algorithm','sqp', 'Display','iter-detailed',...
                        'CheckGradients',true,'SpecifyObjectiveGradient',true);
 lb = [1909,-390;
       2059,-390];
 ub = [1969,-360;
       2059,-360];
 initValue = [1969,-360;
              2059,-360];
 %these values are not important here
 
[x1,~,exitflag] = fmincon(W1,initValue,[],[],[],[], lb1, ub1, [], options);
%where
function [f,gradf] = objfun(X,XQ,YQ,cdf1,cdf2)
p1 = @(Point) interp2(XQ,YQ,cdf1,Point(1,1),Point(1,2),'nearest');
p2 = @(Point) interp2(XQ,YQ,cdf2,Point(2,1),Point(2,2),'nearest');
f =  -1+(1-p1(X))*(1-p2(X));
if nargout > 1
    [gx1, gy1] = gradient(cdf1);
    [gx2, gy2] = gradient(cdf2);
    
    dp1x = @(Point) interp2(XQ,YQ,gx1,Point(1,1),Point(1,2),'nearest');
    dp1y = @(Point) interp2(XQ,YQ,gy1,Point(1,1),Point(1,2),'nearest');            
    dp2x = @(Point) interp2(XQ,YQ,gx2,Point(2,1),Point(2,2),'nearest');
    dp2y = @(Point) interp2(XQ,YQ,gy2,Point(2,1),Point(2,2),'nearest');
    
    gradfx1 = @(P) - dp1x * (1-p2(P)) ;
    gradfy1 = @(P) - dp1y * (1-p2(P)) ;
    gradfx2 = @(P) - (1-p1(P)) * dp2x ;
    gradfy2 = @(P) - (1-p1(P)) * dp2y ;
    gradfx = [gradfx1 ; gradfx2];
    gradfy = [gradfy1 ; gradfy2];
     
    gradf = [gradfx;gradfy];
end
end

I get the following error:

Error using vertcat
Nonscalar arrays of function handles are not allowed; use cell arrays instead.
Error in SBR_Optimisation_tool_1_8_4/objfun (line 2893)
                gradfx = [gradfx1 ; gradfx2];

I am sure I am not calculating the gradient correctly. I simpy cannot use gradient(f) or gradient(p1). That gives other kinds of error. How to calculate the gradient here?

4 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2023-1-11

编辑：Torsten 2023-1-11

在 MATLAB Online 中打开

p1 = @(x,y) interp2(XQ,YQ,cdf1,x,y,'nearest');
p2 = @(x,y) interp2(XQ,YQ,cdf2,x,y,'nearest');

These two functions are not even continuous functions of x and y. So how could you calculate a gradient for them ?

Further, you must supply numerical values for the gradient at x, not a function handle as you do.

Be happy if fmincon works with your non-differentiable function and don't try to make things even worse by supplying non-existing gradients.

Walter Roberson 2023-1-11

在 MATLAB Online 中打开

To get around that particular problem change

    gradfx = [gradfx1 ; gradfx2];
    gradfy = [gradfy1 ; gradfy2];
    gradf = [gradfx;gradfy];

to

    gradfx = {gradfx1 ; gradfx2};
    gradfy = {gradfy1 ; gradfy2};
    gradf = [gradfx;gradfy];

The result will be a 4 x 1 cell array of function handles, which is legal in MATLAB.

MATLAB does not permit you to use [] to put more than one function handle into an array, because if it did, gradfx(1) would be ambiguous about whether you are invoking the array of function handles each on the input [1], or if you were invoking each of the function handles with no input and then indexing the result at the first location, or if you were selecting the first function handle to be returned as a function handle, or if you were selecting the first function handle and invoking it with no input. It therefore forces you to store multiple function handles in cell, as you can then use {} and () to control your meaning.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2023-1-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1892190-find-gradient-of-objective-function-for-fmincon#answer_1146670

编辑：Matt J 2023-1-11

在 MATLAB Online 中打开

Perhaps as follows. Note my restructuring of your array into an Nx2x2 form.

function [f,gradf] = objfun(X,XQ,YQ,cdf1,cdf2)
[x1,y1,x2,y2]=deal(X(:,1,1), X(:,2,1), X(:,1,2), X(:,2,2));
p1 = interp2(XQ,YQ,cdf1,x1,y1,'cubic');
p2 = interp2(XQ,YQ,cdf2,x2,y2,'cubic');
f =  -1+(1-p1)*(1-p2);
 if nargout > 1
      D=1e7*eps(X);
     [dx1,dy1,dx2,dy2]=deal(D(:,1,1), D(:,2,1), D(:,1,2), D(:,2,2));
      
      p1x1=(interp2(XQ,YQ,cdf1,x1+dx1,y1,'cubic') -p1)./dx1;
      
      p1y1=(interp2(XQ,YQ,cdf1,dx1,y1+dy1,'cubic') -p1)./dy1; 
      
      p2x2=(interp2(XQ,YQ,cdf2,x2+dx2,y2,'cubic') -p2)./dx2;
      
      p2y2=(interp2(XQ,YQ,cdf2,dx2,y2+dy2,'cubic') -p2)./dy2;   
      
      gradp1=[p1x1;
              p1y1].*(1-p2);
      
      gradp2=[p2x2;
              p2y2].*(1-p1);
      
      gradf=[gradp1;gradp2];
      
 end
end

2 个评论
显示无隐藏无

Arnab Samaddar-Chaudhuri 2023-1-12

在 MATLAB Online 中打开

Hello Matt,

I have expanded the number of points to 4 to test this code. I get following error both for 2 points and 4 points:

Index in position 2 exceeds array bounds. Index must not exceed 2.
Error in SBR_Optimisation_tool_1_8_4/objfun (line 2900)
    [x1,y1,x2,y2,x3,y3,x4,y4]=deal(X(:,1,1), X(:,2,1), X(:,3,1), X(:,4,1), X(:,1,2), X(:,2,2), X(:,3,2), X(:,4,2));

How did you do the restructuring of the array ?

Also shouldnt it be like this? :

[x1,x2,x3,x4,y1,y2,y3,y4]=deal(X(:,1,1), X(:,2,1), X(:,3,1), X(:,4,1), X(:,1,2), X(:,2,2), X(:,3,2), X(:,4,2));    

I guess it depends how the restructuring is done.

The points are just (x,y) coordinates on 2D space.

BR

Arnab

Arnab Samaddar-Chaudhuri 2023-1-12

在 MATLAB Online 中打开

I modified your code, Matt to this:

[x1,x2,y1,y2]=deal(X(1,1), X(2,1), X(1,2), X(2,2));
[dx1,dx2,dy1,dy2]=deal(D(1,1), D(2,1), D(1,2), D(2,2));
p1y1=(interp2(XQ,YQ,cdf1,x1,y1+dy1,'cubic') -p1)./dy1; % x1,y1+dy1 instead of dx1,y+dy1
% i repeat the above step for p2y2,p3y3, p4y4
gradp1=[p1x1;
        p1y1].*-(1-p2)*(1-p3)*(1-p4); % added minus symbol to all gradients
gradf=[gradp1(1);gradp2(1);gradp3(1);gradp4(1);...
       gradp1(2);gradp2(2);gradp3(2);gradp4(2)]; % changed this

Then the gradientCheck runs perfectly.

Thanks a lot for your help

Best regards

Arnab

请先登录，再进行评论。

Answer 2

Walter Roberson 2023-1-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1892190-find-gradient-of-objective-function-for-fmincon#answer_1146600

Using interp2 with 'nearest' is incompatible with fmincon, fminunc, fminsearch, and most of the other optimizers. In particular it is incompatible with any of the gradient based optimizers, all of which require the the first and second derivatives of the equations are continuous. If you need 'nearest' in the objective function then you must switch to ga() or one of a small number of other optimizers. As none of these other optimizers use gradients then the method of calculation of gradient for them becomes irrelevant.

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

Arnab Samaddar-Chaudhuri 2023-1-11

在 MATLAB Online 中打开

Hello Walter,

Thank you for the answer in one of the previous comments. I normally use interp2 with 'cubic' for fmincon optimization. However since I was more concerned with speed and not with high accuracy, I switched to 'nearest'. At this moment I donot wish to try the optimization with ga() due to time constraints of carrying out new investigative work. I will switch to cubic and try with above these changes suggested by you

gradfx = {gradfx1 ; gradfx2};
gradfy = {gradfy1 ; gradfy2};

请先登录，再进行评论。

Find gradient of objective function for fmincon

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Find gradient of objective function for fmincon

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无

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