最優(yōu)化方法——最小二乘法與梯度下降法

本篇博客主要介紹最小二乘法、梯度下降法的原理與流程，分別使用Matlab、Pycharm分別實現(xiàn)了最小二乘法、不同迭代停止條件的梯度下降法等方法對給定優(yōu)化模型進行求解并進行解之間的誤差分析對比，并進行了一定理論與應用（內附數(shù)據集和python及matlab代碼）。

一、問題

讀取附件“MatrixA_b.mat”文件中的矩陣A和向量b。建立關于矩陣, 向量，未知向量最小二乘優(yōu)化模型： $\min _{x}\|A x-b\|_{2}^{2}$

1）通過最小二乘法的正規(guī)方程，求出優(yōu)化模型的準確解；

2）利用梯度下降法迭代求出模型“近似解”，通過設置迭代停止條件，分析“近似解”與“準確解”之間的誤差。

二、實驗思路綜述

1.實驗工具及算法

本次實驗分別使用Matlab、Pycharm分別實現(xiàn)了最小二乘法、不同迭代停止條件的梯度下降法等方法對給定優(yōu)化模型進行求解并進行解之間的誤差分析對比。

2.實驗數(shù)據

本次實驗使用給定矩陣A（50x40）與向量b（50x1）組成的優(yōu)化模型 $\min _{x}\|A x-b\|_{2}^{2}$ 進行實驗內容的探究，在拓展內容的探究與嘗試中使用了部分網絡數(shù)據集。

3.實驗目標

本次實驗要求使用不同方法對給定優(yōu)化模型（最小二乘問題）進行求解并進行解的誤差分析及對比。此外，本人還在相關理論方面進行了補充，對算法應用進行了實踐。

4.實驗步驟

本次實驗大致流程如表1所示：

表1 實驗3流程

1．實驗思路綜述

2.最小二乘問題的引入

3. 最小二乘法的推導與求解

4. 梯度下降法的推導與求解

5. 不同情況解的分析及誤差對比

6. 理論拓展與應用實踐

三、最小二乘問題引入

1.最小二乘問題樣例

在求解最小二乘問題前，我們需要對其進行定義與數(shù)學模型化，故本部分引入一個二維樣例如圖1所示，一個實際的測量問題如圖2所示：

圖1 二維最小二乘問題樣例

圖2 實際最小二乘問題樣例

分析：對于圖1的問題，無法找到一條直線同時經過A、B、C三點，對于圖2的問題，我們無法求解出一組滿足條件的x1，x2，x3。

2.最小二乘問題解決方案及數(shù)學模型化

最小二乘問題：由于各種誤差，難以求得滿足問題條件的一組解（無法通過現(xiàn)有數(shù)據擬合出一條過所有數(shù)據的線或超平面）的問題。

解決方案：對于最小二乘法，核心的解決方案就是尋找該問題的近似解。并盡可能逼近原問題的目標，使殘差向量r=Ax-b在某種度量下盡可能小。最小二乘問題數(shù)學模型化如圖3所示：

圖3 最小二乘問題模型

3.相關線性代數(shù)知識導入

在后續(xù)需要用不同方法求解最小二乘問題，在此對核心相關線性代數(shù)知識進行一定補充。

3.1 梯度

梯度的本意是一個向量（矢量），表示某一函數(shù)在該點處的方向導數(shù)沿著該方向取得最大值，即函數(shù)在該點處沿著該方向（此梯度的方向）變化最快，變化率最大（為該梯度的模）。

梯度求解樣例如式1所示：

$\nabla f(z)=\left[\begin{array}{c} \frac{\partial f}{\partial z_{1}}(z) \\ \vdots \\ \frac{\partial f}{\partial z_{n}}(z) \end{array}\right]$

式1 梯度求解樣例

3.2 矩陣的逆

當一個矩陣X滿足XA=I時，X被稱為A的左逆，同理可以定義右逆。

矩陣的逆：如果矩陣A存在左逆和右逆，則左逆和右逆一定相等，此時X稱為矩陣的逆（矩陣非奇異），記作A^-1。

逆存在的判斷：對于一個矩陣的逆是否存在，有如表2中所示五種常用方法：

