从Gauss-Newton算法到 LM算法 （详细推导及MATLAB实现、多自变量问题）

已于 2022-02-07 10:25:20 修改
阅读量1.2w 收藏 281
点赞数 51
分类专栏： 算法 文章标签： matlab
于 2019-12-18 22:47:35 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_41869763/article/details/103603089
版权

目录

一、结果回顾

Gauss-Newton算法（高斯牛顿法）用来估计参数，Gauss-Newton算法学习对Gauss-newton算法做了详细的解释，并且使用C++实现，写得很不错。

但是程序其实有微小错误，实际的坐标并不是年代1815-1885,而是1~8，否则 p = A e B t p = Ae^{Bt} p=AeBt拟合时将会迅速增大，也得不到 A = 0.7 A=0.7 A=0.7， B = 0.26 B=0.26 B=0.26 这一结果。 个人觉得原文对程序实现的解释也稍有不足，本文将使用 MATLAB，给出详细实现步骤。

A = 0.7 ， B = 0.26 ； t = 1 − 8 ， y = A e B t A = 0.7，B=0.26；t = 1-8，y = Ae^{Bt} A=0.7，B=0.26；t=1−8，y=AeBt绘图如下（星号为原始数据点，实现为拟合点）

二、Gauss-Newton算法

以下详细说明各个矩阵如何获取，逐步给出计算结果。

2.1 原始数据
时间刻度 t t t12345678
人口数 y y y8.311.014.719.726.735.244.455.9
2.2 数学模型

现需要使用指数模型拟合人口数 p 关于时间刻度 t 的变化。

数学模型为：
y ^ = x 1 e x 2 t (2.1) \hat {y} = x_1 e^{x_2 t} \tag{2.1} y^​=x1​ex2​t(2.1)

现在需要求解合适的参数 x = [ x 1 , x 2 ] x=[x_1, x_2] x=[x1​,x2​]，使得有：
min ⁡ f ( x ) = 1 2 ∑ i = 1 8 r i 2 ( x ) (2.2) \min f(x)=\frac{1}{2} \sum_{i=1}^8 r_i^2(x) \tag{2.2} minf(x)=21​i=1∑8​ri2​(x)(2.2)

r i ( x ) = y i ^ − y i (2.3) r_i(x) = \hat {y_i} - y_i \tag{2.3} ri​(x)=yi​^​−yi​(2.3)

2.3 算法流程

接下来给出 x = [ x 1 , x 2 ] T x=[x_1, x_2]^T x=[x1​,x2​]T 的迭代过程，也即Gauss-Newton算法流程：
1）给定初始点 x 0 x_0 x0​，允许误差 ϵ > 0 \epsilon>0 ϵ>0，置 k = 0 k=0 k=0；
2） f ′ ( x ( k ) ) < ϵ f'(x^{(k)}) < \epsilon f′(x(k))<ϵ，导数接近0即在极小点附近，得到 x x x，算法结束，否则继续；
3）计算 x ( k + 1 ) x^{(k+1)} x(k+1)， x ( k + 1 ) = x k − H − 1 x^{(k+1)}=x^{k}-H^{-1} x(k+1)=xk−H−1 ▽ \bigtriangledown ▽ f f f

2.4 梯度和黑森矩阵

根据定义，可得梯度和黑森矩阵如下：

原函数：

f ( x 1 , x 2 ) = 1 2 ∑ i = 1 8 r i 2 ( x 1 , x 2 ) (2.4) f(x_1, x_2) = \frac{1}{2} \sum_{i=1}^8 r_i^2(x_1, x_2) \tag{2.4} f(x1​,x2​)=21​i=1∑8​ri2​(x1​,x2​)(2.4)

梯度：
▽ f = [ ∂ f ∂ x 1 ∂ f ∂ x 2 ] T = ∑ i = 1 8 [ r i ∂ r i ∂ x 1 r i ∂ r i ∂ x 2 ] T (2.5) \bigtriangledown f= \left [ \frac {\partial f} {\partial x_1}\quad \frac {\partial f} {\partial x_2} \right ]^T = \sum_{i=1}^8 \left [ r_i \frac {\partial r_i} {\partial x_1} \quad r_i \frac {\partial r_i} {\partial x_2} \right ]^T \tag{2.5} ▽f=[∂x1​∂f​∂x2​∂f​]T=i=1∑8​[ri​∂x1​∂ri​​ri​∂x2​∂ri​​]T(2.5)

