什么是多元線性回歸？
在線性回歸分析中，如果有兩個或兩個以上的自變量，就稱為多元線性回歸（multivariable linear regression）。如果我們要預測一套房子的價格，影響房子價格的因素可能包括：面積、卧室數量、層數以及房齡，我們用x1、x2、x3、x4來代表這4個特征（features）。在這里，我們假設房子的價格和這4個特征值都是線性相關的，並且用hθ來表示預測房子的價格，那么有： \[{h_\theta }(x) = {\theta _0} + {\theta _1}{x_1} + {\theta _2}{x_2} + {\theta _3}{x_3} + {\theta _4}{x_4}\]
其中θ為線性回歸中的參數。我們把x、θ分別表示為列向量的形式：
\[\theta = {\left( {{\theta _0},{\theta _1},{\theta _2},{\theta _3},{\theta _4}} \right)^T}\]
\[x = {\left( {1,{x_1},{x_2},{x_3},{x_4}} \right)^T}\]
那么我們可以發現 \[{h_\theta }(x) = {\theta ^T} \bullet x\]
如何實現多元線性回歸？
對於擁有n個特征值、m個數據樣本的數據，我們可以用一個m*(n+1)矩陣的形式表示出來，其中矩陣的每一行（即x(i)）為一個數據樣本，每一列代表一個特征值，xj(i)代表第i個數據樣本中第j個特征值（x0(i)=1）
\[X = \left( {\begin{array}{*{20}{c}}
1&{x_{^1}^{(1)}}&{x_2^{(1)}}& \cdots &{x_n^{(1)}} \\
1&{x_{^1}^{(2)}}&{x_2^{(2)}}& \cdots &{x_n^{(2)}} \\
1&{x_{^1}^{(3)}}&{x_2^{(3)}}& \cdots &{x_n^{(3)}} \\
\vdots & \vdots & \vdots & \cdots & \vdots \\
1&{x_{^1}^{(m)}}&{x_2^{(m)}}& \cdots &{x_n^{(m)}}
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{{x^{(1)}}^T} \\
{{x^{(2)}}^T} \\
{{x^{(3)}}^T} \\
\vdots \\
{{x^{(m)}}^T}
\end{array}} \right)\]
樣本的實際值（對於房價預測來說就是房價）用y來表示：
\[y = {\left( {{y_1},{y_2}, \ldots ,{y_m}} \right)^T}\]
線性回歸參數θ為：
\[\theta = {\left( {{\theta _0},{\theta _1}, \ldots ,{\theta _n}} \right)^T}\]
接下來我們定義代價函數Jθ（cost function）為：
\[{J_\theta }(X) = \frac{1}{{2m}}\sum\limits_{i = 1}^m {{{({h_\theta }({x^{(i)}}) - {y^{(i)}})}^2}} = \frac{1}{{2m}}\sum\limits_{i = 1}^m {{{({\theta ^T} \bullet {x^{(i)}} - {y^{(i)}})}^2}} = \frac{1}{{2m}}\sum\limits_{i = 1}^m {{{({x^{(i)T}} \bullet \theta - {y^{(i)}})}^2}} \]
我們還可以把Jθ利用矩陣表達，獲得更為簡潔的形式
\[{J_\theta }(X) = \frac{1}{{2m}}{(X \bullet \theta - y)^T} \bullet (X \bullet \theta - y)\]
現在，我們只需要讓代價函數Jθ最小，就能得到最優的θ參數。那么，要怎樣才能使Jθ最小呢？有兩個辦法，一個是梯度下降法（gradient descent），一個是標准方程法（norm equation）。
Jθ在樣本數據X確定時，是關於θ的一個函數，而θ是一個（n+1）維的列向量，也就是說Jθ其實是一個（n+1）元函數。學過微積分的人應該知道，要求一個多元函數的極值，可以令多元函數對每個未知元的偏導數為0，從而求出未知元，再代入函數中，便可求出函數的極值。同樣地，我們現在要求Jθ的極值，那么就要先使Jθ關於θ的偏導數，或者說關於θ0、θ1、θ2、...、θn的偏導數都為0。
Jθ對θj的偏導數為：
\[\frac{{\delta J}}{{\delta {\theta _j}}} = \frac{1}{m}\sum\limits_{i = 1}^m {({x^{(i)T}} \bullet \theta - {y^{(i)}})x_j^{(i)}} = \frac{1}{m}\left( {\begin{array}{*{20}{c}}
{x_j^{(1)}}&{x_j^{(2)}}& \cdots &{x_j^{(m)}}
\end{array}} \right) \bullet \left( {X\theta - y} \right)\]
則Jθ對θ的偏導數為：
\[\frac{{\delta J}}{{\delta \theta }} = \left( {\begin{array}{*{20}{c}}
{\frac{{\delta J}}{{\delta {\theta _0}}}} \\
{\frac{{\delta J}}{{\delta {\theta _1}}}} \\
\vdots \\
{\frac{{\delta J}}{{\delta {\theta _n}}}}
\end{array}} \right) = \frac{1}{m}\left( {\begin{array}{*{20}{c}}
1&1& \cdots &1 \\
{x_1^{(1)}}&{x_1^{(2)}}& \cdots &{x_1^{(m)}} \\
\vdots & \vdots &{}& \vdots \\
{x_n^{(1)}}&{x_n^{(2)}}&{}&{x_n^{(m)}}
\end{array}} \right)(X\theta - y) = \frac{1}{m}{X^T}(X\theta - y)\]
標准方程法（Normal Equation）
既然我們知道了Jθ偏導數的表達公式，那么我們只要讓它等於0便可求出我們想要的θ，即令：
\[\begin{gathered}
\frac{{\delta J}}{{\delta \theta }} = \frac{1}{m}{X^T}(X\theta - y) = 0 \hfill \\
X\theta - y = 0 \hfill \\
X\theta = y \hfill \\
{X^T}X\theta = {X^T}y \hfill \\
\theta = {({X^T}X)^{ - 1}}{X^T}y \hfill \\
\end{gathered} \]
注意：在第四步中，兩邊左乘一個XT是為了在左邊獲得一個方陣（XTX），因為只有方陣才有逆矩陣。
在Python中，只需要兩步便可實現這個算法：
1 from numpy.linalg import pinv 2 theta = pinv(X.T.dot(X)).dot(X.T).dot(y)
梯度下降法（Gradient Descent）
標准方程法雖然實現起來很簡單，但是在面對十萬甚至是百萬級別的數據樣本量時，求解一個如此大的矩陣的逆需要花費大量的時間，這時候就需要考慮梯度下降法了。梯度下降法的思想是一點點、逐步靠近極小值，通過多次迭代，不斷更新θ的值，讓代價函數Jθ收斂於極小值。迭代的公式為：
\[\theta = \theta - \alpha \frac{{\delta J}}{{\delta \theta }}\]
我們已經知道了如何求Jθ關於θ的偏導數，只要一直迭代下去，直到Jθ等於0或者足夠的小（比如小於1e-5），最后我們便得到了想要的θ值。
以下我用Python實現了用梯度下降法求θ：
1 import numpy as np 2 3 def partial_derivative(X, theta, y): 4 derivative = X.T.dot(X.dot(theta) - y) / X.shape[0] 5 return derivative 6 7 8 def gradient_descent(X, y, alpha=0.1): 9 theta = np.ones(shape=y.shape, dtype=float) 10 partial_derivative_of_J = partial_derivative(X, theta, y) 11 while any(abs(partial_derivative_of_J) > 1e-5): 12 theta = theta - alpha * partial_derivative_of_J 13 partial_derivative_of_J = partial_derivative(X, theta, y) 14 return theta