作業說明
Exercise 3，Week 4，使用Octave實現圖片中手寫數字 0-9 的識別，采用兩種方式（1）多分類邏輯回歸（2）多分類神經網絡。對比結果。
（1）多分類邏輯回歸：實現 lrCostFunction 計算代價和梯度。實現 OneVsAll 使用 fmincg 函數進行訓練。使用 OneVsAll 里訓練好的 theta 對 X 的數據類型進行預測，得到平均准確率。
（2）多分類神經網絡：兩層 theta 權重值在 ex3weights 里已提供。參數不需要調，只需要在 predict 里進行矩陣計算，即可得到分類結果。
數據集 ：ex3data1.mat 。手寫數字圖片數據，5000個樣例。每張圖片20px * 20px，也就是一共400個特征。數據集X維度為5000 * 400
ex3weights.mat。神經網絡每一層的權重。
文件清單
ex3.m - Octave/MATLAB script that steps you through part 1
ex3 nn.m - Octave/MATLAB script that steps you through part 2
ex3data1.mat - Training set of hand-written digits
ex3weights.mat - Initial weights for the neural network exercise
submit.m - Submission script that sends your solutions to our servers
displayData.m - Function to help visualize the dataset
fmincg.m - Function minimization routine (similar to fminunc)
sigmoid.m - Sigmoid function
[*] lrCostFunction.m - Logistic regression cost function
[*] oneVsAll.m - Train a one-vs-all multi-class classi er
[*] predictOneVsAll.m - Predict using a one-vs-all multi-class classi er
[*] predict.m - Neural network prediction function
* 為必須要完成的
結論
因為Octave里數組下標從1開始。所以這里將分類結果0用10替代。預測結果中的1-10代表圖片數字為1,2,3,4,5,6,7,8,9,0
矩陣運算 tricky 的地方在於維度對應，哪里需要轉置很關鍵。
多分類邏輯回歸模型
一、在數據集X里隨機選擇100個數字，繪制
displayData.m：
function [h, display_array] = displayData(X, example_width) %DISPLAYDATA Display 2D data in a nice grid % [h, display_array] = DISPLAYDATA(X, example_width) displays 2D data % stored in X in a nice grid. It returns the figure handle h and the % displayed array if requested. % Set example_width automatically if not passed in if ~exist('example_width', 'var') || isempty(example_width) example_width = round(sqrt(size(X, 2))); end % Gray Image colormap(gray); % Compute rows, cols [m n] = size(X); example_height = (n / example_width); % Compute number of items to display display_rows = floor(sqrt(m)); display_cols = ceil(m / display_rows); % Between images padding pad = 1; % Setup blank display display_array = - ones(pad + display_rows * (example_height + pad), ... pad + display_cols * (example_width + pad)); % Copy each example into a patch on the display array curr_ex = 1; for j = 1:display_rows for i = 1:display_cols if curr_ex > m, break; end % Copy the patch % Get the max value of the patch max_val = max(abs(X(curr_ex, :))); display_array(pad + (j - 1) * (example_height + pad) + (1:example_height), ... pad + (i - 1) * (example_width + pad) + (1:example_width)) = ... reshape(X(curr_ex, :), example_height, example_width) / max_val; curr_ex = curr_ex + 1; end if curr_ex > m, break; end end % Display Image h = imagesc(display_array, [-1 1]); % Do not show axis axis image off drawnow; end
ex3.m里的調用
load('ex3data1.mat'); % training data stored in arrays X, y m = size(X, 1); % Randomly select 100 data points to display rand_indices = randperm(m); sel = X(rand_indices(1:100), :); displayData(sel);
運行效果如下：
二、代價函數
注意：θ0不參與正則化。
正則化邏輯回歸的代價函數如下，分為三項：
梯度下降算法如下：
lrCostFunction.m：
function [J, grad] = lrCostFunction(theta, X, y, lambda) m = length(y); % number of training examples temp = theta; temp(1) = 0; % because we don't add anything for j = 0 % 第一項 part1 = -y' * log(sigmoid(X * theta)); % 第二項 part2 = (1 - y)' * log(1 - sigmoid(X * theta)); % 正則項 regTerm = lambda / 2 / m * temp' * temp; J = 1 / m * (part1 - part2) + regTerm; % 梯度 grad = 1 / m * X' *((sigmoid(X * theta) - y)) + lambda / m * temp; grad = grad(:); end
正則化之后的代價函數同時返回代價 J 和 grad 梯度。這里和上個星期實現的二分類邏輯回歸代碼是一樣的。
三、訓練oneVsAll多分類模型
week2邏輯回歸分類的目標是區分兩種類型。寫好代價函數之后，調用一次 fminunc 函數，得到 theta、就可以畫出 boundary。
但是多分類中需要訓練 K 個分類器，所以多了一個 oneVsAll.m。其實就是循環 K 次，得到一個 theta矩陣（row 為 K，column為特征值個數）。
function [all_theta] = oneVsAll(X, y, num_labels, lambda) % Some useful variables m = size(X, 1); n = size(X, 2); % You need to return the following variables correctly all_theta = zeros(num_labels, n + 1); % Add ones to the X data matrix X = [ones(m, 1) X]; initial_theta = zeros(n + 1, 1); options = optimset('GradObj', 'on', 'MaxIter', 50); % Run fmincg to obtain the optimal theta % This function will return theta and the cost for i = 1 : num_labels, [theta] = ... fmincg(@(t)(lrCostFunction(t, X, (y == i), lambda)), ... initial_theta, options); all_theta(i,:) = theta'; end;
end
多分類使用的 Octave 函數是 fmincg 不是 fminunc，fmincg更適合參數多的情況。
注意這里用到了 y == i，這個操作將 y 中每個數據和 i 進行比較。得到的矩陣中(這里是向量)相同的為1，不同的為0。
ex3.m 里的調用
lambda = 0.1; [all_theta] = oneVsAll(X, y, num_labels, lambda);
四、使用OneVsAll分類器進行預測
predictOneVsAll.m
function p = predictOneVsAll(all_theta, X) m = size(X, 1); num_labels = size(all_theta, 1); % You need to return the following variables correctly p = zeros(size(X, 1), 1); % Add ones to the X data matrix X = [ones(m, 1) X]; temp = X * all_theta'; [x, p] = max(temp,[],2); end
將 X 與 theta 相乘得到預測結果，得到一個 5000*10 的矩陣，每行對應一個分類結果（只有一個1，其余為0）。
題目要求返回的矩陣 p 是一個 5000*1 的矩陣。每行是 1-K 的數字。使用 max 函數，得到兩個返回值。第一個 x 是一個全1的向量，沒用。第二個是 1 所在的 index，也就是對應的類別。
ex3.m 中的調用：
pred = predictOneVsAll(all_theta, X); fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100);
三層神經網絡
這里g(z) 使用 sigmoid 函數。
神經網絡中，從上到下的每個原點是 feature 特征 x0, x1, x2...，不是實例。計算過程其實就是 feature 一層一層映射的過程。一層轉換之后，feature可能變多、也可能變少。下一層 i+1層 feature 的個數是通過權重矩陣里當前 θ(i) 的 row 行數來控制。
兩層權重 θ 已經在 ex3weights.mat 里給出。從a1映射到a2權重矩陣 θ1為 25 * 401，從a2映射到a3權重矩陣 θ2為10 * 26。因為最后有10個分類。（這意味着運算的時候要注意轉置）
一、使用OneVsAll分類器進行預測
權重已經在文件中給出了，只需要矩陣相乘，預測就好了
predict.m ：
function p = predict(Theta1, Theta2, X) %PREDICT Predict the label of an input given a trained neural network % p = PREDICT(Theta1, Theta2, X) outputs the predicted label ofxgiven the % trained weights of a neural network (Theta1, Theta2) % Useful values m = size(X, 1); num_labels = size(Theta2, 1); % You need to return the following variables correctly p = zeros(size(X, 1), 1); % ====================== YOUR CODE HERE ====================== % Instructions: Complete the following code to make predictions using % your learned neural network. You should set p to a % vector containing labels between 1 to num_labels. % % Hint: The max function might come in useful. In particular, the max % function can also return the index of the max element, for more % information see 'help max'. If your examples are in rows, then, you % can use max(A, [], 2) to obtain the max for each row. %0 % Add ones to thexdata matrix a1 = [ones(m, 1) X]; %5000x401 a2 = sigmoid(a1 * Theta1'); %5000x401乘以401x25得到5000x25。即把401個feature映射到25 a2 = [ones(m, 1) a2]; %5000x26 a3 = sigmoid(a2 * Theta2'); %5000x26乘以26x10得到5000x10。即把26個feature映射到10 [x,p] = max(a3,[],2); %和上面邏輯回歸多分類一樣，最后使用 max 函數獲得分類結果。 % ========================================================================= end
二、預測結果
em3_nn.m 里的調用
%% Setup the parameters you will use for this exercise input_layer_size = 400; % 20x20 Input Images of Digits hidden_layer_size = 25; % 25 hidden units num_labels = 10; % 10 labels, from 1 to 10 % (note that we have mapped "0" to label 10) load('ex3weights.mat'); pred = predict(Theta1, Theta2, X); fprintf('\nTraining Set Accuracy: %f\n', mean(double(pred == y)) * 100); % To give you an idea of the network's output, you can also run % through the examples one at the a time to see what it is predicting. % Randomly permute examples 隨機獲取數字，顯示預測結果 rp = randperm(m); for i = 1:m % Display fprintf('\nDisplaying Example Image\n'); displayData(X(rp(i), :)); pred = predict(Theta1, Theta2, X(rp(i),:)); fprintf('\nNeural Network Prediction: %d (digit %d)\n', pred, mod(pred, 10)); % Pause with quit option s = input('Paused - press enter to continue, q to exit:','s'); if s == 'q' break end end
運行結果：
准確率 Training Set Accuracy: 97.520000
這里我比較疑惑的一點是。代碼中沒對 0 進行特殊處理，為什么能預測成10。后來一想 整個計算過程都是權重矩陣控制的，給的 weights 文件里已經是調整好的結果了，0 就是會被分類成10。
所以神經網絡最主要的還是調參啊，調 θ 矩陣。
三、關於矩陣相乘
矩陣這里轉置總是出問題。仔細思考一下
（1）假設只有一個實例 x，而且它像圖中一樣是一個列向量。那根據神經網絡的公式，最相符的應該是 θ(i) * x。最后的得到的是 K 維列向量。
（2）但實際情況中，訓練數據通常一行一個實例。如果嚴格按圖和公式去計算，(θ(i) * XT)T 就需要將 X 先轉置過來，最后再轉置回去。應該是下面的代碼，但是明顯能看出比上面的方法麻煩：
a1 = [ones(1, m); X']; %401x5000 a2 = sigmoid(Theta1 * a1); %25x401 乘以 401x5000 得到25x5000 a2 = [ones(1, m); a2]; %26x5000 a3 = sigmoid(Theta2 * a2); %10x26 乘以 26x5000 得到10x5000 a3 = a3'; % 5000x10
（3）所以用 X * θ(i)T 替換 (θ(i) * XT)T
a1 = [ones(m, 1) X]; %5000x401 a2 = sigmoid(a1 * Theta1'); %5000x401乘以401x25得到5000x25。即把401個feature映射到25 a2 = [ones(m, 1) a2]; %5000x26 a3 = sigmoid(a2 * Theta2'); %5000x26乘以26x10得到5000x10。即把26個feature映射到10
https://github.com/madoubao/coursera_machine_learning/tree/master/homework/machine-learning-ex3/ex3