Main Content

具有边界约束和带状预条件子的最小化

此示例显示如何使用 fmincon trust-region-reflective 算法求解边界的非线性问题。当问题稀疏且具有解析梯度和已知结构(例如其 Hessian 模式)时,该算法提供了额外的效率。

带梯度的目标函数

对于给定的 n(4 的正倍数),目标函数为

f(x)=1+i=1n|(3-2xi)xi-xi-1-xi+1+1|p+i=1n/2|xi+xi+n/2|p,

其中 p=7/3x0=0xn+1=0本示例末尾tbroyfg 辅助函数实现了目标函数,包括其梯度。

该问题对所有 i 都有边界 -10xi10。使用 n = 800。

n = 800;
lb = -10*ones(n,1);
ub = -lb;

Hessian 模式

Hessian 矩阵的稀疏模式是预先确定的,并存储在文件 tbroyhstr.mat 中。正如您在下面的间谍图中所看到的,该问题的 Hessian 的稀疏结构是带状的。

load tbroyhstr
spy(Hstr)

在此图中,中心条纹本身是一个五带矩阵。下图更清晰地显示了矩阵。

spy(Hstr(1:20,1:20))

问题选项

设置选项以使用 trust-region-reflective 算法。该算法要求您将 SpecifyObjectiveGradient 选项设置为 true

另外,使用 optimoptionsHessPattern 选项设置为 Hstr。如果您没有为具有明显稀疏结构的大型问题设置此选项,则该问题将使用大量内存和计算,因为 fmincon 尝试对 640,000 个非零条目的完整 Hessian 矩阵使用有限差分。

options = optimoptions('fmincon','SpecifyObjectiveGradient',true,'HessPattern',Hstr,...
    'Algorithm','trust-region-reflective');

求解问题

对于奇数索引,将初始点设置为 -1;对于偶数索引,将初始点设置为 +1。

x0 = -ones(n,1);
x0(2:2:n) = 1;

该问题没有线性或非线性约束,因此将这些参数设置为 []

A = [];
b = [];
Aeq = [];
beq = [];
nonlcon = [];

调用 fmincon 以求解问题。

[x,fval,exitflag,output] = ... 
   fmincon(@tbroyfg,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
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.

检查解和求解过程

检查退出标志、目标函数值、一阶最优性测度和求解器迭代次数。

disp(exitflag);
     3
disp(fval)
  270.4790
disp(output.firstorderopt)
    0.0163
disp(output.iterations)
     7

fmincon 不需要经过多次迭代就能得出解。但是该解具有相对较高的一阶最优性测度,这就是退出标志没有首选值 1 的原因。

改进解

尝试使用五带预条件子而不是默认的对角预条件子。使用 optimoptions,将 PrecondBandWidth 选项设置为 2 并再次求解问题。(带宽是上对角线或下对角线的数量,不计算主对角线。)

options.PrecondBandWidth = 2;
[x2,fval2,exitflag2,output2] = ... 
   fmincon(@tbroyfg,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
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.
disp(exitflag2);
     3
disp(fval2)
  270.4790
disp(output2.firstorderopt)
   7.5340e-05
disp(output2.iterations)
     9

退出标志和目标函数值似乎没有改变。然而,随着迭代次数的增加,一阶最优性测度会大幅下降。计算目标函数值的差异。

disp(fval2 - fval)
  -2.9005e-07

目标函数值微小地减少。该解主要改进了一阶最优性测度,而不是目标函数。

辅助函数

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

function [f,grad] = tbroyfg(x,dummy)
%TBROYFG Objective function and its gradients for nonlinear minimization
% with bound constraints and banded preconditioner.
% Documentation example.

%   Copyright 1990-2008 The MathWorks, Inc.

n = length(x);  % n should be a multiple of 4
p = 7/3;
y=zeros(n,1);
i = 2:(n-1);
y(i) = abs((3-2*x(i)).*x(i) - x(i-1) - x(i+1)+1).^p;
y(n) = abs((3-2*x(n)).*x(n) - x(n-1)+1).^p;
y(1) = abs((3-2*x(1)).*x(1) - x(2)+1).^p;
j = 1:(n/2);
z = zeros(length(j),1);
z(j) = abs(x(j) + x(j+n/2)).^p;
f = 1 + sum(y) + sum(z);
%
% Evaluate the gradient.
if nargout > 1
   g = zeros(n,1);
   t = zeros(n,1);
   i = 2:(n-1);
   t(i) = (3-2*x(i)).*x(i) - x(i-1) - x(i+1) + 1;
   g(i) = p*abs(t(i)).^(p-1).*sign(t(i)).*(3-4*x(i));
   g(i-1) = g(i-1) - p*abs(t(i)).^(p-1).*sign(t(i));
   g(i+1) = g(i+1) - p*abs(t(i)).^(p-1).*sign(t(i));
   tt = (3-2*x(n)).*x(n) - x(n-1) + 1;
   g(n) = g(n) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(n));
   g(n-1) = g(n-1) - p*abs(tt).^(p-1).*sign(tt);
   tt = (3-2*x(1)).*x(1)-x(2)+1;
   g(1) = g(1) + p*abs(tt).^(p-1).*sign(tt).*(3-4*x(1));
   g(2) = g(2) - p*abs(tt).^(p-1).*sign(tt);
   j = 1:(n/2);
   t(j) = x(j) + x(j+n/2);
   g(j) = g(j) + p*abs(t(j)).^(p-1).*sign(t(j));
   jj = j + (n/2);
   g(jj) = g(jj) + p*abs(t(j)).^(p-1).*sign(t(j));
   grad = g;
end
end