最近做了好多個資料探勘的小專案，使用並比較了N多演算法，瞭解了很多機器學習的工具，如R語言、Spark機器學習庫、Python、Tensorflow和RapidMiner等等。但是我感覺到自己沒能深入下去，充其量也只是把別人的工具拿來玩玩而已。對演算法本身的優劣及適用範圍不甚了了，更談不上改進優化演算法了。

本著甘當小學生的精神，我最近在網上參加了機器學習牛人Andrew Ng在Coursera上主講的《機器學習》課程（https://www.coursera.org/learn/machine-learning/home/welcome）。這門課程是我見到的最好的入門課程，由淺人深，循序漸進，不需要很深的數學基礎，符合自然的學習規律。剛剛學習了三週，感覺受益頗多。
下面以最簡單的一元線性迴歸為例，簡單介紹一下機器學習的思想及其在Matlab中的實現。

• 1. 基本概念
• 假設函式
  假設模型 (也稱為假設函式，Hypothesis function)，就是根據特徵變數（feature或variable）擬合出目標變數(target variable)的公式或函式。例如在下面的表格中，根據假設函式，即公式 h θ (x) = θ0+ θ 1*x ，我們可以從房屋的面積估算出房屋的總價，這個公式即稱之為假設函式，其中的θ0和θ1稱為引數（parameter）。
  

選擇不同的引數，一元線性迴歸的假設模型也會不同。下圖展示了三個比較簡單的一元線性迴歸模型。

• 代價函式
  代價函式或損失函式（Cost function或Loss function）是用來評價假設模型擬合的精確度。在訓練資料集上，模型的擬合越好，代價函式就越小。在機器學習中，訓練的目的就是要選擇合適的引數（如上圖中的θ0和θ1），讓代價函式達到最小。
  例如，線上性迴歸中，代價函式J(θ0, θ1)的一般定義如下圖所示，可理解為樣本真實輸出值和假設函式估算值之差平方和的平均值。
  

• 梯度下降
  梯度下降（Gradient descent）法是使得代價函式達到最小的經典方法之一。以前我老是用確定性的最小二乘法來求假設模型h θ (x) = θ0+ θ 1*x中的引數θ0和θ1，讓訓練集上的代價函式值最小。其實，代價函式J(θ0, θ1)可以看成兩個變數θ0和θ1組成的函式，將這個函式視覺化的結果如下圖所示。我們學習的目的就是在圖中找到合適的θ0和θ1，讓J(θ0, θ1)接近並達到區域性或全域性最低點（下圖中的紅色圈中的點）。
  