表2 逆矩陣存在判斷常用方法

1.若矩陣行列式不為0，可逆

2.若矩陣的秩為n，可逆

3.若存在一個矩陣B，使AB=BA=I，可逆

4.對于齊次方程AX=0，若方程只有零解，可逆

5.對于非齊次線性方程AX=b，若方程只有特解，可逆

矩陣逆的常用證明框架如圖4所示：

圖4 矩陣逆的證明框架

補充：性質(a)對任意矩陣A都成立，性質(b)對方陣矩陣A都成立。

逆矩陣求解：在編程實現(xiàn)中矩陣求逆一般使用庫函數(shù)，在不同語言中均進行了打包，如matlab中可使用inv（）求逆矩陣，用法詳見：矩陣求逆 - MATLAB inv - MathWorks 中國，用pinv（）求偽逆，用法詳見：Moore-Penrose 偽逆 - MATLAB pinv - MathWorks 中國

3.3 QR分解

QR分解是將一個矩陣A分解成具有標準正交列向量的矩陣Q和上三角矩陣R（對角線元素不為0）的算法。這個分解能夠有效的提高計算機求解線性方程、最小二乘問題、帶約束的最小二乘問題的效率，有效降低計算復雜度，QR分解形式如圖5所示。

圖5 QR分解定義形式

QR分解根據原理分為Gram-Schmidt、Householder、Givens三種實現(xiàn)方法，經本人實驗2探究發(fā)現(xiàn)，對于較稠密矩陣，使用Householder QR分解有較高的效率與穩(wěn)定性。

四、最小二乘法

在本部分對于最小二乘法的定義、數(shù)學模型化、目標求解推導、模型求解做詳解。

1.定義

最小二乘法是一種數(shù)學優(yōu)化技術。它通過最小化誤差的平方和尋找數(shù)據的最佳函數(shù)匹配。利用最小二乘法可以簡便地求得未知的數(shù)據，并使得這些求得的數(shù)據與實際數(shù)據之間誤差的平方和為最小。

最小二乘法還可用于曲線擬合，其他一些優(yōu)化問題也可通過最小化能量或最大化熵用最小二乘法來表達。在誤差估計、不確定度、系統(tǒng)辨識及預測、預報等數(shù)據處理諸多學科領域也得到廣泛應用。

2.數(shù)學模型化

2.1 目標函數(shù)

結合最小二乘問題模型與最小二乘法定義，將最小二乘法數(shù)學模型化，故對于給定的給定A∈R^mxn,b∈R^m,求解x∈R^n讓目標函數(shù)最小，目標函數(shù)如式2所示：

式2 最小二乘法目標函數(shù)

2.2 最小二乘法的解

結合最小二乘法的原理，對于式2的目標函數(shù)求解，得到的x應該滿足式3的條件：

式3 最小二乘法的解的條件

分析：當殘差r=Ax?b=0時，則x是線性方程組Ax=b的解；否則其為誤差最小平方和下方程組的近似解。

2.3 列向量空間的意義

對于滿足最小二乘法目標函數(shù)式2的解x，其列向量空間的意義如圖6所示：

圖6 最小二乘法列向量空間的意義

分析：如圖6所示，Ax∈range(A)中最接近b的向量，r=Ax-b正交(垂直)于值域空間range(A)。

3.目標求解推導

對于最小二乘法的目標函數(shù)式2，我們需要得到滿足式3條件的最優(yōu)解x。由于目標函數(shù)f（x）為可微函數(shù)，故最優(yōu)解x滿足梯度?f(x)=0，如式4所示：

4.正規(guī)方程

由最小二乘法定義與梯度公式求導可知，我們需要找到目標函數(shù)的最優(yōu)解即找梯度?f(x)=0，梯度公式推導后如式8所示，由此定義最小二乘法的正規(guī)方程如式9所示：

$A^{T} A x=A^{T} b$

式9 最小二乘法正規(guī)方程

分析：分析式9中的正規(guī)方程，其等價于?f(x)=0，f(x)= $\|A x-b\|_{2}^{2}$ ，且最小二乘法問題所有解都滿足正規(guī)方程。如果A的列線性無關，則A^TA為非奇異矩陣，此時正規(guī)方程（原問題）有唯一解。

