`fmincon` `interior-point` 算法的黑塞乘法函数

打开实时脚本

本例演示了如何使用基于 fmincon "interior-point" 算法的黑塞矩阵乘法函数。黑塞矩阵乘法函数可以通过直接返回雅可比矩阵的黑塞矩阵与任意向量相乘的结果，从而节省内存和时间，而无需生成完整的黑塞矩阵。黑塞矩阵乘法函数的语法如下：

W = HessMultFcn(x,lambda,v)

x 是当前点，lambda 是拉格朗日乘子结构体，v 是 fmincon 传递给 HessMultFcn 的向量。该函数计算出的输出 W 应为 H*v，其中 H 是拉格朗日量的黑塞矩阵：

$\nabla_{x x}^{2} L (x, λ) = \nabla^{2} f (x) + \sum_{i} λ_{i n e q} (i) \nabla^{2} {i n e q n o n l i n}_{i} (x) + \sum_{i} λ_{e q} (i) \nabla^{2} {e q n o n l i n}_{i} (x) .$

这里， $f (x)$ 是目标函数， $i n e q n o n l i n (x)$ 是当约束条件满足时取非正值的非线性不等式约束函数， $e q n o n l i n (x)$ 是当约束条件满足时取值为零的非线性等式约束函数。

目标函数和约束函数

对于偶数维的情况，使用罗森布罗克函数的多维版本 $n$ ：

$f (x) = \sum_{i = 1}^{n / 2} {(10 (x_{2 i} - (x_{2 i - 1}^{2})))}^{2} + {(1 - x_{2 i - 1})}^{2} .$

$f$ 的梯度为

$\frac{\partial f}{\partial x_{2 i - 1}} = - 400 x_{2 i - 1} x_{2 i} + 400 x_{2 i - 1}^{3}$

$\frac{\partial f}{\partial x_{2 i}} = 200 (x_{2 i} - x_{2 i - 1}^{2}) f o r i = 1 t o n / 2$ .

非线性约束条件是：x 的平方和不超过 $n / 2$ ：

$i n e q n o n l i n (x) = \sum_{i = 1}^{n} x_{i}^{2} - n / 2 .$

该函数的梯度是 $2 x$ 。

将这些函数表示为 MATLAB 程序。

function [f,g] = multirosenbrock(x)
% Get the problem size
n = length(x);  
if mod(n,2) ~= 0
    error('Input vector, x ,must have an even number of components.');
end
% Evaluate the vector function
odds  = (1:2:n)';
evens = (2:2:n)';
F = zeros(n,1);
F(odds,1)  = 1-x(odds);
F(evens,1) = 10.*(x(evens)-x(odds).^2); 
f = sum(F.^2);
if nargout >= 2 % Calculate gradient
    g = zeros(n,1);
    g(evens) = 200*(x(evens)-x(odds).^2);
    g(odds) = -2*(1 - x(odds)) - 400.*(x(evens)-x(odds).^2).*x(odds);
end
end

function [ineqnonlin,eqnonlin,gradineq,gradeq] = cons(x)
eqnonlin = [];
N = length(x);
ineqnonlin = sum(x.^2) - N/2;
if nargout > 2
    gradeq = [];
    gradineq = 2*x;
end
end

黑塞与黑塞乘法函数

目标函数 $f$ 的黑塞矩阵是一个块对角矩阵，其元素为

$\nabla^{2} f (x) = [\begin{array}{cccccc} 2 - 400 x_{2} + 1200 x_{1}^{2} & - 400 x_{1} \\ - 400 x_{1} & 200 \\ 2 - 400 x_{4} + 1200 x_{3}^{2} & - 400 x_{3} \\ - 400 x_{3} & 200 \\ 2 - 400 x_{n} + 1200 x_{n - 1}^{2} & - 400 x_{n - 1} \\ - 400 x_{n - 1} & 200 \end{array}]$ .

黑塞矩阵 $\nabla^{2} i n e q n o n l i n (x)$ 是单位矩阵的两倍。

编写一个函数 HessMultFcn，该函数对向量 $v$ 返回结果 $\nabla^{2} f (x) * v + 2 * l a m b d a . i n e q n o n l i n * v$ 。

function W = HessMultFcn(x,lambda,v)
v = v(:);
n = length(x);
x = x(:);
if mod(n,2) ~= 0
    error('Input vector, x ,must have an even number of components.');
end

m = n/2;
j = (1:m)';
O = 2*j-1;
E = 2*j;

tmp1 = 2 - 400*x(E) + 1200*x(O).^2;
tmp2 = -400*x(O);
HessfTimesV = zeros(n, 1);
HessfTimesV(O) = tmp1.*v(O) + tmp2.*v(E);
HessfTimesV(E) = tmp2.*v(O) + 200*v(E);

% Return Hessian times v
W = HessfTimesV + 2*lambda.ineqnonlin*v;
end

未设置选项时的难度

默认的 fmincon 选项会导致求解器无法求解这个大规模问题。

将变量 $n$ 的数量设置为 1000。

n = 1000;

生成一个具有 n 个分量的伪随机起点 x0。

rng default % For reproducibility
x0 = randn(n,1);