黑森矩阵：

H = [ ∂ 2 f ∂ x 1 2 ∂ 2 f ∂ x 1 x 2 ∂ 2 f ∂ x 2 x 1 ∂ 2 f ∂ x 2 2 ] ≈ ∑ i = 1 8 [ ∂ r i ∂ x 1 ∂ r i ∂ x 1 ∂ r i ∂ x 1 ∂ r i ∂ x 2 ∂ r i ∂ x 2 ∂ r i ∂ x 1 ∂ r i ∂ x 2 ∂ r i ∂ x 2 ] (2.6) H = \left [ \begin {array}{cc} \displaystyle \frac {\partial^2f} {\partial x_1^2} & \displaystyle \frac {\partial^2f} {\partial x_1x_2}\\ \displaystyle \frac {\partial ^2f} {\partial x_2x_1} & \displaystyle \frac {\partial ^2f} {\partial x_2^2} \end {array} \right ] \approx \sum_{i=1}^8 \left [ \begin{array}{cc} \displaystyle \frac {\partial r_i}{\partial x_1} \displaystyle \frac {\partial r_i}{\partial x_1} & \displaystyle \frac {\partial r_i}{\partial x_1} \displaystyle \frac {\partial r_i}{\partial x_2} \\ \displaystyle \frac {\partial r_i}{\partial x_2} \displaystyle \frac {\partial r_i}{\partial x_1} & \displaystyle \frac {\partial r_i}{\partial x_2} \displaystyle \frac {\partial r_i}{\partial x_2} \end {array} \right ] \tag{2.6} H=⎣⎢⎢⎡​∂x12​∂2f​∂x2​x1​∂2f​​∂x1​x2​∂2f​∂x22​∂2f​​⎦⎥⎥⎤​≈i=1∑8​⎣⎢⎡​∂x1​∂ri​​∂x1​∂ri​​∂x2​∂ri​​∂x1​∂ri​​​∂x1​∂ri​​∂x2​∂ri​​∂x2​∂ri​​∂x2​∂ri​​​⎦⎥⎤​(2.6)

现在只剩下 r r r 关于 x 1 、 x 2 x_1、x_2 x1​、x2​ 的偏导，由 r r r 的定义，容易直接根据求导法则得出 r r r 关于 x 1 、 x 2 x_1、 x_2 x1​、x2​ 的偏导（也可使用数值微分，见第2.5小节），如下：

2.5 MATLAB 程序

进行运算时，C++程序需要逐个数据读取，而MATLAB擅长矩阵运算。矩阵运算一方面使得效率提高，运算简洁，但一定程度上也对运算表达能力有更高要求。

%原始数据
t = 1:8;
popu = [8.3 11.0 14.7 19.7 26.7 35.2 44.4 55.9];

%高斯牛顿法
A = 1.0;                    %待求变量1
B = 1.0;                    %待求变量2
k = 0;                      %迭代次数
eps = 1e-3;                 %允许误差
r = A*exp(B*t) - popu;      %误差
r1 = zeros(size(r));        %记录上一次误差

while(1)
    if (sum(abs(r1-r)) < eps) || (k >100)
        A 
        B
        k
        break
    else
        pA = exp(B*t);                      %关于A的偏导数
        pB = A*exp(B*t).*t;                 %关于B的偏导数
        H = [2*sum(pA.*pA) 2*sum(pA.*pB)    
             2*sum(pA.*pB) 2*sum(pB.*pB)];  %黑森矩阵
        df = [sum(r.*pA)  
              sum(r.*pB)];                %梯度
          
        k = k + 1;
        delta = inv(H) * df;                %变化量
        A = A - delta(1);                   %更新待求参数A
        B = B - delta(2);                   %更新待求参数B
        
        r1 = r;                             %记录上一次误差
        r = A*exp(B*t) - popu;              %更新误差
    end
end

运行结果为

A = 7.0001
B = 0.2621
k = 25

