hessian using symbolic variables

Hi there,

I am struggling with obtaining the Hessian of my function. In the first case it works fine but in the second case I can not get it (still symbolic variables X are presented in H). So my question is why this happen? This is small portion of the code where I want to find hessian of my Rtot function with X being x1,x2,x3,x4.

    Rtot= -sum(Rsub(:,i));
    H=hessian(Rtot,X);
    H1 = 0.5.*double(H);

with Rtot:

(6234858085755363*x1*((17865509505801117*x1)/14411518807585587200000 + (17865509505801117*x3)/14411518807585587200000 - 16863388614392257743519/1441151880758558720000))/140737488355328000 + (6069173180546921*x2*((17865509505801117*x1)/7205759403792793600000 + (17865509505801117*x2)/14411518807585587200000 + (17865509505801117*x3)/7205759403792793600000 + (17865509505801117*x4)/14411518807585587200000 - 16867211833426499182557/1441151880758558720000))/140737488355328000

and H1:

1.0e-07 *

      0.5492    0.5346    0.2746         0
      0.5346    0.5346    0.5346    0.2673
      0.2746    0.5346         0         0
           0    0.2673         0         0

This is fine. But then when I want to find hessian of different data set i get this. Rtot:

- (2587079436111693477*x2*(665101922855875194446501149/3153075739026468800000 - (4398046511104*x2)/5912017010674629 - (4398046511104*x3)/5912017010674629 - (4398046511104*x4)/5912017010674629 - (4398046511104*x1)/5912017010674629)^(1/3))/28147497671065600000000 - (4554081285294461427*x1*(1299027180826527113674752/6158351052786071875 - (4398046511104*x3)/5912017010674629 - (4398046511104*x1)/5912017010674629)^(1/3))/56294995342131200000000

and H still have symbolic variables inside which does not allow me to calculate H1. Difference is in Rtot but I still can not find out why it does not work.

Thank you a lot!!!

[Torsten]

In the first case, your function is quadratic in the unknowns. Differentiation twice gives a constant Hessian that does not depend on x1, x2, x3 and x4.

In the second case, you have to substitute values for x1,x2,x3 and x4 to be able to get H1. Use MATLAB's "subs": (H1=double(subs(...))).

Best wishes

Torsten.

[sensation]

Hi Torsten, Thanks for your insight. The code that comes before hessian is:

 %
i=1;
Rsub(:,i) = subs(R(:,i),[c;s0(1,i);P(:,1);F(:,i);ro;eff;g;alpha(i);beta(i)],[C';init_storage(i);price1(1:N,1);inFlow1(1:N,i);ro1;eff1;gg;alpha_par1(i);b_par1(i)]);
Rtot= -sum(Rsub(:,i));
fsym= gradient(Rtot,X);
f= double(subs(fsym,X,zeros(size(X))));
H = hessian(Rtot,X);
H1=0.5.*double(H);

The substitute values youre talking about should be in H = hessian(Rtot,X) right? The X is unknown and the next step would be to do quadprog optimization knowing H and f.

However, if I try to do H1=double(subs(H,X)) I get ofc

"DOUBLE cannot convert the input expression into a double array" since x are unknown. Thanks for any tip!

M

[Torsten]

if you want to use the Hessian in quadprog, use

func_H = matlabFunction(H)

Best wishes

Torsten.

[sensation]

Hi Torsten, Thanks for the tip. I created my Hessian as:

%
     hessObj = hessian(Rtot(:,i),X(:,i)); 
     H=matlabFunction(hessObj);
However, when I use H in my quadprog, H is still not double type. Any clue?

Thanks and cheers

Answer 1

[Alan Weiss]

There are some examples in the documentation that show how to do this:

Using Symbolic Math with Optimization Toolbox Solvers generates function handles.

Symbolic Math Toolbox Calculates Gradients and Hessians generates function files.

Alan Weiss

MATLAB mathematical toolbox documentation

[sensation]

Hi Alan, I had a view on this. My issue is following once I do:

% 
gradf = jacobian(Rtot,X); % column gradf
hessf = jacobian(gradf,X);
qpoptions = optimset('Algorithm','interior-point-convex','Disp','iter');
tic
[x3(:,i),fval3(:,i)] = quadprog(hessf,f,A,b(:,i),Aeq,beq(:,i),LB(:,i),UB(:,i),[],qpoptions);
toc

hessf need to be of type H (double) but I can not convert it to double since there are still x unknowns in hessf. My hessf, instead of being a number able to be converted in double is:

(154325299586858956762382336*x2)/(13653103568166475135677228237359765625*(665101922855875194446501149/3153075739026468800000 - (4398046511104*x2)/5912017010674629 - (4398046511104*x3)/5912017010674629 - (4398046511104*x4)/5912017010674629 - (4398046511104*x1)/5912017010674629)^(5/3)) + (15092307240482528980631552*x1)/(1517011507574052792853025359706640625*(1299027180826527113674752/6158351052786071875 - (4398046511104*x3)/5912017010674629 - etc

thanks for any tip!

[Alan Weiss]

Torsten gave you good advice: use matlabFunction exactly as he said.

Alan Weiss

MATLAB mathematical toolbox documentation