用分别用正规方程, 梯度下降实现多元线性回归. 最后添加L-2正则化项进一步优化实现.

多元线性回归.jpg

线性代数基础

该部分用于回顾一些必要的线性代数基础知识. 并将代码包装成可用的函数方便以后调用.

矩阵的转置:

(A+B)T=AT+BT(AB)T=BTAT(A+B)^T=A^T+B^T\\ (AB)^T=B^TA^T

迭代方法

迭代方法的基石是不动点定理, 完备距离空间中有函数满足:

∣f(x)−f(y)∣<α∣x−y∣|f(x)-f(y)| < \alpha|x-y|

那么ff必然存在不动点, 且对于任意给空间内的初值, 都可以经过有限步骤迭代得到.

所以, 构造一个合理可迭代公式的核心就是构造x=f(x)x=f(x), 再迭代f(x0)f(x_0)

多元线性回归

给定一个样本集XX和一目标集yy,分别有如下形式

X=[x11x12⋯x1mx21x22⋯x2m⋮⋮xn1xn2⋯xnm]X = \left[\begin{matrix} x_{11} & x_{12} & \cdots &x_{1m}\\ x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & &&\vdots\\ x_{n1} & x_{n2} &\cdots & x_{nm}\\ \end{matrix} \right] y=[y1y2⋯yn]Ty = \left[ \begin{matrix} y_1 & y_2 & \cdots &y_n \end{matrix} \right]^T

所谓的线性回归, 实际上就是寻找到某种平面, 使得所有的样本点到这个"平面"的距离之和最短.

