softmax理论及代码解读——UFLDL
前言
看了各種softmax以后迷迷糊糊的,還是研究一下UFLDL的教程穩點。當然還是得參考挺多教程的:UFLDL-softmax 、Softmax的理解與應用 、Logistic 分類器與 softmax分類器 、詳解softmax函數以及相關求導過程 、Exercise:Softmax Regression 。
【UFLDL】網址有所變動:戳這里
其實,在最初比較疑惑的問題是:softmax的回歸和分類到底是不是一個解法,或者它倆是不是一個知識?然后在一些博客中發現,softmax有時候也稱為softmax分類回歸器,而且在機器學習的定義中,回歸問題通常用來預測一個值,如預測房價或者未來的天氣狀況等;分類用于將事物打上一個標簽,通常結果為離散值,但是分類通常建立在回歸上的,最常見的分類方法是邏輯回歸或者叫邏輯分類。詳細請戳回歸(regression)與分類(classification)的區別
還是將logistic回歸和softmax回歸放一起看看。
Logistic回歸
需要注意的是,logistic回歸的中文叫做對數回歸 ,我們經常也稱他為邏輯回歸等
廣義線性模型
上面說了分類通常建立在回歸中,那么這個“建立”方法是什么?周志華老師《機器學習》給出了一個答案說的非常好:“只需要找一個單調可微函數將分類任務的真實標記與線性回歸模型的預測值聯系起來,上公式
y=g?1(wTx+b)y=g^{-1}(\mathbf w^T\mathbf x+b) y=g?1(wTx+b)
這樣得到的一個模型yyy稱為"廣義線性模型(generalized linear model)",其中g(?)g{(\cdot)}g(?)是單調可微函數,稱為聯系函數(link function)。對數廣義函數g(?)=ln(?)g{(\cdot)=ln{(\cdot)}}g(?)=ln(?)是廣義線性函數的一個特例
logistic回歸(對數回歸)
出現原因是,在二分類問題中,線性預測z=wTx+bz=\mathbf w^T\mathbf x+bz=wTx+b的輸出值zzz為實值,轉換為二分類最理想的方法是使用階躍函數
y={0,z<00.5,z=01,z>0y=\begin{cases} 0,\quad z<0\\ 0.5,\quad z=0\\ 1,\quad z>0 \end{cases} y=??????0,z<00.5,z=01,z>0?
但是階躍函數明顯是不可導的,需要采用近似的函數代替,這樣對數幾率函數就出現了(logistic function)
y=11+e?zy=\frac{1}{1+e^{-z}} y=1+e?z1?
它與我們經常說的sigmoid函數一樣。帶入到廣義線性模型中就得到了對數回歸(logistic regression),說是回歸,其實是分類,而且對數回歸的函數也經常被當做樣本x\mathbf xx類別為1的概率:
P(y=1∣x;w,b)=hθ(x)=11+e?(wTx+b)P(y=1|\mathbf x;\mathbf w,b)=h_\theta(\mathbf x)=\frac{1}{1+e^{-(\mathbf w^T\mathbf x+b)}} P(y=1∣x;w,b)=hθ?(x)=1+e?(wTx+b)1?
而且這個函數上下同時乘以e(wTx+b)e^{(w^Tx+b)}e(wTx+b)就能夠變形成
P(y=1∣x;w,b)=e(wTx+b)1+e(wTx+b)P(y=1|\mathbf x;\mathbf w,b)=\frac{e^{(\mathbf w^T\mathbf x+b)}}{1+e^{(\mathbf w^T\mathbf x+b)}} P(y=1∣x;w,b)=1+e(wTx+b)e(wTx+b)?
據此,在UFLDL中定義了對數回歸的損失函數
J(θ)=?1m[∑i=1my(i)log?hθ(x(i))+(1?y(i))log?(1?hθ(x(i)))]J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m y^{(i)}\log h_\theta(x^{(i)})+(1-y^{(i)})\log (1-h_\theta(x^{(i)})) \right] J(θ)=?m1?[i=1∑m?y(i)loghθ?(x(i))+(1?y(i))log(1?hθ?(x(i)))]
更逼格一點的寫法是
J(θ)=?1m[∑i=1m∑j=011{y(i)=j}log?P(y(i)=j∣x(i);θ)]J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=0}^11\{y^{(i)}=j\}\log P(y^{(i)}=j|x^{(i)};\theta) \right] J(θ)=?m1?[i=1∑m?j=0∑1?1{y(i)=j}logP(y(i)=j∣x(i);θ)]
式中的1{?}1\{\cdot\}1{?}稱為示性函數,大括號里面的值為真的時候表達式值取1,反之取0;其中iii表示第幾個樣本,樣本總數為mmm ,因為是有監督訓練訓練更新嘛,那么y(i)y^{(i)}y(i)就是指第iii個樣本的標簽了,由于是二分類問題,所以y(i)∈{0,1}y^{(i)}\in\{0,1\}y(i)∈{0,1} ,更新辦法就是直接用梯度更新就OK了。
Softmax 回歸
梯度求解
由于softmax就是將簡單的二分類擴充為多分類, 所以只需要對logistic回歸代價函數中的兩類問題改為多分類即可,改法就是觀察logistic回歸逼格寫法中的第二個累加和的范圍是從0到1,那么多分類就是從0到類別總數了,這樣就可以得到softmax的代價函數
J(θ)=?1m[∑i=1m∑j=1k1{y(i)=j}log?P(y(i)=j∣x(i);θ)]J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=1}^k1\{y^{(i)}=j\}\log P(y^{(i)}=j|x^{(i)};\theta) \right] J(θ)=?m1?[i=1∑m?j=1∑k?1{y(i)=j}logP(y(i)=j∣x(i);θ)]
由于softmax 輸出的是對每一個類別估算的概率值,為了保證概率和為1,需要進行歸一化,所以條件概率函數,也就是知道樣本輸入和模型參數的時候,這個樣本屬于各個標簽的概率為:
P(y(i)=l∣x(i);θ)=eθlTx(i)∑j=1keθjTx(i)P(y^{(i)}=l|x^{(i)};\theta)=\frac{e^{\theta_l^Tx^{(i)}}}{\sum_{j=1}^ke^{\theta_j^Tx^{(i)}}} P(y(i)=l∣x(i);θ)=∑j=1k?eθjT?x(i)eθlT?x(i)?
需要注意的是,這里的θ\thetaθ代表全部模型參數,大小是k?(n+1)k*(n+1)k?(n+1)
θ=[θ1T,θ2T,?,θkT]T\theta=\left[\theta_1^T,\theta_2^T,\cdots,\theta_k^T \right]^T θ=[θ1T?,θ2T?,?,θkT?]T
其實這樣一看,就感覺應該是每一個類別都有一組模型參數,而每一組模型參數大小就是輸入維度或者權重維度nnn再加上偏置的維度1,跟簡單的神經網絡連接方法一樣,右邊的每一個標簽輸出與左邊的輸入都有權重連接。
然后便可以得到softmax更具體的代價函數
J(θ)=?1m[∑i=1m∑j=1k1{y(i)=j}log?eθlTx(i)∑j=1keθjTx(i)]J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=1}^k1\{y^{(i)}=j\}\log \frac{e^{\theta_l^Tx^{(i)}}}{\sum_{j=1}^ke^{\theta_j^Tx^{(i)}}}\right] J(θ)=?m1?[i=1∑m?j=1∑k?1{y(i)=j}log∑j=1k?eθjT?x(i)eθlT?x(i)?]
插曲:上式是多類分類情況,如果利用sigmoid激活到(0,1)范圍,獲得最后一層值,如果改成多標簽分類的表達式就是:
loss=∑j=1m∑i=1n?yjilog(yji^)?(1?yji)log(1?yji^)loss=\sum_{j=1}^m\sum_{i=1}^n?y_{ji}log(\hat{y_{ji}})?(1?y_{ji})log(1?\hat{y_{ji}}) loss=j=1∑m?i=1∑n??yji?log(yji?^?)?(1?yji?)log(1?yji?^?)
其中y(i)y^{(i)}y(i)是真實標簽,P(y(i))P(y^{(i)})P(y(i))是預測輸出的概率,通常為經過sigmoid激活函數歸一化到(0,1)后的值。
關于多標簽分類與單標簽分類的區別,戳這里
直接令其對權重或者偏置的偏導數為零并不能求出閉式解(解析解),解析解就是有嚴格的公式可以直接根據輸入得到輸出,就像一元二次方程一樣,有確定的求解公式,這一點可以參考百度百科。那么也就只能找近似解了,稱為數值解,采用逼近的方法求解,一般就是梯度下降法或者變種方法了。
說到這個梯度,可以參考Softmax公式推導 ,或者看看本渣的推導,可能比他那個簡單,但正確與否謝謝大家指正。
求導的時候主要注意一點,是針對某一個類別輸出求導,也就是說對應k個θ\thetaθ 中的一個,所以在對某個權重求解偏導?J(θ)θl\frac{\partial J(\theta)}{\theta_l}θl??J(θ)?的時候可以直接對大括號后面一部分求導,即
?log?eθlTx(i)∑j=1keθj?x(i)?θl=∑j=1keθj?x(i)eθlTx(i)??eθlTx(i)∑j=1keθj?x(i)?θl=∑j=1keθj?x(i)eθlTx(i)?eθlTx(i)?x(i)?∑j=1keθj?x(i)?eθlTx(i)?eθlTx(i)?x(i)(∑j=1keθj?x(i))2=x(i)?(∑j=1keθj?x(i)?eθlTx(i))∑j=1keθj?x(i)=x(i)?(1?P(y(l)=1∣x(i),k,θ))\begin{aligned} &\frac{\partial\log \frac{e^{\theta^T_lx^{(i)}}}{\sum_{j=1}^k e^{\theta_j \cdot x^{(i)}}}}{\partial \theta_l}\\ =& \frac{\sum_{j=1}^k e^{\theta_j \cdot x^{(i)}}}{e^{\theta^T_lx^{(i)}}}\cdot\frac{\partial\frac{e^{\theta^T_lx^{(i)}}}{\sum_{j=1}^k e^{\theta_j \cdot x^{(i)}}}}{\partial\theta_l}\\ =& \frac{\sum_{j=1}^k e^{\theta_j \cdot x^{(i)}}}{e^{\theta^T_lx^{(i)}}}\cdot\frac{e^{\theta^T_lx^{(i)}}\cdot x^{(i)}\cdot \sum_{j=1}^k e^{\theta_j \cdot x^{(i)}}-e^{\theta_l^Tx^{(i)}}\cdot e^{\theta_l^Tx^{(i)}}\cdot x^{(i)}}{(\sum_{j=1}^k e^{\theta_j \cdot x^{(i)}})^2}\\ =&\frac{x^{(i)}\cdot(\sum_{j=1}^k e^{\theta_j \cdot x^{(i)}}-e^{\theta_l^Tx^{(i)}})}{\sum_{j=1}^k e^{\theta_j \cdot x^{(i)}}}\\ =&x^{(i)}\cdot (1-P(y^{(l)}=1|x^{(i)},k,\theta)) \end{aligned} ====??θl??log∑j=1k?eθj??x(i)eθlT?x(i)??eθlT?x(i)∑j=1k?eθj??x(i)???θl??∑j=1k?eθj??x(i)eθlT?x(i)??eθlT?x(i)∑j=1k?eθj??x(i)??(∑j=1k?eθj??x(i))2eθlT?x(i)?x(i)?∑j=1k?eθj??x(i)?eθlT?x(i)?eθlT?x(i)?x(i)?∑j=1k?eθj??x(i)x(i)?(∑j=1k?eθj??x(i)?eθlT?x(i))?x(i)?(1?P(y(l)=1∣x(i),k,θ))?
然后帶入到?J(θ)θl\frac{\partial J(\theta)}{\theta_l}θl??J(θ)? 中,就可以得到UFLDL中的那個代價函數
?J(θ)?θj=?1m[∑i=1mx(i)?(1{y(i)=j}?P(y(i)=j∣x(i);θ))]\frac{\partial J(\theta)}{\partial\theta_j}=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-P(y^{(i)}=j|x^{(i)};\theta))\right] ?θj??J(θ)?=?m1?[i=1∑m?x(i)?(1{y(i)=j}?P(y(i)=j∣x(i);θ))]
梯度下降法的時候就直接使用
θj=θj?α?J(θ)?θj\theta_j=\theta_j-\alpha\frac{\partial J(\theta)}{\partial \theta_j} θj?=θj??α?θj??J(θ)?
權重衰減
先看UFLDL中的一個式子
P(y(i)=j∣x(i);θ)=e(θj?ψ)T?x(i)∑l=1ke(θj?ψ)T?x(i)=e(θj)T?x(i)∑l=1ke(θj)T?x(i)\begin{aligned} P(y^{(i)}=j|x^{(i)};\theta)&=\frac{e^{(\theta_j-\psi)^T\cdot x^{(i)}}}{\sum_{l=1}^k e^{(\theta_j-\psi)^T\cdot x^{(i)}}}\\ &=\frac{e^{(\theta_j)^T\cdot x^{(i)}}}{\sum_{l=1}^k e^{(\theta_j)^T\cdot x^{(i)}}}\\ \end{aligned} P(y(i)=j∣x(i);θ)?=∑l=1k?e(θj??ψ)T?x(i)e(θj??ψ)T?x(i)?=∑l=1k?e(θj?)T?x(i)e(θj?)T?x(i)??
可以看出權重減去一個向量與不減去這個向量得到的結果一樣,說明這個求得的權重參數是冗余的,也就是說最小化代價函數J(θ)J(\theta)J(θ)的解不唯一,如果θ\thetaθ是解,那么θ?ψ\theta-\psiθ?ψ也是一個最優解。這就引出了“權重衰減”,用于解決softmax回歸的參數冗余帶來的數值問題。
新的代價函數為
J(θ)=?1m[∑i=1m∑j=1k1{y(i)=j}log?P(y(i)=j∣x(i);θ)]+12λ∑i=1k∑j=0nθij2J(\theta)=-\frac{1}{m}\left[\sum_{i=1}^m\sum_{j=1}^k1\{y^{(i)}=j\}\log P(y^{(i)}=j|x^{(i)};\theta) \right]+\frac{1}{2}\lambda\sum_{i=1}^k\sum_{j=0}^n\theta_{ij}^2 J(θ)=?m1?[i=1∑m?j=1∑k?1{y(i)=j}logP(y(i)=j∣x(i);θ)]+21?λi=1∑k?j=0∑n?θij2?
這樣加入了權重衰減以后,能夠保證代價函數J(θ)J(\theta)J(θ)是嚴格的凸函數,并且有唯一解。
?J(θ)?θj=?1m[∑i=1mx(i)?(1{y(i)=j}?P(y(i)=j∣x(i);θ))]+λθj\frac{\partial J(\theta)}{\partial\theta_j}=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-P(y^{(i)}=j|x^{(i)};\theta))\right]+\lambda\theta_j ?θj??J(θ)?=?m1?[i=1∑m?x(i)?(1{y(i)=j}?P(y(i)=j∣x(i);θ))]+λθj?
何時使用softmax?
當所有的類別互斥的時候就選擇softmax分類器,比如:古典音樂、鄉村音樂、搖滾樂和爵士樂;但是如果是人聲音樂、舞曲、影視原聲、流行歌曲,那就最好還是用四個二元logistic分類器了,因為一首歌曲可以來源于影視原聲,同時也包含人聲。
代碼解讀
參考UFLDL上的順序編程,文末會給出全部程序。
(1) 四準備
- 手寫數字數據文件
- 讀取數據的MATLAB程序
- 需要補全的softmax分類器程序
- computeNumericalGradient源碼
(2) 初始化參數和讀數據
需要初始化的參數是輸入大小、總類別數目、權重衰減、初始模型參數θ\thetaθ
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 % Randomly initialise theta theta = 0.005 * randn(numClasses * inputSize, 1); theta = reshape(theta, numClasses, inputSize);讀數據有提供實例,主要注意的是將標簽0改為10,讀取數據的時候已經自動歸一化到
addpath ../data/mnist images = loadMNISTImages('train-images-idx3-ubyte'); labels = loadMNISTLabels('train-labels-idx1-ubyte'); labels(labels==0) = 10; % Remap 0 to 10(3) 損失函數的書寫
在官方代碼中,損失函數并未補全,但是提供了實現技巧:
- 計算真實值矩陣(ground truth matrix )
其實就是每個樣本的示性函數,matlab中有自帶的函數可以做到這一點,利用sparse可以將對應樣本的對應標簽標注出來,比如labels是60000×160000\times160000×1的標簽向量,標記了60000個樣本分屬標簽,numCases是總樣本數目60000,那么sparse(labels,1:numCases,1)表示的就是一個60000行的稀疏矩陣,最后一行是(8,60000) 1意思是第60000個樣本對應第八個標簽;為了將稀疏矩陣轉換為10×6000010\times 6000010×60000的單熱度(one hot)編碼,那就需要對稀疏矩陣調用一個full函數,其實計算真實值矩陣就是一句話
groundTruth = full(sparse(labels, 1:numCases, 1));- 計算梯度函數
?J(θ)?θj=?1m[∑i=1mx(i)?(1{y(i)=j}?P(y(i)=j∣x(i);θ))]+λθj=?1m[∑i=1mx(i)?(1{y(i)=j}?e(θj)T?x(i)∑l=1ke(θj)T?x(i)]+λθj\begin{aligned} \frac{\partial J(\theta)}{\partial\theta_j}&=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-P(y^{(i)}=j|x^{(i)};\theta))\right]+\lambda\theta_j\\ &=-\frac{1}{m}\left[\sum_{i=1}^mx^{(i)}\cdot (1\{y^{(i)}=j\}-\frac{e^{(\theta_j)^T\cdot x^{(i)}}}{\sum_{l=1}^k e^{(\theta_j)^T\cdot x^{(i)}}}\right]+\lambda\theta_j \end{aligned} ?θj??J(θ)??=?m1?[i=1∑m?x(i)?(1{y(i)=j}?P(y(i)=j∣x(i);θ))]+λθj?=?m1?[i=1∑m?x(i)?(1{y(i)=j}?∑l=1k?e(θj?)T?x(i)e(θj?)T?x(i)?]+λθj??
式子是分別對連接到每一個輸出標簽的權重求梯度,這里直接用矩陣把所有的權重梯度一次性求出來,先看看參數θ\thetaθ是10×78410\times78410×784的矩陣,data是784×60000784\times 60000784×60000的矩陣,第一維是每張圖片的大小28×2828\times2828×28被拉長為一個向量,groundTruth是10×6000010\times6000010×60000的矩陣,梯度求解三個步驟:求分母,求分數項,求這個損失函數梯度(可以不用人工計算梯度,方法看下一個小步驟)
M = exp(theta*data); M = bsxfun(@rdivide, M, sum(M)); thetagrad = (-1/size(data,2)).*(((groundTruth-(M))*data') + (lambda*theta));順便也可以計算一下損失函數
cost = (-1/size(data,2)).*sum(sum(groundTruth.*log(M))); cost = cost + (lambda/2).*sum(sum(theta.^2));為了后續方便,為這一步驟定義一個函數稱為
function [cost, grad] = softmaxCost(theta, numClasses, inputSize, lambda, data, labels)(4) matlab自動梯度優化
隨后可能與理論部分有點區別了,正常求解是直接用當前權重減去步長和梯度的乘積就是新的模型參數也就是θj=θi?α?θi\theta_j=\theta_i-\alpha\nabla\theta_iθj?=θi??α?θi?,但是呢UFLDL中提供了另一種快速解法,直接使用MATLAB的優化函數minFunc,關于此函數的解釋可以戳這里 ,所以我們可以利用此函數對上面定義的損失函數的最小值,而且minFunc提供了幾個很好地求梯度方法,這里使用大規模優化算法(LBFGS)代替上面人工計算的損失函數的梯度。
if ~exist('options', 'var')options = struct; end if ~isfield(options, 'maxIter')options.maxIter = 400; end % initialize parameters theta = 0.005 * randn(numClasses * inputSize, 1); % Use minFunc to minimize the function % addpath minFunc/ options.Method = 'lbfgs'; % Here, we use L-BFGS to optimize our cost% function. Generally, for minFunc to work, you% need a function pointer with two outputs: the% function value and the gradient. In our problem,% softmaxCost.m satisfies this. options.display = 'on'; % options.displayImage = false; [softmaxOptTheta, cost] = minFunc( @(p) softmaxCost(p, ...numClasses, inputSize, lambda, ...inputData, labels), ... theta, options); % Fold softmaxOptTheta into a nicer format softmaxModel.optTheta = reshape(softmaxOptTheta, numClasses, inputSize); softmaxModel.inputSize = inputSize; softmaxModel.numClasses = numClasses;這里優化的時候會自動優化(3)中輸出項的第一個cost,而非第二項grad
這樣就得到了模型參數了,可以進入測試階段了
(5)測試階段
先讀數據
images = loadMNISTImages('t10k-images-idx3-ubyte'); labels = loadMNISTLabels('t10k-labels-idx1-ubyte'); labels(labels==0) = 10; % Remap 0 to 10 inputData = images;預測標簽
M = exp(theta*inputdata); M = bsxfun(@rdivide, M, sum(M)); [temp pred] = max(M);計算準確率
acc = mean(labels(:) == pred(:)); fprintf('Accuracy: %0.3f%%\n', acc * 100);logistic回歸和softmax關系
其實logistic回歸就是具有兩種類別的softmax, 直接由softmax推導看看(注意softmax具有參數冗余特性, 靠的是權重衰解決的):
1eθ1x+eθ2x[eθ1?xeθ2?x]=1e0?x+e(θ2?θ1)x[e0?xe(θ2?θ1)?x]=11+eθ′?x[1eθ′?x]\begin{aligned} &\frac{1}{e^{\theta_1 x}+e^{\theta_2 x}} \begin{bmatrix} e^{\theta_1 * x } \\ e^{\theta_2* x} \end{bmatrix} \\ =&\frac{1}{e^{0* x}+e^{(\theta_2-\theta_1) x}} \begin{bmatrix} e^{0* x } \\ e^{(\theta_2-\theta_1)* x} \end{bmatrix} \\ =&\frac{1}{1+e^{\theta '*x}} \begin{bmatrix} 1\\ e^{\theta'*x} \end{bmatrix} \end{aligned} ==?eθ1?x+eθ2?x1?[eθ1??xeθ2??x?]e0?x+e(θ2??θ1?)x1?[e0?xe(θ2??θ1?)?x?]1+eθ′?x1?[1eθ′?x?]?
結語
以上便是整個流程了,可運行代碼下載戳此處
給各位大佬們提個建議哈,大部分代碼不是能夠一次性運行成功的,如果是計算機專業,強烈建議自行分析整個程序,嘿嘿,加油。若有錯誤,多多指正
本文已經同步到微信公眾號中,公眾號與本博客將持續同步更新運動捕捉、機器學習、深度學習、計算機視覺算法,敬請關注
總結
以上是生活随笔為你收集整理的softmax理论及代码解读——UFLDL的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: IDEA中配置Project Struc
- 下一篇: 金山词霸2007免费下载