[线性模型] 线性回归（Linear Regression）原理及MATLAB实现

线性回归解决的问题

“线性回归” 试图学得一个通过属性的线性组合来进行预测的函数，以尽可能准确地预测实值输出标记，一般形式为

\[f(\boldsymbol{x})=\boldsymbol{w}^T\boldsymbol{x}+b \tag 1 \]

其中 \(\boldsymbol{x}\) 表示一组属性，长度为 \(n\) 的列向量. \(\boldsymbol{w}=(w_1;w_2;w_3;...;w_n)\) 表示一组参数，长度为 \(n\) 的列向量，每个 \(w_i\) 可以理解为属性 \(x_i\) 的 “重要程度”. \(b\) 也为参数，可以理解为对 \(\boldsymbol{w}^T\boldsymbol{x}\) 的值进行修正. \(f(\boldsymbol{x})\) 表示通过函数算得的预测值.

比如，对同一个地区的父亲和儿子的身高进行记录，把所有父亲的身高作为 \(\boldsymbol{x}\)，所有儿子的身高作为 \(\boldsymbol{y}\) ，通过一些数学方法，解出一组 \(\boldsymbol{w}\) 和 \(b\) 使得对于每一个输入 \(x_i\)，尽可能满足

\[f(x_i) = wx_i+b\approx y_i \tag 2 \]

之后，对于任何一个不在 \(\boldsymbol{x}\) 中的父亲身高，都可以使用函数预测出其对应的儿子身高.

显然，该模型预测的准确程度取决于训练集在坐标轴上的分布，如果大致呈线性的话，使用线性回归可以得到较好的结果. 因此，在一开始，可以输出训练集的坐标图，粗略地看一下，再决定是否使用线性回归训练。

\(\boldsymbol{w}\) 和 \(b\) 的求解过程

常用的解决思路有两种：

1. 使用最小二乘法进行参数估计，求得参数的闭式解；
2. 给定一个合适的学习率，使用梯度下降法，求出参数的数值解。

最小二乘法

首先，要定一个优化目标，考虑式 \((2)\)，回归任务中常使用均方误差作为损失函数，并使其最小化，可以表示为

\[\begin{align} E_{(w,b)} &= \sum\limits_{i=1}^{m}(f(x_i)-y_i)^2 \\ & = \sum\limits_{i=1}^{m}(y_i-wx_i-b)^2 \tag 3\\ (w^*, b^*) &= \arg\min\limits_{(w,b)}E_{(w,b)} \tag 4\\ \end{align} \]

该公式可以理解为，使所有样本到直线上的欧式距离和最小. 求解 \(w\) 和 \(b\) 使 \(E_{(w,b)}\) 最小化的过程，可以使用概率论中参数估计的方法. 将 \(E_{(w,b)}\) 对 \(w\) 和 \(b\) 分别求导，得：

\[\begin{align} \frac{\partial E_{(w,b)}}{\partial w} &= 2\left( \sum\limits_{i=1}^{m}x_i^2- \sum\limits_{i=1}^{m}(y_i-b)x_i\right) \tag 5\\ \frac{\partial E_{(w,b)}}{\partial b} &= 2\left(mb- \sum\limits_{i=1}^{m}(y_i-wx_i)\right) \tag 6\\ \end{align} \]

此处 \(E_{(w,b)}\) 为关于 \(w\) 和 \(b\) 的凸函数，因此当式 \((5)(6)\) 都为 \(0\) 时，得到参数 \(w\) 和 \(b\) 的最优的闭式解

\[\begin{align} w &= \frac{\sum\limits_{i=1}^{m}y_i(x_i-\bar x)}{\sum\limits_{i=1}^{m}x_i^2 - \frac{1}{m}\left(\sum\limits_{i=1}^{m}x_i\right)^2}\ = \frac{\sum\limits_{i=1}^{m}(x_i-\bar x)(y_i-\bar y)}{\sum\limits_{i=1}^{m}(x_i-\bar x)^2} \tag 7\\ b &= \frac{1}{m}\sum\limits_{i=1}^{m}(y_i-wx_i) \tag 8\\ \end{align} \]

