线性回归

本文本文所述内容，参考了机器学习实战：基于Scikit-Learn和TensorFlow一书。
从上往下，本文的代码可完整运行。

线性模型

线性模型就是对输入特征加权求和，再加上一个
称为偏置项（也称为截距项）的常数。

线性回归模型预测
$<mover> y \land </mover> = θ_{0} + θ_{1} x_{1} + θ_{2} x_{2} + \dots + θ_{n} x_{n}$
向量化
$<mover> y \land </mover> = h_{θ} (X) = θ^{T} \cdot X$
线性回归模型的MSE成本函数
$M S T E (X, h_{θ}) = \frac{1}{m} <munderover> \sum i = 1 m </munderover> {(θ^{T} \cdot X^{(i)} - y^{(i)})}^{2}$
也记作MSE(θ)。
标准方程
为了得到是成本函数最小的θ值：
$<mover> θ \land </mover> = (X^{T} \cdot X)^{- 1} \cdot X^{T} \cdot y$

import numpy as np

X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

X_b = np.c_[np.ones((100, 1)), X]       # # add x0 = 1 to each instance
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print(theta_best)   # [[3.99431695][3.18877861]] 可用于预测
import matplotlib.pyplot as plt

plt.scatter(X, y, )
plt.show()

from sklearn.linear_model import LinearRegression

lin_reg = LinearRegression()
lin_reg.fit(X, y)
print(lin_reg.intercept_, lin_reg.coef_)   # 偏置项，权重
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
print(lin_reg.predict(X_new))
#[4.01548197] [[2.95851536]]
#[[4.01548197] [9.93251269]]

计算复杂度
标准方程求逆的矩阵XT·X，是一个n×n矩阵（n是特征数量）。对这种矩阵求逆的计算复杂度通常为O（n^2.4 ）到O（n^3）之间（取决于计算实现）。优点是：相对于训练集中的实例数量（O（m））来说，方程是线性的，所以能够有效地处理大量的训练集，只要内存足够。线性回归模型一经训练（不论是标准方程还是其他算法），预测就非常快速：因为计算复杂度相对于想要预测的实例数量
和特征数量来说，都是线性的。

梯度下降

梯度下降的两个挑战：陷入局部最优而非全局最优，陷入高原区。
线性回归模型MSE成本函数是凸函数，意味着只存在全局最优，且斜率不会陡峭变化。也就是，只要时间足够长，即便是乱走，也可以趋近全局最小值。
左图使用了特征缩放，右图特征值无缩放。训练模型也就是搜寻使成本函数（在训练集上）最小化的参数组合。

成本函数偏导数
$\frac{\partial}{\partial θ_{j}} M S E (θ) = \frac{2}{m} <munderover> \sum i = 1 m </munderover> (θ^{T} \cdot X^{(i)} - y^{(i)}) x_{j}^{(i)}$
成本函数梯度向量
梯度下降步长
$θ (n e x t s t e p) = θ - η \nabla_{θ} M S E (θ)$

eta = 0.1
n_iterations = 1000
m = 100

theta = np.random.randn(2, 1)
for iteration in range(n_iterations):
    gradients = 2 / m * X_b.T.dot(X_b.dot(theta) - y)
    theta = theta - eta * gradients

print(theta)  # [[3.90352217] [2.98694858]]

收敛率： $O (\frac{1}{迭代次数})$
随机梯度下降
在每一步的训练中随机选择一个实例，基于单个实例来计算梯度。有算法的随机性质，比批量梯度下降要不规则的多，成本函数将不再是缓缓降低直到抵达最小值，而是不断上上下下，但是从整体来看，还是在慢慢下降。随着时间推移，最终会非常接近最小值，但是即使它到达了最小值，依旧还会持续反弹，永远不会停止。
随机性的好处在于可以避免局部最优，但无法定位出最小值。要解决这个问题，可以逐步降低学习率。开始时步长较大（快速进展和逃离局部最小值），然后减小，让算法尽量靠近全局最小值。这个过程叫模拟退火。确定每个迭代学习率的函数叫作学习计划。

def learning_schedule(t):
    return t0/(t+t1)

theta=np.random.randn(2,1)      # random
for epoch in range(n_epochs):
    for i in range(m):
        random_index=np.random.randint(m)
        xi=X_b[random_index:random_index+1]
        yi=y[random_index:random_index+1]
        gradients=2*xi.T.dot(xi.dot(theta)-yi)
        eta=learning_schedule(epoch*m+i)
        theta=theta-eta*gradients
print(theta)  # [[3.77856352] [3.19706037]]

也可以使用scikit-learn来实现，

