什么是多元线性回归？
在线性回归分析中，如果有两个或两个以上的自变量，就称为多元线性回归（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)=0）
\[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