本页对应的英文页面已更新,但尚未翻译。 若要查看最新内容,请点击此处访问英文页面。

编写标量目标函数

函数文件

标量目标函数文件接受一个输入,例如 x,并返回一个实数标量输出,例如 f。输入 x 可以是标量、向量或矩阵。一个函数文件可以返回多个输出(请参阅包括梯度和 Hessian 矩阵)。

例如,假设您的目标是包含三个变量(x、y 和 z)的函数:

f(x) = 3*(x – y)4 + 4*(x + z)2 / (1 + x2 + y2 + z2) + cosh(x – 1) + tanh(y + z)。

  1. 将此函数编写为接受向量 xin = [x;y;z] 并返回 f 的文件:

    function f = myObjective(xin)
    f = 3*(xin(1)-xin(2))^4 + 4*(xin(1)+xin(3))^2/(1+norm(xin)^2) ...
        + cosh(xin(1)-1) + tanh(xin(2)+xin(3));
  2. 将其以名为 myObjective.m 的文件保存到 MATLAB® 路径上的文件夹中。

  3. 检查函数是否能够正确进行计算:

    myObjective([1;2;3])
    
    ans =
        9.2666

有关如何包含额外参数的信息,请参阅传递额外参数。有关函数文件的更复杂示例,请参阅Minimization with Gradient and Hessian Sparsity PatternMinimization with Bound Constraints and Banded Preconditioner

局部函数和嵌套函数

函数可以作为局部函数 (MATLAB)或嵌套函数 (MATLAB)存在于其他文件中。使用局部函数或嵌套函数可以减少保存的不同文件的数量。使用嵌套函数还允许您访问额外的参数,如嵌套函数中所示。

例如,假设您要在非线性约束中所述的 ellipseparabola.m 约束条件下对函数文件中所述的 myObjective.m 目标函数进行最小化。您不必编写两个文件 myObjective.mellipseparabola.m,而只需编写一个将这两个函数作为局部函数包含在其中的文件:

function [x fval] = callObjConstr(x0,options)
% Using a local function for just one file

if nargin < 2
    options = optimoptions('fmincon','Algorithm','interior-point');
end

[x fval] = fmincon(@myObjective,x0,[],[],[],[],[],[], ...
    @ellipseparabola,options);

function f = myObjective(xin)
f = 3*(xin(1)-xin(2))^4 + 4*(xin(1)+xin(3))^2/(1+sum(xin.^2)) ...
    + cosh(xin(1)-1) + tanh(xin(2)+xin(3));

function [c,ceq] = ellipseparabola(x)
c(1) = (x(1)^2)/9 + (x(2)^2)/4 - 1;
c(2) = x(1)^2 - x(2) - 1;
ceq = [];

从点 [1;1;1] 开始求解受约束的求最小化问题:

[x fval] = callObjConstr(ones(3,1))

Local minimum found that satisfies the constraints.

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

x =
    1.1835
    0.8345
   -1.6439

fval =
    0.5383

匿名函数目标

使用匿名函数编写简单的目标函数。有关匿名函数的详细信息,请参阅What Are Anonymous Functions? (MATLAB)。Rosenbrock 函数非常简单,可以写成匿名函数:

anonrosen = @(x)(100*(x(2) - x(1)^2)^2 + (1-x(1))^2);
检查 anonrosen[-1 2] 处的计算是否正确:
anonrosen([-1 2])

ans =
   104
使用 fminunc 最小化 anonrosen 会生成以下结果:
options = optimoptions(@fminunc,'Algorithm','quasi-newton');
[x fval] = fminunc(anonrosen,[-1;2],options)

Local minimum found.

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

x =
    1.0000
    1.0000

fval =
  1.2266e-10

包括梯度和 Hessian 矩阵

为求解器提供导数