對上圖中梯度下降的一個直觀的解釋是：“比如我們在一座大山上的某處位置，由於我們不知道怎麼下山，於是決定走一步算一步，也就是在每走到一個位置的時候，求解當前位置的梯度，沿著梯度的負方向，也就是當前最陡峭的位置向下走一步，然後繼續求解當前位置梯度，向這一步所在位置沿著最陡峭最易下山的位置走一步。這樣一步步的走下去，一直走到覺得我們已經到了山腳。當然這樣走下去，有可能我們不能走到山腳，而是到了某一個區域性的山峰低處。”(

從幾何上來說，在代價函式J(θ0, θ1)的曲面上，任意給定一個起始點（θ0, θ1），我們可沿著梯度下降的方向，即對代價函式J(θ0, θ1) 中的θ0和θ1分別求偏導數的方向，按照一定的步長（a, learning rate，學習效率）進行不斷迭代，直到θ0, θ1收斂，也就是達到上圖中的區域性或全域性最低點。對一個簡單的線性迴歸，其演算法描述如下：

• 2. 簡單實現
  在Matlab或Octave中，我們可以通過損失函式和梯度下降的原理來簡單實現一元線性迴歸。

假設原始資料如下, x表示城市的人口，y為城市的利潤。我們可以通過線性迴歸，用某個城市的人口來預測該城市的利潤。詳細程式碼如下：

%% 在Matlab或者在Octave中實現梯度下降法，進行一元線性迴歸
%% 初始化
clear ; close all; clc;
%%載入資料
% x表示城市的人口數量
x = [6.1101; 5.5277; 8.5186; 7.0032; 5.8598; 8.3829; 7.4764; 8.5781; 6.4862; 5.0546; 5.7107; 14.1640; 5.7340; 8.4084; 5.6407; 5.3794; 6.3654; 5.1301; 6.4296; 7.0708; 6.1891; 20.2700; 5.4901; 6.3261; 5.5649; 18.9450; 12.8280; 10.9570; 13.1760; 22.2030; 5.2524; 6.5894; 9.2482; 5.8918; 8.2111; 7.9334; 8.0959; 5.6063; 12.8360; 6.3534; 5.4069; 6.8825; 11.7080; 5.7737; 7.8247; 7.0931; 5.0702; 5.8014; 11.7000; 5.5416; 7.5402; 5.3077; 7.4239; 7.6031; 6.3328; 6.3589; 6.2742; 5.6397; 9.3102; 9.4536; 8.8254; 5.1793; 21.2790; 14.9080; 18.9590; 7.2182; 8.2951; 10.2360; 5.4994; 20.3410; 10.1360; 7.3345; 6.0062; 7.2259; 5.0269; 6.5479; 7.5386; 5.0365; 10.2740; 5.1077; 5.7292; 5.1884; 6.3557; 9.7687; 6.5159; 8.5172; 9.1802; 6.0020; 5.5204; 5.0594; 5.7077; 7.6366; 5.8707; 5.3054; 8.2934; 13.3940; 5.4369];
% y代表城市的利潤
y = [ 17.5920;9.1302;13.6620;11.8540;6.8233;11.8860;4.3483;12.0000;6.5987;3.8166;3.2522;15.5050;3.1551;7.2258;0.7162;3.5129;5.3048;0.5608;3.6518;5.3893;3.1386;21.7670;4.2630;5.1875;3.0825;22.6380;13.5010;7.0467;14.6920;24.1470;-.2200;5.9966;12.1340;1.8495;6.5426;4.5623;4.1164;3.3928;10.1170;5.4974;0.5566;3.9115;5.3854;2.4406;6.7318;1.0463;5.1337;1.8440;8.0043;1.0179;6.7504;1.8396;4.2885;4.9981;1.4233;-1.4211;2.4756;4.6042;3.9624;5.4141;5.1694;-0.7428;17.9290;12.0540;17.0540;4.8852;5.7442;7.7754;1.0173;20.9920;6.6799;4.0259;1.2784;3.3411;-.6807;0.2968;3.8845;5.7014;6.7526;2.0576;0.4795;0.2042;0.6786;7.5435;5.3436;4.2415;6.7981;0.9270;0.1520;2.8214;1.8451;4.2959;7.2029;1.9869;0.1445;9.0551;0.6170];
% m是訓練樣本的個數
m = length(y);
%% 對x, y進行作圖
% 資料點以一個十字交叉的符號表示，大小設定為10
plot(x, y, 'rx', 'MarkerSize', 10);
% 設定X軸
xlabel('人口');
% 設定Y軸
ylabel('利潤');

然後，我們可在當前工作目錄下新建一個computeCost.m檔案，其中只包含一個computeCost函式，來計算代價函式的值，程式碼如下：

function J = computeCost(X, y, theta)
%COMPUTECOST Compute cost for linear regression
%   J = COMPUTECOST(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly 
J = 0;

hypothesis = X * theta;
J = sum((hypothesis - y) .^ 2);
J = J/(2*m);

% =========================================================================
end

同樣的，我們可以可在當前工作目錄下新建另外一個gradientDescent.m檔案，其中只包含一個gradientDescent函式，按照迭代次數和設定的alpha（學習效率）來進行梯度下降迭代。最後，我們還可視化了線性迴歸的結果。具體程式碼如下：

function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
%GRADIENTDESCENT Performs gradient descent to learn theta
%   theta = GRADIENTDESCENT(X, y, theta, alpha, num_iters) updates theta by 
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

temp = zeros(2, 1);

for iter = 1:num_iters

    hypothesis = X * theta -  y;
    n = length(theta);
    for j = 1:n
        theta(j) = theta(j) - alpha/m * (hypothesis' * X(:,j));
    end

    % ============================================================
    % Save the cost J in every iteration    
    J_history(iter) = computeCost(X, y, theta);

    %fprintf('the cost is %f\n', J_history(iter));

end

end

寫完這兩個函式後，我們就可以計算代價函式的值，並用梯度下降迭代法進行線性迴歸了。具體見程式碼：

%% =================== 計算代價函式並進行梯度下降迭代 ===================
% 在x向量前面加上一列由m個1組成的向量，構成一個X矩陣
X = [ones(m, 1), x]; % Add a column of ones to x
% 初始化theta引數
theta = zeros(2, 1); % initialize fitting parameters

% 梯度下降的迭代次數為3000次，學習率步長alpha為0.01；如果不收斂，可減少步長。
iterations = 3000;
alpha = 0.01;
fprintf('\n測試代價函式 ...\n')
% 計算初始代價函式的值
J = computeCost(X, y, theta);
fprintf('With theta = [0 ; 0]\nCost computed = %f\n', J);
fprintf('Expected cost value (approx) 32.07\n');

fprintf('\n開始梯度下降迭代 ...\n')
% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent:\n');
fprintf('%f\n', theta);
fprintf('Expected theta values (approx)\n');
fprintf(' -3.6303\n  1.1664\n\n');

% Plot the linear fit
hold on; % keep previous plot visible
plot(X(:,2), X*theta, '-')
legend('Training data', 'Linear regression')
hold off % don't overlay any more plots on this figure

% 根據擬合的線性迴歸曲線，對人口為35,000和70,000的兩個城市的利潤進行預測。
predict1 = [1, 3.5] *theta;
fprintf('For population = 35,000, we predict a profit of %f\n',...
    predict1*10000);
predict2 = [1, 7] * theta;
fprintf('For population = 70,000, we predict a profit of %f\n',...
    predict2*10000);

fprintf('Program paused. Press enter to continue.\n');


%% ============= 視覺化 J(theta_0, theta_1) =============
fprintf('Visualizing J(theta_0, theta_1) ...\n')

% Grid over which we will calculate J
theta0_vals = linspace(-10, 10, 100);
theta1_vals = linspace(-1, 4, 100);

% initialize J_vals to a matrix of 0's
J_vals = zeros(length(theta0_vals), length(theta1_vals));

% Fill out J_vals
for i = 1:length(theta0_vals)
    for j = 1:length(theta1_vals)
      t = [theta0_vals(i); theta1_vals(j)];
      J_vals(i,j) = computeCost(X, y, t);
    end
end


% Because of the way meshgrids work in the surf command, we need to
% transpose J_vals before calling surf, or else the axes will be flipped
J_vals = J_vals';
% Surface plot
figure;
surf(theta0_vals, theta1_vals, J_vals)
xlabel('\theta_0'); ylabel('\theta_1');

% Contour plot
figure;
% Plot J_vals as 15 contours spaced logarithmically between 0.01 and 100
contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20))
xlabel('\theta_0'); ylabel('\theta_1');
hold on;
plot(theta(1), theta(2), 'rx', 'MarkerSize', 10, 'LineWidth', 2);

視覺化的結果顯示如下：

Matlab控制檯的輸出結果為：
測試代價函式 …
With theta = [0 ; 0]
Cost computed = 32.030658
Expected cost value (approx) 32.07

開始梯度下降迭代 …
Theta found by gradient descent:
-3.795429
1.184909
Expected theta values (approx)
-3.6303
1.1664

迭代3000次後的代價函式值=4.366376
For population = 35,000, we predict a profit of 3517.530646
For population = 70,000, we predict a profit of 44989.347534
Program paused. Press enter to continue.
Visualizing J(theta_0, theta_1) …

根據控制檯結果，我們可以知道，在迭代3000次後，θ0和θ1的值大約為-3.6303和1.1664。用最小二乘法，我們會得到θ0和θ1的精確結果為3.8128和1.1867，可見要得到更精確的結果，我們還需要更多次的迭代。