机器学习算法实践-标准与局部加权线性回归

前言

最近开始总结学习回归相关的东东了，与分类的目标变量是标称型不同，回归是对连续型数据进预测。当然还是从最简单的线性回归开始，本文主要介绍无偏差的标准线性回归和有偏局部加权线性回归的理论基础以及相应的Python实现。

标准线性回归

标准线性回归的理论知识很简单，我们既可以写出它的标量表达式也可以写成矩阵的形式，其中矩阵的形式也可以通过投影矩阵进行推到得到。本部分就对标准线性回归的表达式进行下简单的推导。

给定一组数据其中包括特征矩阵$X$, 目标变量向量$y$:
$$
y = \left[ \begin{matrix}
y_1 \\
y_2 \\
: \\
y_m
\end{matrix} \right]
$$

$$
X = \left[ \begin{matrix}
1 & x_{11} & x_{12} & … & x_{1n} \\
1 & x_{21} & x_{22} & … & x_{2n} \\
: & : & : & … & : & \\
1 & x_{n1} & x_{n2} & … & x_{nn} \\
\end{matrix} \right]
$$
其中$X$第一列为截距项，我们做线性回归是为了得到一个最优回归系数向量$w$使得当我们给定一个$x$能够通过$y=xw$预测$y$的值。其中

$$
w = \left[ \begin{matrix}
w_0 \\
w_1 \\
w_2 \\
: \\
w_n
\end{matrix} \right]
$$

最小二乘法获取回归系数

那么怎样的$w$才是最优的呢？在标准线性回归中我们需要找到是误差最小的$w$, 即预测的$y$值与真实的$y$值之间的差值，为了避免简单累加造成的正负差值相互抵消，这里采用了平方误差:
$$
f(w) = \sum_{i=1}^{m} (y_i - x_{i}^{T}w)^2
$$

对于上述式子$f(w)$可以通过梯度下降等方法得到最优解。但是使用矩阵表示将会是的求解和程序更为简单:
$$
f(w) = (y - Xw)^{T}(y - Xw)
$$
将$f(w)$对$w$求导可得:
$$
\frac{\partial f(w)}{\partial w} = -2X^{T}(y - Xw)
$$
使其等于0，便可得到:
$$
X^{T}(y - Xw) = 0 \\
X^{T}y = X^{T}Xw \\
\hat{w} = (X^{T}X)^{-1}X^{T}y \\
$$

通过投影矩阵获取回归系数

除了通过最小平方差的方法推导得到$w$的表达式，我们还可以通过投影矩阵(Projection Matrix)来得到。

我们知道如果我们能够求得一个$w$使得$Xw = y$肯定是最好的，但是实际情况中$y$一般并不在矩阵$X$的列空间中，也就是此方程无解，于是我们希望通过将向量$y$投影到$X$的列空间中得到投影矩阵$p$, 然后求解$Xw = p$来获取一个最接近的一个解, 矩阵$X$的投影矩阵形式为
$$P = X(X^{T}X)^{-1}X^{T}$$
于是得到$y$在$X$列空间的投影为
$$P\cdot y = X(X^{T}X)^{-1}X^{T}y$$
此时方程$Xw = p$是有解的，得到最接近$Xw = y$的解为:
$$
\hat{w} = (X^{T}X)^{-1}X^{T}y
$$

标准线性回归的Python实现

通过矩阵形式我么可以很方便的通过Numpy的接口进行矩阵运算获取线性回归系数向量$\hat{w}$, 实现如下:

1
2
3
4
5
6

def std_linreg(X, Y):
    xTx = X.T*X
    if np.linalg.det(xTx) == 0:
        print('xTx is a singular matrix')
        return
    return xTx.I*X.T*Y