對于正規(guī)方程的求解一般有三種方法，分別為直接求解正規(guī)方程組求解、通過Gram矩陣求解與QR分解求解，后兩種方法實現(xiàn)流程詳解如下：

4.1 通過Gram矩陣求解正規(guī)方程

通過Gram矩陣求解正規(guī)方程一般流程如表3所示：

Tips:經過四舍五入之后，Gram矩陣為奇異矩陣。

4.2 通過QR分解求解正規(guī)方程

方法②比方法①更穩(wěn)定，因為它避免構造Gram矩陣，通過QR分解求解正規(guī)方程一般流程如表4所示：

5.編程實踐

根據實驗任務1）的要求，本部分將編程實踐通過最小二乘法的正規(guī)方程，求出給定數(shù)據優(yōu)化模型的準確解。

5.1 QR分解

對實驗給定矩陣A與向量b進行導入，并對A進行QR分解（Householder），算法流程如表5所示：

在代碼實現(xiàn)上，可使用matlab庫函數(shù)[Q,R] = qr(A)使用Householder進行QR分解，用法解析可見：QR 分解 - MATLAB qr - MathWorks 中國，也可自己構建QR分解函數(shù)，分解與穩(wěn)定性分析代碼可見：最優(yōu)化方法——QR分解_@李憶如的博客-CSDN博客

5.2 求最優(yōu)解

在得到給定矩陣A分解出的Q、R矩陣（不同QR分解得到不同矩陣，需要轉換）后，根據公式10對Q、R、b進行求最優(yōu)解并編程實現(xiàn)（逆矩陣可用inv（）函數(shù)求），最終得到最優(yōu)解x_least并保存供后續(xù)對比，代碼如下：

x_least=inv(R)*Q'*b; %精確解

五、梯度下降法

除了最小二乘法，梯度下降法也常用于最優(yōu)化問題最優(yōu)解的逼近，尤其是對于R^mxn列向量線性相關或n非常大的情況，本部分對于梯度下降法法的定義、數(shù)學模型化、目標求解推導、模型求解做詳解。

1.定義

梯度下降法是一個一階最優(yōu)化算法。要使用梯度下降法找到一個函數(shù)的局部極小值，必須向函數(shù)上當前點對應梯度（或者是近似梯度）的反方向的規(guī)定步長距離點進行迭代搜索。即梯度下降法求解目標問題最優(yōu)解的過程為：x1,x2,…,xk→x，其中xk是第k步迭代，期望更新xk+1，滿足f(xk+1)<f(xk),核心原理如圖12所示：

圖12 梯度下降法核心原理

2.目標函數(shù)推導

3.操作與算法流程

對于優(yōu)化問題 $\min _{x \in \mathbb{R}^{n}} \frac{1}{2}\|A x-b\|_{2}^{2}, A \in \mathbb{R}^{m \times n}, b \in \mathbb{R}^{m}$ ，根據梯度下降法原理與目標求解推導總結其算法流程如表6所示：

其中，迭代停止條件一般有設置迭代次數(shù)與相鄰迭代解之間的“相對接近程度”兩種。

4.編程實踐

根據實驗任務2）的要求，本部分將編程實踐通過梯度下降法求出給定優(yōu)化模型的“近似解”，核心代碼如下：

%%梯度下降法
min=0.01;
x=zeros(40,1);
for k = 1:30 %或指定迭代次數(shù)
f(1,k)=0.5*norm(A*x-b,2)^2; % 目標函數(shù)值
p = A'*(A*x-b);
a = norm(p,2)^2 / norm(A*p,2)^2;
y = x - a * p; %y為x（k+1）
temp(1,k) = norm((x-y),2)/norm(x,2); %迭代解間的相對接近程度
error(1,k) = norm((x_least - x),2); %誤差迭代
% 
%   if norm((x-y),2)/norm(x,2) < min
%       break
%   end
x = y; %迭代
end

4.1 迭代次數(shù)

本部分以迭代次數(shù)作為迭代停止條件，為探究最優(yōu)迭代次數(shù)，應觀察分析不同迭代次數(shù)對于目標求解的影響（目標函數(shù)值的變化）。本次實驗中，迭代次數(shù)與目標函數(shù)值之間的關系如圖14所示：

