BP神經網路實現（R語言+Python）

04-25

看吳恩達老師機器學習課程的反向傳播神經網路（BP-Network）一節時，網上有個《十一行Python代碼實現一個神經網路》小例子，非常小巧玲瓏，不過和吳恩達老師實現過程不太一樣（歡迎一起討論），於是又找到了一篇《R語言13行代碼實現神經網路》，正是想要的思路，和大家分享！

神經網路建模任務描述：

已知輸入為Input，輸出為Output，進行BP-Network建模，即權重矩陣 $Theta$ 的求解。

Python BP神經網路實現代碼：（一層隱藏層）

import numpy as npX = np.array([ [0,0,1],[1,1,1],[1,0,1],[0,1,1] ]) #輸入值，4*3矩陣y = np.array([[0,1,1,0]]).T #輸出值，4*1矩陣syn1 = 2*np.random.random((3,4)) - 1 # 輸入層到隱藏層的3*4 權重陣，初始值隨機產生syn2 = 2*np.random.random((4,1)) - 1 # 隱藏層到輸出層的4*1 權重陣，初始值隨機產生for j in range(20000): #迭代 a1 = 1/(1+np.exp(-np.dot(X,syn1))) #隱藏層經權重和sigmod函數轉換後的值，4*4 矩陣 a2 = 1/(1+np.exp(-np.dot(a1,syn2))) #輸出層經權重和sigmod函數轉換後的值，4*1 矩陣 error_2 = a2-y #輸出層誤差，4*1 矩陣 error_1 = np.dot(error_2,syn2.T)*(a1*(1-a1)) #隱藏層誤差，4*4 矩陣 syn2 = syn2 - np.dot(a1,error_2) #應用梯度下降法求輸出層權重陣 syn1 = syn1 - np.dot(X.T,error_1) #應用梯度下降法求隱藏層權重陣print(a2)

Python 運行結果：

[[ 3.76271664e-05] [ 9.99943981e-01] [ 9.99943983e-01] [ 5.19326828e-05]]

R 語言 BP神經網路實現代碼：（一層隱藏層）

c <- c(0,0,1,1,1,1,1,0,1,0,1,1)x <- matrix(c,nrow = 4,ncol = 3,byrow = T) #輸入值，4*3矩陣y <- matrix(c(0,1,1,0),ncol = 1) #輸出值，4*1矩陣#set.seed(111)syn1 <- matrix(rnorm(12),nrow = 3,byrow = T) # 輸入層到隱藏層的3*4 權重陣，初始值隨機產生syn2 <- matrix(rnorm(4),nrow = 4,byrow = T) # 隱藏層到輸出層的4*1 權重陣，初始值隨機產生cost <- c(NA) #損失函數，畫圖用for (i in 1:10000){ a1 <- 1/(1+exp(-(x %*% syn1))) #隱藏層經權重和sigmod函數轉換後的值，4*4 矩陣 a2 <- 1/(1+exp(-(a1%*% syn2))) #輸出層經權重和sigmod函數轉換後的值，4*1 矩陣 error_2 <- a2 - y #輸出層誤差，4*1 矩陣 error_1 <- error_2 %*% t(syn2) * a1 * (1-a1 ) #隱藏層誤差，4*4 矩陣 syn2 <- syn2 - 0.06 * a1 %*% error_2 #應用梯度下降法求輸出層權重陣，0.06為學習率 syn1 <- syn1 - 0.06 * t(x) %*% error_1 #應用梯度下降法求隱藏層權重陣，0.06為學習率 cost[i] <- error_2 # 存儲輸出層誤差為損失函數}a1a2plot(cost)

R 語言運行結果：

> a1 [,1] [,2] [,3] [,4][1,] 0.088806742 0.071804692 0.95895670 1.334214e-02[2,] 0.002789701 0.001092767 0.99987193 8.641429e-01[3,] 0.736807129 0.817886530 0.09021681 3.025079e-05[4,] 0.074378684 0.059717055 0.97070441 1.403001e-02> a2 [,1][1,] 0.003723701[2,] 0.999205564[3,] 0.997794973[4,] 0.003076961

R 語言損失函數迭代圖：

先讓我們來看看神經網路是怎麼回事

網路圖

圖中，有四行輸入觀測變數，變數的特徵數為3，即 $X[4 imes3]$ :

> x [,1] [,2] [,3][1,] 0 0 1[2,] 1 1 1[3,] 1 0 1[4,] 0 1 1

輸出變數為 $Y[4 imes1]$ :

> y [,1][1,] 0[2,] 1[3,] 1[4,] 0

該網路共有三層，第一層 $a^{0}$ 層為輸入層，即 $X$ ；第二層 $a^{1}$ 層為隱藏層；第三層 $a^{2}$ 層為輸出層。

網路最上層的兩個值為1的圈為偏置單元，所謂偏置單元，可以這樣理解：

比如擬合一個線性模型 $h_{ heta}(x)= heta_{1}x+ heta_{0}$ ,此時為了表示成與 $heta$ 向量或矩陣同維度的形式，輸入矩陣 $x$ 實際為[1, $x$ ]的形式，而這個1就是偏置單元。

