使用密集、结构化的 Hessian 矩阵进行二次最小化

打开实时脚本

quadprog trust-region-reflective 方法可以求解 Hessian 密集但结构化的大型问题。对于这些问题，quadprog 不会直接使用 Hessian H*Y 计算 H，就像它对稀疏 H 的信赖区域反射问题所做的那样，因为形成 H 会占用大量内存。相反，您必须为 quadprog 提供一个函数，该函数给定矩阵 Y 和有关 H 的信息，计算 W = H*Y。

在此示例中，Hessian 矩阵 H 具有结构 H = B + A*A'，其中 B 是一个稀疏的 512×512 对称矩阵，而 A 是一个由多个密集列组成的 512×10 稀疏矩阵。由于 H 较为密集，为了避免直接使用 H 时可能发生的内存占用过多，示例提供了一个 Hessian 乘法函数 qpbox4mult。此函数在传递矩阵 Y 时，使用稀疏矩阵 A 和 B 来计算 Hessian 矩阵乘积 W = H*Y = (B + A*A')*Y。

在此示例的第一部分中，需要将矩阵 A 和 B 提供给 Hessian 乘法函数 qpbox4mult。您可以将一个矩阵作为第一个参量传递给 quadprog，然后将其传递给 Hessian 乘法函数。您可以使用嵌套函数来提供第二个矩阵的值。

示例的第二部分展示了如何收紧 TolPCG 容差来补偿近似预条件子而不是精确的 H 矩阵。

步骤 1：决定将 H 的哪部分传递给 `quadprog` 作为第一个参量。

A 或 B 都可以作为 quadprog 的第一个参量传递。该示例选择传递 B 作为第一个参量，因为这样会产生更好的预条件子（请参阅预条件）。

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

现在，定义一个函数 runqpbox4，它：

包含一个嵌套函数 qpbox4mult，该函数使用 A 和 B 来计算 Hessian 矩阵乘积 W，其中 W = H*Y = (B + A*A')*Y。嵌套函数必须具有以下形式，其中前两个参量 Hinfo 和 Y 是必需的：

W = qpbox4mult(Hinfo,Y,...)

从 qpbox4.mat 加载问题参数。
使用 optimoptions 将 HessianMultiplyFcn 选项设置为指向 qpbox4mult 的函数句柄。
以 quadprog 作为第一个参量调用 B。

嵌套函数 qpbox4mult 的第一个参量必须与传递给 quadprog 的第一个参量相同，在本例中为矩阵 B。

qpbox4mult 的第二个参量是矩阵 Y (W = H*Y 的)。因为 quadprog 期望 Y 用于形成 Hessian 矩阵乘积，所以 Y 始终是一个具有 n 行的矩阵，其中 n 是问题中的维数。Y 中的列数可能有所不同。函数 qpbox4mult 是嵌套的，因此矩阵 A 的值来自外部函数。runqpbox4.m 文件出现在本示例的末尾。

步骤 3：调用具有起点的二次最小化例程。

要调用 runqpbox4 中包含的二次最小化例程，请输入

[fval,exitflag,output] = runqpbox4;

Local minimum possible.

quadprog stopped because the relative change in function value is less than the sqrt of the function tolerance, the rate of change in the function value is slow, and no negative curvature was detected.

然后显示 fval、exitflag、output.iterations 和 output.cgiterations 的值。

fval,exitflag,output.iterations, output.cgiterations

fval = 
-1.0538e+03

exitflag = 
3

ans = 
18

ans = 
30

经过 18 次迭代，共 30 次 PCG 迭代后，函数值降低为

fval

fval = 
-1.0538e+03

一阶最优解为

output.firstorderopt

ans = 
0.0043

预条件

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

由于预条件子比明确提供 H 时更加近似，因此可能需要将 TolPCG 参数调整为稍小的值。此示例与前一个示例相同，但将 TolPCG 从默认的 0.1 降低到 0.01。runqpbox4prec 函数列在本示例的末尾。

[fval,exitflag,output] = runqpbox4prec;

Local minimum possible.

quadprog stopped because the relative change in function value is less than the sqrt of the function tolerance, the rate of change in the function value is slow, and no negative curvature was detected.

经过 18 次迭代和 50 次 PCG 迭代后，函数值具有五位有效数字相同的值

fval

fval = 
-1.0538e+03

并且一阶最优性本质上是相同的：

output.firstorderopt

ans = 
0.0028

注意：过多减少 TolPCG 会显著增加 PCG 迭代次数。

辅助函数

此代码创建 runqpbox4 辅助函数。

function [fval, exitflag, output, x] = runqpbox4
%RUNQPBOX4 demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds.

problem = load('qpbox4'); % Get xstart, u, l, B, A, f
xstart = problem.xstart; u = problem.u; l = problem.l;
B = problem.B; A = problem.A; f = problem.f;
mtxmpy = @qpbox4mult; % function handle to qpbox4mult nested function

% Choose algorithm and the HessianMultiplyFcn option
options = optimoptions(@quadprog,'Algorithm','trust-region-reflective','HessianMultiplyFcn',mtxmpy);

% Pass B to qpbox4mult via the H argument. Also, B will be used in
% computing a preconditioner for PCG.
[x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options);

    function W = qpbox4mult(B,Y)
        %QPBOX4MULT Hessian matrix product with dense structured Hessian.
        %   W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where
        %   INPUT:
        %       B - sparse square matrix (512 by 512)
        %       Y - vector (or matrix) to be multiplied by B + A'*A.
        %   VARIABLES from outer function runqpbox4:
        %       A - sparse matrix with 512 rows and 10 columns.
        %
        %   OUTPUT:
        %       W - The product (B + A*A')*Y.
        %

        % Order multiplies to avoid forming A*A',
        %   which is large and dense
        W = B*Y + A*(A'*Y);
    end

end

以下代码创建 runqpbox4prec 辅助函数。

function [fval, exitflag, output, x] = runqpbox4prec
%RUNQPBOX4PREC demonstrates 'HessianMultiplyFcn' option for QUADPROG with bounds.

problem = load('qpbox4'); % Get xstart, u, l, B, A, f
xstart = problem.xstart; u = problem.u; l = problem.l;
B = problem.B; A = problem.A; f = problem.f;
mtxmpy = @qpbox4mult; % function handle to qpbox4mult nested function

% Choose algorithm, the HessianMultiplyFcn option, and override the TolPCG option
options = optimoptions(@quadprog,'Algorithm','trust-region-reflective',...
    'HessianMultiplyFcn',mtxmpy,'TolPCG',0.01);

% Pass B to qpbox4mult via the H argument. Also, B will be used in
% computing a preconditioner for PCG.
% A is passed as an additional argument after 'options'
[x, fval, exitflag, output] = quadprog(B,f,[],[],[],[],l,u,xstart,options);

    function W = qpbox4mult(B,Y)
        %QPBOX4MULT Hessian matrix product with dense structured Hessian.
        %   W = qpbox4mult(B,Y) computes W = (B + A*A')*Y where
        %   INPUT:
        %       B - sparse square matrix (512 by 512)
        %       Y - vector (or matrix) to be multiplied by B + A'*A.
        %   VARIABLES from outer function runqpbox4prec:
        %       A - sparse matrix with 512 rows and 10 columns.
        %
        %   OUTPUT:
        %       W - The product (B + A*A')*Y.
        %

        % Order multiplies to avoid forming A*A',
        %   which is large and dense
        W = B*Y + A*(A'*Y);
    end

end

另请参阅

quadprog

主题

雅可比乘法函数与线性最小二乘法