Solve many similar fmincon problems with a fully Vectorized Objective function and Analytical gradient + hessian

I have many such problems:

Basically, is the location in my state-space. The objective function and constraints are almost identical at every point s, they just differ by some constant (for instance, the distance to the boundary ).

The function f (resp. ) is fully vectorized, meaning that I can perform all s evaluations at once, even when x is a vector, i.e.

Also, the evaluation of f and g is very cheap.

For now, to optimize using fmincon, I have to solve this problem for each s separately because fmincon can only handle scalar valued functions. parfor allows for some efficiency gains, but it is overall still quite slow as each problem is separated and sent to a different core.

This is my (stylized) code for now:

parfor ss = 1:NN
    
    min_ss = STATEmin(ss);
    max_ss = STATEmax(ss);
    Uss    = U(ss);
    Vss    = V(ss);

    % Objective function
    F   = @(X)  fobjective(X,Uss) ;
    DF  = @(X) [ f1grad(X,Uss)      ; f2grad(X,Uss)      ; f3grad(X,Uss)      ; f4grad(X,Uss) ] ; %4 partial derivatives

    D2F = @(X) [ f1_1hess(X(1),Uss) ; f1_2hess(X(1),Uss) ; f1_3hess(X(1),Uss) ; f1_4hess(X(1),Uss) ;
                 f2_1hess(X(2),Uss) ; f2_2hess(X(2),Uss) ; f2_3hess(X(2),Uss) ; f2_4hess(X(2),Uss) ;
                 f3_1hess(X(3),Uss) ; f3_2hess(X(3),Uss) ; f3_3hess(X(3),Uss) ; f3_4hess(X(3),Uss) ;
                 f4_1hess(X(4),Uss) ; f4_2hess(X(4),Uss) ; f4_3hess(X(4),Uss) ; f4_4hess(X(4),Uss) ] ; %4-by-4 (sparse)

    % non-linear inequality constraints (<=0)
    G  = @(X) [ g(X,Vss) - max_ss , min_ss   - g(X,Vss) ] ; %2M constraints

    DG = @(X) [ g1grad(X,Vss) ;
                g2grad(X,Vss) ;
                g3grad(X,Vss) ;
                g4grad(X,Vss) ]; %4-by-2M


    % non-linear equality constraints (==0)
    H  = @(X) [];
    DH = @(X) [];

    Objective = @(X) deal(F(X),DF(X));  
    ConsGrad  = @(X) deal(G(X),H(X),DG(X),DH(X));
    Hessian   = @(X,lmbd) D2F(X)  ] ;

    % % NB: can also specify Hessian of non-linear inequality constraints
    % Hessian   = @(X,lmbd) D2F(X) ...
    %     +  lmbd.ineqnonlin(1)*[ g1hess(X,W_ss)] ...
    %     ...
    %     +  lmbd.ineqnonlin(2M)*[ g2Mhess(X,W_ss) ] ; %each line is 4-by-4

    X0      = Xold(ss);
    Ax      = [];
    b       = [];
    Aeq     = [];
    beq     = [];
    lb      = [ Xmin(ss) ] -5 ;
    ub      = [ Xmax(ss) ] +5 ;

    options_fmin = optimoptions('fmincon','Algorithm','interior-point','SpecifyObjectiveGradient',true, ...
                    'SpecifyConstraintGradient',true','HessianFcn',Hessian,'Display','none', ...
                    'OptimalityTolerance',1e-6,'MaxIterations',1e3,'ConstraintTolerance',1e-6);

    [Xsol,fsol,exitflag,~] = fmincon(Objective,X0,Ax,b,Aeq,beq,lb,ub,ConsGrad,options_fmin);

    Xnew(ss)     = Xsol;
    f_mod(ss)    = fsol;
    flag_mod(ss) = exitflag;

end

What I would like to do is to run fmincon but for all s at once. It would be much faster as f and g are fully vectorized. However, the problem is that there are independent problems so f returns a vector of dimension and fmincon only works for scalar valued functions. I have tried to minimize the norm but it is a dirty trick and it doesn't work well in practise (the results are quite different from the concatenated results of each independent problem).

Any idea how to perform these optimization problem in a vectorized fashion, and not in a parralelized fashion?

I have seen that ga or MultiStart could be some options but I don't know how to implement these. Could someone please provide an example?

Thanks a lot! I would really appreciate the help :)

Answer 1

[Matt J]

To vectorize, you must,

(1) Sum the elements of f, not take its norm.

(2) Vectorize your constraints and their gradient evaluations.

(3) Your HessianFcn needs to implement a sparse block diagonal matrix multiplication, where each block corresponds to a state.

However, I would be rather surprised if vectorization turns out to be the faster strategy (or at least a hybrid of vectorization and parfor might be best).

One thing I can't understand from your code is how X0 is chosen, as Xold is never defined or updated anywhere. If the problem data varies gradually as a function of the state, I would expect you to use the solution from the previous state in the loop as the initial x0 at for the current state: X0=Xsol.

[Francois Miguet]

Hi, thanks for the answer!

Apologies, norm indeed doesn't make sense as you would change the objective.

In the above, there is an outer loop that is not shown. At every iteration k, I need to solve an optimization problem at every point s of my state space. It speeds things a bit to give as initial condition for iteration k the solution from step .

Xold = zeros(NN,n);

for k=1:maxit
    
    parfor ss = 1:NN
        ... # code shown above with Xold as initial guess
        ... # get Xnew      
    end
    
    if someotherobj(Xnew) < tol
        break;
    end
    
    Xold = Xnew; 
end