圖14 迭代次數(shù)與目標函數(shù)值間的關系

分析：由圖14可見，迭代次數(shù)在30之后，隨迭代次數(shù)增大，目標函數(shù)趨于穩(wěn)定，故在本實驗中，選取迭代次數(shù)為30作為停止條件為較優(yōu)選擇。

本部分以相鄰迭代解之間的“相對接近程度”作為迭代停止條件，本實驗以公式： $\left\|x^{k}-x^{k+1}\right\|_{2} /\left\|x^{k}\right\|_{2}$ 為例作為評估標準。為探究最優(yōu)閾值，應觀察分析不同閾值對于目標求解的影響（目標函數(shù)值的變化）。本次實驗中，迭代次數(shù)與相鄰迭代解之間的“相對接近程度”之間的關系如圖15所示：

圖15 迭代次數(shù)與相鄰迭代解之間的“相對接近程度”之間的關系

分析：結合圖14，由圖15可見，隨迭代次數(shù)增加，相鄰迭代解之間的“相對接近程度”波動下降。經統(tǒng)計分析后，對本實驗，我選擇相鄰迭代解之間的“相對接近程度”的閾值為0.01,作為梯度下降法的終止條件。

5.不同情況解的分析及誤差對比

不同方法、語言求解最小二乘問題得到結果、效率都有所不同，本部分進行對比分析。

5.1 不同算法分析

最小二乘法與梯度下降法得到的解及對應效率是不同的，結合兩種算法的原理與流程分析可做解釋。

對于最小二乘法，核心就是求偏導，然后使偏導為0，得到理論上的“準確解”。其最后一步解方程組，計算量相對較大。

而對于梯度下降法，可以看作是更簡單的一種求最小二乘法最后一步解方程的方法，本質上是在以梯度的方向和步長向目標“準確解”迭代逼近的算法。誤差存在于梯度下降會有一個初始解，這個解往往與“準確解”的距離較遠，所以每一次迭代的步長的方向和長度都是盡量“減小”誤差，但是得到最后的解還是會與“準確解”存在一定的誤差。

總的來說，最小二乘法可以得到全局最優(yōu)的閉式解，梯度下降法是通過迭代更新來逐步進行的參數(shù)優(yōu)化方法，最終結果為局部最優(yōu)。

5.2 誤差分析

本部分對兩種不同迭代停止條件的梯度下降法求出的“近似解”與最小二乘法得到的“準確解”進行對比，然后用 $\text { distance error }=\| x_{-} \text {least }-x \|_{2}$ 做誤差分析。其中，迭代解與準確解的誤差如圖16所示，近似解與準確解之間的誤差如表7所示：

圖16 迭代解與準確解的誤差關系

分析：由圖16可見，初始化x(0)=0，以 $\| x_{-} \text {least }-x \|_{2}$ 度量誤差的情況下，梯度下降法求得的迭代解與最小二乘法求出的準確解之間的誤差隨迭代次數(shù)增加而減少，由0次迭代時誤差為2.0007，到100次迭代時誤差降為0.752。

表7 本實驗近似解與準確解之間的誤差

分析：由表7可見，初始化x(0)=0，以 $\| x_{-} \text {least }-x \|_{2}$ 度量誤差的情況下，本實驗使用的兩種梯度下降法（迭代次數(shù)=30停止， $\left\|x^{k}-x^{k+1}\right\|_{2} /\left\|x^{k}\right\|_{2}$ <0.01停止）得到的近似解與最小二乘法得到的精確解（閉式最優(yōu)）的誤差分別為1.33798527與1.491785332。而兩種梯度下降法得到的目標函數(shù)值與最小二乘法得到的函數(shù)值誤差分別為0.084586155與0.1191079252。

5.3 效率對比

為探究不同方法的效率對比，本部分用上述提到的三種方法分別針對實驗給定的最小二乘問題求解，每種方法運行20次，運行時間數(shù)據匯總如表8所示，效率對比如圖17所示：

表8 不同方法求解最小二乘問題平均運行時間匯總

圖17 不同方法求解最小二乘問題效率對比