因此迭代13次达到满足条件，也与文章最开始的结果一致。

2.6 微分的数值运算

在很多实际问题中，往往需要计算微分、偏微分等，本例就运用了偏微分。但是并不是所有的数学模型都像本例这样能够给出微分的解析式，很多时候需要数值运算。
回到第一次接触导数时，导数的定义为

f ′ ( x ) = lim ⁡ Δ x → 0 f ( x + Δ x ) − f ( x − Δ x ) 2 Δ x f'(x)=\lim\limits_{\Delta x \to 0} \frac {f(x+\Delta x) - f(x-\Delta x)} {2 \Delta x} f′(x)=Δx→0lim​2Δxf(x+Δx)−f(x−Δx)​

这里两边都取 Δ \Delta Δ x x x，有助于提高精度，偏微分与之类似。
本例中

r i = x 1 e x 2 t i − y i r_i = x_1 e^{x_2 t_i} - y_i ri​=x1​ex2​ti​−yi​

偏导即为：

∂ r i ∂ x 1 = r ( x 1 + Δ x 1 , x 2 ) − r ( x 1 − Δ x 1 , x 2 ) 2 Δ x 1 = ( x 1 + Δ x 1 ) e x 2 t i − ( x 1 − Δ x 1 ) e x 2 t i 2 Δ x 1 \frac {\partial r_i}{\partial x_1} = \frac {r(x_1+\Delta x_1, x_2) - r(x_1-\Delta x_1, x_2)} {2\Delta x_1} = \frac { (x_1+\Delta x_1)e^{x_2t_i} - (x_1-\Delta x_1) e^{x_2t_i}} {2\Delta x_1} ∂x1​∂ri​​=2Δx1​r(x1​+Δx1​,x2​)−r(x1​−Δx1​,x2​)​=2Δx1​(x1​+Δx1​)ex2​ti​−(x1​−Δx1​)ex2​ti​​

∂ r i ∂ x 2 = r ( x 1 , x 2 + Δ x 2 ) − r ( x 1 , x 2 − Δ x 2 ) 2 Δ x 2 = x 1 e ( x 2 + Δ x 2 ) t i − x 1 e ( x 2 − Δ x 2 ) t i 2 Δ x 2 \frac {\partial r_i}{\partial x_2} = \frac {r(x_1, x_2 + \Delta x_2) - r(x_1, x_2-\Delta x_2)} {2\Delta x_2} = \frac {x_1 e^{(x_2 + \Delta x_2) t_i} - x_1 e^{(x_2 - \Delta x_2)t_i}} {2\Delta x_2} ∂x2​∂ri​​=2Δx2​r(x1​,x2​+Δx2​)−r(x1​,x2​−Δx2​)​=2Δx2​x1​e(x2​+Δx2​)ti​−x1​e(x2​−Δx2​)ti​​

再取 Δ x 1 = Δ x 2 = 0.0001 \Delta x_1 = \Delta x_2 =0.0001 Δx1​=Δx2​=0.0001 例如，一个很小的值，我们就能够得出关于 x 1 、 x 2 x_1、x_2 x1​、x2​的偏导数。值得注意的是，我们其实也没必要过多地化简，因为很多时候是很复杂的，例如此时 x 2 x_2 x2​ 的偏导数就比较难化简，何不充分发挥计算机的计算能力呢？

修改程序如下，程序中使用 A A A 代替 x 1 x_1 x1​， B B B 代替 x 2 x_2 x2​：

%原始数据
t = 1:8;
popu = [8.3 11.0 14.7 19.7 26.7 35.2 44.4 55.9];

