机器学习基础之 Linear Regression

• 1. Linear Regression

  • 1.1 正规方程

  • 1.2 梯度下降

• 2. Ridge Regression

• 3. python实现

• 4. 波士顿房价预测

  • 4.1 加载数据

  • 4.2 线性回归拟合

  • 4.3 模型预测

1. Linear Regression

Linear Regression 属于线性模型，模型假设输入变量 和输出目标值 之间存在线性关系，通常可以表示成如下形式：

下图为一个示例，图中的红色直线为最佳拟合直线。基于给定的数据点 (训练样本)，我们尝试拟合出一条直线，使得数据点到直线的误差平方和最小。

给定 个训练样本，Linear Regression 模型的损失函数定义为：

写成矩阵形式如下：

矩阵 中的每一行表示一个训练样本的特征，通过最小化损失函数 ，求解得到模型参数 。

1.1 正规方程

损失函数对参数 的梯度如下：

令上式为 0，可得参数 的解析解：

1.2 梯度下降

使用梯度下降法求解模型参数 的流程如下：

• 初始化输入：特征矩阵 ，真实输出目标 ，初始权重 ，学习率 ，最大迭代次数 。

• 更新参数：

• 重复步骤 2，直至达到最大迭代次数 。

梯度下降正规方程需要选择学习率 不需要需要多次迭代一次运算得出当特征数量  较大时也能较好适用需要计算 ，当特征数量  较大时运算代价大，因为矩阵逆的计算时间复杂度为 。通常来说，当  时还是可以接受的适用于各种模型只适用于线性回归，不适用于逻辑回归等其他模型

2. Ridge Regression

在线性回归的解析解中需要计算 ，如果训练样本的特征数比样本数还多，即 时， 不是满秩矩阵，此时不可逆。为了解决这个问题，统计学家引入了岭回归的概念。岭回归的损失函数中加入了一个 正则项 (不包含截距项 )：

写成矩阵形式如下：

其中 ，， 中不包含截距项 。

岭回归最早用来处理特征数多于样本数的情况，现在也用于在参数估计中加入偏差，从而得到更好的估计，同时也可以解决多重共线性的问题。岭回归是一种有偏估计。

岭回归损失函数对参数 的梯度如下：

矩阵形式下的损失函数对参数向量 的梯度可以表示为：

矩阵 为 对角矩阵，去除第一行和第一列后为 单位矩阵。

岭回归的模型参数也可以通过正规方程或梯度下降两种方法求解。正规方程解析解为：

使用梯度下降法的求解过程与线性回归类似，参数更新公式如下：

下图从直观上解释了岭回归解与线性回归最小二乘解的区别：

3. python实现

如下代码 (linear_regression.py) 基于 python 实现了一个 LinearRegression 类，支持使用正规方程和梯度下降两种方法求解一般线性回归模型和岭回归模型。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np

def normalize(features):
    """
    Normalizes input features X.

Parameters
    ----------
    features: numpy.ndarray, feature matrix of data points.

Returns
    -------
    features_normalized: numpy.ndarray, normalized feature matrix.
    features_mean: the average value for each feature.
    features_deviation: the standard deviation for each feature.
    """
    # Copy original array to prevent it from changes.
    features_normalized = np.copy(features).astype(float)

# Get average values for each feature (column) in X.
    features_mean = np.mean(features, 0)

# Calculate the standard deviation for each feature.
    features_deviation = np.std(features, 0)

# Subtract mean values from each feature (column) of every example (row)
    # to make all features be spread around zero.
    if features.shape[0] > 1:
        features_normalized -= features_mean

# Normalize each feature values so that all features are close to [-1:1] boundaries.
    # Also prevent division by zero error.
    features_deviation[features_deviation == 0] = 1
    features_normalized /= features_deviation

return features_normalized, features_mean, features_deviation


class LinearRegression:
    def __init__(self, method="normal", reg_lambda=0, fit_intercept=True):
        """
        A simple linear regression algorithm using two methods: 
            a. normal equation
            b. gradient descent
        Optionally, this class can perform ridge regression, if the reg_lambda parameter
        is set to a number greater than zero.

Parameters
        ----------
        method: str, "normal": normal equation, "gradient": gradient descent.
        reg_lambda: int, regularization constant for ridge regression.
        fit_intercept: bool, whether or not to include an intercept term.
        """
        assert(method in ["normal", "gradient"])
        self.method = method
        self.normalize_data = None

assert(reg_lambda >= 0)
        self.reg_lambda = reg_lambda

self.fit_intercept = fit_intercept
        self.theta = None

def fit(self, X, y, lr=0.01, tol=1e-4, max_iters=500, normalize_data=True):
        """
        Fit the regression coefficients

Parameters
        ----------
        X: numpy.ndarray, matrix of data points.
        y: numpy.ndarray, the measured data for each point in X. 
        lr: float, the gradient descent learning rate.
        tol: float, tolerance for stopping criteria.
        max_iters: int, the maximum number of iterations to run the gradient descent solver.
        """
        if self.method == "normal":
            self.normal_equation(X, y)
        else:
            self.normalize_data = normalize_data
            if normalize_data:
                X = normalize(X)[0]
            self.gradient_descent(X, y, lr, tol, max_iters)

