MATLAB中LSTM算法实例_bresenham直线算法

大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE稳定放心使用

一、结果回顾

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

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

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

二、Gauss-Newton算法

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

2.1 原始数据

2.2 数学模型

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

数学模型为：
$$\hat{y}=x_1{e}^{x_{2}t}\text{\hspace{1em}(2.1)}$$

现在需要求解合适的参数 $x=[x_1,x_2]$，使得有：
$$\min f(x)=\frac{1}{2}\sum _{i=1}^8{r}_i^2(x)\text{\hspace{1em}(2.2)}$$

$$r_i(x)=\widehat{y_i}-y_i\text{\hspace{1em}(2.3)}$$

2.3 算法流程

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

2.4 梯度和黑森矩阵

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

原函数：

$$f(x_1,x_2)=\frac{1}{2}\sum _{i=1}^8{r}_i^2(x_1,x_2)\text{\hspace{1em}(2.4)}$$

梯度：
$$▽f={\left[\frac{\partial f}{\partial x_1}\frac{\partial f}{\partial x_2}\right]}^T=\sum _{i=1}^8{\left[r_i\frac{\partial r_i}{\partial x_1}{r}_i\frac{\partial r_i}{\partial x_2}\right]}^T\text{\hspace{1em}(2.5)}$$

黑森矩阵：

$$H=\left[\begin{array}{cc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1{x}_2}\\ \frac{{\partial }^{2}f}{\partial x_2{x}_1}& \frac{{\partial }^{2}f}{\partial x_2^2}\end{array}\right]\approx \sum _{i=1}^8\left[\begin{array}{cc}\frac{\partial r_i}{\partial x_1}\frac{\partial r_i}{\partial x_1}& \frac{\partial r_i}{\partial x_1}\frac{\partial r_i}{\partial x_2}\\ \frac{\partial r_i}{\partial x_2}\frac{\partial r_i}{\partial x_1}& \frac{\partial r_i}{\partial x_2}\frac{\partial r_i}{\partial x_2}\end{array}\right]\text{\hspace{1em}(2.6)}$$

现在只剩下 $r$ 关于 $x_1、x_2$ 的偏导，由 $r$ 的定义，容易直接根据求导法则得出 $r$ 关于 $x_1、x_2$ 的偏导（也可使用数值微分，见第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^{\prime }(x)=\underset{\Delta x\to 0}{\lim }\frac{f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}$$

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

$$r_i=x_1{e}^{x_2{t}_i}-y_i$$

偏导即为：

$$\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_2{t}_i}-(x_1-\Delta x_1)e^{x_2{t}_i}}{2\Delta x_1}$$

$$\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}$$

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

修改程序如下，程序中使用$A$ 代替 $x_1$，$B$ 代替 $x_2$：

%原始数据
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$，允许误差 $ϵ>0$，置$k=0$；
2）$f^{\prime }(x^{(k)})<ϵ$，导数接近0即在极小点附近，得到 $x$，算法结束，否则继续；
3）计算 $x^{(k+1)}$，$x^{(k+1)}=x^k-(H+$ $\alpha$ $I{)}^{-1}$ $▽$$f$ （高斯牛顿法为$x^{(k+1)}=x^k-H^{-1}$ $▽$$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									

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

四、多变量问题

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

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

数学模型为：
$$\hat{y}=x_1{e}^{x_2{t}_1}+x_3{e}^{x_4{t}_2}$$

现在需要求解合适的参数 $x=[x_1,x_2,x_3,x_4]$，使得有：
$$\begin{gathered} \min f(x)=\frac{1}{2}\sum _{i=1}^8{r}_i^2(x) \\ {r}_i(x)=\widehat{y_i}-y_i \end{gathered}$$

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

类比式(2.5)，梯度为：

$$▽f={\left[\frac{\partial f}{\partial x_1}\frac{\partial f}{\partial x_2}\frac{\partial f}{\partial x_3}\frac{\partial f}{\partial x_4}\right]}^T=\sum _{i=1}^8{\left[r_i\frac{\partial r_i}{\partial x_1}{r}_i\frac{\partial r_i}{\partial x_2}{r}_i\frac{\partial r_i}{\partial x_3}{r}_i\frac{\partial r_i}{\partial x_4}\right]}^T\text{\hspace{1em}(4.1)}$$

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

$$H=\left[\begin{array}{ccc}\frac{{\partial }^{2}f}{x_1^2}& \cdots & \frac{{\partial }^{2}f}{\partial x_1{x}_4}\\ ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{x_4{x}_1}& \cdots & \frac{{\partial }^{2}f}{\partial x_4^2}\end{array}\right]\approx \sum _{i=1}^8\left[\begin{array}{ccc}\frac{\partial r_i}{\partial x_1}\frac{\partial r_i}{\partial x_1}& \cdots & \frac{\partial r_i}{\partial x_1}\frac{\partial r_i}{\partial x_4}\\ ⋮& \ddots & ⋮\\ \frac{\partial r_i}{\partial x_4}& \cdots & \frac{\partial r_i}{\partial x_4}\frac{\partial r_i}{\partial x_4}\end{array}\right]\text{\hspace{1em}(2.6)}$$

公式有了，写程序，程序中使用 $A,B,C,D$ 代替 $x_1,x_2,x_3,x_4$ ：

%原始数据
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算法，并给出数值计算偏导数方法。