What are suggestions if I need to write a function to calculate Hessian Matrix, and there is non-polynomina part in the constraints?

[Xinwei Shen]

Can any other function in Yalmip help me to do that?

I knew hessian(Obj, x) can help a little bit. But, 

• What are the orders of all the variables that I defined? So that I can put it in the x in the right order. e.g. x=[x1; x2; ...; xn].
• How to establish the Lagrangian?

[Johan Löfberg]

Not clear what you are asking about. YALMIP does not supply Hessians to the solver (only gradients) and you cannot change this fact despite knowing how to compute the hessian

[Xinwei Shen]

Dear Prof. 

I built up the model with Yalmip, but I was trying to use fmincon and interior-point method with newton-raphson, so I need a complete function to compute Hessian matrix of the Lagragian.

See

Hessian for fmincon interior-point algorithm 

at

.http://cn.mathworks.com/help/optim/ug/writing-scalar-objective-functions.html#bu2w6a9-1

My questions are:

• What are the orders of all the variables that I defined (with Yalmip) ? So that I can put it in the x in the right order. e.g. x=[x1; x2; ...; xn].
• How to establish the Lagrangian (effciently with functions in Yalmip)?

[Johan Löfberg]

If you are calling fmincon via YALMIP, YALMIP automatically supplies fmincon with function values and gradients. The Hessian is not computed and fmincon thus uses quasi-newton methods.

 

If you want to use fmincon and supply hessians, you will have to work completely outside of YALMIP.

[Xinwei Shen]

OK I see. 

Thank you so much

[Xinwei Shen]

btw, does IPOPT support solving problem by interio point method with newton-raphson and hessian?

[Johan Löfberg]

Same thing there. As YALMIP doesn't compute hessians, ipopt will estimate them using BFGS and similiar

[Johan Löfberg]

misread your question. Yes, ipopt iimplements a (quasi-)newton-based interior-point algorithm

[Xinwei Shen]

Do you mean that IPOPT don't compute exact hessian?

[Xinwei Shen]

But sir, according to

https://www.coin-or.org/Ipopt/documentation/node53.html#SECTION0001113010000000000000

Quasi-Newton

hessian_approximation:

Indicates what Hessian information is to be used. 
This determines which kind of information for the Hessian of the Lagrangian function is used by the algorithm. The default value for this string option is "exact". 
Possible values:

• exact: Use second derivatives provided by the NLP.
• limited-memory: Perform a limited-memory quasi-Newton approximation

I think choosing "exact" here means that using exact hessian, is that true?

[Johan Löfberg]

solvers never compute hessians (or gradients or function values). These are supplied from outside (unless you are writing a special solver for a particular problem)

IPOPT will use the Hessian if it is supplied. If not, it will estimate it using quasi-newton methods. Default approximation method appears to be limited-memory BFGS

>> ops = sdpsettings;

>> ops.ipopt.hessian_approximation

ans =

limited-memory

[Johan Löfberg]

It USES the exact hessian in that case, it does not compute it.It is provided by the code that interfaces IPOPT.

[Xinwei Shen]

Hi Prof. 

I was wondering if you can tell me why do I have to work outside Yalmip to include a function to calculate Hessian?

I see the options 

ops=sdpsettings;
ops.fmincon

...
                   HessFcn: []
                   Hessian: []
                  HessMult: []

I think it is possible.

The only problem is that how to correspond those variables and constraits with "x" and "lambda" in hessian = hessianfcn(x,lambda) at http://cn.mathworks.com/help/optim/ug/writing-scalar-objective-functions.html#bu2xbye-1

[Johan Löfberg]

If you use YALMIP to model your problem, and then solve ipopt via yalmip to solve that problem, you will never get that hessian to be sent to ipopt. YALMIP automatically computes function values and gradients directly from your model completely automatically, and that is all the information that is sent to the solvers.

OK, perhaps you intend to do some mix where you create a separate function to compute the hessian completely outside yalmip, but then I have to question why you are using YALMIP at alll? In that case, you must have a pretty simple model (simple as in straightforward indexing  and no high-level operators etc) Why not create functions for the function values and gradients also, and get rid of the head-ache of trying to glue together yalmip stuff (where a lot of stuff can happen under the hood) and manual stuff

The variables in the solver are typically sorted according to the order that they where created. Note though that some operators might introduce additional variables introduced internally to normalize expressions and simplify gradient computations etc,  that you as a user never see or use. For instance, the model 

sdpvar x
optimize([],exp(1+x))

will have two variables and a model with an equality constraint. The first one will correspond to x, and the second will correspond to a variable y which is constrained to be equal to 1+x

[Xinwei Shen]

As what you've mentioned, I am just trying to avoid the complexity of indexing the variables and modeling the constraints, so I'd like to remaining the model with Yalmip and try to supply Hessian with another function that I defined.

Now that you've say that, could you please offer me some suggestions about how to index the variables in X, and perhaps, about producing the A, b, Aeq, beq (AX<=b, AeqX=beq)?

Maybe I can produce these matrixs and vectors by export() to other solvers? If I can make sure that the orders of variables are corresponding to X.

[Johan Löfberg]

sdpvar x y z w
[model,internal] = export([x + y + w] == 1,x^2)

internal.used_variables
  used_variables: [1 2 4]

getvariables(w)

ans =

     4

[Xinwei Shen]

after half-day labor work, I succeed in using export() to export my YALMIP model to CPLEX and get the Aineq, bineq, Aeq, beq (AineqX<=bineq, AeqX=beq), even though I only use these matrix and vectors in fmincon and solve the problem.

However,

(1) you must move all the nonlinear constraints into a self-defined function, because CPLEX cannot handle them. 

(2) it's necessary to check if all the variables are used in remaining linear constraints. If not,  Aineq, bineq, Aeq, beq could be in a wrong size since not all variables are used.

Above codes are helpful for indexing all the variables.

[Xinwei Shen]

Thank you so much.

[Johan Löfberg]

to ensure a complete map, simply add a redundant linear constraint which involves all your variables