%高斯牛顿法
A = 1.0;                    %待求变量1
B = 1.0;                    %待求变量2
k = 0;                      %迭代次数
eps = 1e-3;                 %允许误差
r = A*exp(B*t) - popu;      %误差
r1 = zeros(size(r));        %记录上一次误差
dA = 1e-4;                  %求导时变化率
dB = 1e-4;                  %求导时变化率
while(1)
    if (sum(abs(r1-r)) < eps) || (k >100)
        A 
        B
        k
        break
    else
        pA = (((A+dA)*exp(B*t)) - (A-dA)*exp(B*t)) / (2*dA);%关于A的偏导数
        pB = (A*exp((B+dB)*t) - A*exp((B-dB)*t)) / (2*dB) ; %关于B的偏导数
        
        H = [2*sum(pA.*pA) 2*sum(pA.*pB)    
             2*sum(pA.*pB) 2*sum(pB.*pB)];  %黑森矩阵
        df = [sum(r.*pA)  
              sum(r.*pB)];                %梯度
          
        k = k + 1;
        delta = inv(H) * df;                %变化量
        A = A - delta(1);                   %更新待求参数A
        B = B - delta(2);                   %更新待求参数B
        
        r1 = r;                             %记录上一次误差
        r = A*exp(B*t) - popu;              %更新误差
    end
end

代码主要修改了pA，pB的计算方式，结果跟上面一致

A = 7.0001
B = 0.2621
k = 25

三、LM算法

3.1 算法原理

LM算法结合高斯牛顿法和梯度下降法，虽然听起来很复杂的样子，实际上很简单，就是在高斯牛顿法的基础上添加一个梯度下降因子，算法如下：

1）给定初始点 x 0 x_0 x0​，允许误差 ϵ > 0 \epsilon>0 ϵ>0，置 k = 0 k=0 k=0；
2） f ′ ( x ( k ) ) < ϵ f'(x^{(k)}) < \epsilon f′(x(k))<ϵ，导数接近0即在极小点附近，得到 x x x，算法结束，否则继续；
3）计算 x ( k + 1 ) x^{(k+1)} x(k+1)， x ( k + 1 ) = x k − ( H + x^{(k+1)}=x^{k}-(H+ x(k+1)=xk−(H+ α \alpha α I ) − 1 I)^{-1} I)−1 ▽ \bigtriangledown ▽ f f f （高斯牛顿法为 x ( k + 1 ) = x k − H − 1 x^{(k+1)}=x^{k}-H^{-1} x(k+1)=xk−H−1 ▽ \bigtriangledown ▽ f f f）

3.2 MATLAB程序
% 原始数据
t = 1:8;
popu = [8.3 11.0 14.7 19.7 26.7 35.2 44.4 55.9];

% LM法
A = 1.0;                    %待求变量1
B = 1.0;                    %待求变量2
k = 0;                      %迭代次数
eps = 1e-3;                 %允许误差
r = A*exp(B*t) - popu;      %误差
r1 = zeros(size(r));        %记录上一次误差
dA = 1e-4;                  %求导时变化率
dB = 1e-4;                  %求导时变化率
alpha = 10;					%LM法梯度下降因子 =====此行添加=====
while(1)
    if (sum(abs(r1-r)) < eps) || (k >100)
        A 
        B
        k
        break
    else
        pA = (((A+dA)*exp(B*t)) - (A-dA)*exp(B*t)) / (2*dA);%关于A的偏导数
        pB = (A*exp((B+dB)*t) - A*exp((B-dB)*t)) / (2*dB) ; %关于B的偏导数
        
        H = [2*sum(pA.*pA) 2*sum(pA.*pB)    
             2*sum(pA.*pB) 2*sum(pB.*pB)];  %黑森矩阵
        df = [sum(r.*pA)  
              sum(r.*pB)];                 	%梯度
          
        k = k + 1;
        delta = (inv(H + alpha*eye(2))) * df;	%变化量 =====此行修改=====
        A = A - delta(1);                   %更新待求参数A
        B = B - delta(2);                   %更新待求参数B
        
        r1 = r;                             %记录上一次误差
        r = A*exp(B*t) - popu;              %更新误差
    end
end

A = 6.9999
B = 0.2621
k = 33

结果基本一致，但比高斯牛顿法还多运行了几代。

四、多变量问题

先说明一下，多变量时结果并不如预期所想，存在一些问题！

由于其他博友问到多自变量时怎么办这个问题，思考后觉得可行，做此补充。方便起见，直接把上面的时间换成两个吧：

时间刻度 t 1 t_1 t1​12345678
时间刻度 t 2 t_2 t2​1.11.21.31.41.51.61.71.8
人口数 y y y8.311.014.719.726.735.244.455.9