使用默认选项解决该问题。

[x1,fval1,eflag1,output1] = fmincon(@multirosenbrock,x0,[],[],[],[],[],[],@cons)

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+03.

x1 = 1000×1

    0.9958
    1.9081
   -0.2722
    1.5397
    0.2169
   -1.2915
   -0.5438
    0.3822
   -4.4181
    4.7854
   -2.0833
    3.4741
    0.7097
    0.0370
    0.6668
      ⋮

fval1 = 
6.3932e+04

eflag1 = 
0

output1 = struct with fields:
         iterations: 2
          funcCount: 3010
    constrviolation: 312.2817
           stepsize: 7.9055
          algorithm: 'interior-point'
      firstorderopt: 2.6050e+04
       cgiterations: 2
            message: 'Solver stopped prematurely.↵↵fmincon stopped because it exceeded the function evaluation limit,↵options.MaxFunctionEvaluations = 3.000000e+03.'
       bestfeasible: []

为了尝试获得更好的解，请将 SpecifyObjectiveGradient 和 SpecifyConstraintGradient 选项设置为 true，然后重试。记录求解器的运行时间。

opts1 = optimoptions("fmincon",...
    SpecifyObjectiveGradient=true,...
    SpecifyConstraintGradient=true);
tic
[x1,fval1,eflag1,output1] = fmincon(@multirosenbrock,x0,[],[],[],[],[],[],@cons,opts1)

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+03.

x1 = 1000×1

    0.7890
    0.6218
    0.7836
    0.6132
    0.7863
    0.6175
    0.7863
    0.6174
    0.7848
    0.6152
    0.7847
    0.6151
    0.7866
    0.6179
    0.7866
      ⋮

fval1 = 
22.8431

eflag1 = 
0

output1 = struct with fields:
         iterations: 651
          funcCount: 3000
    constrviolation: 0
           stepsize: 0.0148
          algorithm: 'interior-point'
      firstorderopt: 0.2428
       cgiterations: 6657
            message: 'Solver stopped prematurely.↵↵fmincon stopped because it exceeded the function evaluation limit,↵options.MaxFunctionEvaluations = 3.000000e+03.'
       bestfeasible: [1×1 struct]

toc

Elapsed time is 78.720656 seconds.

求解器无法得出精确解，且耗时极长。

使用黑塞矩阵乘法函数求解带约束的非线性问题

要使用黑塞矩阵乘法函数来求解问题，必须使用适当的选项。

options = optimoptions("fmincon",...
    Algorithm="interior-point",...
    SpecifyObjectiveGradient=true,...
    SpecifyConstraintGradient=true,...
    SubproblemAlgorithm="cg",...
    HessianMultiplyFcn=@HessMultFcn);

在满足非线性约束的条件下，使 multirosenbrock 函数值最小。

tic
[x,fval,eflag,output] = fmincon(@multirosenbrock,x0,[],[],[],[],[],[],@cons,options)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

fval = 
22.8374

eflag = 
1

output = struct with fields:
         iterations: 28
          funcCount: 870
    constrviolation: 0
           stepsize: 3.8640e-07
          algorithm: 'interior-point'
      firstorderopt: 2.0000e-08
       cgiterations: 788
            message: 'Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 2.000017e-08,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'
       bestfeasible: [1×1 struct]

tmult = toc

tmult = 
0.2275

fmincon 仅用约 30 次迭代和 900 次函数计算就解决了这个包含 1000 个变量的问题。

与 `LBFGS` 黑塞矩阵比较

fmincon "lbfgs" 黑塞矩阵近似提供了一种节省内存的替代方法，无需创建黑塞矩阵乘法函数。

optslbfgs = optimoptions("fmincon",SpecifyObjectiveGradient=true,...
    SpecifyConstraintGradient=true,HessianApproximation="lbfgs");
tic
[x2,fval2,eflag2,output2] = fmincon(@multirosenbrock,x0,[],[],[],[],[],[],@cons,optslbfgs)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

x2 = 1000×1

    0.7864
    0.6177
    0.7864
    0.6177
    0.7864
    0.6177
    0.7864
    0.6177
    0.7864
    0.6177
    0.7864
    0.6177
    0.7864
    0.6177
    0.7864
      ⋮

fval2 = 
22.8374

eflag2 = 
1

output2 = struct with fields:
         iterations: 206
          funcCount: 287
    constrviolation: 0
           stepsize: 7.2407e-08
          algorithm: 'interior-point'
      firstorderopt: 8.0768e-07
       cgiterations: 28
            message: 'Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 8.076780e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.'
       bestfeasible: [1×1 struct]

tlbfgs = toc

tlbfgs = 
0.6444

与黑塞矩阵乘法函数相比，在此情况下，LBFGS 黑塞矩阵近似法：

基本上返回相同的点和函数值
需要更多的迭代次数
需要更少的函数计算次数
具有更大的最终一阶优化度量
完成所需的时间大约是原来的三倍

根据具体问题不同，这些比较结果可能会有所变化。

另请参阅

fmincon

主题

使用稠密结构 Hessian 和线性等式进行最小化

fmincon interior-point 算法的黑塞乘法函数