下面來看看前向傳播的概念：

$a^{0}=X_{[4 imes3]};$

第一層為輸入層，不做處理，上標 j = 0為層數；

$z_{[4 imes1]}^{1}=X_{[4 imes3]} imes Theta_{[3 imes1]}^{1};$

第二層有3個單元（先不考慮偏置單元），這裡以其中某一單元為例，注意這裡為單獨的一個單元，相當於網路圖中藍色箭頭指向的單元， $z$ 為輸入經過了權重 $Theta^{1}$ 的變換；（一定要搞清楚這裡矩陣相乘的維度）

對於第二層的所有單元，有： $z_{[4 imes4]}^{1}=X_{[4 imes3]} imes Theta_{[3 imes4]}^{1};$

這裡矩陣的維度發生了變化， $Theta^{j}$ 代表從 $a^{0}$ 層映射到第 $a^{1}$ 層時的權重矩陣，其維度為以要映射層的單元數量為行數，以前一層的單元數加一（加一為偏置單元）為列數的矩陣。例如本例中的尺寸為3*4。

$a_{[4 imes4]}^{1}=sigmod(z_{[4 imes4]}^{1});$

就是如下代碼部分：

a1 <- 1/(1+exp(-(x %*% syn1))) #隱藏層經權重和sigmod函數轉換後的值，4*4 矩陣a2 <- 1/(1+exp(-(a1%*% syn2))) #輸出層經權重和sigmod函數轉換後的值，4*1 矩陣

$sigmod()$ 函數為將連續值轉換成 $(0,1)$ 區間的值，具有如下形狀：

$sigmod()$ 函數公式： $sigmod(x)=frac{1}{1+e^{-x}};$

$sigmod()$ 函數的導數公式： $S(x)^{}=frac{e^{-x}}{(1+e^{-x})^{2}}=S(x)cdot(1-S(x));$

因此，注意求 $z$ 時是矩陣之間的乘（Python中是np.dot乘，R中是%*%乘），而 $sigmod()$ 函數變換是對矩陣 $z$ 中各個數值的變換，在程序里就是矩陣中各個元素與實數的乘（*乘）！

以上即為前向傳播的概念，其實就是兩個步驟，一步是輸入乘以權重矩陣，並經過sigmod函數變化後作為下一層的輸入，直到得到最後一層的 $a^{j}$ ,即輸出值。

再來看一下反向傳播的概念：

反向網路圖

圖中，我們從最後一層，即輸出層計算輸出誤差 $e^{2}=a^{2}-y$ ,從輸出層向左看，即下一層是隱藏層 $a^{1}$ 層，而 $a^{1}$ 層的誤差 $e^{1}$ 可以看成是前一層誤差 $e^{2}$ 的函數（反向傳播即指的是誤差的反向函數），之前計算的權重陣 $Theta^{2}$ 依然是此處的權重陣，不過在此基礎上另加了一項 $a^1$ 的導數（ $a^{1}=a^1 imes(1-a^1)$ ）,即：

$e_{[4 imes4]}^1=e_{[4 imes1]}^2 imes (Theta_{[4 imes1]}^2).T* a^{1}$ # T是轉置的意思，此處*是數值直接相乘的意思；

代碼部分為：

error_2 <- a2 - y #輸出層誤差，4*1 矩陣error_1 <- error_2 %*% t(syn2) * a1 * (1-a1 ) #隱藏層誤差，4*4 矩陣

如果有多個隱藏層同樣求法（公式此處不做證明），注意輸入沒有誤差。

由此，我們求出了除了輸入層以後的各層誤差（各層誤差與各層的 $a^j$ 同維度），那麼求這個誤差有什麼用呢？下面就是神經網路迭代使損失函數最小的關鍵環節：

我們先寫梯度下降函數：

$Theta=Theta-alpha cdotpartialfrac{J(Theta)}{Theta}$ # $J( heta)$ 為損失函數，可以簡單理解為模型輸出與實際輸出的差值； $alpha$ 為學習率；

即，我們想要損失函數最小時的 $Theta$ ,那我們就沿著損失函數對 $Theta$ 求偏導的方向來更新 $Theta$ ,學習率可以嘗試著選取（文章開頭的程序里選的0.06）。

而經數學推導有： $partialfrac{J(Theta)}{Theta}=a^j imes e^{j+1}$

代碼部分：

syn2 <- syn2 - 0.06 * a1 %*% error_2 #應用梯度下降法求輸出層權重陣，0.06為學習率syn1 <- syn1 - 0.06 * t(x) %*% error_1 #應用梯度下降法求隱藏層權重陣，0.06為學習率

這樣利用前面由前向傳播求出的 $a^j$ ,由反向傳播求出的 $e^{j+1}$ ,我們就可以應用梯度下降演算法來更新 $Theta$ （ $Theta$ 陣的初始值可取一個（-1，1）之間的數）了！

讓我們來重新看一下文章開頭部分的代碼，是不是有一種恍然大明白的感覺！現在快動手來嘗試一下吧~~

（如有疏漏，歡迎指正！）