分析：由表8與圖17可見，無論是哪種梯度下降法，平均運行時間均低于最小二乘法。

結合正確性與效率分析，最小二乘法雖然能求出相對準確的解，但需要更長的運行時間，故在面對給定的問題時，應該有選擇性的根據問題的性質選擇兩種方法中的一個。

具體來說，最小二乘法中需要計算矩陣的逆，這是相當耗費時間的，而且求逆也會存在數(shù)值不穩(wěn)定的情況，因而這樣的計算方法在應用中有時不值得提倡。

相比之下，梯度下降法雖然有一些弊端，迭代的次數(shù)可能也比較高，但是相對來說計算量并不大.而且，在最小二乘法這個問題上，收斂性有保證。故在大數(shù)據量的時候，反而是梯度下降法（其實應該是其他一些更好的迭代方法）更加值得被使用。

6.不同語言與平臺對求解的影響

為探究不同語言與平臺對最小二乘問題求解的影響，分別將最小二乘法、兩種梯度下降法在Pycharm2021中使用Python重構，具體代碼詳見附件。

分別使用maatlab與python實現(xiàn)的三種方法對實驗給定矩陣A（50x40）與向量b（50x1）進行求解，每個平臺的各個方法均進行20次求平均運行時間，數(shù)據匯總如表8所示，效果對比如圖18所示：

表8 不同語言、不同方法求解最小二乘問題的平均運行時間

圖18 不同語言、不同方法求解最小二乘問題的效率對比

分析：由表8與圖18可見，在不同方法求解最小二乘問題中，matlab的運行時間均略低于Python，效率較高。

六、理論補充與應用拓展

對于最小二乘法與梯度下降法，除了本實驗中對于矩陣向量構成的優(yōu)化模型求解，在其他方面上也有廣泛的應用，在本部分做簡單嘗試與實踐。

1.最小二乘法

1.1 線性回歸定義與算法步驟

線性回歸及其詳細應用可見：機器學習——LR（線性回歸）、LRC（線性回歸分類）與人臉識別

回歸與線性回歸：回歸分析是指一種預測性的建模技術，主要是研究自變量和因變量的關系。線性回歸為最基礎的一種回歸算法。用線（面）等模型對于現(xiàn)有相對線性的數(shù)據進行較小損失的擬合，并使擬合出的模型可較好預測數(shù)據，一般算法流程如表9所示：

表9 線性回歸算法流程

輸入: 數(shù)據集

過程：

1、變量的篩選與控制

2、對正態(tài)性分布的數(shù)據做散點圖與相關分析

3、通過最小化損失函數(shù)來確定參數(shù)，得到（擬合）回歸方程

4、不斷檢驗模型，優(yōu)化參數(shù)，得到最優(yōu)的回歸方程

5、使用回歸方程做預測

輸出：回歸方程

1.2 最小二乘法的應用

根據線性回歸定義與表9所示，在線性回歸問題中，常常使用最小二乘法來擬合數(shù)據，包括但不限于基于正規(guī)方程的解去擬合直線或超平面，預測數(shù)據。在本部分以一個實際樣例探究最小二乘法在線性回歸中的應用。

問題描述：探究學生成績與學生學習時間的關系

線性回歸實現(xiàn)：將學習時間作為變量，成績作為預測值，建立回歸方程，并用最小二乘法最小化損失函數(shù)，得到回歸方程并驗證，驗證后用其預測。核心代碼如下：

%%最小二乘法應用
x=[23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40];
y=[41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0];
figure
plot(x,y,'r*') %作散點圖(制定橫縱坐標)
xlabel('x(學生學習時間)','fontsize',12)
ylabel('y(學生成績)','fontsize',12)
set(gca,'linewidth',2)
%采用最小二乘擬合
Lxx=sum((x-mean(x)).^2);
Lxy=sum((x-mean(x)).*(y-mean(y)));
b1=Lxy/Lxx;
b0=mean(y)-b1*mean(x);
y1=b1*x+b0; %線性方程用于預測和擬合
hold on
plot(x,y1,'linewidth',2);
m2=LinearModel.fit(x,y); %函數(shù)進行線性回歸

