非精确线搜索(Zhang & Hager)的 BB 步长梯度下降法
考虑无约束优化问题：
其中目标函数
是可微的。
对于可微的目标函数 ，梯度下降法通过使用如下重复如下迭代格式：
求解 的最小值，其中
为第
步的步长。
令 ,
，定义两种 BB 步长：
和
。在经验上，采用 BB步长的梯度法往往可以取得比固定步长更好的结果。
初始化和迭代准备
输入信息： 为迭代的初始值（可以为一个矩阵），函数 fun 依次返回给定点的函数和梯度值，以及给出参数的结构体 opt 。
输出信息： 为求得的最优点，
为该点处梯度，out 为一个包含其它信息的结构体。
- out.nfe ：函数 fun 调用的次数
- out.fval0 ：初始点函数值
- out.msg ：标记是否达到收敛
- out.nrmG ：退出时该点处梯度的F-范数
- out.itr ：迭代步数。
varargin：可变长度输入参数列表
function [x, g, out]= fminGBB(x, fun, opts, varargin)
非精确线搜索(Zhang & Hager)的 BB 步长梯度下降法
从输入的结构体中读取参数或采取默认参数。
- opts.xtol ：针对自变量的停机准则
- opts.gtol ：针对梯度的停机准则
- opts.ftol ：针对函数值的停机准则
- opts.tau ：默认步长（第一步或 BB 步长失效时采用默认步长）
- opts.rhols ：线搜索准则中下降量参数
- opts.eta ：步长衰减率
- opts.gamma ： Zhang & Hager 非单调线索准则参数
- opts.maxit ：最大迭代步数
- opts.record ：输出的详细程度，当为 0 时不进行输出，大于等于 1 时输出每一步信息，为10 时额外记录每一步的函数值
if ~isfield(opts, 'xtol'); opts.xtol = 0; end
if ~isfield(opts, 'gtol'); opts.gtol = 1e-6; end
if ~isfield(opts, 'ftol'); opts.ftol = 0; end
if ~isfield(opts, 'tau'); opts.tau = 1e-3; end
if ~isfield(opts, 'rhols'); opts.rhols = 1e-4; end
if ~isfield(opts, 'eta'); opts.eta = 0.2; end
if ~isfield(opts, 'gamma'); opts.gamma = 0.85; end
if ~isfield(opts, 'maxit'); opts.maxit = 1000; end
if ~isfield(opts, 'record'); opts.record = 0; end
复制参数。
xtol = opts.xtol;
ftol = opts.ftol;
gtol = opts.gtol;
maxit = opts.maxit;
rhols = opts.rhols;
eta = opts.eta;
gamma = opts.gamma;
record = opts.record;
初始化。计算初始点 处的函数值和梯度。
[n,p] = size(x);
[f,g] = feval(fun, x, varargin{:});
out.nfe = 1; % 函数 fun 调用的次数
out.fval0 = f; % 初始点函数值
nrmG = norm(g, 'fro'); % F范数==2范数
Zhang & Hager 非单调线索准则参数。
Q = 1; % 初始值
Cval = f; % 初始值
tau = opts.tau; % 默认步长（第一步或 BB 步长失效时采用默认步长）
当 record 大于等于 1 时为每一步的输出打印表头，每一列分别为：当前迭代步、当前步长、当前步函数值、梯度范数、相对 变化量、相对函数值变化量、线搜索次数。
if (record >= 1)
fprintf('----------- fminBB ----------- \n');
fprintf('%4s \t %6s \t %8s \t %5s \t %7s \t %7s \t %3s \n', ...
'Iter', 'tau', 'f(X)', 'nrmG', 'XDiff', 'FDiff', 'ls-Iter');
end
当 record 为 10 时，额外记录每一步的函数值。
out.fvec
if record == 10; out.fvec = f; end
初始化求解结果为超过最大迭代次数（即在达到最大迭代步数过程中没有满足任何针对自变量、函数值、梯度的停机准则）
out.msg = 'exceed max iteration';
迭代主循环
以 maxit 为最大迭代次数。在每一步迭代开始时复制前一步的 、函数值和梯度。
for itr = 1 : maxit
xp = x;
fp = f;
gp = g;
非精确线搜索。初始化线搜索次数 nls = 1。满足线搜索准则(Zhang & Hager) 或进行超过 10 次步长衰减后退出线搜索，否则以
的比例对步长进行衰减。
nls = 1;
while 1
x = xp - tau*gp;
[f,g] = feval(fun, x, varargin{:});
out.nfe = out.nfe + 1; % 函数 fun 调用的次数
if f <= Cval - tau*rhols*nrmG^2 || nls >= 10 % d=-nrmG
break
end
tau = eta*tau; % 步长更新，eta衰减率
nls = nls+1;
end
当 record 为 10 时，记录每一步迭代的函数值。
if record == 10
out.fvec = [out.fvec; f];
end
nrmG 为 处的2范数，XDiff 表示
与上一步迭代
之前的相对变化， FDiff 表示函数值的相对变化，当 record 大于等于 1 时将这些信息输出。
XDiff = norm(s,‘fro’)/sqrt(n);
FDiff = abs(fp-f)/(abs(fp)+1);
nrmG = norm(g, 'fro');
s = x - xp;
XDiff = norm(s,'fro')/sqrt(n);
FDiff = abs(fp-f)/(abs(fp)+1);
if (record >= 1)
fprintf('%4d \t %3.2e \t %7.6e \t %3.2e \t %3.2e \t %3.2e \t %2d\n', ...
itr, tau, f, nrmG, XDiff, FDiff, nls);
end
判断是否收敛，当下列停机准则至少一个被满足时停止迭代，并记录 out.msg 为收敛：
(1)相对的 和函数值的相对变化量均小于给定的阈值；
(2)当前梯度的范数小于给定的阈值。
if ( XDiff < xtol && FDiff < ftol ) || nrmG < gtol
out.msg = 'converge';
break;
end
BB 步长的计算，以 BB 步长作为线搜索的初始步长。令 ,
，这里在偶数与奇数步分别对应
和
两个 BB 步长，当
时使用默认步长。
y = g - gp;
sy = abs(iprod(s,y)); % 辅助函数->矩阵形式的内积
tau = opts.tau;
if sy > 0
if mod(itr,2)==0
tau = abs(sum(sum(s.*s)))/sy;
else
tau = sy/abs(sum(sum(y.*y)));
end
限定在 中。
tau = max(min(tau, 1e20), 1e-20);
end
计算 (Zhang & Hager) 线搜索准则中的递推常数，其满足 ，序列
满足
。
Qp = Q;
Q = gamma*Qp + 1;
Cval = (gamma*Qp*Cval + f)/Q;
end
记录输出。
out.nrmG = nrmG; % 退出时该点处梯度的F-范数
out.fval = f;
out.itr = itr; % 迭代步数。
end
辅助函数：矩阵形式的内积。
为了在复数域上良定义，取实部。
function a = iprod(x,y)
a = real(sum(sum(x.*y)));
end
实例：Tikhonov 正则化模型用于图片去噪
对于真实图片 和带噪声的图片
（其中
是高斯白噪声）。Tikhonov 正则化模型为：
其中 ,
分别表示
在水平和竖直方向上的向前差分，
为正则化系数。上述优化问题的目标函数中，第二项要求恢复的
有较好的光滑性，以达到去噪的目的。注意到上述目标函数是可微的，我们利用结合BB步长和非精确搜索的梯度下降对其进行求解。
图片和参数准备
设定随机种子。
clear;
seed = 97006855;
ss = RandStream('mt19937ar','Seed',seed); % 使用指定的伪随机数生成器算法创建随机数流
RandStream.setGlobalStream(ss); % 设置全局随机数流
载入未加噪的原图作为参考，记录为 u0。
u = load ('tower.mat'); % 加载图片数据
u = u.B1;
u = double(u);
[m,n] = size(u);
u0 = u;
生成加噪的图片，噪声 的每个元素服从独立的高斯分布
，并对每个像素进行归一化处理（将像素值转化到[0,1]区间内）。注意到 MATLAB 的 imshow函数（当第二个参数设定为空矩阵时），能够自动将矩阵中最小的元素对应到黑色，将最大的元素对应为白色。
u = u + 20*randn(m,n);
maxu = max(u(:)); minu = min(u(:));
u = (u - minu)/(maxu - minu); % 归一化
参数设定，以一个结构体提供各参数，分别表示 ，梯度和函数值的停机标准，输出的详细程度，和最大迭代次数。
opts = struct();
opts.xtol = 1e-8; %
opts.gtol = 1e-6; %
opts.ftol = 1e-16; %
opts.record = 0; %
opts.maxit = 200; % 最大迭代次数
求解正则化优化问题
分别取正则化系数为 和
，利用带BB 步长的梯度下降求解对应的优化问题，见<fminGBB.html 带BB步长线搜索的梯度法> 。
lambda = 0.5;
fun = @(x) TV(x,u,lambda);
[x1,~,out1] = fminGBB(u,fun,opts);
lambda = 2;
fun = @(x) TV(x,u,lambda);
[x2,~,out2] = fminGBB(u,fun,opts);
结果可视化，将不同正则化系数的去噪结果以图片形式展示。
subplot(2,2,1);
imshow(u0,[]);
title('原图')
subplot(2,2,2);
imshow(u,[]);
title('高斯噪声')
subplot(2,2,3);
imshow(x1,[]);
title('\lambda = 0.5')
subplot(2,2,4);
imshow(x2,[]);
title('\lambda = 2')
print(gcf,'-depsc','tv.eps')
Tikhonov 正则化模型的目标函数值和梯度计算
该无约束优化问题的目标函数为：
function [f,g] = TV(x,y,lambda)
分别表示带噪声图片和正则化参数， f，g 表示在 x 点处的目标函数值和梯度。
第一项 用于控制去噪后的图片
和带噪声的图片
之间的距离。
f = .5*norm(x - y, 'fro')^2; % F范数
计算两个方向上的离散差分：,
。
离散的拉普拉斯算子 d2x : 。
ip1 = min(i+1,m); jp1 = min(j+1,n);
im1 = max(i-1,1); jm1 = max(j-1,1);
[m,n] = size(y); % 实则为u
dx = zeros(m,n);
dy = zeros(m,n);
d2x = zeros(m,n);
for i = 1:m
for j = 1:n
% 为了确保i、j不超过范围，保证求解
ip1 = min(i+1,m); jp1 = min(j+1,n);
im1 = max(i-1,1); jm1 = max(j-1,1);
dx(i,j) = x(ip1,j) - x(i,j);
dy(i,j) = x(i,jp1) - x(i,j);
d2x(i,j) = x(ip1,j) + x(im1,j) + x(i,jp1) + x(i,jm1) - 4*x(i,j);
end
end
计算目标函数的第二项（Tikhonov 正则化）并与第一项合并得到当前点处的目标函数值
f = f + lambda * (norm(dx,'fro')^2 + norm(dy,'fro')^2);
目标函数的梯度可以解析地写出：
g = x - y - 2*lambda*d2x;
end
结果分析
首先针对图片去噪的效果进行分析。我们发现利用 Tikhonov 正则化模型可以有效地去除图片中的噪声。
当正则化系数 增大时，去噪的效果逐渐增强，但是图片中的物体边界也逐渐模糊。
同时我们也对带BB 步长的梯度下降法在其中的表现进行分析：在这两个问题中 BB
步长的梯度下降法都以非常迅速地速度收敛到了最优值。当最终收敛时，我们看到梯度的范数 |nrmG|
已经很小，这表明算法有较好的收敛性。同时注意到，虽然我们采用了回退法的线搜索方法，
但是在上面的应用中 BB 步长总是自然地满足了线搜索准则的要求，因此没有进行额外的步长衰减
（每一步的步长试探次数 |ls-Iter| 均为1）。
最优化：建模、算法与理论（刘浩洋、户将、李勇锋、文再文编著）
文章出处登录后可见！