这个所谓的距离, 也就是L-2范数, 它定义了向量的长度, 也称为欧几里得距离. 在后文的[正则化项](#### 范数)中会有更加详细的解释, 此处先做大致了解.

∥x∥2=∑i=1nxi2\| x\|_2 = \sqrt{\sum^n_{i=1} x_i^2}

用法实例:

∥x∥2=((∑i=1Nx(i)2)1/2)2;∇∥x∥22=2x\|x\|^2 = ((\sum^N_{i=1}x^{(i)^2})^{1/2})^2;\nabla\|x\|^2_2 = 2x

而我们要寻找的"平面", 实际上就是超平面, 数学意义上的超平面是要经过原点的, 但是直接使用会造成一些误差. 例如在一元线性回归中, 回归函数y=k∗x+by=k*x+b如果去掉截距项, 那么误差大概率是增加的. 但是我们可以改写数据集XX, 使其的维度增加一维:

Xˉ=[1x11x12⋯x1m1x21x22⋯x2m⋮⋮⋮1xn1xn2⋯xnm]\bar X = \left[\begin{matrix} 1 &x_{11} & x_{12} & \cdots &x_{1m}\\ 1&x_{21} & x_{22} & \cdots & x_{2m}\\ \vdots&\vdots & &&\vdots\\ 1&x_{n1} & x_{n2} &\cdots & x_{nm}\\ \end{matrix} \right]

这样, 我们寻找的超平面就可以写为y=θXˉy = \theta\bar X, 其中

θ=[θ0θ1⋯θn]T\boldsymbol{\theta} = \left[ \begin{matrix} \theta_0 & \theta_1 &\cdots &\theta_n \end{matrix} \right]^T

注1: 所有的字母默认都为矩阵(向量), 除法公式中特别标注. 如上式中θ\boldsymbol\theta为矩阵, 其余为标量参数.

注2: 实际上, 大多数文档解释中所谓的超平面 是某一个超平面的仿射集, 但是一般就不做严格区分了. 具体了解可以参阅凸优化 相关书籍.

而损失函数就可以使用范数的形式定义:

L(θ)=12m∥XTθ−y^∥2\mathscr L(\theta) =\frac1{2m}\|X^T\theta -\hat y\|^2

其中, mm为为XX的行数.

观察损失函数, 其实它在某种程度上十分类似于二次函数, 我们可以将其当做二次函数来思考. 要求一个二次函数的最小值, 我们只需要对其求导, 其导数值为零的点就是极小值点, 而二次函数的极小点就是最小值点. 事实上也是同样的做法, 损失函数对θ\theta求梯度, 梯度为零的点一定就是全局最小的"点". 当然, 在高维的情况下, 能够解出最小值的点很可能是一个解集而非单个的向量.

因此, 损失函数对θ\theta求梯度, 得到:

∇θ(L(θ))=1mX(XTθ−y)\nabla_\theta(\mathscr{L}(\theta))= \frac1{m}X(X^T\theta - y)

令其等于零, 反解θ\theta得到:

θ=(XTX)−1XTy\theta = (X^TX)^{-1} X^Ty

这种方法得到的就是解析解, 显然, 这需要XX是非奇异的. 当样本数量很少的时候, 不能使用该方法. 但是这并不意味着我们就无法求出损失函数的最小值, 可以考虑梯度下降法(Gradient Descent).

如[迭代方法](# 迭代方法)中所言, 构造一个迭代函数, 用新计算出来的值覆盖原来的值, 即:

θj+1:=θj−α∇θL(θj)\theta_{j+1} := \theta_j - \alpha\nabla_\theta\mathscr{L}(\theta_j)

其中, α\alpha称为学习率.

在设置了合适的学习率后, 如果不断进行迭代, 那么最终梯度会区趋近于零, θ\theta也就收敛到了最优解附近.

如此一来, 我们只要求解得出梯度, 就能够得到最优解了(很可能不是唯一解)\

核心代码如下:

while ∣∇θL(θ)∣>εθ=θ−α∇θL(θ)\begin{align} &\text{while}\ |\nabla_\theta\mathscr{L}(\theta)|>\varepsilon\\ &\quad\quad\theta=\theta - \alpha \nabla_\theta\mathscr{L}(\theta) \\ \end{align}

使用梯度下降算法实现数值解程序详解:

引入库:

import numpy as np
import time #比运行时间

具体函数实现:

def cal_gradient(theta=[], X=[],y=[]):
    # 计算梯度
    rows , columns = X.shape
    '''
    梯度根据公式
    gradient = (1/n) * X * ( X * theta - y)
    '''
    diff = np.dot(X,theta) - y
    return (1./rows) * np.dot(np.transpose(X),diff)

def MLR_gradient_descent(X_0=[],y = [], tol = 1e-5,alpha = 0.01):
       rows , columns = X_0.shape
    X_add_one  = np.ones((rows,1))
    X = np.append(X_add_one,X_0,axis=1) # 给X_0增加一列1

    # theta 目标矩阵
    theta = np.random.rand(columns+1,1)
    # 某点处的梯度
    gradient = cal_gradient(theta,X,y)

    while not np.all(np.absolute(gradient) <= tol):
        theta = theta - alpha * gradient
        gradient = cal_gradient(theta,X,y)

    return theta

if __name__ == '__main__':

    X_0 = data_in.data_in_X()
    y = data_in.data_in_Y() # 数据读取

    rows, columns = X_0.shape

    X_add_one = np.ones((rows, 1))
    X = np.append(X_add_one, X_0, axis=1)  # 给X_0增加一列1
    begin_time = time.time()

    theta = MLR_gradient_descent(X_0,y)

    end_time = time.time()

    print(np.dot(X,theta)) # 验证结果是否正确
    print('运算时间:',end_time - begin_time.'s')

正则化项

用以上方法进行拟合的时候, 会出现各种各样的问题. 这种拟合方法对数据十分敏感, 可能出现过拟合, 共线性等情况. 而且, 当数据量稀少, 而特性繁多时, 满足最小的θ\theta浩瀚如海. 于是, 我们引入正则化项, 企图减少模型的复杂度. 让结果更加通用和稳定. 但是在使用正则化项之前, 我们不妨先了解一下范数.

范数

范数可以看做是距离概念的推广. 在平面几何中, 已知两点坐标(x1,y1),(x2,y2)(x_1,y_1),(x_2,y_2)我们很容易就知道两点的距离坐标公式为:

d=(x1−x2)2+(y1−y2)2d = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2}

即使放到三维空间中, 已知两个点分别为(x1,y1,z1),(x2,y2,z2)(x_1,y_1,z_1),(x_2,y_2,z_2), 那么它们的距离就是

d=(x1−x1)2+(y1−y2)2+(z1−z2)2d = \sqrt{(x_1-x_1)^2+(y_1-y_2)^2+(z_1-z_2)^2}

那么四维空间呢? 五维甚至更高维? 虽然我们无法想象出这些空间的形状, 但是不难猜想, 这些空间的"距离"和二维三维一样有着一种统一的定义. 它的形式就是:

d=∑i=1n(Ai−Bi)2d = \sqrt{\sum^n_{i=1}(A_i - B_i)^2}

其中A,BA,B都是nn维向量. 不难验证, 一维情况下也是成立的, 说明这种猜想十分合理. 事实上, 只要是满足如下三个条件, 就可以称为是一种"距离"了. 而此处所谓的猜想, 其实也仅仅是L-2范数, 它从属于L-P范数. L-P范数不是某个范数, 而是一组范数, 满足形式:

Lp=∥X∥p=∑i=1nxippX=(x1,x2,⋯ ,xn)\begin{aligned} L_p = \|X\|_p& = \sqrt[p]{\sum^n_{i=1}x_i^p}\\ X = (x_1,&x_2,\cdots,x_n) \end{aligned}

注意, p=0p=0的时候其实并不满足距离的性质(三角不等式), 故不是真正的范数.

当p=2p=2的时候, 就是我们所说的L-2范数啦, 常常也被称为欧几里得范数(Euclidean norm)

在正则化中, 除了使用L-2范数, L-1范数也是常常被用到的, 称为曼哈顿距离(Manhattan Distance)

L-2正则化

回归分析中, 所有的特征都是为了建立求解的系数使用的. 但是显然, 直接使用就会遇到刚才所说的种种问题. 我们引入正则化项来减少这些无意义的特性的影响.

于是, 加上正则化项后, 损失函数变为:

L(θ)=12∑i=1N(f(x(i))−y(i))2+λ2∥θ∥22\mathscr{L}(\theta) = \frac{1}{2}\sum^N_{i=1}(f(x^{(i)}) - y^{(i)})^2 + \frac\lambda2\|\theta\|^2_2

写成范数的形式就是:

Lθ=12∥XTθ−Y∥22+λ2∥θ∥22\mathscr{L}\theta =\frac1{2} \|X^T\theta-Y\|^2_2 + \frac\lambda 2\|\theta\|^2_2

参数λ\lambda控制这模型的复杂度, λ\lambda越小, 那么越接近与最原始的线性回归, 也越容易出现过拟合等问题.

接下来对损失函数求梯度, 得到:

∇θL(θ)=(XTX+λI2)θ−XTY\nabla_{\theta}\mathscr{L}(\theta)=(X^TX + \frac{\lambda I}{2})\theta - X^TY

令其为零, 反解θ\theta得到:

θ=(XTX+λI2)−1XTY\theta=(X^TX+\frac{\lambda I}{2})^{-1}X^TY

当XTXX^TX为非满秩矩阵的时候(行列式为零), 加上一个λ2I\frac{\lambda}2I就可以保证该矩阵可逆.

不过要注意的是, 上述的单位矩阵 II, 其实并不是一个严格的单位矩阵, 它是一个(n+1)×(n+1)(n+1)\times(n+1),形如:

[000⋯0010⋯0001⋯0⋮⋮000⋯1]\left[ \begin{matrix} 0 & 0 &0\cdots &0\\ 0 & 1 &0\cdots &0\\ 0 & 0 &1\cdots &0\\ \vdots & & & \vdots\\ 0 & 0 &0 \cdots&1\\ \end{matrix} \right]

原因也只是第一项系数不需要做缩放.

如果使用梯度下降迭代, 那么在原先的迭代公式

θ:=θ−α⋅∇θLθ\theta := \theta- \alpha \cdot \nabla_\theta \mathscr{L}\theta

上增加正则化项,

θ:=θ−α⋅∇θLθ−αλ⋅θ\theta := \theta- \alpha \cdot \nabla_\theta \mathscr{L}\theta-\alpha \lambda \cdot \theta

然而L2正则化比较难收敛到零,代码中需要添加额外的限制:

iter = 500
Lambda = 0.5
while (not np.all(np.absolute(gradient) <= tol)) and iter > 0:
    theta = theta - alpha * gradient - (1./rows) *alpha*Lambda*theta
    gradient = cal_gradient(theta,X,y)
    iter -= 1

L-1正则化

L-1正则化和L-2正则化相比, 只是改变了正则项的形式. 其损失函数如下:

L(θ)=12∥XTθ−Y∥22+λ∥θ∥1\mathscr{L}(\theta)=\frac12\|X^T\theta-Y\|^2_2+\lambda\|\theta\|_1

但是当我们企图求梯度的时候就会发现, L-1范数是不可微的, 但是存在次梯度:

∇∥X∥1=sign(X)\nabla \|X\|_1 = \text{sign}(X)

因此, 损失函数的梯度可以写为:

∇θL(θ)=λ×sign(θ)+∇θJ(θ)\nabla_{\theta}\mathscr{L}(\theta) = \lambda\times\text{sign}(\theta) + \nabla_{\theta}\mathscr{J}(\theta)

其中J(θ)=12∥XTθ−Y∥\mathscr{J}(\theta)=\frac{1}{2}\|X^T\theta-Y\|

L-1正则化实现只需要略微改动求解析解时的代码即可.

equ1 = np.linalg.inv(np.dot(X_t, X))

改为:

I = Lambda/2 * np.identity(max(rows+1,columns+1))
I[0][0] = 0
equ1 = np.linalg.inv(np.dot(X_t,X) + I )

L1 和 L2 正则化比较

先说结论, L1正则化使拟合趋于稀疏, 而L2正则化使拟合趋于平滑.

观察L1正则化的方法, 对于某一个值而言, 如果它大于零, 那么λ∥θ∥1\lambda \|\theta\|_1就会使其变小; 反之如果小于零, 那么λ∥θ∥1\lambda \|\theta\|_1就会使其变大. 经过若干次迭代后, 它使得更多的特征为零, 也就让拟合结果更加稀疏了.

而L2正则化每次迭代都会将自身乘上一个小于1的权值, 再减去下降的梯度. 一个比较不贴合的解释就是L2正则化使所有的参数都减小, 但是对突出的特征会有更大的惩罚, 因此最后结果就是让拟合变得平滑.