数学模型为：
y ^ = x 1 e x 2 t 1 + x 3 e x 4 t 2 \hat {y} = x_1 e^{x_2 t_1} + x_3e^{x_4 t_2} y^​=x1​ex2​t1​+x3​ex4​t2​

现在需要求解合适的参数 x = [ x 1 , x 2 , x 3 , x 4 ] x=[x_1, x_2, x_3, x_4] x=[x1​,x2​,x3​,x4​]，使得有：
min ⁡ f ( x ) = 1 2 ∑ i = 1 8 r i 2 ( x ) r i ( x ) = y i ^ − y i \min f(x)=\frac{1}{2} \sum_{i=1}^8 r_i^2(x) \\ r_i(x) = \hat {y_i} - y_i minf(x)=21​i=1∑8​ri2​(x)ri​(x)=yi​^​−yi​

由于梯度和雅可比矩阵都使用一阶偏导即可计算，因此需要分别计算： ∂ y ^ ∂ x 1 \displaystyle \frac {\partial \hat y}{\partial x_1} ∂x1​∂y^​​， ∂ y ^ ∂ x 2 \displaystyle \frac {\partial \hat y}{\partial x_2} ∂x2​∂y^​​， ∂ y ^ ∂ x 3 \displaystyle \frac {\partial \hat y}{\partial x_3} ∂x3​∂y^​​， ∂ y ^ ∂ x 4 \displaystyle \frac {\partial \hat y}{\partial x_4} ∂x4​∂y^​​ 。此处导数可以直接运用公式计算，也可利用数值方式计算，程序将采用数值方式计算。

类比式(2.5)，梯度为：

▽ f = [ ∂ f ∂ x 1 ∂ f ∂ x 2 ∂ f ∂ x 3 ∂ f ∂ x 4 ] T = ∑ i = 1 8 [ r i ∂ r i ∂ x 1 r i ∂ r i ∂ x 2 r i ∂ r i ∂ x 3 r i ∂ r i ∂ x 4 ] T (4.1) \bigtriangledown f= \left [ \frac {\partial f} {\partial x_1}\quad \frac {\partial f} {\partial x_2} \quad \frac {\partial f} {\partial x_3} \quad \frac {\partial f} {\partial x_4} \right]^T = \sum_{i=1}^8 \left [ r_i \frac {\partial r_i} {\partial x_1} \quad r_i \frac {\partial r_i} {\partial x_2} \quad r_i \frac {\partial r_i} {\partial x_3} \quad r_i \frac {\partial r_i} {\partial x_4}\right ]^T \tag{4.1} ▽f=[∂x1​∂f​∂x2​∂f​∂x3​∂f​∂x4​∂f​]T=i=1∑8​[ri​∂x1​∂ri​​ri​∂x2​∂ri​​ri​∂x3​∂ri​​ri​∂x4​∂ri​​]T(4.1)

类比式(2.6)，黑森矩阵为：

H = [ ∂ 2 f x 1 2 ⋯ ∂ 2 f ∂ x 1 x 4 ⋮ ⋱ ⋮ ∂ 2 f x 4 x 1 ⋯ ∂ 2 f ∂ x 4 2 ] ≈ ∑ i = 1 8 [ ∂ r i ∂ x 1 ∂ r i ∂ x 1 ⋯ ∂ r i ∂ x 1 ∂ r i ∂ x 4 ⋮ ⋱ ⋮ ∂ r i ∂ x 4 ⋯ ∂ r i ∂ x 4 ∂ r i ∂ x 4 ] (2.6) H = \left [ \begin {array}{ccc} \displaystyle \frac {\partial^2f} {x_1^2}& \cdots &\displaystyle \frac {\partial^2f} {\partial x_1x_4}\\ \vdots & \ddots & \vdots \\ \displaystyle \frac {\partial^2f} {x_4x_1}& \cdots &\displaystyle \frac {\partial^2f} {\partial x_4^2} \end {array} \right ] \approx \sum_{i=1}^8 \left [ \begin{array}{ccc} \displaystyle \frac {\partial r_i}{\partial x_1} \frac {\partial r_i}{\partial x_1} & \cdots & \displaystyle \frac {\partial r_i}{\partial x_1} \frac {\partial r_i}{\partial x_4} \\ \vdots & \ddots & \vdots \\ \displaystyle \frac {\partial r_i}{\partial x_4} & \cdots & \displaystyle \frac {\partial r_i}{\partial x_4} \displaystyle \frac {\partial r_i}{\partial x_4} \end {array} \right ] \tag{2.6} H=⎣⎢⎢⎢⎢⎢⎡​x12​∂2f​⋮x4​x1​∂2f​​⋯⋱⋯​∂x1​x4​∂2f​⋮∂x42​∂2f​​⎦⎥⎥⎥⎥⎥⎤​≈i=1∑8​⎣⎢⎢⎢⎢⎡​∂x1​∂ri​​∂x1​∂ri​​⋮∂x4​∂ri​​​⋯⋱⋯​∂x1​∂ri​​∂x4​∂ri​​⋮∂x4​∂ri​​∂x4​∂ri​​​⎦⎥⎥⎥⎥⎤​(2.6)

