|On this page…|
The fmincon interior-point and trust-region-reflective algorithms, and the fminunc trust-region algorithm can solve problems where the Hessian is dense but structured. For these problems, fmincon and fminunc do not compute H*Y with the Hessian H directly, because forming H would be memory-intensive. Instead, you must provide fmincon or fminunc with a function that, given a matrix Y and information about H, computes W = H*Y.
In this example, the objective function is nonlinear and linear equalities exist so fmincon is used. The description applies to the trust-region reflective algorithm; the fminunc trust-region algorithm is similar. For the interior-point algorithm, see the 'HessMult' option in Hessian. The objective function has the structure
where V is a 1000-by-2 matrix. The Hessian of f is dense, but the Hessian of is sparse. If the Hessian of is , then H, the Hessian of f, is
To avoid excessive memory usage that could happen by working with H directly, the example provides a Hessian multiply function, hmfleq1. This function, when passed a matrix Y, uses sparse matrices Hinfo, which corresponds to , and V to compute the Hessian matrix product
W = H*Y = (Hinfo - V*V')*Y
In this example, the Hessian multiply function needs and V to compute the Hessian matrix product. V is a constant, so you can capture V in a function handle to an anonymous function.
However, is not a constant and must be computed at the current x. You can do this by computing in the objective function and returning as Hinfo in the third output argument. By using optimoptions to set the 'Hessian' options to 'on', fmincon knows to get the Hinfo value from the objective function and pass it to the Hessian multiply function hmfleq1.
Now, define a function hmfleq1 that uses Hinfo, which is computed in brownvv, and V, which you can capture in a function handle to an anonymous function, to compute the Hessian matrix product W where W = H*Y = (Hinfo - V*V')*Y. This function must have the form
W = hmfleq1(Hinfo,Y)
The first argument must be the same as the third argument returned by the objective function brownvv. The second argument to the Hessian multiply function is the matrix Y (of W = H*Y).
Because fmincon expects the second argument Y to be used to form the Hessian matrix product, Y is always a matrix with n rows where n is the number of dimensions in the problem. The number of columns in Y can vary. Finally, you can use a function handle to an anonymous function to capture V, so V can be the third argument to 'hmfleqq'.
function W = hmfleq1(Hinfo,Y,V); %HMFLEQ1 Hessian-matrix product function for BROWNVV objective. % W = hmfleq1(Hinfo,Y,V) computes W = (Hinfo-V*V')*Y % where Hinfo is a sparse matrix computed by BROWNVV % and V is a 2 column matrix. W = Hinfo*Y - V*(V'*Y);
Load the problem parameter, V, and the sparse equality constraint matrices, Aeq and beq, from fleq1.mat, which is available in the optimdemos folder. Use optimoptions to set the GradObj and Hessian options to 'on' and to set the HessMult option to a function handle that points to hmfleq1. Call fmincon with objective function brownvv and with V as an additional parameter:
function [fval, exitflag, output, x] = runfleq1 % RUNFLEQ1 demonstrates 'HessMult' option for FMINCON with linear % equalities. problem = load('fleq1'); % Get V, Aeq, beq V = problem.V; Aeq = problem.Aeq; beq = problem.beq; n = 1000; % problem dimension xstart = -ones(n,1); xstart(2:2:n,1) = ones(length(2:2:n),1); % starting point options = optimoptions(@fmincon,'Algorithm','trust-region-reflective','GradObj','on', ... 'Hessian','user-supplied','HessMult',@(Hinfo,Y)hmfleq1(Hinfo,Y,V),'Display','iter', ... 'TolFun',1e-9); [x,fval,exitflag,output] = fmincon(@(x)brownvv(x,V),xstart,,,Aeq,beq,,, ... ,options);
To run the preceding code, enter
[fval,exitflag,output,x] = runfleq1;
Because the iterative display was set using optimoptions, this command generates the following iterative display:
Norm of First-order Iteration f(x) step optimality CG-iterations 0 1997.07 916 1 1072.57 6.31716 465 1 2 480.247 8.19711 201 2 3 136.982 10.3039 78.1 2 4 44.416 9.04685 16.7 2 5 44.416 100 16.7 2 6 44.416 25 16.7 0 7 -9.05631 6.25 52.9 0 8 -317.437 12.5 91.7 1 9 -405.381 12.5 1.11e+003 1 10 -451.161 3.125 327 4 11 -482.688 0.78125 303 5 12 -547.427 1.5625 187 5 13 -610.42 1.5625 251 7 14 -711.522 1.5625 143 3 15 -802.98 3.125 165 3 16 -820.431 1.13329 32.9 3 17 -822.996 0.492813 7.61 2 18 -823.236 0.223154 1.68 3 19 -823.245 0.056205 0.529 3 20 -823.246 0.0150139 0.0342 5 21 -823.246 0.00479085 0.0152 7 22 -823.246 0.00353697 0.00828 9 23 -823.246 0.000884242 0.005 9 24 -823.246 0.0012715 0.00125 9 25 -823.246 0.000317876 0.0025 9 Local minimum possible. fmincon stopped because the final change in function value relative to its initial value is less than the selected value of the function tolerance.
Convergence is rapid for a problem of this size with the PCG iteration cost increasing modestly as the optimization progresses. Feasibility of the equality constraints is maintained at the solution.
problem = load('fleq1'); % Get V, Aeq, beq V = problem.V; Aeq = problem.Aeq; beq = problem.beq; norm(Aeq*x-beq,inf) ans = 2.3093e-14
In this example, fmincon cannot use H to compute a preconditioner because H only exists implicitly. Instead of H, fmincon uses Hinfo, the third argument returned by brownvv, to compute a preconditioner. Hinfo is a good choice because it is the same size as H and approximates H to some degree. If Hinfo were not the same size as H, fmincon would compute a preconditioner based on some diagonal scaling matrices determined from the algorithm. Typically, this would not perform as well.