设 \(\boldsymbol x_d = (x_1-\bar x, x_2-\bar x, ..., x_m-\bar x), \boldsymbol y_d= (y_1-\bar y, y_2-\bar y, ..., y_m-\bar y)\), 式 \((7)\) 可以改写为

\[w = \frac{\boldsymbol x_d^T \boldsymbol y_d}{\boldsymbol x_d^T \boldsymbol x_d} \tag 9\\ \]

这些表示是为了与 ”多元线性回归“ 中 \(\boldsymbol w\) 的解对应. 接下来，我们讨论求解一开始的线性模型，即属性 \(\boldsymbol x\) 为一个列向量，而不是单一的数，大多数情况下线性规划处理的也是这样的问题。

求解 “多元线性回归” 中参数的过程与上述类似，仍旧采用最小二乘法。令 \(\hat{\boldsymbol w} = (\boldsymbol w;b), \boldsymbol y=(y_1;y_2;...;y_m)\)，将所有的属性向量构造成一个 \(m(n+1)\) 大小的矩阵 \(\boldsymbol X\), 表示为

\[\boldsymbol X = \begin{pmatrix}        x_{11} & x_{12} & \cdots & x_{1n} &1\\        x_{21} & x_{22} & \cdots & x_{2n} &1\\       \vdots & \vdots & \ddots & \vdots & \vdots\\        x_{m1} & x_{m2} & \cdots & x_{mn} &1 \\        \end{pmatrix} \]

类似式 \((3)(4)\) 此时均方误差可以表示为

\[\begin{align} E_{\hat{\boldsymbol w}} &= (\boldsymbol y - \boldsymbol X\hat{\boldsymbol w})^T(\boldsymbol y - \boldsymbol X\hat{\boldsymbol w}) \tag {11}\\ \hat{\boldsymbol w}^* &= \arg\min\limits_{\hat{\boldsymbol w}^*}E_{\hat{\boldsymbol w}} \tag {10} \end{align} \]

将 \(E_{\hat{\boldsymbol w}}\) 对 \(\hat{\boldsymbol w}\) 求偏导，得

\[\begin{align} \frac{\partial E_{\hat{\boldsymbol w}}}{\partial \hat{\boldsymbol w}} &= \frac{\partial \boldsymbol y^T\boldsymbol y}{\partial \hat{\boldsymbol w}} +\frac{\partial \boldsymbol y^T\boldsymbol X \hat{\boldsymbol w}}{\partial \hat{\boldsymbol w}} + \frac{\partial \hat{\boldsymbol w}^T \boldsymbol X^T \boldsymbol y}{\partial \hat{\boldsymbol w}} + \frac{\partial \hat{\boldsymbol w}^T \boldsymbol X^T \boldsymbol X \hat{\boldsymbol w}}{\partial \hat{\boldsymbol w}}\\ & = 0-\boldsymbol X^T \boldsymbol y -\boldsymbol X^T \boldsymbol y +(\boldsymbol X^T \boldsymbol X + \boldsymbol X^T \boldsymbol X) \hat{\boldsymbol w}\\ & = 2\boldsymbol X^T (\boldsymbol X\hat{\boldsymbol w}-\boldsymbol y) \tag {11}\\ \end{align} \]

令其等于 \(0\), 当 \(\boldsymbol X^T \boldsymbol X\) 为满值矩阵或正定矩阵时，解出 \(\hat{\boldsymbol w}\) 为