from sklearn.linear_model import SGDRegressor
sgd_reg = SGDRegressor(n_iter=50, penalty=None, eta0=0.1)   # 50 轮，没有使用正则化
sgd_reg.fit(X, y.ravel())

print(sgd_reg.intercept_, sgd_reg.coef_)    # [3.77573813] [3.24655236]

小批量梯度下降
训练基于一小部分随机实例集。可以从矩阵运算的硬件优化中获得显著性能提升（GPU）。在参数空间的层面前进过程不向SGD那样不稳定，最终会比SGD更接近最小值；但另一方面可能更难重局部最小值逃脱。
线性回归算法比较

多项式回归

m = 100
X = 6 * np.random.rand(m, 1) - 3
X = np.sort(X, axis=0)  # axis=0 表示按列排序
y = 0.5 * X ** 2 + X + 2 + np.random.randn(m, 1)

from sklearn.preprocessing import PolynomialFeatures

poly_features = PolynomialFeatures(degree=2, include_bias=False)  # x^2
X_poly = poly_features.fit_transform(X)
print(X[0], X_poly[0])  # [-2.99993479] [-2.99993479 8.99960877] 新加入特征

lin_reg = LinearRegression()
lin_reg.fit(X_poly, y)
print(lin_reg.intercept_, lin_reg.coef_)  # [1.81698516] [[1.06065537 0.57368645]] 差不多的
y_pred = lin_reg.predict(X_poly)
plt.scatter(X, y)
plt.plot(X, y_pred, "g-", label="predict")
plt.legend()
plt.show()

当存在多个特征时，多项式回归能够发现特征和特征之间的关系。在特征a,b，阶数degree=3时，会添加特征征a^2 ， a^3、b^2和b^3,ab、a^2b及ab^2。

学习曲线

曲线绘制的是模型在训练集和验证集上，关于“训练集大小”的性能函数。

def plot_learning_curves(model, X, y):
    X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)
    train_errors, val_errors = [], []
    for m in range(1, len(X_train)):
        model.fit(X_train[:m], y_train[:m])
        y_train_pred = model.predict(X_train[:m])
        y_val_pred = model.predict(X_val)
        train_errors.append(mean_squared_error(y_train_pred, y_train[:m]))
        val_errors.append(mean_squared_error(y_val_pred, y_val))
    plt.plot(np.sqrt(train_errors), 'r-+', linewidth=2, label="train")
    plt.plot(np.sqrt(val_errors), 'b--', linewidth=3, label="val")
    plt.lenged()
    plt.show()


lin_reg = LinearRegression()
plot_learning_curves(lin_reg, X, y)

由图中可以看出，模型拟合不足。两条曲线均达到高地，非常接近，而且相当高。这时，需要使用更复杂的模型或者更好的特征。在相同的数据集上，使用10阶多项式模型：

from sklearn.pipeline import Pipeline

polynomial_regression = Pipeline([
    ('poly_features', PolynomialFeatures(degree=10, include_bias=False)),
    ('sgd_reg', LinearRegression()),
])

plot_learning_curves(polynomial_regression, X, y)

有图中曲线可以看出，过度拟合。可以提供更多的训练数据，知道验证误差接近训练误差。

偏差/方差的权衡

模型的泛化误差可以被表示为以下三个误差之和：

偏差：错误的假设。比如假设数据是线性的，而实际上是二次的。高偏差模型最有可能对训练数据拟合不足。
方差：于模型对训练数据的微小变化过度敏感导致的。具有高自由度的模型（例如高阶多项式模型）很可能也有高方差，所以很容易对训练数据过度拟合。
不可避免的误差：数据本身的噪声所致。减少这部分误差的唯一方法就是清理数据（例如修复数据源，如损坏的传感器，或者是检测
并移除异常值）。

增加模型的复杂度会显著提升模型的方差，减少偏差。反之亦然，所以要做出权衡。

正则线性模型

减少过拟合的方法就是对模型加以约束（正则化）。比如将多项式模型正则化的简单方法就是降低多项式的阶数。
对线性模型来说，正则化通常通过约束模型的权重来实现。有以下三中方法：

岭回归（Ridge Regression，吉洪诺夫正则化）
在成本函数中添加 $α \sum_{i = 1}^{n} θ_{i}^{2}$ 正则项。学习算法不仅需要拟合数据，还需要让模型权重最小。超参数α控制的是对模型进行正则化的程度。如果α=0，则岭回归就是线性模型。如果α非常大，那么所有的权重都将非常接近于
零，结果是一条穿过数据平均值的水平线。
岭回归成本函数：
$J (θ) = M S E (θ) + \frac{1}{2} α \sum_{i = 1}^{n} θ_{i}^{2}$
这里偏置项θ0没有正则化。如果我们将w定义为特征权重的向量（θ1到θn），那么正则项即等于1/2（||w||^ 2）2其中||w||^2为权重向量的l2范数。
在执行岭回归之前，必须对数据进行缩放（例如使StandardScaler），因为它对输入特征的大小非常敏感。
闭式解的岭回归：
$<mover> θ \land </mover> = (X^{T} \cdot X + α A)^{- 1} \cdot X^{T} \cdot y$