公式有了，写程序，程序中使用 A , B , C , D A, B, C, D A,B,C,D 代替 x 1 , x 2 , x 3 , x 4 x_1, x_2, x_3, x_4 x1​,x2​,x3​,x4​ ：

%原始数据
t1 = 1:8;
t2 = 1.1:0.1:1.8;
popu = [8.3 11.0 14.7 19.7 26.7 35.2 44.4 55.9];

%高斯牛顿法
A = 1.0;                    %待求变量1
B = 1.0;                    %待求变量2
C = 1.0;
D = 1.0;

k = 0;                      %迭代次数
eps = 1e-3;                 %允许误差（相比之前，降低了一个数量级精度）
r = A*exp(B*t1)+ C*exp(D*t2) - popu;      %误差
r1 = zeros(size(r));        %记录上一次误差
dA = 1e-4;                  %求导时变化率
dB = 1e-4;                  %求导时变化率
dC = 1e-4;
dD = 1e-4;

alpha = 5;                 % 梯度下降因子

while(1)
    if (sum(abs(r1-r)) < eps) || (k >50)
        A 
        B
        C
        D
        k
        break
    else
        % 计算偏导
        pA = (((A+dA)*exp(B*t1)) - (A-dA)*exp(B*t1)) / (2*dA);%关于A的偏导数
        pB = (A*exp((B+dB)*t1) - A*exp((B-dB)*t1)) / (2*dB) ; %关于B的偏导数
        pC = (((C+dC)*exp(D*t2)) - (C-dC)*exp(D*t2)) / (2*dC);%关于C的偏导数
        pD = (C*exp((D+dD)*t2) - C*exp((D-dD)*t2)) / (2*dD) ; %关于D的偏导数
        
        % 黑森矩阵
        H = [2*sum(pA.*pA) 2*sum(pA.*pB) 2*sum(pA.*pC)  2*sum(pA.*pD)   
             2*sum(pB.*pA) 2*sum(pB.*pB) 2*sum(pB.*pC)  2*sum(pB.*pD)
             2*sum(pC.*pA) 2*sum(pC.*pB) 2*sum(pC.*pC)  2*sum(pC.*pD)
             2*sum(pD.*pA) 2*sum(pD.*pB) 2*sum(pD.*pC)  2*sum(pD.*pD)];  
         
        % 梯度
        df = [sum(r.*pA)  
              sum(r.*pB)
              sum(r.*pC)
              sum(r.*pD)];                
          
        k = k + 1;
        delta =  (inv(H + alpha*eye(4))) * df;                %变化量
        A = A - delta(1);                   %更新待求参数A
        B = B - delta(2);                   %更新待求参数B
        C = C - delta(3);
        D = D - delta(4);
        
        r1 = r;                             %记录上一次误差
        r = A*exp(B*t1)+ C*exp(D*t2) - popu;              %更新误差
    end
end

结果如下：

A = -0.2168
B = 0.5926
C = 0.2473
D = 3.2152
k = 51

注意，不使用LM算法只使用Gauss-Newton法会出现矩阵奇异的问题，难以继续。

总结

本文在Gauss-Newton算法学习的基础上给出具体公式及推导，使用 MATLAB 详细地实现Gauss-Newton算法和LM算法，并给出数值计算偏导数方法。

