Generating Hessian of the Lagrangian with dynamic number of nonlinear constraints in Fmincon

Hi,

(I posted this in stack exchange 2 days back but didn't get a response. Hope it works here).

I'm using interior point algorithm for solving a nonlinear optimization problem and want to provide Hessian of the Lagrangian as part of fmincon to speed up the process (running couple of thousand variations for different sets of parameters)

The workflow that I'm designing has the (nonlinear) objective function, form of linear and nonlinear constraints each in separate function files. I use a master script which has the fmincon. Since i'm running different variations of the problems, sometimes there may be zero nonlinear equality & zero nonlinear inequality constraints and in other scenarios, nonzero number of nonlinear equality/inequality constraints. In short the number of nonlinear constraints are dynamic and stored in the 2 variables "NumOfNonLinEqConstr" and "NumOfNonLinInEqConstr"

I have taken a look at Generating Hessian using Symbolic toolbox and few other web-pages but cannot see an example where the Hessian of the Lagrangian is constructed for dynamic number of constraints. In the referred matlab webpage example, like in one variation I tried replacing 10 with NumOfNonLinInEqConstr bu it doesn't work as matlabFunction does not work on cell data type. Can anybody provide a working example for constructing hessian of the lagrangian with dynamic number of nonlinear constraints

Relevant code lines from the matlab link provided below: My question is how could I avoid hard-coding 10 in the below code (and also avoid manual addition of the product of hessian of the constraint and lambda.ineqnonlin 10 times)?

    hessc = cell(1, 10);
    for i = 1:10
        hessc{i} = jacobian(gradc(:,i),x);
    end

    for i = 1:10
        ii = num2str(i);
        thename = ['hessc',ii];
        filename = [currdir,thename,'.m'];
        matlabFunction(hessc{i},'file',filename,'vars',{x});
    end

    function H = hessfinal(X,lambda)
    %
    % Call the function hessenergy to start
    H = hessenergy(X);

    % Add the Lagrange multipliers * the constraint Hessians
    H = H + hessc1(X) * lambda.ineqnonlin(1);
    H = H + hessc2(X) * lambda.ineqnonlin(2);
    H = H + hessc3(X) * lambda.ineqnonlin(3);
    H = H + hessc4(X) * lambda.ineqnonlin(4);
    H = H + hessc5(X) * lambda.ineqnonlin(5);
    H = H + hessc6(X) * lambda.ineqnonlin(6);
    H = H + hessc7(X) * lambda.ineqnonlin(7);
    H = H + hessc8(X) * lambda.ineqnonlin(8);
    H = H + hessc9(X) * lambda.ineqnonlin(9);
    H = H + hessc10(X) * lambda.ineqnonlin(10);

    end

Answer 1

[Matt J]

The number of nonlinear equality and inequality constraints can be recovered from the lambda structure passed to hessfinal, and could look something like the following,

    function H = hessfinal(X,lambda)

      lam_ineq=[lambda.ineqnonlin].';
      lam_eq=[lambda.eqnonlin].';

      M_ineq=length(lam_ineq); %Same as "NumOfNonLinInEqConstr" 
      M_eq=length(lam_eq);     %Same as "NumOfNonLinEqConstr"

      for i=M_ineq:-1:1

       H_ineq(:,:,i)=...  %Hessians of all the inequality constraints

      end 

      for i=M_eq:-1:1

       H_eq(:,:,i)=...  %Hessians of all the equality constraints

      end

     H= reshape(H_ineq,[],M_ineq)*lambda_ineq + ...
        reshape(H_eq,[],M_eq)*lambda_eq;

     H= reshape(H,numel(X),[]) +Hessian_objective(X);

    end

[Hari]

Thanks a lot Matt. The main code snippets provided by you have been immensely useful and with some modification it works for me perfectly. In the spirit of trying to learn why you followed a certain approach can you please clarify following:-

1. Why have you used transpose for lam_ineq and lam_eq

2. Why is the "for i =" has the start/endcounter in reverse format?

Also some of your codes lines doesnt work for me: -

A. If there are zero constraints, the algorithm fails

B. The code "reshape(H_ineq,[],M_ineq)*lambda_ineq + ... reshape(H_eq,[],M_eq)*lambda_eq" errors out saying matrix dimensions do not agree. The column result of reshape and row result of lambda do not agree. Infact my question is, is it right to multiply entire lambda_eq or lamda_ineq with corresponding H_eq or H_ineq as we need to do the muliplication for each constraint seperately and then add them up, right?

In view of above, please see the modifications I made..the below code seems to work but I would love to hear your feedback and any suggestions for improvement.

(one more question: In my code below, I'm forced to specify global statement for incorporating the Hessian of individual constraint functions as well as Hessian of Objective function. Is there a way global can be avoided?)

function H = hessfinal(X,lambda)
global NonLinInEqConstrHessianF NonLinEqConstrHessianF DCRObjHessianF qvalue; %qvalue is an additional parameter associated with objective function;

lam_ineq=[lambda.ineqnonlin].';
lam_eq=[lambda.eqnonlin].';
M_ineq=length(lam_ineq); %Same as "NumOfNonLinInEqConstr" 
M_eq=length(lam_eq);     %Same as "NumOfNonLinEqConstr"

ProdHAndLambda_ineq = zeros(numel(X),numel(X)); %for the case when M_ineq is zero;%prod stands for product of H and corresponding lamda;
ProdHAndLambda_eq = zeros(numel(X),numel(X)); %for the case when M_eq is zero;

for CnterNonLinInEqConstr=M_ineq:-1:1
    H_ineq(:,:,CnterNonLinInEqConstr)= NonLinInEqConstrHessianF{CnterNonLinInEqConstr}(X);%Hessians of all the equality constraints
    ProdHAndLambda_ineq = ProdHAndLambda_ineq + H_ineq(:,:,CnterNonLinInEqConstr)*lam_ineq(CnterNonLinInEqConstr);
end

for CnterNonLinEqConstr=M_eq:-1:1
    H_eq(:,:,CnterNonLinEqConstr)= NonLinEqConstrHessianF{CnterNonLinEqConstr}(X);%Hessians of all the equality constraints
    ProdHAndLambda_eq = ProdHAndLambda_eq + H_eq(:,:,CnterNonLinEqConstr)*lam_eq(CnterNonLinEqConstr);
end     

H = ProdHAndLambda_eq + ProdHAndLambda_ineq + DCRObjHessianF(X,qvalue);

% H= reshape(H_ineq,[],M_ineq)*lam_ineq + ...
% reshape(H_eq,[],M_eq)*lam_eq;
% H= reshape(H,numel(X),[]) +DCRObjHessianF(X);    
end

Thanks again, Hari

[Matt J]

Why have you used transpose for lam_ineq and lam_eq

That looks likes like a mistake, now that I think of it. I wanted to ensure that they are column vectors so

 lam_eq=reshape(lambda.eqnonlin,[],1)
 lam_ineq=reshape(lambda.ineqnonlin,[],1)

may be more appropriate.

Why is the "for i =" has the start/endcounter in reverse format?

It's an alternative to pre-allocating the arrays H_eq and H_ineq. If you loop forward, then H_eq and H_ineq will grow in size throughout the loop (unless you pre-allocate them) resulting in inefficient repeated memory allocation.

If there are zero constraints, the algorithm fails

No idea what "zero constraints" means. No idea what "algorithm fails" means. If you mean that you sometimes have no nonlinear constraints and get errors that H_eq or H_ineq doesn't exist, then it's probably because of looping backwards instead of explicitly pre-allocating H_eq, H_ineq.

The code "reshape(H_ineq,[],M_ineq)*lambda_ineq + ... reshape(H_eq,[],M_eq)*lambda_eq" errors out saying matrix dimensions do not agree.

I'd have to know the dimensions of all matrix quantities in that line to know what's wrong. It might be because lambda_eq,lambda_ineq are not column vectors and my correction above might fix this.

In any case, the whole idea of that line was to do the multiply-and-sum with the Lagrange multipliers as a matrix-vector multiplication, which of course MATLAB does very fast. However, since the number of constraints is small, it's probably just as fast to apply the lambda weights inside the for-loops as you have done in your code. Except that in that case it is unnecessary and less memory efficient to build and store the 3D arrays H_eq and H_ineq. You could just as well do,

 for CnterNonLinInEqConstr=1:M_ineq

 ProdHAndLambda_ineq = ProdHAndLambda_ineq +lam_ineq(CnterNonLinInEqConstr) * NonLinInEqConstrHessianF{CnterNonLinInEqConstr}(X);

 end

and similarly for the equality constraints.

I'm forced to specify global statement for incorporating the Hessian of individual constraint functions as well as Hessian of Objective function. Is there a way global can be avoided?

See here for better techniques for passing fixed parameters, like qvalue, to functions

http://www.mathworks.com/help/optim/ug/passing-extra-parameters.html

Incidentally, the variables NonLinInEqConstrHessianF, NonLinEqConstrHessianF, DCRObjHessianF look like external function handles of some sort. No idea why you wouldn't just write all the code for the Hessian computation inside hessfinal().

[Hari]

Thank you for the clarifications regarding "for" loop and zero # of constraints. I will try the reshape approach for multiplying hessian of constraint and lambda and get back in case of ques.

You are right about removing H_eq/H_ineq.

The best part for me was the matlab link on passing extra parameters. I have managed to avoid global in all my scripts/functions. I used the anonymous function approach. One question. Are there any known speed differences between anonymous versus nested function approach?

[Hari]

One additional ques regarding usage of anonymous function for passing extra parameter.

My (limited) analysis indicated that for cases when I can create a malabFunction with the function saved in specific location using 'file' option and if for that same scenario if I use anonymous function, then the named matlabFunction is faster than anonymous func approach. If it is true, my question is, can the functionality of anonymous function used to pass extra parameters be replicated by non-anonymous function using matlabFunction 'file' approach? (I couldn't work out the syntax unless I use the subs function within matlabFunction definition)

[Matt J]

I'm not really sure I understand the comparison that was done. As best I can tell, you have 2 implementations of the same function(s):

1. As a symbolic expression generated with the Symbolic Math Toolbox.
2. As an mfile, generated from 1. using matlabFunction using optimization flags.

I would definitely expect the optimized mfile implementation to run faster than the symbolic implenetation.

If you were to then wrap either 1. or 2. in an anonymous function, for the purpose of passing extra parameters, I wouldn't expect either implementation's performance to change very much. That's assuming the computations done inside the function are considerably more time-consuming than the overhead needed to call the function itself. Otherwise, if your objective function and gradient calculations are trivially fast, you're at a point where pursuing speed optimization doesn't really make much sense.