对于 fminconfminunc,您可以在目标函数中包含梯度。一般情况下,当包含梯度时,求解器更稳健,速度也会有所提升。请参阅 包含导数的好处。要同时包括二阶导数(Hessian 矩阵),请参阅包含 Hessian 矩阵

下表显示哪些算法可以使用梯度和 Hessian 矩阵。

求解器算法梯度Hessian 矩阵
fminconactive-set可选
interior-point可选可选(请参阅适用于 fmincon 内点算法的 Hessian 矩阵
sqp可选
trust-region-reflective必需可选(请参阅适用于 fminunc 信赖域或 fmincon 信赖域反射算法的 Hessian 矩阵
fminuncquasi-newton可选
trust-region必需可选(请参阅适用于 fminunc 信赖域或 fmincon 信赖域反射算法的 Hessian 矩阵

如何包含梯度

  1. 编写返回以下内容的代码:

    • 目标函数(标量)作为第一个输出

    • 梯度(向量)作为第二个输出

  2. 使用 optimoptionsSpecifyObjectiveGradient 选项设置为 true。适用时,也将 SpecifyConstraintGradient 选项设置为 true

  3. (可选)检查梯度函数是否与有限差分逼近匹配。请参阅 Checking Validity of Gradients or Jacobians

提示

为了获得最大的灵活性,请编写条件化代码。条件化意味着函数输出的数目可以变化,如以下示例中所示。条件化代码不会因 SpecifyObjectiveGradient 选项的值而出错。非条件化代码要求您适当地设置选项。

例如,假设有以下 Rosenbrock 函数

f(x)=100(x2x12)2+(1x1)2,

基于求解器求解受约束的非线性问题中对该函数进行了说明和绘图。f(x) 的梯度为

f(x)=[400(x2x12)x12(1x1)200(x2x12)],

rosentwo 是条件化函数,可返回求解器所需的任何内容:

function [f,g] = rosentwo(x)
% Calculate objective f
f = 100*(x(2) - x(1)^2)^2 + (1-x(1))^2;

if nargout > 1 % gradient required
    g = [-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1));
        200*(x(2)-x(1)^2)];
    
end

nargout 检查调用函数指定的参数个数。请参阅 Find Number of Function Arguments (MATLAB)。

为无约束优化设计的 fminunc 求解器允许您最小化 Rosenbrock 函数。通过设置 optionsfminunc 使用梯度和 Hessian 矩阵:

options = optimoptions(@fminunc,'Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true);

[-1;2] 开始运行 fminunc

[x fval] = fminunc(@rosentwo,[-1;2],options)
Local minimum found.

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

x =
    1.0000
    1.0000

fval =
  1.9886e-17

如果您有 Symbolic Math Toolbox™ 许可证,您可以自动计算梯度和 Hessian 矩阵,如 Symbolic Math Toolbox™ Calculates Gradients and Hessians中所述。

包含 Hessian 矩阵

您可以使用 fmincon 'trust-region-reflective''interior-point' 算法以及 fminunc 'trust-region' 算法来包含二阶导数。根据信息的类型和算法,可通过几种方法来包括 Hessian 矩阵信息。

您还必须包含梯度(将 SpecifyObjectiveGradient 设置为 true,如果适用,还必须将 SpecifyConstraintGradient 设置为 true)以便包含 Hessian 矩阵。

适用于 fminunc 信赖域或 fmincon 信赖域反射算法的 Hessian 矩阵.  这些算法要么没有约束,要么只有边界或线性等式约束。因此,Hessian 矩阵是由目标函数的二阶导数组成的矩阵。

包含 Hessian 矩阵作为目标函数的第三个输出。例如,Rosenbrock 函数的 Hessian 矩阵 H(x) 如下(请参阅如何包含梯度

H(x)=[1200x12400x2+2400x1400x1200].

在目标中包含此 Hessian 矩阵:

function [f, g, H] = rosenboth(x)
% Calculate objective f
f = 100*(x(2) - x(1)^2)^2 + (1-x(1))^2;

if nargout > 1 % gradient required
    g = [-400*(x(2)-x(1)^2)*x(1)-2*(1-x(1));
        200*(x(2)-x(1)^2)];
    
    if nargout > 2 % Hessian required
        H = [1200*x(1)^2-400*x(2)+2, -400*x(1);
            -400*x(1), 200];  
    end

end

HessianFcn 设置为 'objective'。例如,

options = optimoptions('fminunc','Algorithm','trust-region',...
    'SpecifyObjectiveGradient',true,'HessianFcn','objective');

适用于 fmincon 内点算法的 Hessian 矩阵.  该 Hessian 矩阵是拉格朗日函数的 Hessian 矩阵,其中 L(x,λ) 是

L(x,λ)=f(x)+λg,igi(x)+λh,ihi(x).

g 和 h 是向量函数,分别表示所有不等式和等式约束(指有界、线性和非线性约束),因此最小化问题的公式为

minxf(x) subject to g(x)0, h(x)=0.

有关详细信息,请参阅 受约束的最优性理论。拉格朗日方程的 Hessian 矩阵公式为

xx2L(x,λ)=2f(x)+λg,i2gi(x)+λh,i2hi(x).(1)

要包含 Hessian 矩阵,请用以下语法编写函数

hessian = hessianfcn(x,lambda)

hessian 是 n×n 矩阵(稀疏或稠密),其中 n 是变量的数目。如果 hessian 很大并且非零项相对较少,请将 hessian 表示为稀疏矩阵,以节省运行时间和内存。lambda 是具有与非线性约束相关联的拉格朗日乘数向量的结构体:

lambda.ineqnonlin
lambda.eqnonlin

fmincon 计算结构体 lambda,并将其传递给您的 Hessian 函数。hessianfcn 必须计算公式 1 中的总和。通过设置以下选项表明您提供了 Hessian 函数:

options = optimoptions('fmincon','Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
    'HessianFcn',@hessianfcn);

例如,要针对 Rosenbrock 函数包含 Hessian 矩阵(约束为单位圆盘 x12+x221),请注意约束函数 g(x)=x12+x2210 具有梯度和二阶导数矩阵

g(x)=[2x12x2]Hg(x)=[2002].

将 Hessian 函数写成

function Hout = hessianfcn(x,lambda)
% Hessian of objective
H = [1200*x(1)^2-400*x(2)+2, -400*x(1);
            -400*x(1), 200];
% Hessian of nonlinear inequality constraint
Hg = 2*eye(2);
Hout = H + lambda.ineqnonlin*Hg;

hessianfcn 保存到 MATLAB 路径。为了完成此示例,包含梯度的约束函数为

function [c,ceq,gc,gceq] = unitdisk2(x)
c = x(1)^2 + x(2)^2 - 1;
ceq = [ ];

if nargout > 2
    gc = [2*x(1);2*x(2)];
    gceq = [];
end

求解包括梯度和 Hessian 矩阵的问题。

fun = @rosenboth;
nonlcon = @unitdisk2;
x0 = [-1;2];
options = optimoptions('fmincon','Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
    'HessianFcn',@hessianfcn);
[x,fval,exitflag,output] = fmincon(fun,x0,[],[],[],[],[],[],@unitdisk2,options);

有关使用内点 Hessian 矩阵的其他示例,请参阅使用解析 Hessian 函数的 fmincon 内点算法Symbolic Math Toolbox™ Calculates Gradients and Hessians

Hessian 矩阵乘法函数.  fmincon interior-pointtrust-region-reflective 算法都允许您提供 Hessian 乘法函数,而无需提供完整的 Hessian 函数。此函数给出 Hessian 矩阵乘以向量的乘积结果,而不直接计算 Hessian 矩阵。这可以节省内存。SubproblemAlgorithm 选项必须为 'cg',Hessian 乘法函数才能发挥作用;这是 trust-region-reflective 的默认值。

这两种算法的语法不同。

  • 对于 interior-point 算法,语法如下

    W = HessMultFcn(x,lambda,v);

    结果 W 应为 H*v 的积,其中 H 是在 x 处拉格朗日方程的 Hessian 矩阵(请参阅公式 1),lambda 是拉格朗日乘数(由 fmincon 计算),v 是大小为 n×1 的向量。按如下所示设置选项:

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

    提供函数 HessMultFcn,该函数返回 n×1 向量,其中 n 是 x 的维数。HessianMultiplyFcn 选项使您能够传递 Hessian 矩阵乘以向量的结果,而无需计算 Hessian 矩阵。

  • trust-region-reflective 算法不涉及 lambda

    W = HessMultFcn(H,v);

    结果 W = H*vfmincon 传递 H 作为在目标函数的第三个输出中返回的值(请参阅 适用于 fminunc 信赖域或 fmincon 信赖域反射算法的 Hessian 矩阵)。fmincon 还传递 v,即包含 n 个行的向量或矩阵。v 中的列数可能会有所不同,因此请写入 HessMultFcn 以接受任意数量的列。H 不一定是 Hessian 矩阵;它可以是任何让您能够计算 W = H*v 的项。

    按如下所示设置选项:

    options = optimoptions('fmincon','Algorithm','trust-region-reflective',... 
        'SpecifyObjectiveGradient',true,'HessianMultiplyFcn',@HessMultFcn);

    有关将 Hessian 乘法函数和 trust-region-reflective 算法结合使用的示例,请参阅Minimization with Dense Structured Hessian, Linear Equalities

包含导数的好处

如果您不提供梯度,求解器会通过有限差分来估计梯度。如果提供梯度,则求解器不需要执行这种有限差分估计,从而可以节省时间且结果更准确,但对于复杂导数来说,有限差分估计可能更快。此外,求解器使用逼近的 Hessian 矩阵,它可能与真实的 Hessian 矩阵相差甚远。提供 Hessian 矩阵可以用更少的迭代次数来生成解。例如,请参阅

对于约束问题,提供梯度还有另一个优势。求解器可以到达一个点 x,这说明 x 是可行点,但对于此处的 xx 附近的有限差分总是导致不可行点。进一步假设在不可行点处的目标函数返回复数输出、InfNaN 或错误。在这种情况下,求解器可能会失败或过早停止。提供梯度可允许求解器继续。要获得此优势,您可能还需要包括非线性约束函数的梯度,并将 SpecifyConstraintGradient 选项设置为 true。请参阅 非线性约束

为内点 fmincon 选择输入 Hessian 矩阵逼近方法

fmincon interior-point 算法提供了许多选项用来选择输入 Hessian 矩阵逼近方法。有关语法详细信息,请参阅 Hessian as an Input。以下是选项以及对其相对特征的估计。

Hessian 矩阵相对内存使用量相对效率
'bfgs'(默认值)高(适用于大型问题)
'lbfgs'低到中等中等
'fin-diff-grads'中等
'HessianMultiplyFcn'低(取决于您的代码)中等
'HessianFcn'不确定(取决于您的代码)High (不确定(取决于您的代码)

请使用默认的 'bfgs' Hessian 矩阵,除非您

'lbfgs' 只有中等效率的原因是需要进行双倍计算。它采用计算成本相对高昂的 Sherman-Morrison 更新。由于 'lbfgs' 有内存使用限制,得到的迭代步可能不够准确。

'fin-diff-grads'HessianMultiplyFcn 只有中等效率的原因是它们使用共轭梯度法。它们能准确估计目标函数的 Hessian 矩阵,但它们不生成最准确的迭代步。有关详细信息,请参阅 fmincon Interior Point Algorithm,及其关于求解公式 36 的 LDL 方法和共轭梯度法的讨论。

相关主题