什么是多元线性回归？
在线性回归分析中，如果有两个或两个以上的自变量，就称为多元线性回归（multivariable linear regression）。如果我们要预测一套房子的价格，影响房子价格的因素可能包括：面积、卧室数量、层数以及房龄，我们用x1、x2、x3、x4来代表这4个特征（features）。在这里，我们假设房子的价格和这4个特征值都是线性相关的，并且用hθ来表示预测房子的价格，那么有： [{h_ heta }(x) = { heta _0} + { heta _1}{x_1} + { heta _2}{x_2} + { heta _3}{x_3} + { heta _4}{x_4}]
其中θ为线性回归中的参数。我们把x、θ分别表示为列向量的形式：
[ heta = {left( {{ heta _0},{ heta _1},{ heta _2},{ heta _3},{ heta _4}} ight)^T}]
[x = {left( {1,{x_1},{x_2},{x_3},{x_4}} ight)^T}]
那么我们可以发现 [{h_ heta }(x) = { heta ^T} ullet x]
如何实现多元线性回归？
对于拥有n个特征值、m个数据样本的数据，我们可以用一个m*(n+1)矩阵的形式表示出来，其中矩阵的每一行（即x(i)）为一个数据样本，每一列代表一个特征值，xj(i)代表第i个数据样本中第j个特征值（x0(i)=1）
[X = left( {egin{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}}
ight) = left( {egin{array}{*{20}{c}}
{{x^{(1)}}^T} \
{{x^{(2)}}^T} \
{{x^{(3)}}^T} \
vdots \
{{x^{(m)}}^T}
end{array}}
ight)]
样本的实际值（对于房价预测来说就是房价）用y来表示：
[y = {left( {{y_1},{y_2}, ldots ,{y_m}} ight)^T}]
线性回归参数θ为：
[ heta = {left( {{ heta _0},{ heta _1}, ldots ,{ heta _n}} ight)^T}]
接下来我们定义代价函数Jθ（cost function）为：
[{J_ heta }(X) = frac{1}{{2m}}sumlimits_{i = 1}^m {{{({h_ heta }({x^{(i)}}) - {y^{(i)}})}^2}} = frac{1}{{2m}}sumlimits_{i = 1}^m {{{({ heta ^T} ullet {x^{(i)}} - {y^{(i)}})}^2}} = frac{1}{{2m}}sumlimits_{i = 1}^m {{{({x^{(i)T}} ullet heta - {y^{(i)}})}^2}} ]
我们还可以把Jθ利用矩阵表达，获得更为简洁的形式
[{J_ heta }(X) = frac{1}{{2m}}{(X ullet heta - y)^T} ullet (X ullet heta - 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 { heta _j}}} = frac{1}{m}sumlimits_{i = 1}^m {({x^{(i)T}} ullet heta - {y^{(i)}})x_j^{(i)}} = frac{1}{m}left( {egin{array}{*{20}{c}}
{x_j^{(1)}}&{x_j^{(2)}}& cdots &{x_j^{(m)}}
end{array}}
ight) ullet left( {X heta - y}
ight)]
则Jθ对θ的偏导数为：
[frac{{delta J}}{{delta heta }} = left( {egin{array}{*{20}{c}}
{frac{{delta J}}{{delta { heta _0}}}} \
{frac{{delta J}}{{delta { heta _1}}}} \
vdots \
{frac{{delta J}}{{delta { heta _n}}}}
end{array}}
ight) = frac{1}{m}left( {egin{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}}
ight)(X heta - y) = frac{1}{m}{X^T}(X heta - y)]
标准方程法（Normal Equation）
既然我们知道了Jθ偏导数的表达公式，那么我们只要让它等于0便可求出我们想要的θ，即令：
[egin{gathered}
frac{{delta J}}{{delta heta }} = frac{1}{m}{X^T}(X heta - y) = 0 hfill \
X heta - y = 0 hfill \
X heta = y hfill \
{X^T}X heta = {X^T}y hfill \
heta = {({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θ收敛于极小值。迭代的公式为：
[ heta = heta - alpha frac{{delta J}}{{delta heta }}]
我们已经知道了如何求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