1. 程式人生 > >deep learning Softmax分類器(L-BFGS,CG,SD)

deep learning Softmax分類器(L-BFGS,CG,SD)

前言:

      在最優化計算方法中,我已經講到了機器學習常用的一些引數優化的方法,如梯度法,共軛梯度法,牛頓法,擬牛頓法,在《最優化計算方法》板塊,我都用迴歸分析比較了這些引數優化的方法,從現在開始,我將把這些引數優化的方法用來訓練分類器。

  在本節中,我講介紹softmax分類器,該分類器是在logistic迴歸模型在多分類問題上的推廣。在多分類問題中,分類標籤y可以取兩個以上的值。softmax分類器對於諸如MNIST手寫數字分類等問題上有很好的效果,該問題就是識別不同的單個數字和影象(可以應用於車牌識別、目標檢測等方向),softmax迴歸時有監督,後面我會介紹它在deep learning

中的應用,對於以後我介紹的分類器我都將用到MNIST

同時我會講Steepest Descent(SD)、CG方法、L-BFGS這些調參方法對識別效果的影響。

softmax簡介  

我們的訓練集有m個已標記的樣本構成:{(x(1),y(1)),...,(x(m),y(m))},輸入特徵是n+1維的,其中x_0代表截距項,也就是通常所說的偏置。在softmax迴歸中,我們解決的是多分類的問題,類標y可以去k個不停的值。因此,對於訓練集{(x(1),y(1)),...,(x(m),y(m))},我們有,分類的時候0替換為10。在MNIST數字識別的任務中,我們有k=10個不同的類別。

  對於給定的測試輸入x,我們想用假設函式針對每一類別j估算出概率值p(y=j|x)。也就是說,我們想估計x的每一種分類結果出現的概率。因此,我們的假設函式將要輸入一個k維的向量(向量元素的總和為1)來表示這k個估計的概率值。具體地說,我們的假設函式

形式如下:

                 

其中\theta_1, \theta_2, \ldots, \theta_k \in \Re^{n+1}是模型的引數。請注意\frac{1}{ \sum_{j=1}^{k}{e^{ \theta_j^T x^{(i)} }} }這一項對概率分佈進行歸一化,使得所有的概率之和為1。

為了方便起見,我們同樣的使用符號θ來表示全部的模型引數。在實現softmax迴歸時,將θ用一個k×(n+1)的矩陣來表示會很方便的,該矩陣是將\theta_1, \theta_2, \ldots, \theta_k 按列羅列起來的,如下所示:

'

代價函式

現在我們來介紹softmax分類器中的代價函式。在下面的公式中,1{·}是示性函式,其取值規則為:1{值為真的表示式}=1,1{值為假的表示式} =0.舉例來說明1{1+1=2}=1,1{1+1=3}=0,我們的代價函式可以表示如下:

熟悉logistic迴歸的讀者都知道,softmax的損失函式和

logistic的損失函式非常類似,只是在softmax損失函式中對類標記的k個可能值進行累加。注意在softmax迴歸中將x分類為類別j的概率為:

p(y^{(i)} = j | x^{(i)} ; \theta) = \frac{e^{\theta_j^T x^{(i)}}}{\sum_{l=1}^k e^{ \theta_l^T x^{(i)}} }.

可以用如下圖例進行表示,類似於神經網路,只是少了一層隱含層。


    輸出最大值概率值則屬於相應的類,比如說,p1最大,就屬於第一類,這個網路在後續的CNNs中會用到。

    對於J(θ)的最小化問題,目前還沒有閉式的解法,因此,我們使用迭代優化演算法(SD,CG,L-BFGS)。經過求導,我們得到梯度公式如下:

\begin{align}\nabla_{\theta_j} J(\theta) = - \frac{1}{m} \sum_{i=1}^{m}{ \left[ x^{(i)} \left( 1\{ y^{(i)} = j\}  - p(y^{(i)} = j | x^{(i)}; \theta) \right) \right]  }\end{align}

有了上面的偏導數公式之後,我們就可以將它帶入到這些引數優化的演算法中去,進而來優化J(θ)。例如,在梯度下降的標準實現中,每一次迭代需要進行如下更新:\textstyle \theta_j := \theta_j - \alpha \nabla_{\theta_j} J(\theta)(\textstyle j=1,\ldots,k)。

softmax迴歸模型引數化的特點

Softmax分類器有一個不尋常的特點:它有一個“冗餘”的引數集。為了便於闡述這一特點,假設我們從引數向量 \textstyle \theta_j 中減去了向量 \textstyle \psi,這時,每一個 \textstyle \theta_j 都變成了 \textstyle \theta_j - \psi(\textstyle j=1, \ldots, k)。此時假設函式變成了以下的式子:

\begin{align}p(y^{(i)} = j | x^{(i)} ; \theta)&= \frac{e^{(\theta_j-\psi)^T x^{(i)}}}{\sum_{l=1}^k e^{ (\theta_l-\psi)^T x^{(i)}}}  \\&= \frac{e^{\theta_j^T x^{(i)}} e^{-\psi^Tx^{(i)}}}{\sum_{l=1}^k e^{\theta_l^T x^{(i)}} e^{-\psi^Tx^{(i)}}} \\&= \frac{e^{\theta_j^T x^{(i)}}}{\sum_{l=1}^k e^{ \theta_l^T x^{(i)}}}.\end{align}

   換句話說,從 \textstyle \theta_j 中減去 \textstyle \psi 完全不影響假設函式的預測結果!這表明前面的softmax分類器中存在冗餘的引數。更正式一點來說,Softmax分類器被過度引數化了。對於任意一個用於擬合數據的假設函式,可以求出多組引數值,這些引數得到的是完全相同的假設函式 \textstyle h_\theta

   進一步而言,如果引數 \textstyle (\theta_1, \theta_2,\ldots, \theta_k) 是代價函式 \textstyle J(\theta) 的極小值點,那麼 \textstyle (\theta_1 - \psi, \theta_2 - \psi,\ldots,\theta_k - \psi) 同樣也是它的極小值點,其中 \textstyle \psi 可以為任意向量。因此使 \textstyle J(\theta) 最小化的解不是唯一的。(有趣的是,由於 \textstyle J(\theta) 仍然是一個凸函式,因此梯度下降時不會遇到區域性最優解的問題。但是 Hessian 矩陣是奇異的/不可逆的,這會直接導致採用牛頓法優化就遇到數值計算的問題)

   注意,當 \textstyle \psi = \theta_1 時,我們總是可以將 \textstyle \theta_1替換為\textstyle \theta_1 - \psi = \vec{0}(即替換為全零向量),並且這種變換不會影響假設函式。因此我們可以去掉引數向量 \textstyle \theta_1 (或者其他 \textstyle \theta_j 中的任意一個)而不影響假設函式的表達能力。實際上,與其優化全部的 \textstyle k\times(n+1) 個引數 \textstyle (\theta_1, \theta_2,\ldots, \theta_k) (其中 \textstyle \theta_j \in \Re^{n+1}),我們可以令 \textstyle \theta_1 =\vec{0},只優化剩餘的 \textstyle (k-1)\times(n+1) 個引數,這樣演算法依然能夠正常工作。

   在實際應用中,為了使演算法實現更簡單清楚,往往保留所有引數 \textstyle (\theta_1, \theta_2,\ldots, \theta_n),而不任意地將某一引數設定為 0。但此時我們需要對代價函式做一個改動:加入權重衰減。權重衰減可以解決 softmax分類的引數冗餘所帶來的數值問題。

權重衰減

我們通過新增一個權重衰減項 \textstyle \frac{\lambda}{2} \sum_{i=1}^k \sum_{j=0}^{n} \theta_{ij}^2 來修改代價函式,這個衰減項會懲罰過大的引數值,現在我們的代價函式變為:

\begin{align}J(\theta) = - \frac{1}{m} \left[ \sum_{i=1}^{m} \sum_{j=1}^{k} 1\left\{y^{(i)} = j\right\} \log \frac{e^{\theta_j^T x^{(i)}}}{\sum_{l=1}^k e^{ \theta_l^T x^{(i)} }}  \right]              + \frac{\lambda}{2} \sum_{i=1}^k \sum_{j=0}^n \theta_{ij}^2\end{align}


有了這個權重衰減項以後 (\textstyle \lambda > 0),代價函式就變成了嚴格的凸函式,這樣就可以保證得到唯一的解了。 此時的 Hessian矩陣變為可逆矩陣,並且因為\textstyle J(\theta)是凸函式,梯度下降法和L-BFGS等演算法可以保證收斂到全域性最優解。

為了使用優化演算法,我們需要求得這個新函式 \textstyle J(\theta) 的導數,如下:

\begin{align}\nabla_{\theta_j} J(\theta) = - \frac{1}{m} \sum_{i=1}^{m}{ \left[ x^{(i)} ( 1\{ y^{(i)} = j\}  - p(y^{(i)} = j | x^{(i)}; \theta) ) \right]  } + \lambda \theta_j\end{align}

通過最小化 \textstyle J(\theta),我們就能實現一個可用的softmax 分類器。

下面將給出實驗結果
圖1 損失函式隨迭代次數變化的曲線     其中softmax-CG-DY的CG-DY方法中β是我們老師的老師發明的,所以CG方法我用了DY方法來β,具體的方法見最優化計算方法板塊。
圖2 識別率隨著迭代次數變化的曲線

   下面只給出softmax-L-BFGS的程式碼。

%% STEP 0: Initialise constants and parameters
%
%  Here we define and initialise some constants which allow your code
%  to be used more generally on any arbitrary input.
%  We also initialise some parameters used for tuning the model.

