神經網絡與機器學習
第5章 隨機梯度下降法-BP的起源
神經網絡的訓練有很多方法，以數值優化為基礎的隨機梯度學習算法能夠處理大規模的數據集合，它也是后面多層神經網絡后向傳播算法的基礎。
隨機梯度下降是以均方誤差為目標函數的近似最速下降算法，該算法被廣泛用於自適應信號處理領域，是Widrow與學生Hoff提出的，他們當時叫做自適應線性神經元（線性組合濾波器）和最小均方學習。
§5.1 自適應濾波器
我們先學習一下什么是濾波。濾波、平滑、預測都屬於信號估計，三種方式的區別如下面圖所示。濾波是給出信號$s(K)$在當前時刻$n$的估計$\widehat{s}(k)=f[x(k),\cdots,x(0)]$，$k$時刻前的數據可以用。
平滑是0在$k$到時刻數據中，$k$時刻數據缺少了，那么根據$k^'$前后時刻的數據估計$k^'$時刻的數據，稱為平滑$\widehat{s}(k)=f[x(k),\cdots,x(k^'+1),x(k^'-1),\cdots,x(0)]$。數學里叫插值。
預測是依據從0到$k$時刻數據，得到並沒有觀測到$k+\tau$時刻的數據信息，給出一個估計$\widehat{s}(k+\tau)=f[x(k),\cdots,x(0)]$。
線性濾波器：是構造輸入數據的函數$f[x(k),\cdots,x(0)]$，如果$f$為線性函數，比如
\[\widehat{s}(k)=\frac{x(k)}{4}+\frac{x(k-1)}{2}+\frac{x(k-2)}{4}\]
這個濾波器的結構非常簡單，就是當前值為3個時刻觀測值的加權平均值，線性系數[1/4，1/2，1/4]稱為濾波器系數，這類線性函數是線性濾波器，否則是非線性濾波器$f$。維納濾波和卡爾曼濾波都是線性濾波，維納濾波適用於平穩的輸入，而卡爾曼濾波是時變的最優線性濾波，適用於非穩態的輸入。
自適應濾波：按照遞歸算法自動調節濾波器系數，平穩輸入情況下收斂到維納濾波器，非平穩情況下，如果變化足夠緩慢，也能慢慢跟蹤輸入信號的變化，自適應濾波在通訊、地震、雷達、醫學等方面應用廣泛，主要可以完成系統辨識，均衡、預測、噪聲對消這四個信號處理任務，應用十分廣泛，並且由此衍生的隨機梯度下降法對於機器學習也影響很大。
自適應濾波器是一類特殊的感知神經元，因為自適應濾波器由抽頭延遲線、線性神經元組成的，如
從圖5.1可以看出自適應橫向濾波器就是一種特殊的線性神經元（感知機），激活函數選擇的是線性函數$y=x$，我們可以將神經元$k$時刻的輸出寫為
\[a(k)=w^{\mathrm{T}}p+b=\sum_{i=1}^Rw_ip_i+b=\sum_{i=1}^Rw_ix(k-i+1)+b\]
如果增廣向量
\[p\rightarrow \begin{bmatrix} p\\ 1 \end{bmatrix},w\rightarrow \begin{bmatrix} w\\ b \end{bmatrix}\] 5.3
自適應濾波器或者線性神經元可以簡化為
$a(k)=\bm{w}^{\mathrm{T}}\bm{p}$ 5.4
自適應濾波器也和神經網絡很像，有訓練和測試兩個過程，在訓練過程中，已知期望響應$t(k)$，通過對權向量的設計，使從神經元輸出在某種意義下最佳逼近期望響應。在工作過程中，權向量由訓練過程得到，對未知輸入（通常是隨機信號）濾波得到響應信號。
性能曲面函數
誤差信號
$e(k)=t(k)-a(k)=t(k)-\bm{W}^{\mathrm{T}}\bm{p}$ 5.5
其平方稱為瞬時平方誤差
$e^2(k)=t^2(k)-2t(k)\bm{w}^{\mathrm{T}}\bm{p}+\bm{w}^{\mathrm{T}}\bm{p}\bm{p}^{\mathrm{T}}\bm{w}$ 5.6
設$p(k)$、$t(k)$是平穩的隨機過程樣本-隨機信號，滿足一定的分布，那么我們希望均方誤差
$E[e^2]=E[t^2]-2\bm{w}^{\mathrm{T}}E[t\bm{p}]+\bm{w}^{\mathrm{T}}E[\bm{p}\bm{p}^{\mathrm{T}}]\bm{w}$ 5.7
最小，這里輸入數據向量$p$的自相關矩陣
\[R=E[pp^{\mathrm{T}}]\] 5.8
輸入和期望信號的互相關向量為
$d=E[t\bm{p}]$ 5.9
那么可以進一步將均方誤差$E[\epsilon^2(k)]$寫成如下形式
\[\xi =E[e^2(k)]=E[t^2(k)]-2dw+w^{\mathrm{T}}Rw\] 5.10
可以看出均方誤差是權向量$w$的二次函數，其二階導數-Hessian矩陣
\[\triangledown ^2\xi =\frac{\partial ^2\xi}{\partial w^2} =2R\succeq 0\]
$R\succeq 0$表示$R$矩陣非負定，我們一般遇到$R\succ 0$正定，因此均方誤差是權向量$w$的二次凸函數，必定有最小解，那么通過一次導數等於零，解出最優權向量$w^o$
\[\triangledown \xi =\frac{\partial \xi}{\partial w} =2Rw- 2d=0\Rightarrow w^o=R^{-1}d\] 5.11
那么最優權向量$w^o$代入公式（3），最小均方誤差$\xi_{min}$寫為
\[\xi_{min}=E[t^2(k)]-d^{\mathrm{T}}w^o=E[t^2(k)]-d^{\mathrm{T}}R^{-1}d\] 5.12
我們也一般利用第4章學的最速下降法、牛頓法等迭代求解最優權向量$w^o$。
§5.2 隨機梯度下降法
我們現實中遇到這樣的問題，自相關矩陣$R$得不到，我們不知道隨機輸入信號的分布，我們只知道樣本集合
\[\left \{ p_1,t_1 \right \},\left \{ p_2,t_2 \right \},\cdots,\left \{ p_Q,t_Q \right \}\]
下標$q=1,2,\cdots,Q$是樣本數量。最速下降法、牛頓法等迭代求解方法都無法進行，Widrow提出了一種近似最速下降算法，將均方誤差$\xi=E[e^2(k)]$被瞬時平方誤差代替
\[\widetilde{\xi}=e^2(k)\]
那么每次迭代計算尋優時，都計算梯度估計值
\[\widetilde{\triangledown\xi}=\triangledown e^2(k)=2e(k)\triangledown e(k)=2e(k)\frac{\partial e(k)}{\partial w}\]
由於$e(k)=t(k)-a(k)=t(k)-\bm{W}^{\mathrm{T}}\bm{p}$，所以上式可以寫為
$\widetilde{\triangledown\xi}=\triangledown e^2(k)=-2e(k)\bm{p}$
可以看出梯度近似值只是誤差和輸入向量的乘積而已，所以依次放入樣本${ p_1,t_1},{ p_2,t_2 },\cdots,{ p_Q,t_Q }$，恆定的學習率$\alpha$，則依次更新權向量
\[w(q+1)=w(q)-\alpha \widetilde{\triangledown\xi}=w(q)+2\alpha e(k)p\]
這樣每次樣本輸入，神經網絡權值進行更新，這種學習方式又稱為在線學習，最小均方算法，$delta$規則，隨機梯度下降法，或者Widrow-Hoff學習。
按照第4章最優化中的二次函數分析，理論分析學習率$\alpha<\frac{2}{\lambda_{max}}$與自相關矩陣$R$最大特征值$\lambda_{max}$相關，但是現實中不知道自相關矩陣，因此學習率經常是試錯法選擇。理論分析學習可以在自適應濾波理論中詳細學習，這里不再一一講解。
這里列舉兩個例子來看隨機梯度下降法的應用。一是利用隨機梯度下降法訓練感知機，另外一個是自適應噪聲對消。
例子5.1 感知機的訓練
我們有兩個樣本分別表示橘子和蘋果，各自的期望目標輸出為
\[\left \{ p_1=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_1=-1 \right \},\left \{ p_2=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_2=-1 \right \}\]
我們利用隨機梯度下降法，依次輸入${p_1,t_1},{p_2,t_2}$，${p_1,t_1},{p_2,t_2}$，...,因為只有兩個樣本，所以可以組成這樣的在線訓練序列，${p_1,t_1},{p_2,t_2}$的順序可以隨機輸入。初始化，我們選擇這樣的感知機
\[a=w^{\mathrm{T}}p=[w_1,w_2,w_3]\begin{bmatrix} p_1\\ p_2\\ \vdots\\ p_R \end{bmatrix}\]
$f$取線性函數$f(x)=x$ ，偏置為$b=0$，初始的權向量為$w(0)=[0,0,0]^{\mathrm{T}}$，學習率$\alpha=0.2$，那么我們輸入${ p_1=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_1=-1 }$，得到
\[a(0)=w^{\mathrm{T}}(0)p(0)=w^{\mathrm{T}}(0)p_1=[0,0,0]\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix}=0\]
誤差為
\[e(0)=t(0)-a(0)=t_1-a(0)=-1\]
那么按照隨機梯度下降方法更新權
\[w(1)=w(0)+2\alpha e(0)p_1=\begin{bmatrix} 0\\ 0\\ 0 \end{bmatrix}+2\times 0.2\times (-1)\times \begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix}=\begin{bmatrix} -0.4\\ 0.4\\ 0.4 \end{bmatrix}\]
我們輸入${ p_2=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_2=-1 }$，得到
\[a(1)=w^{\mathrm{T}}(1)p(1)=w^{\mathrm{T}}(1)p_2=[-0.4,0.4,0.4]\begin{bmatrix} 1\\ 1\\ -1 \end{bmatrix}=-0.4\]
誤差為
\[e(1)=t(1)-a(1)=t_2-a(0)=1.4\]
那么按照隨機梯度下降方法更新權
\[w(2)=w(1)+2\alpha e(1)p_2=\begin{bmatrix} -0.4\\ 0.4\\ 0.4 \end{bmatrix}+2\times 0.2\times (1.4)\times \begin{bmatrix} 1\\ 1\\ -1 \end{bmatrix}=\begin{bmatrix} 0.16\\ 0.96\\ -0.16 \end{bmatrix}\]
我們再次輸入${p_1,t_1}$，得到
\[a(2)=w^{\mathrm{T}}(2)p(2)=w^{\mathrm{T}}(2)p_1=-0.64\]
\[e(2)=t(2)-a(2)=t_1-a(2)=-0.36\]
\[w(3)=w(2)+2\alpha e(2)p_1=\begin{bmatrix} 0.016\\ 1.104\\ -0.016 \end{bmatrix}\]
這樣訓練下去，我們會得到
\[w(\infty)=\begin{bmatrix} 0\\ 1\\ 0 \end{bmatrix}\]
這比前面第3章得到感知機決策邊界更加好，因為前面的決策邊界模式一旦划分正確則停止，而現在的隨機梯度下降法是使得決策邊界盡量原來樣本模式，使得近似均方誤差最小，兩種決策邊界雖然都處於樣本模式之間，但是隨機梯度法得到決策邊界在測試中是非常實用的。
例子5.2 自適應噪聲對消
圖5.3 自適應噪聲對消
任務：設噪聲源$u(k)$將語音等信號$s(k)$污染，我們不知道噪聲強度多少，我們想在語音源附近也測量同源噪聲$\beta(k)$，然后用自適應算法預測噪聲$u(k)$，然后對消噪聲，得到誤差就是相對純凈的信號發送出去。
設信號是隨機相位信號
\[s(k)=cos(\frac{\pi k}{16}+\varphi )\]
污染信號的噪聲是
\[u(k)=0.9u(k)+\eta(k)\]
$\eta(k)$是標准正態分布的高斯白噪,$\sigma_{\eta}^2=1$，而另外一個渠道得到的噪聲和$u(k)$同源，但是不同的信道
\[\beta(k)=-0.75\beta(k)+\eta(k)\]
用隨機梯度下降算法進行自適應濾波，不需要任何關於隨機輸入的二階矩的知識。主程序如下：
%System parameters
omg0 = pi/16; a1 = 0.9; a2 = -0.75; var_eta = 1;
% Generate Desired and Input signals
order=20;
N = 1000;%data length
n = [0:N-1]';
phi = 2*pi*rand(1,1);
s = cos(omg0*n+phi);
eta = randn(N,1);
u = filter(1,[1,-a1],eta);
beta = filter(1,[1,-a2],eta);
t = s + u;
mu = 0.005;%step
w0 = zeros(order,1);
[w,e,yhat] = firlms(beta,t,mu,order,w0);
plot(e,'b-');
hold on
plot(s,'r--')
| function [w,e,yhat] = firlms(x,t,mu,L,w0) % Computes filter coefficients w and the corresponding error in % e, given signal x, desired signal y, step size mu, filter order L, and % the initial coefficient vector w0. N = length(x); % Length of x(n) & number of iterations x = reshape(x,N,1); % x as a column vector w = zeros(L,N); % Initialization of w X = zeros(L,N); % Data matrix xx = [zeros(L-1,1);x]; % Prepend L-1 zeros for i = 1:L X(i,:) = xx(L-i+1:end-i+1)'; end e = zeros(1,N); yhat = zeros(N,1); %--LMS Algorithm for FIR filter - real-valued case % Algorithm initialization yhat(1) = w0'*X(:,1); e(1) = t(1) - yhat(1); w(:,1) = w0 + 2*mu*X(:,1)*e(1); % Iterations for n = 2:N yhat(n) = w(:,n-1)'*X(:,n); e(n) = t(n) - yhat(n); w(:,n) = w(:,n-1) + 2*mu*X(:,n)*e(n); end |
一開始濾波器並不能消除噪聲，隨着時間推移，慢慢濾波器系數（神經元權向量收斂），對消了噪聲，余下相對純凈的信號。這種系統在飛機、戰場、大風等惡劣天氣中個人通信等等場景中經常遇到。盡管飛行器發動機、風嘯等噪聲和人的語音同時進入聲音發射器（麥克風），但是相對清晰的語音被發送過去，當然如果遇到與語音頻率特別接近的聲音是無法消除的。
圖5.4 自適應噪聲對消示例
例子5.3（*選講）：二階自回歸AR(2)過程
\[x(n)-0.96x(n_1)+0.9x(n-2)=\eta(n)\]
$\eta(n)$是高斯白噪聲，均值為0，方差為$\sigma_{\eta}^2$，我們利用一個二階濾波器進行系統辨識
\[y(n)=w_1x(n-1)+w_2x(n-2)\]
%AR過程
N = 501; mu= 0.01;
a1 = -0.950; a2 = 0.9; varx = 1; r0 = varx;
vareta = ((1-a2)*((1+a2)^2-a1^2)/(1+a2))*varx;
r1 = -a1/(1+a2)*r0; r2 = (-a2+a1^2/(1+a2))*r0;
lam1 = (1-a1/(1+a2))*varx; lam2 = (1+a1/(1+a2))*varx;
XR = lam1/lam2;d = [r1;r2];%條件數
R = toeplitz([r0,r1,r2]);%自相關陣
w = sqrt(vareta)*randn(N,1);
x = filter(1,[1 a1 a2],w);%傳遞函數形式
%觀測隨機梯度算法的性能
[w,e] = lplms(x,0.01,2,[0 0]');
subplot(121)
plot(w(1,:),w(2,:));
hold on
[w1,w2]=meshgrid(-1:0.01:3,-2:0.01:1);
z=r0+r0*w1.^2+r0*w2.^2+2*r1*w1.*w2-2*r1*w1-2*r2*w2;%MSE
contour(w1,w2,z,50);
subplot(122)
plot(w(1,:),'b--*');
hold on
plot(w(2,:),'r-^')
plot([0 501],[0.95 0.95],'--')
plot([0 501],[-0.9 -0.9],'--')
| function [w,e] = lplms(x,mu,L,w0) % Computes filter coefficients in w and the corresponding error in % e, given signal x, desired signal d, step size mu, filter order L, and % the initialcoefficient vector w0. N = length(x); % Length of x(n) x = reshape(x,N,1); % x as a column vector d = x(2:end); % Generate desired d(n) N = length(d); % Number of iterations c = zeros(L,N); % Initialization of c X = zeros(L,N); % Data matrix xx = [zeros(L-1,1);x]; for i = 1:L X(i,:) = xx(L-i+1:end-i)'; end e = zeros(1,N); %--LMS Algorithm % Algorithm initialization yhat = w0'*X(:,1); e(1) = d(1) - yhat; w(:,1) = 2*mu*X(:,1)*e(1); % Iterations for n = 2:N yhat = w(:,n-1)'*X(:,n); e(n) = d(n) - yhat; w(:,n) = w(:,n-1) + 2*mu*X(:,n)*e(n); end w = [w0,w]; e = [x(1),e]; |
圖5.5 隨機梯度下降法中函數值和權向量收斂學習曲線
\[y(n)=w_1x(n-1)+w_2x(n-2)\]
權系數$w_1$和$w_2$相對"隨機"地收斂到-0.95和0.9，達到系統辨識的目的。
作業：利用隨機梯度下降算法訓練四分類判定決策界面問題。
圖3.2 四類數據的幾何表示
訓練數據為8個目標向量，分別是
\[\left \{ (p_1=\begin{bmatrix} 1\\ 1 \end{bmatrix},t_1=\begin{bmatrix} -1\\ -1 \end{bmatrix}),(p_2=\begin{bmatrix} 1\\ 2 \end{bmatrix},t_2=\begin{bmatrix} -1\\ -1 \end{bmatrix}) \right \}\]
類1標記為○
\[\left \{ (p_3=\begin{bmatrix} 2\\ -1 \end{bmatrix},t_3=\begin{bmatrix} -1\\ 1 \end{bmatrix}),(p_4=\begin{bmatrix} 2\\ 0 \end{bmatrix},t_4=\begin{bmatrix} -1\\ 1 \end{bmatrix}) \right \}\]
類2標記為□
\[\left \{ (p_5=\begin{bmatrix} -1\\ 2 \end{bmatrix},t_5=\begin{bmatrix} 1\\ -1 \end{bmatrix}),(p_6=\begin{bmatrix} -2\\ 1 \end{bmatrix},t_6=\begin{bmatrix} 1\\ -1 \end{bmatrix}) \right \}\]
類3標記為●
\[\left \{ (p_7=\begin{bmatrix} -1\\ -1 \end{bmatrix},t_7=\begin{bmatrix} 1\\ 1 \end{bmatrix}),(p_8=\begin{bmatrix} -2\\ -2 \end{bmatrix},t_8=\begin{bmatrix} 1\\ 1 \end{bmatrix}) \right \}\]
類4標記為■
神經元選為帶偏置的
\[a=Wp+b=\begin{bmatrix} W & b \end{bmatrix}\begin{bmatrix} p\\ 1 \end{bmatrix}\]
設初值
\[W(0)=\begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix},b(0)=\begin{bmatrix} 1\\ 1 \end{bmatrix}\]
或者增廣
\[W(0)=\begin{bmatrix} 1 & 0 & 1\\ 0 & 1 & 1 \end{bmatrix}\]
那么利用隨機梯度下降算法40次，並畫出收斂的決策邊界。
| 2021年3月27日星期六 clc,clear |