數(shù)據與擬合結果如圖20所示：

圖20 數(shù)據與擬合結果

分析：由圖20可見，可看出模型對數(shù)據擬合較好，預測相對線性。如需預測不在圖中的數(shù)據，只需將對應學習時間作為x代入回歸模型（方程）中即可。驗證了最小二乘法在線性回歸應用中的正確性。

2.梯度下降法

2.1 BP神經網絡

BP神經網絡及其應用詳見：機器學習——深度神經網絡實踐（FCN、CNN、BP

BP神經網絡是一種簡單的神經網絡，核心思路是模仿人的大腦工作原理，構造的一個數(shù)學模型，它的仿生結構如圖21所示：

圖21 BP神經網絡拓撲圖

其中，BP神經網絡的結構包含三層，最靠前的是輸入層，中間是隱層（可以有多個隱層，每層隱層可以有多個神經元），最后是輸出層，一般工作流程如表10所示：

表10 BP神經網絡流程

輸入: 數(shù)據集

過程：

1.輸入層負責接收輸入，在輸入層接收到輸入后，每個輸入神經元會把值加權傳遞到各個隱層神經元。

2.各個隱神經元接收到輸入神經元傳遞過來的值后，與自身的基礎閾值b匯總求和，經過一個激活函數(shù)（通常激活函數(shù)是tansig函數(shù)），然后加權傳給輸出層。

3.輸出神經元把各個隱神經元傳過來的值與自身閾值b求和（求和后也可以再經過一層轉換），即是輸出值。

輸出：對應結果

2.2 梯度下降法的應用

根據神經網絡定義與表10所示，對于相關算法，參數(shù)更新是重要步驟。對于BP神經網絡而言，常用梯度下降法去更新參數(shù)，即通過反向傳播計算不同參數(shù)的梯度，再用梯度進行參數(shù)的優(yōu)化。在本部分以一個實際樣例探究梯度下降在BP神經網絡中的應用。

問題描述：鳶尾花數(shù)據的分類（根據鳶尾花的四種特征屬性去分三類）

BP神經網絡實現(xiàn)：本次實踐我選擇了四層BP神經網絡，第一層為輸入層，第二層和第三層都為中間層，第四層為輸出層。

輸入層為四個神經元（每類特征屬性都能參與計算），輸出層為三個神經元（分別對應三個類別的概率大?。?，兩個中間層根據經驗確定為二十五個神經元。

將不同層神經元進行全鏈接，中間的鏈接即為權重w，除了輸入層，其它層神經元都賦予一個偏置b以及激活函數(shù)f，另外給最后輸出結果一個評判誤差的損失函數(shù)。

權重w和偏置b通過隨機數(shù)生成，中間層激活函數(shù)設置為relu函數(shù)，輸出層激活函數(shù)設置為softmax函數(shù)（用來分類），損失函數(shù)設置為交叉熵誤差（因為分類時用到了獨熱編碼，因此適合用交叉熵誤差）。

參數(shù)更新的方法設置為隨機梯度下降法。即通過反向傳播計算不同參數(shù)的梯度，再用梯度進行參數(shù)的優(yōu)化。代碼如下：

# 訓練集：鳶尾花150*50%
# 網絡結構:輸入層(4)+中間層(25)+中間層(25)+輸出層(3)
# 中間層激活函數(shù):relu,輸出層激活函數(shù):softmax
# 損失函數(shù):交叉熵誤差
# 隨機梯度下降
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
# 鳶尾花數(shù)據讀入
iris_data = datasets.load_iris()
input_data = iris_data.data
correct = iris_data.target
n_data = len(correct)
# 對數(shù)據進行預處理
# 標準化
ave_input = np.average(input_data, axis=0)
std_input = np.std(input_data, axis=0)
input_data = (input_data - ave_input) / std_input
print(input_data)
# 標簽轉化為獨熱編碼
correct_data = np.zeros((n_data, 3))
for i in range(n_data):
    correct_data[i, correct[i]] = 1.0
