上一個部落格已經講了softmax理論部分，接下來我們就來做個實驗，我們有一些手寫字型圖檔（28*28），訓練樣本（train-images.idx3-ubyte裡面的圖像對應train-labels.idx1-ubyte）和測試樣本（t10k-images.idx3-ubyte裡面的圖檔對應t10k-labels.idx1-ubyte），我們用訓練樣本訓練softmax模型，測試樣本用來做測試。資料和下面講解的程式下載下傳位址：（這裡）
我們首先展示下我們訓練樣本部分的圖檔和label：
1: images = loadMNISTImages('train-images.idx3-ubyte');%得到的images是一個784*60000的矩陣，意思是每一列是一幅28*28的圖像展成了一列， 2: %一共有60000幅圖像。 3: labels = loadMNISTLabels('train-labels.idx1-ubyte'); 4: display_network(images(:,1:100)); % Show the first 100 images 5: disp(labels(1:10)); 圖檔展示： 部分label：
我們發現lable是從0~9的，為了便于操作，我們最好轉化成1~10（後面轉換）。
下面我們進行訓練，首先我們定義一些softmax模型常量：
1: inputSize = 28 * 28; % Size of input vector (MNIST images are 28x28) 2: inputSize =inputSize +1;% softmx的輸入還要加上一維（x0=1），也是θj向量的次元 3: numClasses = 10; % Number of classes (MNIST images fall into 10 classes) 4: lambda = 1e-4; % Weight decay parameter 導入訓練樣本資料
1: images = loadMNISTImages('train-images.idx3-ubyte');%得到的images是一個784*60000的矩陣，意思是每一列是一 2: %幅28*28的圖像展成了一列，一共有60000幅圖像。 3: labels = loadMNISTLabels('train-labels.idx1-ubyte'); 4: labels(labels==0) = 10; % 因為這裡類别是1,2..k，從0開始的，是以這裡把labels中的0映射成10 5: 6: inputData = images; 7: inputData = [ones(1,60000); inputData];%每個樣本都要增加一個x0=1 初始化模型參數：
1: theta = 0.005 * randn(inputSize*numClasses, 1); 接下來也是最重要的一步就是：給定模型參數的情況下，求訓練樣本的softmax的cost function和梯度，即
1: [cost, grad] = softmax_regression_vec(theta,inputData ,labels,lambda ); 接下來我們就要寫softmax_regression_vec函數：
1: function [f,g] = softmax_regression_vec(theta, X,y,lambda ) 2: %下面的n和inputSize指資料有多少維（包括新加的x0=1這一維），也是θj向量的次元 3: %這裡y是1,2....到k，從1開始的 4: m=size(X,2);%X每一列是一個樣本，m是指有m個樣本 5: n=size(X,1); %n指代的前面說了 6: theta=reshape(theta, n, []); %也就是把theta設定成這樣矩陣：有inputSize行也就是n行，每一列是一個θj，有k列。這樣的θ矩陣跟前面理論部分的θ矩陣不一樣，存在
%轉置關系，為什麼這樣呢？這樣這樣的話在後面的reshape和矩陣A（:）這樣的操作，友善，都是按列進行的，還原也友善。是以隻好程式中出現的θ矩陣都是這樣的，k列，跟理論部分的相反。 7: % initialize objective value and gradient. 8: f = 0; 9: g = zeros(size(theta)); 10: h = theta'*X;%h是k行m列的矩陣，見圖1. 圖1
1: a = exp(h); 2: p = bsxfun(@rdivide,a,sum(a)); % sum(a)是一個行向量，每個元素是a矩陣的每一列的和。然後運用bsxfun(@rdivide，，） 3: %是a矩陣的第i列的每個元素除以 sum(a)向量的第i個元素。得到的p矩陣大小和圖1一樣，每個元素如圖2. 圖2
1: c = log(p); %然後我們取log的對數，c矩陣大小和圖1一樣，每個元素如圖3 要注意取以e為底的對數，如果是其他的最後結果也正确，但是在下面梯度驗證 %部分會有一個倍數關系，而不是相等。
圖3
1: i = sub2ind(size(c), y',1:size(c,2)); %y',1:size(c,2)這兩個向量必須同時是行向量或列向量 2: %因為我們接下來每一個樣本xi對應的yi是幾，就去找到p的每一列中，所對應的第幾個元素就是要找的，如圖4.首先使用sub2ind 3: %sub2ind: 在matlab中矩陣是按一列一列的存儲的，比如A=[1 2 3;4 5 6] 4: %那麼A（2）=4，A(3)=2...而這個函數作用就是比如 sub2ind（size（A）,2,1）就是傳回A的第2行第一列的元素存儲的下标，因為 5: %A（2）=4，是以存儲的下标是2，是以這裡傳回2.這裡sub2ind（size（A）,2,1）的2,1也可以換成向量[a1,a2..],[b1,b2..]但是注意 6: %這兩個向量必須同時是行向量或列向量，而不能一個是行向量一個是列向量。是以傳回的 7: %第一個元素是A的第a1行第b1列的元素存儲的下标，傳回的第,二個元素是A的第a2行第b2列的元素存儲的下标...i是一個向量，c（i）得到的 8: %向量的每一個元素就是p中每一列你前面要找的的元素。 圖4
1: values = c(i); 2: f = -(1/m)*sum(values)+ lambda/2 * sum(theta(:) .^ 2); %這個就是cost function 最後求梯度：
1: d = full(sparse(1:m,y,1)); %d為一個稀疏矩陣，有m行k列（k是類别的個數），這個矩陣的（1，y（1））、（2，y（2）） 2: %....(m,y(m))位置都是1。 3: g = (1/m)*X*(p'-d)+ lambda * theta; %這個g和theta矩陣的結構一樣。 4: g=g(:); % 再還原成向量的形式，這裡（：）和reshape都是按列進行的，是以裡面位置并沒有改變。 解釋下這幾行：
我們想求梯度矩陣g，這裡的g和θ=[θ1,θ2,…,θk]矩陣大小size一樣（跟部落格中的θ矩陣存在轉置關系，之所有代碼中這麼做，是因為這樣再把參數矩陣轉成一個向量或轉回去利用g（：）或reshape函數按列比較方面），是n行k列的矩陣。n是θj或一個樣本xi（包括截距1這一維）的次元大小，k是類别個數。m是樣本個數。
我們已經從上一篇部落格知道我們g（u，v）的大小通過
公式求得。
我們知道
就是前面d矩陣（有m行k列） i 行v列的值，
是前面p矩陣 （有m列k行）i 列 v行的值。
我們想用矢量程式設計來求g矩陣：
我們有樣本X（代碼中每一列是一個樣本，也即X為n行m列），那麼g = (1/m).*X*(p'-d)即是。比如，X的第 i 行乘以（p'-d）的第 j 列就是X（i，j）的值。（正是這種行向量乘以列向量是對應元素相乘再相加就完成了公式裡的Σ，這也是矢量程式設計的核心）
其實cost function不含正則項部分我們還可以用-1/m*d（：）*c（：）來實作。
ok,現在我們這個函數寫完了，我們想驗證下，我們寫的這個求導數或着說梯度的這個公式正确不正确，我們還是用之前部落格提到的用求導公式來驗證，因為你求softmax模型某個參數的導數跟你輸入的資料是什麼、多少都沒有關系，是以我們這有用一些簡單的随意寫得資料和label，然後随意取一個參數進行驗證是不是正确，這些程式在前面已經有了，就不進行講解了。
1: % DEBUG = true; % Set DEBUG to true when debugging. 2: DEBUG = false; 3: if DEBUG 4: inputSize = 9; 5: inputData = randn(8, 100); 7: inputData = [ones(1,100);inputData]; 8: labels = randi(10, 100, 1);%從[1,100]中随機生成一個100*1的列向量 10: end 12: % Randomly initialise theta 13: theta = 0.005 * randn(inputSize*numClasses, 1); 1: [cost, grad] = softmax_regression_vec(theta,inputData ,labels,lambda ); 2: 3: if DEBUG 4: numGrad = computeNumericalGradient( @(theta) softmax_regression_vec(theta,inputData ,labels,lambda) ,theta); 5: 6: % Use this to visually compare the gradients side by side 7: disp([numGrad grad]); 8: 9: % Compare numerically computed gradients with those computed analytically 10: diff = norm(numGrad-grad)/norm(numGrad+grad); 11: disp(diff); 12: % The difference should be small. 13: % In our implementation, these values are usually less than 1e-7. 14: 15: % When your gradients are correct, congratulations! 16: end 然後我們看看檢驗結果（部分對比）：
我們發現最終結果為 9.0915e-10，誤差已經足夠小了，說明我們的寫得函數正确。
接下來，我們要通過疊代算法來最小化我們的cost function，我們仍然使用之前的minfunc 程式包來優化：
1: options.maxIter = 100; %疊代100次 2: softmaxModel = softmaxTrain(inputSize, numClasses, lambda, ... 3: inputData, labels, options); 4: %得到的softmaxModel是一個結構體。softmaxModel.optTheta是一個n行k列的參數矩陣，每一列是一個θj 1: function [softmaxModel] = softmaxTrain(inputSize, numClasses, lambda, inputData, labels, options) 2: % options (optional): options 3: % options.maxIter: number of iterations to train for 4: if ~exist('options', 'var') 5: options = struct; 6: end 7: if ~isfield(options, 'maxIter') 8: options.maxIter = 400; 9: end 10: % initialize parameters 11: theta = 0.005 * randn(numClasses* inputSize, 1); 12: 13: % Use minFunc to minimize the function 14: addpath minFunc/ 15: options.Method = 'lbfgs'; % Here, we use L-BFGS to optimize our cost 16: % function. Generally, for minFunc to work, you 17: % need a function pointer with two outputs: the 18: % function value and the gradient. In our problem, 19: % softmaxCost.m satisfies this. 20: minFuncOptions.display = 'on'; 21: 22: [softmaxOptTheta, cost] = minFunc( @(theta) softmax_regression_vec(theta,inputData,labels,lambda),theta, options); 23: 24: % Fold softmaxOptTheta into a nicer format 25: softmaxModel.optTheta = reshape(softmaxOptTheta,inputSize, numClasses ); 26: softmaxModel.inputSize = inputSize; 27: softmaxModel.numClasses = numClasses; 28: 29: end 然後看下結果：
我們發現cost function 的值 收斂在0.2599附近。
接下來來預測我們的預測樣本，看看識别的正确率是多少：
1: images = loadMNISTImages('t10k-images.idx3-ubyte'); 2: labels = loadMNISTLabels('t10k-labels.idx1-ubyte'); 3: labels(labels==0) = 10; % Remap 0 to 10 4: 5: inputData = images; 6: inputData = [ones(1,size(inputData,2)); inputData];%每個樣本都要增加一個x0=1 7: 8: [pred] = softmaxPredict(softmaxModel, inputData); 9: acc = mean(labels(:) == pred(:)); 10: fprintf('Accuracy: %0.3f%%\n', acc * 100); 11: 12: function [pred] = softmaxPredict(softmaxModel, data) 13: 14: theta = softmaxModel.optTheta; %theta是k列，n行的矩陣 15: pred = zeros(1, size(data, 2)); 16: [~,pred]= max(theta'*data);%theta'*data這個矩陣如圖5，某一個樣本softmax最大值與這個矩陣某一列最大 17: %值是等價的，因為每一列除以同一個分母和不除是一樣的，并且exp（.）是增函數，是以隻求裡面的最大值即可。 圖5
最後看一下我們的正确識别率是多少：
（完）
參考：http://deeplearning.stanford.edu/wiki/index.php/Exercise:Softmax_Regression
轉載于:https://www.cnblogs.com/happylion/p/4225830.html