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

用于低内存的 Hessian 乘法函数

fmincon interior-point 和 trust-region-reflective 算法以及 fminunc trust-region 算法可以求解 Hessian 密集但结构化的问题。对于这些问题，fmincon 和 fminunc 不会直接使用 Hessian H 计算 H*Y，因为形成 H 会占用大量内存。相反，您必须为 fmincon 或 fminunc 提供一个函数，该函数给定矩阵 Y 和有关 H 的信息，计算 W = H*Y。

在这个示例中，目标函数是非线性的，并且存在线性等式，因此使用 fmincon。该描述适用于 trust-region reflective 算法；fminunc trust-region 算法类似。对于 interior-point 算法，请参阅黑塞矩阵乘法函数中的 HessianMultiplyFcn 选项。目标函数的结构如下

$f (x) = f_{}^{ˆ} (x) - \frac{1}{2} x^{T} V V^{T} x,$

其中 $V$ 是一个 1000×2 的矩阵。 $f$ 的 Hessian 是密集的，但 $f_{}^{ˆ}$ 的 Hessian 是稀疏的。如果 $f_{}^{ˆ}$ 的 Hessian 是 $H_{}^{ˆ}$ ，那么 $f$ 的 Hessian $H$ 就是

$H = H_{}^{ˆ} - V V^{T} .$

为了避免直接使用 $H$ 可能发生的过多内存占用，示例提供了一个 Hessian 乘法函数 hmfleq1。此函数在传递矩阵 Y 时，使用稀疏矩阵 Hinfo（对应于 $H_{}^{ˆ}$ ）和 V 来计算 Hessian 矩阵乘积