from sklearn.linear_model import Ridge

ridge_reg = Ridge(alpha=1, solver='cholesky')   #André-Louis Cholesky的矩阵因式分解法，闭式解
ridge_reg.fit(X, y)
print(ridge_reg.predict([[1.5]]))    #[[4.95870051]]

sgd_reg = SGDRegressor(penalty='l2')     # l2范数正则
sgd_reg.fit(X, y)
print(sgd_reg.predict([[1.5]]))     # [3.7999064]

套索回归（Lasso Regression）
最小绝对收缩和选择算子回归（Least Absolute Shrinkage and Selection Operator Regression，简称Lasso回归，或套索回归）。与岭回归一样，它也是向成本函数增加一个正则项，但是它增加的是权重向量的l1范数，而不是l2范数的平方的一半。
Lasso回归成本函数：
$J (θ) = M S E (θ) + α \sum_{i = 1}^{n} ∣ θ_{i} ∣$
Lasso回归的重要特点是倾向于完全消除最不重要的特征的权重（即设置为0）。如在下图中，α＝10-7看起来像是一个二次的，快要接近于线性：因为所以的高阶多项式特征权重等于0。也即，Lasso回归会自动执行特征选择并输出一个稀疏矩阵（只有很少的特征有非零权重）。
左上图中，背景轮廓（椭圆）表示未正则化MSE成本函数，白色圆点表示该成本函数批梯度下降的路径，前景轮廓（菱形）表示l1正发函数，黄色三角形表示该惩罚函数下，批量梯度下降的路径（α→∞）。此路径下先达到达θ1＝0，然后一路沿轴滚动，直到θ2＝0。右上图中，背景轮廓表示同样的成本函数加上一个α＝0.5的l1惩罚函数。全局最小值位于θ2＝0轴上。批量梯度下降先是到达了θ2
＝0，再沿轴滚动到全局最小值。底部的两张图与上图的含义相同，但是把l1换成了l2惩罚函数。可以看出，正则化后的最小值虽然比未正则化的最小值更接近于θ＝0，但是权重并没有被完全消除。

在Lasso成本函数下，BGD最后的路线似乎在轴上不断上下反弹，这是因为当θ2＝0时，斜率突变。需要降低学习率来保证全局最小值收敛。
当时θi＝0（i＝1，2，…，n），Lasso成本函数是不可微的，但是，当任意θi＝0时，如果使用次梯度向量g作为替代，依旧可以让梯度下降正常运转。公式如下：

from sklearn.linear_model import Lasso

lassp_reg = Lasso(alpha=0.1)
lassp_reg.fit(X, y)
print(lassp_reg.predict([[1.5]]))   # [5.03894038]

弹性网络（Elastic Net）
其正则项就是岭回归和Lasso回归的正则项的混合，混合比例通过r来控制。当r＝0时，弹性网络即等同于岭回归，而当r＝1时，即相当于Lasso回归。弹性网络成本函数：
$J (θ) = M S E (θ + r α \sum_{i = 1}^{n} ∣ θ_{i} ∣ + \frac{1 - r}{2} α \sum_{i = 1}^{n} θ_{i}^{2}$
如何选择上面三种方法呢？默认情况下，可以选用岭回归，但若是实际用到的特征只有少数几个，则要更倾向于Lasso回归或是弹性网络，因为它们会将无用特征的权重降为0。弹性网络一般由于Lasso回归。

from sklearn.linear_model import ElasticNet

elastic_net = ElasticNet(alpha=0.1, l1_ratio=0.5)
elastic_net.fit(X, y)
print(elastic_net.predict([[1.5]]))

早期停止法：验证误差达到最小值时停止训练，避免过度拟合。

from sklearn.base import clone

sgd_reg = SGDRegressor(n_iter=1, warm_start=True, penalty=None)   # warm_start为True会从停下来的地方开始训练
minimum_val_error = float('inf')
best_epoch = None
best_model = None
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2)
for epoch in range(1000):
    sgd_reg.fit(X_train, y_train)
    y_val_pred = sgd_reg.predict(X_val)
    val_error = mean_squared_error(y_val_pred, y_val)
    if val_error < minimum_val_error:
        minimum_val_error = val_error
        best_epoch = epoch
        best_model = clone(sgd_reg)