print(correct_data)
# 切分訓練集和測試集
index = np.arange(n_data)
index_train = index[index % 2 == 0]
index_test = index[index % 2 != 0]
input_train = input_data[index_train, :]
input_test = input_data[index_test, :]
correct_train = correct_data[index_train, :]
corre_test = correct_data[index_test, :]
n_train = input_train.shape[0]
n_test = input_test.shape[0]
# 設置參數(shù)
n_in = 4
n_mid = 10
n_out = 3
wb_width = 0.1
eta = 0.1
epoch = 100
batch_size = 8
interval = 100
# 實現(xiàn)網絡層
class Baselayer:
    def __init__(self, n_upper, n):
        self.w = wb_width * np.random.randn(n_upper, n)
        self.b = wb_width * np.random.randn(n)
    def updata(self, eta):
        self.w = self.w - eta * self.grad_w
        self.b = self.b - eta * self.grad_b
class MiddleLayer(Baselayer):
    def forward(self, x):
        self.x = x
        self.u = np.dot(x, self.w) + self.b
        self.y = np.where(self.u <= 0, 0, self.u)  # relu函數(shù)
    def backward(self, grad_y):
        delta = grad_y * np.where(self.u <= 0, 0, 1.0)  # relu函數(shù)的求導--!!
        self.grad_w = np.dot(self.x.T, delta)
        self.grad_b = np.sum(delta, axis=0)
        self.grad_x = np.dot(delta, self.w.T)
class OutputLayer(Baselayer):
    def forward(self, x):
        self.x = x
        u = np.dot(x, self.w) + self.b
        self.y = np.exp(u) / np.sum(np.exp(u), axis=1, keepdims=True)  # SoftMax函數(shù)
    def backward(self, t):
        delta = self.y - t
        self.grad_w = np.dot(self.x.T, delta)
        self.grad_b = np.sum(delta, axis=0)
        self.grad_x = np.dot(delta, self.w.T)
# 實例化
middle_layer_1 = MiddleLayer(n_in, n_mid)
middle_layer_2 = MiddleLayer(n_mid, n_mid)
output_layer = OutputLayer(n_mid, n_out)
# 定義函數(shù)
def forward_propagation(x):
    middle_layer_1.forward(x)
    middle_layer_2.forward(middle_layer_1.y)
    output_layer.forward(middle_layer_2.y)
def back_propagation(t):
    output_layer.backward(t)
    middle_layer_2.backward(output_layer.grad_x)
    middle_layer_1.backward(middle_layer_2.grad_x)
def update_wb():
    middle_layer_1.updata(eta)
    middle_layer_2.updata(eta)
    output_layer.updata(eta)
def get_error(t, batch_size):
    return -np.sum(t * np.log(output_layer.y + 1e-7)) / batch_size
train_error_x = []
train_error_y = []
test_error_x = []
test_error_y = []
# 學習過程
n_batch = n_train // batch_size
for i in range(epoch):
    # 統(tǒng)計誤差
    forward_propagation(input_train)
    error_train = get_error(correct_train, n_train)
    forward_propagation(input_test)
    error_test = get_error(corre_test, n_test)
    train_error_x.append(i)
    train_error_y.append(error_train)
    test_error_x.append(i)
    test_error_y.append(error_test)
    index_random = np.arange(n_train)
    np.random.shuffle(index_random)
    for j in range(n_batch):
        mb_index = index_random[j * batch_size:(j + 1) * batch_size]
        x = input_train[mb_index, :]
        t = correct_train[mb_index, :]
        forward_propagation(x)
        back_propagation(t)
        update_wb()
plt.plot(train_error_x, train_error_y, label="Train")
plt.plot(test_error_x, test_error_y, label="Test")
plt.legend()
plt.xlabel("epoch")
plt.ylabel("error")
plt.show()

不同梯度下降方法更新參數(shù)下的分類結果與epoch的關系如圖22所示：

圖22 不同梯度下降方法更新參數(shù)下的分類結果與epoch的關系

分析：由圖22可見，無論哪種梯度下降更新參數(shù)，隨著epoch增加，訓練集與測試集的誤差均會減小并呈現(xiàn)較相似趨勢。但對于隨機梯度下降而言，波動與誤差較大。而自適應梯度下降更新BP神經網絡的參數(shù)較為穩(wěn)定，且兩數(shù)據集擬合效果好（分類效果好，誤差?。?。