W = H*Y = (Hinfo - V*V')*Y.

在此示例中，Hessian 乘法函数需要 $H_{}^{ˆ}$ 和 V 来计算 Hessian 矩阵乘积。V 是一个常量，因此您可以在匿名函数的函数句柄中捕获 V。

但是， $H_{}^{ˆ}$ 不是常量，必须在当前 x 进行计算。您可以通过在目标函数中计算 $H_{}^{ˆ}$ 并在第三个输出参量中将 $H_{}^{ˆ}$ 作为 Hinfo 返回来实现此目的。通过使用 optimoptions 将 HessianMultiplyFcn 选项设置为此处列出的 hmfleq1 Hessian 乘法函数的句柄，fmincon 知道从目标函数中获取 Hinfo 值并将其传递给 Hessian 乘法函数 hmfleq1。

步骤 1：编写一个文件 brownvv.m，计算目标函数、梯度和 Hessian 的稀疏部分。

该示例将 brownvv 传递给 fmincon 作为目标函数。brownvv.m 文件很长，这里就不列出了。您可以使用命令查看代码

type brownvv

因为 brownvv 计算梯度以及目标函数，所以示例（步骤 3）使用 optimoptions 将 SpecifyObjectiveGradient 选项设置为 true。

步骤 2：编写一个函数，给定矩阵 Y 来计算 H 的 Hessian 矩阵乘积。

现在，定义一个函数 hmfleq1，它使用 Hinfo（在 brownvv 中计算）和 V（可以在匿名函数的函数句柄中捕获）来计算 Hessian 矩阵乘积 W 其中 W = H*Y = (Hinfo - V*V')*Y。该函数必须具有以下形式

W = hmfleq1(Hinfo,Y)

第一个参量必须与目标函数 brownvv 返回的第三个参量相同。Hessian 乘法函数的第二个参量是矩阵 Y（W = H*Y 的）。

因为 fmincon 期望第二个参量 Y 用于形成 Hessian 矩阵乘积，所以 Y 始终是一个具有 n 行的矩阵，其中 n 是问题中的维数。Y 中的列数可能有所不同。最后，您可以使用匿名函数的函数句柄来捕获 V，因此 V 可以作为 hmfleq1 的第三个参量。该技术在传递额外参数中有更详细的描述。这是文件 hmfleq1.m 的列表。

type hmfleq1.m

function W = hmfleq1(Hinfo,Y,V)
%HMFLEQ1 Hessian-matrix product function for BROWNVV objective.
% Documentation example.
% W = hmfbx4(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.

%   Copyright 1984-2008 The MathWorks, Inc.
  
W = Hinfo*Y - V*(V'*Y);

步骤 3：调用具有起点和线性等式约束的非线性最小化例程。

从 fleq1.mat 加载问题参数 V 和稀疏等式约束矩阵 Aeq 和 beq，运行此示例时可用。使用 optimoptions 将 SpecifyObjectiveGradient 选项设置为 true ，并将 HessianMultiplyFcn 选项设置为指向 hmfleq1 的函数句柄。使用目标函数 brownvv 调用 fmincon，并以 V 作为附加参数。查看运行该过程的代码：

type runfleq1.m

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',...
    'SpecifyObjectiveGradient',true, ...
    'HessianMultiplyFcn',@(Hinfo,Y)hmfleq1(Hinfo,Y,V),...
    'Display','iter',...
    'OptimalityTolerance',1e-9,...
    'FunctionTolerance',1e-9);
[x,fval,exitflag,output] = fmincon(@(x)brownvv(x,V),xstart,[],[],Aeq,beq,[],[], ...
                                    [],options);

要运行此代码，请输入

[fval,exitflag,output,x] = runfleq1;

                                Norm of      First-order 
 Iteration        f(x)          step          optimality   CG-iterations
     0            2297.63                      1.41e+03                
     1            1084.59         6.3903            578           1
     2            1084.59            100            578           3
     3            1084.59             25            578           0
     4            1084.59           6.25            578           0
     5            1047.61         1.5625            240           0
     6            761.592          3.125           62.4           2
     7            761.592           6.25           62.4           4
     8            746.478         1.5625            163           0
     9            546.578          3.125           84.1           2
    10            274.311           6.25           26.9           2
    11            55.6193        11.6597             40           2
    12            55.6193             25             40           3
    13            22.2964           6.25           26.3           0
    14            -49.516           6.25             78           1
    15           -93.2772         1.5625             68           1
    16           -207.204          3.125           86.5           1
    17           -434.162           6.25           70.7           1
    18           -681.359           6.25           43.7           2
    19           -681.359           6.25           43.7           4
    20           -698.041         1.5625            191           0
    21           -723.959          3.125            256           7
    22            -751.33        0.78125            154           3
    23           -793.974         1.5625           24.4           3
    24           -820.831        2.51937           6.11           3
    25           -823.069       0.562132           2.87           3
    26           -823.237       0.196753          0.486           3
    27           -823.245      0.0621202          0.386           3
    28           -823.246      0.0199951           0.11           6
    29           -823.246     0.00731333         0.0404           7
    30           -823.246     0.00505883         0.0185           8
    31           -823.246     0.00126471        0.00268           9
    32           -823.246     0.00149326        0.00521           9
    33           -823.246    0.000373314        0.00091           9

Local minimum possible.

fmincon stopped because the final change in function value relative to 
its initial value is less than the value of the function tolerance.

对于这种规模的问题来说，收敛很快，并且随着优化的进展，PCG 迭代成本会适度增加。在解中保持等式约束的可行性。

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

预条件

在这个示例中，fmincon 不能使用 H 来计算预条件子，因为 H 仅隐式存在。fmincon 不使用 H，而是使用 Hinfo（由 brownvv 返回的第三个参量）来计算预条件子。Hinfo 是一个不错的选择，因为它的大小与 H 相同，并且在某种程度上近似于 H。如果 Hinfo 与 H 的大小不一样，fmincon 将根据算法确定的一些对角缩放矩阵计算预条件子。通常情况下，这不会有很好的效果。

另请参阅

fmincon | fminunc