接上一篇：https://www.guyuehome.com/44909
一、牛顿法（Newton-Raphson）
1.1 介绍
牛顿法（Newton’s method）是由数学家Isaac Newton提出的，也被称为二阶梯度法，通过使用目标函数的一阶和二阶导数信息来逐步逼近最优解，该方法在目标函数具有二阶连续导数的情况下收敛较快。
牛顿法的原理基于泰勒级数展开，使用局部线性逼近来逼近函数的根或极小值点
优点：
- 收敛速度快：牛顿法通常具有较快的收敛速度，特别是在目标函数局部具有良好的二次收敛性的情况下。它能够更快地逼近最优解，尤其对于复杂的非线性问题效果显著；
- 更精确的逼近：相比梯度下降法等一阶优化算法，牛顿法利用目标函数的局部二次近似来进行参数更新，因此更能精确地逼近最优解；
- 鲁棒性：在目标函数局部具有二次收敛性的情况下，牛顿法通常表现较稳健，能够较好地收敛到最优解。
不足：
- 计算复杂度高：牛顿法需要计算目标函数的一阶导数（梯度）和二阶导数（海森矩阵），特别是在参数较多的情况下，计算海森矩阵的代价较高，导致计算复杂度增加。
- 对初始点敏感：牛顿法对于初始参数的选择较为敏感，不同的初始点可能会导致不同的收敛结果，甚至可能发散。在应用中，需要进行仔细的初始点选择或采取一些初始化策略来提高稳定性。
- 需要海森矩阵可逆：牛顿法需要计算海森矩阵的逆，因此海森矩阵必须是可逆的。在某些情况下，海森矩阵可能不可逆，导致无法应用牛顿法。
1.2 牛顿法的基本步骤：
就是加入了海森矩阵计算
- 定义损失函数：首先，我们需要定义一个损失函数，它度量模型预测与实际观测之间的差异。在优化问题中，我们希望最小化损失函数。一般损失函数被定义为误差的最小二乘形式。
- 初始化参数：选择初始参数作为优化的起始点。
- 计算损失函数的梯度（Jacobian矩阵）和Hessian矩阵：计算损失函数关于参数的一阶导数（梯度）和二阶导数（Hessian矩阵）。梯度（Jacobian矩阵）告诉我们损失函数在当前参数点的变化方向和速率，而Hessian矩阵提供了更多关于目标函数曲率的信息。
- 更新参数：使用Hessian矩阵的逆乘以梯度向量来获得参数的改进方向（即delta）。然后将delta添加到当前参数，得到新的参数值。delta = - inv(H) * J; theta = theta + delta; 这一步是牛顿法与梯度下降法的主要区别，牛顿法使用了二阶导数信息，可以更快地收敛到最优解。
- 重复迭代：重复步骤2和步骤3，直到满足停止条件，例如达到一定的迭代次数或参数的变化足够小。
1.3 matlab例子
如下：
clear;clc;
%% 真实值
a = 0.5;
b = -0.5;
c = 1.0;
% 横坐标取值：1-100每隔两个数取一个数，并除以100
x = (1:2:100)./100;
% 产生待拟合数据
nosie = 0.005*randn(1,size(x,2)); % size(x,1) 指矩阵的行数；size(x,2) 就是矩阵的列数；
y = a*x.^2 + b*x + c + nosie;
% 显示
% scatter(x,y,100,'black','.');
%% 牛顿法（Newton-Raphson）、二阶梯度法
% step0:定义最大迭代次数
max_iter = 100;
% step1:定义损失函数，初始化loss值
loss_function = @(val) mean((val(1) * x.^2 + val(2) * x + val(3)) - y).^2;
loss = zeros(1, max_iter);
% step2:定义雅克比矩阵（一阶梯度）和海森矩阵（二阶梯度）
Jacobian_function = @(val) [mean(2 * x.^2 .* (val(1) * x.^2 + val(2) * x + val(3) - y)); % 返回一个列向量（与优化参数theta对应）
mean(2 * x .* (val(1) * x.^2 + val(2)* x + val(3) - y));
mean(2 * (val(1) * x.^2 + val(2) * x + val(3) - y))];
hessian_function = @(val) [ mean(2 * x.^4) mean(2 * x.^3) mean(2 * x.^2); % 返回一个矩阵（与优化参数theta对应）
mean(2 * x.^3) mean(2 * x.^2) mean(2 * x);
mean(2 * x.^2) mean(2 * x) 2];
% step3:初始化参数（theta）、收敛容差（epsilon）
theta = [0.35; -0.25; 1.5]; % 列向量
epsilon = 1e-5;
% 使用牛顿法进行优化
for iter = 1:max_iter
% step4:计算梯度和Hessian矩阵
gradient = Jacobian_function(theta);
hessian = hessian_function(theta);
% step5:更新参数
delta = - hessian\gradient; % inv(hessian) * gradient;
theta = theta + delta;
% step6:计算损失函数的值
loss(iter) = loss_function(theta);
% step7:根据损失函数值以及迭代次数判断是否退出拟合
if norm(Jacobian_function(theta)) < epsilon
fprintf('结束时迭代次数： %d, 结束梯度值： %.8f \n', iter, norm(Jacobian_function(theta)));
break;
end
% ...
end
% 输出拟合的参数和损失函数的值
fprintf('拟合的参数：a = %.4f, b = %.4f, c = %.4f\n', theta(1), theta(2), theta(3));
%% 绘制原始数据点和拟合的曲线
figure;
scatter(x, y, 100, 'black', '.');
hold on;
plot(x, a * x.^2 + b * x + c, 'r', 'LineWidth', 2);
xlabel('X');
ylabel('Y');
title('Newton-Raphson: Curve Fitting');
legend(' raw data', ' Fit curve');
grid on;
% 绘制loss变化曲线，这里为了方便观察只看前50个loss数据
figure;
plot(loss(1:50));
axis([-inf, inf,0,0.1])
title('Newton-Raphson: Loss');
grid on;
它只需要一步就拟合了
二、高斯牛顿法（Gauss-Newton method）
2.1 介绍
高斯牛顿法（Gauss-Newton method）是梯度下降法和牛顿法的一种折衷方法，综合了两者的优点和不足。
对比牛顿法， GN用误差函数f的 J^T.J 代替牛顿法中的二阶 Hession 矩阵计算。
注意这里，J和H不是损失函数的，是误差函数的
高斯牛顿法特别适用于拟合非线性模型到实际数据中。
参考https://zhuanlan.zhihu.com/p/42383070/?utm_id=0
优点：
- 收敛速度较快：高斯牛顿法相对于梯度下降法来说，通常具有较快的收敛速度。它利用目标函数的局部二次近似来进行参数更新，因此能够更快地逼近最优解；
- 不需要计算海森矩阵的逆：与传统的牛顿法相比，高斯牛顿法不需要计算目标函数的海森矩阵的逆，而是使用Jacobian矩阵的伪逆（或广义逆）J^T.J 来近似 H，从而减少了计算复杂度；
- 鲁棒性较好：相比于传统的牛顿法，高斯牛顿法对于初始参数的选择较为鲁棒，不太敏感，一般能够比较稳定地收敛到最优解。
不足：
- 可能发散：高斯牛顿法在目标函数局部具有二次收敛性的情况下收敛速度较快，但在目标函数的局部线性特性较差的地方，可能出现发散的情况；
- 对初始点的选择：虽然相对于传统的牛顿法鲁棒性较好，但高斯牛顿法仍然对初始参数的选择有一定的依赖，需要仔细选择初始点以保证收敛性。
- 受噪声影响：高斯牛顿法对于噪声敏感，如果目标函数受到较大的噪声影响，可能导致拟合结果不够稳定。
- 会出现 J^T.J 为奇异的情况，可能是奇异矩阵或者病态的，此时会导致稳定性很差，算法不收敛。（要求J列满秩）
2.2 高斯牛顿法的基本步骤：
用 J^T.J 替代海森矩阵的计算
需要注意这里的f不是损失函数，而是误差函数。
- 定义损失函数：首先，我们需要定义一个损失函数，它度量模型预测与实际观测之间的差异。在优化问题中，我们希望最小化损失函数。一般损失函数被定义为误差的最小二乘形式。
- 初始化参数：选择初始参数作为优化的起始点。
- 计算 Jacobian 矩阵：计算误差函数f关于参数 θ 的 Jacobian 矩阵 J。
- 计算 H = JT.J ，B = JT.f ，delta = - inv(H) * B，theta = theta + delta。
- 重复迭代：重复步骤2、步骤3、步骤4，直到满足停止条件，例如达到一定的迭代次数或参数的变化足够小。
2.3 matlab例子
如下：
clear; clc;
%% 真实值
a = 0.5;
b = -0.5;
c = 1.0;
% 横坐标取值：1-100每隔两个数取一个数，并除以100
x = (1:2:100)./100;
% 产生待拟合数据
noise = 0.005*randn(1, size(x, 2)); % size(x, 1) 指矩阵的行数；size(x, 2) 就是矩阵的列数；
y = a*x.^2 + b*x + c + noise;
%% 高斯-牛顿法（Gauss-Newton）
% step0:定义最大迭代次数
max_iter = 100;
% step1:定义损失函数，初始化loss值
loss_function = @(val) mean((val(1) * x.^2 + val(2) * x + val(3)) - y).^2;
loss = zeros(1, max_iter);
% step2:定义误差函数、误差函数雅克比矩阵
error_function = @(val) val(1) * x.^2 + val(2) * x + val(3) - y;
Jacobian_function = @(val) [x.^2; x; ones(size(x))]';
% step3:初始化参数（theta）、收敛容差（epsilon）
theta = [0.35; -0.25; 1.5]; % 列向量
epsilon = 1e-5;
% 使用高斯-牛顿法进行优化
for iter = 1:max_iter
% step4:计算雅克比矩阵
J = Jacobian_function(theta);
% step5:计算残差向量
residuals = error_function(theta);
% step6:计算近似的Hessian矩阵和梯度
hessian_approx = J' * J;
gradient = J' * residuals';
% step7:更新参数
delta = - hessian_approx \ gradient; % 等价于inv(hessian_approx) * gradient
theta = theta + delta;
% step8:计算损失函数的值
loss(iter) = loss_function(theta);
% step9:根据损失函数值以及迭代次数判断是否退出拟合
if norm(Jacobian_function(theta)' * error_function(theta)') < epsilon
fprintf('结束时迭代次数： %d, 结束梯度值： %.8f \n', iter, norm(gradient));
break;
end
end
% 输出拟合的参数和损失函数的值
fprintf('拟合的参数：a = %.4f, b = %.4f, c = %.4f\n', theta(1), theta(2), theta(3));
%% 绘制原始数据点和拟合的曲线
figure;
scatter(x, y, 100, 'black', '.');
hold on;
plot(x, theta(1) * x.^2 + theta(2) * x + theta(3), 'r', 'LineWidth', 2);
xlabel('X');
ylabel('Y');
title('Gauss-Newton: Curve Fitting');
legend('raw data', 'Fit curve');
grid on;
% 绘制loss变化曲线，这里为了方便观察只看前50个loss数据
figure;
plot(loss(1:50));
axis([-inf, inf, 0, 0.1])
title('Gauss-Newton: Loss');
grid on;
这里它也是一步就完成了优化
三、Levenberg-Marquardt
3.1 介绍
主要为了克服高斯牛顿法 J^T.J 为奇异的情况。
LM在求解过程中引入了阻尼因子，在一定程度上修正了高斯牛顿法的缺点，它比高斯牛顿法更加鲁棒，以牺牲一定的收敛速度为代价。
LM是被应用最广泛的无约束优化方法，曲线拟合、机器学习中的参数优化、图像配准等。
针对高斯牛顿的 J^T.J . △x = - J^T . f
改为：(J^T.J + μ.I) . △x = - J^T . f
其中μ被称为阻尼因子。实际是对误差函数被忽略的二次项以上进行补偿。
μ首先要大于0，以保证类似于海森矩阵需要正定的要求。
μ取比较大时，则△x ≈ - J^T . f / μ ≈ - F’^T /μ ，F’这里表示损失函数的导数，接近牛顿法
μ取比较小时，则表示依旧忽略误差函数的二次项以上，接近高斯牛顿法
这里μ的取值就非常关键，要限制其范围。
μ的初值选取，一般在 1e-8 到 1e-3 之间。
LM采用的搜索方法是信赖域（Trust Region）方法，和梯度下降、牛顿法、高斯牛顿法采用的线性搜索（Line Search）方法不同。
信赖域的主要思想是：固定搜索区域（防止 △x 过大），寻找区域内的最优点，在信赖区域里，认为近似是有效的，出了这个区域，近似会出问题。
在使用Levenberg-Marquart时，先设置一个比较小的μ值，当发现目标函数反而增大时，将μ增大使用梯度下降法快速寻找，然后再将μ减小使用牛顿法进行寻找。
为了合理更新阻尼因子μ，引入比例因子 ρ 等于（误差函数实际差值）/（误差函数近似差值）
若 ρ ≤0.25 ，说明步子迈得太大了，应缩小信赖域半径，则令 μ(t+1) = μ(t)/2
若 0.25 < ρ < 0.75 ，说明这一步迈出去之后，处于“可信赖”和“不可信赖”之间，可以维持当前的信赖域半径，令 μ(t+1) = μ(t)
若 ρ > 0.75 说明这一步近似比较准确，可以尝试扩大信赖域半径，则令 μ(t+1) = μ(t)*2
若 ρ <= 0 说明函数值是向着上升而非下降的趋势变化了，也就是优化错了方向，这时不应该朝这个方向更新参数了， 则 x(t+1) = x(t)。
若 ρ > 0 则更新参数 x(t+1) = x(t) + delta
（这里x为被优化参数，与我代码里的θ一个意思）
优点：
- 稳定性好：由于结合了梯度下降法和牛顿法，LM算法在靠近解时表现出牛顿法的快速收敛特性，而在离解较远时具有梯度下降法的稳定性。
- 自适应调整：通过调整阻尼因子μ，LM算法可以在不同阶段动态选择步长，使得算法在收敛过程中更加灵活。
- 高效：在处理非线性最小二乘问题时，LM算法能够有效地利用误差的二阶信息（近似海森矩阵），从而加快收敛速度。速度上：牛顿法 > 高斯牛顿法 > LM法 > 一阶梯度法
不足：
- 依赖初始值：算法的表现很大程度上依赖于初始参数值。若初始值选择不当，可能会导致收敛到局部最优解，或者收敛速度较慢。
- 计算量大：每次迭代需要计算雅克比矩阵和近似海森矩阵，计算量较大，对于大规模问题或者高维数据，计算代价较高。
- 参数调优复杂：需要设置阻尼因子μ的初始值和调整策略，如果选择不当，可能影响收敛速度和结果的准确性。
- 对噪声敏感：在数据中存在较大噪声的情况下，LM算法可能会出现过拟合，导致拟合效果不佳。
3.2 LM法的基本步骤：
大致与高斯牛顿法一致，主要在于调整阻尼因子μ
步骤中的部分参数，0.25、0.75、* 2、/2之类的，不一定严格，可以自己调整。
- 定义损失函数：首先，我们需要定义一个损失函数，它度量模型预测与实际观测之间的差异。在优化问题中，我们希望最小化损失函数。一般损失函数被定义为误差的最小二乘形式。
- 初始化参数：选择初始参数作为优化的起始点。
- 计算 Jacobian 矩阵：计算误差函数f关于参数 θ 的 Jacobian 矩阵 J。
- 计算 H = J^T.J + μ.I，B = JT.f ，delta = - inv(H) * B。
- 计算比例因子 ρ = (误差函数实际差值）/（误差函数近似差值）
- 调节阻尼因子，若 ρ ≤0.25 则 μ(t+1) = μ(t)/2；若 0.25 < ρ < 0.75 则 μ(t+1) = μ(t)；若 ρ > 0.75 则 μ(t+1) = μ(t)* 2
- 调节参数，若 ρ <= 0 则 θ(t+1) = θ(t)，若 ρ > 0 则 θ(t+1) = θ(t) + delta
- 重复迭代：重复步骤2、步骤3、步骤4，直到满足停止条件，例如达到一定的迭代次数或参数的变化足够小。
3.3 matlab例子
如下：
clear; clc;
%% 真实值
a = 0.5;
b = -0.5;
c = 1.0;
% 横坐标取值：1-100每隔两个数取一个数，并除以100
x = (1:2:100)./100;
% 产生待拟合数据
noise = 0.005*randn(1, size(x, 2)); % size(x, 1) 指矩阵的行数；size(x, 2) 就是矩阵的列数；
y = a*x.^2 + b*x + c + noise;
%% Levenberg-Marquardt方法
% step0:定义最大迭代次数
max_iter = 100;
% step1:定义损失函数，初始化loss值
loss_function = @(val) mean((val(1) * x.^2 + val(2) * x + val(3) - y).^2);
loss = zeros(1, max_iter);
% step2:定义误差函数和雅克比矩阵
error_function = @(val) val(1) * x.^2 + val(2) * x + val(3) - y;
Jacobian_function = @(val) [x.^2; x; ones(size(x))]';
% step3:初始化参数（theta）、收敛容差（epsilon）和阻尼因子（lambda）
theta = [0.35; -0.25; 1.5]; % 列向量
theta_new = theta;
epsilon = 1e-5; % 收敛容差
lambda = 1e-3; % 初始阻尼因子
Delta_max = 2; % 最大信赖域半径
Delta = 1; % 初始信赖域半径
% 使用Levenberg-Marquardt方法进行优化
for iter = 1:max_iter
% step4:计算雅克比矩阵
J = Jacobian_function(theta);
% step5:计算残差向量
residuals = error_function(theta);
% step6:计算LM更新
H = J' * J;
H_lm = H + lambda * eye(size(H));
gradient = J' * residuals';
delta = - H_lm \ gradient; %- inv(H_lm) * gradient
theta_new = theta + delta;
% % 信赖域调整策略
% 计算模型和实际误差的匹配程度
rho = (loss_function(theta) - loss_function(theta_new)) / (delta' * (lambda*delta - gradient) / 2);
% 判断是否接受新的参数值并调整信赖域半径和阻尼因子
if rho > 0
if rho > 0.75
lambda = lambda * 2;
elseif rho < 0.25
lambda = lambda / 2;
else
% lambda = lambda;
end
% step7:更新参数
theta = theta_new;
else
lambda = lambda / 2;
% step7:更新参数
% theta = theta;
end
% % % % %
% % % % 更简便方法
% % 判断是否接受新的参数值，调整阻尼参数
% if loss_function(theta) > loss_function(theta_new)
% lambda = lambda * 2;
% % step7:更新参数
% theta = theta_new;
% else
% lambda = lambda / 2;
% % step7:更新参数
% % theta = theta;
% end
% % % % % % %
% step8:计算损失函数的值
loss(iter) = loss_function(theta);
% step10:判断是否收敛
if norm(Jacobian_function(theta)' * error_function(theta)') < epsilon
fprintf('结束时迭代次数： %d, 结束梯度值： %.8f \n', iter, norm(delta));
break;
end
end
% 输出拟合的参数
fprintf('拟合的参数：a = %.4f, b = %.4f, c = %.4f\n', theta(1), theta(2), theta(3));
%% 绘制原始数据点和拟合的曲线
figure;
scatter(x, y, 100, 'black', '.');
hold on;
plot(x, theta(1) * x.^2 + theta(2) * x + theta(3), 'r', 'LineWidth', 2);
xlabel('X');
ylabel('Y');
title('Levenberg-Marquardt: Curve Fitting');
legend('raw data', 'Fit curve');
grid on;
% 绘制loss变化曲线，这里为了方便观察只看前50个loss数据
figure;
plot(loss(1:50));
axis([-inf, inf, 0, 0.1])
title('Levenberg-Marquardt: Loss');
grid on;
其比高斯牛顿慢一步
评论（0）
您还未登录，请登录后发表或查看评论