Preconditioning in fmincon

I am trying to find the minimum of a nonlinear, multivariate function using fmincon in MATLAB. Currently, my set of options is

 options = optimset('Algorithm', 'interior-point', 'Hessian',{'lbfgs', 5}, 'Display', 'iter','MaxIter', 75, 'MaxFunEvals', 12000, 'TolCon', 1e-10, 'TolFun', 1e-10);

I would like to precondition the Hessian matrix, but I can't figure out how to do so using the current command and options set. Any advice or direction on this matter would be great.

[Mariano]

I have a similar question.

I am using fmincon with the algorithm trust-region-reflective and the option HessianMultiplyFcn, so that the quadratic subproblems that appear in the process are solved internally by the preconditioned conjugate gradient method.

If I have understood correcty the documentation, fmincon somehow builds a preconditioner by itsef, but for my problem it is not very effective.

I would like to know if there is a way to pass a specific preconditioner. I have in mind a diagonal matrix D.

Thanks,

Mariano

Answer 1

[Catalytic]

I would like to know if there is a way to pass a specific preconditioner. I have in mind a diagonal matrix D.

Unless I am mistaken, preconditioning is equivalent to making a change of variables x=D*y in the optimization problem. So, you could just do -

fun=@(y)wrapper(y,D,fun);

x = D*fmincon(fun,D\x0,A*D, b, Aeq*D,beq,D\lb,D\ub);

function varargout=wrapper(y,D,fun)

  [varargout{1:nargout}]=fun(D*y);
  
  if nargout>1
      varargout{2}=D'*varargout{2}; %scale the gradient
  end
  
  if nargout>2
     varargout{3}= D'*varargout{3}*D; %scale the Hessian
  end

end

If you have nonlinear constraints, you would need a similar wrapper for the nonlcon argument as well.

[Mariano]

Thanks for the answer.

It works fine.

Just one minor comment. I think that for the lower and upper bounds one has to use D\lb and D\ub.

[Catalytic]

Yes, I think you're right.

Answer 2

[Matt J]

EDIT: You cannot pass a preconditioner on its own, nor would you want to if the true Hessian can be computed. However, using either the HessianFcn or HessianMultiplyFcn options, you can return a matrix of the form D*H*D' to simulate the effect of preconditioing. In other words, you have to take on the responsibility of computing the entirety of what you want the Hessian approximation to be.

[Mariano]

Thanks for your answer, Matt. I have tried it, but I think that this idea cannot work.

To my best knowledge, the algorithms inside fmincon compute a direction d, apply some rules to compute a stepsize r, and then iterate x(n+1) =x(n) + r*d

To find the direction d, the idea is to solve a the linear system of the form

Hd=-g

where H is (a submatrix of) the Hessian, and g is (a subvector) of the gradient.

If I pass D*H*D' as Hessian, the system

(D*H*D')d = -g

has a different solution, and we have wrong descend directions. The result is that each linear system is solved in fewer iterations, but the resulting direction is not correct,

I have also tried passing inv(D)*H for the hessian and inv(D)*g for the gradient. In this way I should get the correct directions d. The problem now is that inv(D)*g is not the gradient of the objective functional, and hence I get a wrong solution at the end.

Finally, an option that works is passing D*H*inv(D) for the Hessian and g for the gradient, but the system has the same condition number as the original one, so it is useless.

All the best,

Mariano

[Matt J]

The true Hessian is the ideal preconditioner, so if the true Hessian can be computed, there is no point to preconditioning artificially. However, if computing the true Hessian is too burdensome, an approximate Hessian can work in the trust-region algorithm, as shown in the example below.

Q=rand(4); Q=Q*Q';
Qapprox=diag(diag(Q));

x0=rand(4,1);
e=ones(4,1);

tol=1e-10;
opts=optimoptions('fmincon','Algorithm','trust-region-reflective', ...
                            'SpecifyObjectiveGradient', true, ...
                            'HessianFcn',"objective",'StepTol',0, ...
                            'FunctionTol',0,'OptimalityTol',tol, ...
                            'MaxFunEvals',inf,'MaxIter',1e5);

%% Use true Hessian
fun=@(x)objFcn(x,Q);
[x,fval,ef,stats]=fmincon(fun,x0,[],[],[],[],-5*e,+5*e,[],opts)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 4×1

1.0000 2.0000 3.0000 4.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 5.8628e-20

ef = 1

stats = struct with fields:
         iterations: 1940
          funcCount: 1941
           stepsize: 4.1789e-11
       cgiterations: 3879
      firstorderopt: 9.8305e-11
          algorithm: 'trust-region-reflective'
            message: 'Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The first-order optimality measure, 9.830509e-11, ↵is less than options.OptimalityTolerance = 1.000000e-10, and no negative/zero↵curvature is detected in the trust-region model.'
    constrviolation: 0
       bestfeasible: [1×1 struct]

%% Use approximate Hessian
fun=@(x)objFcn(x,Q,Qapprox);
[x,fval,ef,stats]=fmincon(fun,x0,[],[],[],[],-5*e,+5*e,[],opts)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 4×1

1.0000 2.0000 3.0000 4.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

fval = 1.0224e-18

ef = 1

stats = struct with fields:
         iterations: 17243
          funcCount: 17244
           stepsize: 2.0751e-11
       cgiterations: 14714
      firstorderopt: 9.1189e-11
          algorithm: 'trust-region-reflective'
            message: 'Local minimum found.↵↵Optimization completed because the size of the gradient is less than↵the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The first-order optimality measure, 9.118858e-11, ↵is less than options.OptimalityTolerance = 1.000000e-10, and no negative/zero↵curvature is detected in the trust-region model.'
    constrviolation: 0
       bestfeasible: [1×1 struct]

function [f,g,H]=objFcn(x,Q,Qapprox)
arguments
  x (:,1);
  Q; Qapprox=Q;
end

  dx=(x-[1;2;3;4]);
  f=dx'*Q*dx/2;

  if nargout>1
       g=Q*dx;
       H=Qapprox;
  end

end

[Mariano]

Thank you again for your detailed response. I see what you mean, and it is a good hint.