上一個博客已經講了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