通过对现有数据进行标准线性回归并绘制回归直线得到如下图(完整代码和数据见: https://github.com/PytLab/MLBox/tree/master/linear_regression)

1

w = [3.00774324, 1.69532264]

相关系数(Correlation Coefficient)计算

如何判断获得的模型预测能力的好坏呢？我们需要计算模型计算得到的$y$的值向量与实际$y$值向量的匹配程度, 也就是计算相关系数Correlation Coefficient。

相关系数的计算公式:
$$
r(X, Y) = \frac{Cov(X, Y)}{\sigma_{X}\sigma_{Y}}
$$
也就是两个数据序列的协方差并除上各自的标准差，本质上就是一种剔除了两个变量量纲影响、标准化后的特殊协方差。

而协方差便是衡量两个变量变化趋势是否相似的一种方法，是同向变化(同时变大或变小)还是反向变化(一个变大一个变小), 同向或者反向的程度如何，计算公式如下:

$$
Cov(X, Y) = E[(X - E[X])(Y - E[Y])] = E[XY] - E[X]E[Y]
$$

通过公式可以看出，如果对于向量中的每个$x, y$同时大于或同时小于各自的期望值，协方差为正，相反则为负。可见如果协方差越大相似程度就越高，协方差越小相似程度就越小。也可以看到如果$X, Y$相同，协方差就是方差，也就是方差是一种特殊情况下的协方差。

关于协方差与相关系数的通俗解释可以参考知乎上的回答:如何通俗易懂地解释「协方差」与「相关系数」的概念？

虽然Numpy中有计算协方差的接口numpy.corrcoef，是分别对两两向量进行比较并计算协方差，得到协方差矩阵。为了练习，我还是稍微自己计算了下协方差并只计算两列不同数据之间的相关系数:

1
2
3
4

def get_corrcoef(X, Y):
    # X Y 的协方差
    cov = np.mean(X*Y) - np.mean(X)*np.mean(Y)
    return cov/(np.var(X)*np.var(Y))**0.5

通过对上面得到的线性回归模型得到的预测的值与实际的值进行相关系数计算可以得到相关系数为

1

Correlation coeffient: 0.9864735622335125

局部加权线性回归(Local Weighted Linear Regression)

上面的数据点是通过公式$y = 3 + 1.7x + 0.1sin(30x)$添加噪声生成的数据，而标准的线性回归是一种无偏差估计，在计算所有点的时候都是无偏差的计算误差并通过优化方法优化误差，如果针对不同的点能够对误差进行调整便可以一定程度上避免标准线性回归带来的欠拟合现象。

也就是引入偏差来降低预测的均方误差，本部分总结下局部加权线性回归的方法。当我们获取某个$x$的预测值的时候，我们需要计算回归系数$w$，但是如果针对样本中的数据，距离$x$越近我们就给个越大的权重，如果距离越远就给一个小的权重，这样就会使得针对$x$的预测值$y_{predict}$能够更贴合样本数据。

当我们需要对数据点$x$相应的目标值进行预测的时候，我们需要给样本中的每个点赋予一个权重值$w_i$(为了区分权重和回归系数，在这里用$\theta$表示回归系数，$w$表示权重), 那么平方误差的表达式就变成:

$$
f(\theta) = \sum_{i=1}^{m}w_{i}(y_i - x_{i}^{T}\theta)^{2}
$$

通过矩阵可以表示成:
$$
f(\theta) = (y - X\theta)^{T}W(y - X\theta)
$$
$f(\theta)$对$\theta$求导等于0得到:
$$
\frac{\partial f(w)}{\partial w} = -2X^{T}W(y - X\theta) = 0 \\
X^{T}Wy = X^{T}WX\theta \\
\theta = (X^{T}WX)^{-1}X^{T}Wy
$$

通过上面的公式，对于任意给定的未知数据可以计算出对应的回归系数$\theta$，并得到相应的预测值$y_{predict}$, 其中$W$是一个对角矩阵，对角线上的元素$w_{ii}$对应样本点$x_i$的权重值。

使用高斯核来赋值权重

那么权重的表达式又是怎样的呢，我们需要距离给定$x$的样本点的权重越高，LWRL使用核来对附近的点赋予更高的权重，最常用的是高斯核函数，对应表达式如下:
$$
w_{ii} = exp(\frac{\vert x_i - x \vert}{-2k^2})
$$
通过公式可以看到如果$x_i$距离$x$的距离越小，$w_{ii}$就会越大，其中参数$k$决定了权重的大小。k越大权重的差距就越小，k越小权重的差距就很大，仅有局部的点参与进回归系数的求取，其他距离较远的权重都趋近于零。如果k去进入无穷大，所有的权重都趋近于1，$W$也就近似等于单位矩阵，局部加权线性回归变成标准的无偏差线性回归，会造成欠拟合的现象；当k很小的时候，距离较远的样本点无法参与回归参数的求取，会造成过拟合的现象。

LWLR的Python实现

本部分对局部加权线性回归进行Python实现，对于给定数据求取相应回归系数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

def lwlr(x, X, Y, k):
    ''' 局部加权线性回归，给定一个点，获取相应权重矩阵并返回回归系数
    '''
    m = X.shape[0]

    # 创建针对x的权重矩阵
    W = np.matrix(np.zeros((m, m)))
    for i in range(m):
        xi = np.array(X[i][0])
        x = np.array(x)
        W[i, i] = exp((np.linalg.norm(x - xi))/(-2*k**2))

    # 获取此点相应的回归系数

    xWx = X.T*W*X
    if np.linalg.det(xWx) == 0:
        print('xWx is a singular matrix')
        return
    w = xWx.I*X.T*W*Y

    return w

我们对上部分使用的数据进行回归并绘制回归曲线:

当k = 0.5, 基本上就是无偏差的标准线性回归

1

Correlation coefficient: 0.9869292425124014

当k = 0.1, 可以较好的反应数据的潜在规律

1

Correlation coefficient: 0.997890402610583

当k = 0.03, 拟合的曲线较多的考虑了噪声数据导致过拟合的现象

1

Correlation coefficient: 0.9995975498285312

总结

本文总结了标准线性回归以及局部加权线性回归的基础知识，并对两张回归方式给与了Python的实现。可见局部加权线性回归在取得适当的$k$，便可以较好的发现数据的内在潜质，但是局部加权线性回归有个缺点就是类似kNN一样，每计算一个点的预测值就需要利用所有数据样本进行计算，如果数据量很大，计算量会是一个问题。

参考

• 《Machine Learning in Action》
• 如何通俗易懂地解释「协方差」与「相关系数」的概念？