def normal_equation(self, X, y):
        """
        Fit the regression coefficients via maximum likelihood.
        """
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]
            A = self.reg_lambda * np.eye(X.shape[1])
            A[0, 0] = 0
        else:
            A = self.reg_lambda * np.eye(X.shape[1])

pseudo_inverse = np.linalg.inv(X.T @ X + A) @ X.T
        self.theta = pseudo_inverse @ y

def gradient_descent(self, X, y, lr=0.01, tol=1e-4, max_iters=500):
        """
        Fit the regression coefficients via gradient descent.
        """
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]

# Calculate the number of training examples.       
        N = X.shape[0]

# Initialization of model parameters.
        self.theta = np.zeros(X.shape[1], )

# Calculate regularization parameter.
        reg_factor = 1 - lr * self.reg_lambda / N

for k in range(max_iters):
            # The difference between predictions and actual values for all m examples.
            delta = X @ self.theta - y

# Create theta shortcut.
            theta = self.theta

# Gradient descent
            theta = theta * reg_factor - lr * (X.T @ delta) / N

# We should not regularize the parameter theta_zero in ridge regression.
            if self.fit_intercept and self.reg_lambda:
                theta[0] = self.theta[0] - lr * np.sum(delta) / N

# Judge whether the stopping criteria is reached.
            if np.linalg.norm(theta - self.theta) < np.sqrt(tol):
                return

# Update model parameters.
            self.theta = theta

def predict(self, X):
        """
        Calculate y_i for each data point in points.

Parameters
        ----------
        X: numpy.ndarray, the data points to calculate with.

Returns
        -------
        y_pred: numpy.ndarray, the model predictions for the points in X.
        """
        if self.normalize_data:
            X = normalize(X)[0]

if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]

return X @ self.theta

4. 波士顿房价预测

4.1 加载数据

from sklearn import datasets
data = datasets.load_boston()
X = data.data
y = data.target
print(X.shape)
print(y.shape)

(506, 13)
(506,)

波士顿房价数据集中一共有 506 个样本点，每个样本包含 13 个特征。

from sklearn.model_selection import train_test_split
train_X, test_X, train_y, test_y = train_test_split(X, y, test_size=0.2)
print(train_X.shape)
print(train_y.shape)
print(test_X.shape)
print(test_y.shape)

(404, 13)
(404,)
(102, 13)
(102,)

进行训练集和测试集划分，训练集中包含 404 个样本点，测试集中包含 102 个样本点。

4.2 线性回归拟合

from sklearn import linear_model
from linear_regression import LinearRegression
sklearn_model = linear_model.LinearRegression(fit_intercept=True)
our_model = LinearRegression(method="normal", fit_intercept=True)
sklearn_model.fit(train_X, train_y)
our_model.fit(train_X, train_y)
print("*************************************************")
print("sklearn 线性回归模型的参数计算为：")
print(sklearn_model.intercept_)
print(sklearn_model.coef_)
print("*************************************************")
print("我们实现的线性回归模型的参数计算为(正规方程方法)：")
print(our_model.theta[0])
print(our_model.theta[1:])

*************************************************
sklearn 线性回归模型的参数计算为：
36.74046795544322
[-1.05325231e-01  4.37831965e-02  4.11878526e-02  3.25713789e+00
 -1.64746096e+01  3.91802570e+00  9.24703495e-04 -1.38576364e+00
  2.86059585e-01 -1.18340743e-02 -1.00610172e+00  6.56588250e-03
 -5.38512638e-01]
*************************************************
我们实现的线性回归模型的参数计算为(正规方程方法)：
36.74046795540561
[-1.05325231e-01  4.37831965e-02  4.11878526e-02  3.25713789e+00
 -1.64746096e+01  3.91802570e+00  9.24703495e-04 -1.38576364e+00
  2.86059585e-01 -1.18340743e-02 -1.00610172e+00  6.56588250e-03
 -5.38512638e-01]

可以看出，在正规方程求解方法下，我们实现的线性回归模型参数计算结果与 sklearn 完全一致。在梯度下降方法中，由于首先需要对训练数据进行标准化，因而计算出的模型参数和 sklearn 没有直接可比性，我们在测试集中比较模型的预测结果。

4.3 模型预测

our_model = LinearRegression(method="gradient", fit_intercept=True)
our_model.fit(train_X, train_y)
print("*************************************************")
print("sklearn 线性回归模型预测结果：")
print(sklearn_model.predict(test_X)[0:10])
print("*************************************************")
print("我们实现的线性回归模型预测结果(梯度下降方法)：")
print(our_model.predict(test_X)[0:10])
print("*************************************************")
print("真实房价结果：")
print(test_y[0:10])

*************************************************
sklearn 线性回归模型预测结果：
[17.18336482 27.30896182 16.88648576 25.33209922 34.01943058 37.60117887
 22.54285594 13.37048284 26.78622285 20.52525983]
*************************************************
我们实现的线性回归模型预测结果(梯度下降方法)：
[15.85690714 24.54861086 16.00965351 23.71145578 30.75878755 34.44602422
 20.81995609 14.36823878 25.71795901 20.43610554]
*************************************************
真实房价结果：
[19.1 22.6 10.2 29.8 27.  37.6 23.2 18.5 22.  21.7]

可以看出，在梯度下降求解方法下，我们实现的线性回归模型预测结果与 sklearn 预测结果以及真实房价结果较为接近，趋势变化一致。