inputSize = 28 * 28; % Size of input vector (MNIST images are 28x28)
numClasses = 10;     % Number of classes (MNIST images fall into 10 classes)
lambda = 1e-4; % Weight decay parameter
itera_num=200;
Jtheta = zeros(itera_num, 1);
a=1;roi=0.5;c=0.6;m=10;
%%======================================================================
%% STEP 1: Load data
images = loadMNISTImages('train-images.idx3-ubyte');
labels = loadMNISTLabels('train-labels.idx1-ubyte');
labels(labels==0) = 10; % Remap 0 to 10
inputData = images;
DEBUG = false;
if DEBUG
    inputSize = 8;
    inputData = randn(8, 100);
    labels = randi(10, 100, 1);
end

theta = 0.005 * randn(numClasses * inputSize, 1);%輸入的是一個列向量
% Randomly initialise theta
theta = reshape(theta, numClasses, inputSize);%將輸入的引數列向量變成一個矩陣

numCases = size(inputData, 2);%輸入樣本的個數
groundTruth = full(sparse(labels, 1:numCases, 1));%這裡sparse是生成一個稀疏矩陣,該矩陣中的值都是第三個值1
%稀疏矩陣的小標由labels和1:numCases對應值構成
cost = 0;

thetagrad = zeros(numClasses, inputSize);
%% STEP 2: Implement softmaxCost

p = weight(theta,inputData);
cost(1) = -1/numCases * groundTruth(:)' * log(p(:)) + lambda/2 * sum(theta(:) .^ 2);
thetagrad = -1/numCases * (groundTruth - p) * inputData' + lambda * theta;
B=eye(numClasses);
H=-inv(B);
d1=H*thetagrad;
theta_new=theta+a*d1;
theta_old=theta;
fprintf('%10s %10s %15s %15s %15s','Iteration','cost','Accuracy');
fprintf('\n');
%% Training
for i=2:itera_num %計算出某個學習速率alpha下迭代itera_num次數後的引數
    a=0.5;
    theta_new=reshape(theta_new, numClasses,inputSize);
    theta_old=reshape(theta_old,numClasses,inputSize);
    p=weight(theta_new,inputData);
    Mp=weight(theta_old,inputData);
    cost(i)=-1/numCases * groundTruth(:)' * log(p(:)) + lambda/2 * sum(theta_new(:) .^ 2);
    thetagrad_new = -1/numCases * (groundTruth - p) * inputData' + lambda * theta_new;
    thetagrad_old = -1/numCases * (groundTruth - Mp) * inputData' + lambda * theta_old;
    thetagrad_new=reshape(thetagrad_new,numClasses*inputSize,1);
    thetagrad_old=reshape(thetagrad_old,numClasses*inputSize,1);
    theta_new=reshape(theta_new,numClasses*inputSize,1);
    theta_old=reshape(theta_old,numClasses*inputSize,1);
    M(:,i-1)=thetagrad_new-thetagrad_old;
    BB(:,i-1)=theta_new-theta_old;
    roiJ(i-1)=1/(M(:,i-1)'*BB(:,i-1));
    gamma=(BB(:,i-1)'*M(:,i-1))/(M(:,i-1)'*M(:,i-1));
    HK=gamma*eye(inputSize*numClasses);
    r=lbfgsloop(i,m,HK,BB,M,roiJ,thetagrad_new);
    d=-r;
    d=reshape(d,numClasses,inputSize);
    theta_new=reshape(theta_new,numClasses,inputSize);
    theta_old=theta_new;
    theta_new = theta_new + a*d;
    %% Test the Accuracy
    images = loadMNISTImages('t10k-images.idx3-ubyte');
    labels = loadMNISTLabels('t10k-labels.idx1-ubyte');
    labels(labels==0) = 10;
    inputDatatest = images;
    pred = zeros(1, size(inputData, 2));
    [nop,pred]=max(theta_new*inputDatatest);
    acc(i-1) = mean(labels(:) == pred(:));
    acc(i-1)=acc(i-1) * 100;
    %%
    fprintf('%5d     %13.4e      %0.3f%%\n',i,cost(i),acc(i-1));
end
plot(0:199, cost(1:200),'b--','LineWidth', 3);
figure
plot(1:199,acc(1:199),'b--','LineWidth',3);
%%

legend('softmax-SD','softmax-L-BFGS');
xlabel('Number of iterations')
ylabel('Cost function')

Weight函式

function Q=weight(theta,Data)
    M = bsxfun(@minus,theta*Data,max(theta*Data, [], 1));
    M = exp(M);
    M= bsxfun(@rdivide,M, sum(M));
    Q=M;
end
 


==============================================================================================

下節講sparse autoencoder softmax分類器,敬請期待!

==============================================================================================

                                                                                   懷柔風光