\[\hat{\boldsymbol w} = (\boldsymbol X^T \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol y \tag{12} \]

与式 \((9)\) 拥有类似的结构. 当 \(\boldsymbol X^T \boldsymbol X\) 不满秩时，通常引入正则化项，决定其中一组 \(\hat{\boldsymbol w}\) 作为最优解.

梯度下降

对于式 \((11)\) 可以采用一种迭代求解的方法，即梯度下降法. 首先随机给定 \(\hat{\boldsymbol w}\) 的一组初始解，可以是零向量. 然后，让 \(\hat{\boldsymbol w}\) 沿着梯度下降最快的方向递减，当偏导数为 \(0\) 时，得到 \(\hat{\boldsymbol w}\) 的最优解. 其中 \(\hat{\boldsymbol w}\) 的迭代式为

\[\hat{\boldsymbol w} \leftarrow \hat{\boldsymbol w} - \alpha\frac{\partial E_{\hat{\boldsymbol w}}}{\partial \hat{\boldsymbol w}}\tag {13} \]

其中 \(\alpha\) 表示为 \(\hat{\boldsymbol w}\)沿着梯度下降最快方向递减的速率，也即学习率. 将式 \((11)\) 代入式 \((13)\) 并令 \(\alpha \leftarrow 2\alpha\)，得

\[\begin{align} \hat{\boldsymbol w} \leftarrow \hat{\boldsymbol w} - \alpha\boldsymbol X^T (\boldsymbol X\hat{\boldsymbol w}-\boldsymbol y)\tag {14}\\ \hat{\boldsymbol w}_j \leftarrow \hat{\boldsymbol w}_j - \alpha\boldsymbol x_j^T (\boldsymbol X\hat{\boldsymbol w}-\boldsymbol y)\tag {15}\\ \end{align} \]

每一次迭代使用训练集的所有数据的算法称为批量梯度下降法，一次迭代的时间复杂度为 \(O(nm)\). 因此，当训练集样本数和为维数都很大时，得出一组解需要大量的时间。针对这样的情况，引入随机梯度下降法.

每次更新 \(\hat{\boldsymbol w}\) 仅使用一个训练样本，然后判断是否收敛. 当有几万条以上的训练样本时，采用随机梯度下降会快很多，其伪代码如下

do until Converge 
    foreach x in X 
        foreach j in w 
            w[j] = w[j] - a*x[j]*(x[i]*w-y)

判断收敛方法一般由两种：

1. \(\hat{\boldsymbol w}\) 中每个参数的变化均小于一个阈值；
2. \(E_{\hat{\boldsymbol w}}\) 的变化小于一个阈值;

使用随机梯度下降，虽然速度快，但是收敛性能可能会不太好. 综合批量梯度下降法和随机梯度下降法，可每次迭代随机从训练集中挑选一定数量的样本训练，即小批量梯度下降法.

MATLAB 实现

代码采用批量梯度下降法，使用随机生成在坐标系中线性分布的 \(100\) 条训练样本.

% 生成随机训练样本，大致为 y=0.7x+12;
kb = [0.7, 12];
x = zeros(100,2);
y = zeros(100,1);

for i = 1:100
    x(i,1) = randi(1000,1);
    x(i,2) = 1;
end
for i = 1:100
    y(i) = kb(1)*x(i)+kb(2) + randi(200,1)-100;
end
hold on;

% w: 权重
% X: 训练集属性
% alpha: 学习率
% eps: 阈值
% m: 训练集样本数
% n: 训练集维数
function w = BDG(X, y, alpha, eps)
    [m,n] = size(X);
    X = [X ones(m,1)];
    n = n+1;
    w = zeros(n,1);
    prew = zeros(n,1);
    while(1)
        flag = 0;
        w(:) = prew(:) - alpha*X(:,:)'*(X(:,:)*prew(:)-y(:));
        for j = 1:n
            if abs(w(j)-prew(j))>eps
                flag = 1;
                break;
            end
       end
       prew = w; 
       if flag==0
           break;
       end
    end
end

% 测试并输出图像
alpha = 0.00000001;
eps = 0.01;
w = BDG(x, y, alpha, eps);
p = zeros(1000,1);
q = zeros(1000,1);
for j = 1:1000
    p(j,1) = j;
    p(j,2) = j;
    q(j) = w(1)*p(j,1)+w(2);
end
plot(x(:,1), y, '*');
hold on;
plot(p(:,1),  q, '-r');

直线如图所示

参考文献

1. 周志华《机器学习》
2. 南瓜书 https://github.com/datawhalechina/pumpkin-book
3. 线性回归与梯度下降算法 https://www.cnblogs.com/eczhou/p/3951861.html