大强强小强强
关注
  • 51
    点赞
  • 踩
  • 281
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 25
    评论
专栏目录
matlab编写的LM优化算法
05-23
使用matlab编写的LM优化算法，基于imm3215，不是matlab内部函数。FK.m是原函数，JFK是一阶Jaccobi导数
lm算法的matlab程序.rar
04-29
lm算法的matlab程序 解压后即可，还不错
lm算法的matlab实现
12-14
该算法是自己编的一个lm算法，是基于MATLAB的程序
机器学习训练算法十(列文伯格-马夸尔特法(LM 法))
黎明鱼儿
01-09 3365
本文详细介绍了列文伯格-马夸尔特算法（LM算法），这是一种用于解决非线性最小二乘问题的优化算法。文章从算法的介绍、数学原理、阻尼因子更新策略、列文伯格方法与马夸尔特方法的区别，到Matlab程序实现，全面阐述了LM算法的理论和应用。LM算法结合了高斯牛顿法的快速局部收敛性和梯度下降法的全局搜索能力，通过引入阻尼因子实现算法的鲁棒性。
【Matlab】LM迭代估计法
最新发布
咸鱼漫谈的博客
06-18 2479
其中 y(i) 为第 i 组数据的真实值或理论真实值，fi(θ) 为第 i 组数据的预测值，目的就是找到一组 θ，使得代价函数达到最小值。7.更改估计参数之后，再次使用 n 组采样数据，重复 6 步骤，如果代价函数值变化小于设定的迭代精度，或达到最大迭代次数，则停止迭代，返回跌代后的参数向量 θ。6.迭代更新参数向量，雅可比矩阵(有一些简单的情况雅可比矩阵只与采样值有关，就可以只计算一次)，这里将收敛补偿设置为1，即α = 1。迭代次数为 1000 次，采样 100 组 x，y 数据。
LM(Levenberg-Marquadrdt )算法在MATLAB中的实现及实例
热门推荐
waitingwinter的博客
05-15 1万+
LM算法简介 考虑如下非线性最小二乘问题min⁡f(x)=12∥r(x)∥2, \min \quad f(x) = \frac{1}{2}\Vert r(x)\Vert^2, minf(x)=21​∥r(x)∥2, 其中残差向量 r:Rn→Rr:\mathbb{R}^n \to \mathbb{R}r:Rn→R 为 r(x)=(r1(x),r2(x),⋯ ,rm(x))T.r(x) = \big(r_1(x),r_2(x),\cdots,r_m(x)\big)^T. r(x)=(r1​(x),r2​(x),
MATLAB--高斯牛顿法
老邹的博客
05-16 2600
利用高斯牛顿法求解非线性问题，包含了利用LDLT分解法求解Ax=b的方程: clear all; clc %%演示牛顿高斯法--非线性优化 %真实值 ar = 1.0; br = 2.0; cr = 1.0; x = 1:2:100; x = x./100; nosie = 0.9*randn(1,size(x,2));%产生高斯噪音 y_real = exp(ar*x.^2+br*x+cr)+nosie; %估计初始值 ae = 2.0; be = 4.0; ce = 3.0; itr = 20;%迭
LM算法-Matlab/Opencv
11-25
LM算法介绍：牛顿法需要时刻计算H矩阵，即二阶导数信息，是很麻烦的一件事情，LM算法的核心思想就是用雅可比矩阵（易计算）代替H矩阵的计算，使得优化效率有了提升，文档包括了Matlab/Opencv实现案例
LM拟合优化算法的MATLAB实现
12-11
Levenberg-Marquardt Method线性拟合优化算法的实现。给予MATLAB平台，.m文件
算法_lm算法_lm方程组_matlab_小波方程_python_
10-03
几种求解非线性方程组的算法，如Legendre小波函数，LM算法等
matlablm算法代码-OptimisationSA-LM:编写程序来演示模拟退火和LevenbergMarquardt的工作方式
05-26
matlab lm算法代码这个小项目的完成主要是为了更好地了解Levenberg Marquardt和Simulated Annealing两者如何工作，并且他们也变得更加精通C＃.NET。 它们是MATLAB代码的改编。 对于模拟退火： 对于Levenberg-Marquardt：第7页的lm.m 对于LM Matrix操作，使用Math.Net.Numerics使事情变得更快。 但是，稍后写出矩阵求逆和乘法将很有趣。 Program.cs执行以下操作： -使用已知参数实现功能 -采样 -向这些样本添加随机噪声 -使用这两种方案尝试再次查找参数。 不幸的是，由于需要改进内存管理，结果不能令人满意。 该代码仍可以（并且将会）得到改进。 但是，从第一个出版物起，我很高兴能够理解每种算法的步骤。
LM算法的讲解
04-29
LM法概论
LM算法+PPT讲解
06-04
MATLAB中LM算法的实现+PPT讲解，希望对大家有帮助！LM算法非线性最小二乘法
LM-BP神经网络的Matlab程序源代码
10-24
LM-BP神经网络的源代码，有注释，容易读懂
matlablm算法代码-Levenberg-Marquardt-algorithm:在此存储库中，我将通过重写代码并手动计算Jacobian
05-26
matlab lm算法代码
gauss-newton:Gauss-Newton算法的简单Python实现
05-12
高斯牛顿 ... datasets.py - 非线性回归问题。 gaussnewton.py简单的非线性最小二乘问题求解器。 graph.py图形生成脚本。 img/ -由graph.py生成的graph.py 。 要求 Python 2.7 NumPy 象征 Matplotlib
Gauss-FFT算法：用于地球物理重磁场建模的Gauss-FFT算法-matlab开发
05-29
在MATLAB环境中实现Gauss-FFT算法，程序员可以利用其内置的FFT函数和其他数值工具，如矩阵运算和优化算法，来简化代码编写和提高计算效率。MATLAB的优势在于其友好的编程界面和丰富的科学计算库，使得研究人员和...
Image Filtering by Gauss-Seidel：使用 Gauss-Seidel 算法过滤图像的程序。-matlab开发
05-30
Gauss Seidel 算法包括边缘鲁棒性实现，使用 beta 变量实现。 该程序会拍摄一张图像，将高斯型噪声添加到其中，然后对其进行滤波，从而计算出相对于原始图像的均方误差，并计算出必须始终收敛的图像能量函数。
如何编写gauss-newton matlab
05-13
Gauss-Newton算法是一种非线性最小二乘问题的优化算法，可以用于拟合非线性模型。在Matlab中，可以使用以下步骤编写Gauss-Newton算法： 1. 定义非线性模型。例如，假设我们要拟合一个二次函数模型y=ax^2+bx+c，可以定义一个函数： function [y] = myfun(x,coeffs) y = coeffs(1)*x.^2 + coeffs(2)*x + coeffs(3); end 其中，x是自变量，coeffs是模型参数。 2. 定义观测数据。例如，我们有10个数据点，可以定义一个向量： x = linspace(0, 1, 10); % 自变量 y_data = 2*x.^2 + 0.5*x + 1 + 0.1*randn(1,10); % 观测数据 其中，y_data是带有噪声的观测数据。 3. 定义Gauss-Newton算法。可以使用以下代码： function [coeffs] = gauss_newton(x, y_data, coeffs_0) max_iter = 100; % 最大迭代次数 tol = 1e-6; % 收敛条件 coeffs = coeffs_0; for i = 1:max_iter % 计算残差 r = y_data - myfun(x, coeffs); % 计算雅可比矩阵 J = [x.^2, x, ones(size(x))]; % 计算增量 delta = (J'*J)\(J'*r'); % 更新参数 coeffs = coeffs + delta'; % 判断收敛 if norm(delta) < tol break end end end 其中，coeffs_0是初始参数值。 4. 调用Gauss-Newton算法并输出结果。可以使用以下代码： coeffs_0 = [1, 1, 1]; % 初始参数值 coeffs = gauss_newton(x, y_data, coeffs_0); % 调用Gauss-Newton算法 y_fit = myfun(x, coeffs); % 计算拟合结果 plot(x, y_data, 'o', x, y_fit, '-') % 绘制观测数据和拟合结果 注意，Gauss-Newton算法可能会陷入局部最优解，因此需要多次运行算法以获得更好的结果。
评论 25
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

大强强小强强

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值