线性回归

补遗: 机器学习手记系列 3 线性回归样例程序

距离 2012 的两三年后(这篇的草稿时间)又过了两三年,这个补遗看起来也烂尾了 -.-

之前在机器学习手记系列 3: 线性回归和最小二乘法后面留了个问题, 也给了结果, 但是当时说好的程序代码并没给出来, 那个手记系列的坑感觉填不上了, 但是已经刨好的小坑还是填上吧

现在已经有很多深度学习框架和教程来教这个,自己也忘得差不多了,就不班门弄斧裸写。推荐看一下 动手学深度学习 http://zh.gluon.ai/index.html,Deep Learning 领域大神 李沐 等人在维护(我能凑不要脸的蹭热度说下这是前百度同事我们还一起吃饭打牌来着么)。刨的小坑就按 线性回归的从零开始实现 http://zh.gluon.ai/chapter_deep-learning-basics/linear-regression-scratch.html 里的做法来实现

先重复下问题

如下式子里不同的阿拉伯数字只是一个符号, 实际表示的可能是其他数字 967621 = 3 797321 = 1 378581 = 4 422151 = 0 535951 = 1 335771 = 0

根据上述式子, 判断下式等于?
565441 = ?

这题的脑筋急转弯版本答案是看每个数字有几个圈,就代表几,这样 1/2/3/4/5/7 都是 0 个圈,6/9 是 1 个圈,8 是 2 个圈,所以最后 565441 里面只有 6 有 1 个圈,答案为 1

按 gluon 上的教程我们也来走一遍,装环境什么的就看 gluon 了,先引入要用的包

from mxnet import autograd, nd

真正做线性回归是没法只用这么一点数据来模拟的,所以我们要先根据真实值来构造一些数据(这里跟 gluon 不一样的是我没有 bias 因子 b,后面也请一并注意)

num_inputs = 9          # 特征数,当前问题里的变量数 1-9
num_examples = 1000     # 样例数,我们会随机生成多少份样例来学习
true_w = nd.array([0, 0, 0, 0, 0, 1, 0, 2, 1])  # 真实值
features = nd.random.normal(scale=1, shape=(num_examples, num_inputs))  # 随机生成数据集
labels = nd.dot(features, true_w)                                   # 数据集对应的结果

初始化模型参数并创建梯度

w = nd.random.normal(scale=0.01, shape=(9, 1))
w.attach_grad()

定义模型,我们就是做的矩阵乘法

def linreg(X, w):
    return nd.dot(X, w)

定义损失函数,用平方损失

def squared_loss(y_hat, y):
    return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2

定义优化算法,用小批量随机梯度下降(因为我们只用了一个大参数 w,所以还是比 gluon 的样例简单)

def sgd(param, lr, batch_size):
    param[:] = param - lr * param.grad / batch_size

训练,取步长 lr 为 0.01,轮次为 1000 轮

def train():
    lr = 0.01
    num_epochs = 1000
    net = linreg
    loss = squared_loss

    for epoch in range(num_epochs):
        with autograd.record():
            l = loss(net(features, w), labels)
        l.backward()
        sgd(w, lr, labels.size)
        train_l = loss(net(features, w), labels)
        if epoch % 100 == 99:
            print("epoch {}, loss {}, w {}".format(epoch + 1, train_l.mean().asnumpy(), w))

验证下结果看看

if __name__ == "__main__":
    train()
    test = nd.array([1, 0, 0, 2, 2, 1, 0, 0, 0])    # 测试集,565441
    print(nd.dot(test, w))

随便跑了一次输出如下,注意模型里每个值的科学计数法的指数

epoch 1000, loss [  5.72006487e-09], w
[[ -6.20802666e-06]
 [  1.62000088e-05]
 [ -1.03610901e-05]
 [  7.82768348e-06]
 [  2.59973749e-05]
 [  9.99964714e-01]
 [  1.86312645e-05]
 [  1.99990368e+00]
 [  1.00001490e+00]]
<NDArray 9x1 @cpu(0)>

[ 1.00002611]
<NDArray 1 @cpu(0)>

忽略精度问题,可以认为符合真实结果

全部代码详见 https://gist.github.com/whusnoopy/af0aa6fd276ace8a7c4d483e586e936d

机器学习手记系列 3: 线性回归和最小二乘法

好几个月没再继续, 挖坑不填是不对的, 还是继续写手记.

线性回归

线性回归一般用来学习一维自变量和一维因变量之间的线性关系. 如果存在一维自变量 x, 同时还有一维因变量 y = f(x), 如果有一堆对不同 x 下 y 的观测值, 即 <x, y> 的观测对, 且如果 x, y 之间存在较明显的线性关系, 可以用 f(x) = a*x + b 这样的方程表示, 则可以用线性回归的方法学习出 a 和 b 的值, 同时估计这个拟合方法的误差 r.

扩展一下, 线性回归也可以指 y = a1*x1 + a2*x2 + … + ak*xk + b 这样多维自变量和一维因变量之间的线性关系 (多项式里自变量的最高幂次都是 1), 同样也可以用回归的方式学出来里面不同的系数和常数项的值. 只有一维自变量的称之为一元线性回归, 否则是多元线性回归.

判断线性回归好坏, 一般就用平方误差和来描述, 其表达为 (f(x1)-y1)2 + (f(x2)-y2)2 + … + (f(xk)-yk)2), 此值如果为 0 则说明自变量和因变量存在完全的线性关系, 否则是近似线性, 越小近似的越好. 这个东西看着有没有一点面熟? 其实就是机器学习手记系列 2: 离线效果评估里最后提到的 MSE 的非平均版本.

解方程法

假如输入样本存在绝对的线性关系, 即最后的误差为 0, 则问题变为解二元一次方程 (一元线性回归里的系数加常数) 或 N+1 元一次方程 (多元线性回归里 N 个自变量的系数加常数). 这没什么好说的, 直接对原输入直接求解就行了, 类似计算方法这样的课本上有的是解法, 做梯度下降或牛顿法乃至矩阵分解都可以解.

梯度下降法

考虑到绝大部分情况下不存在绝对的线性关系, 则问题可以变成怎么求平方误差和的最小值点. 如果是一元线性回归, 目标函数变为 g(a, b) = (f(x1)-y1)2 + (f(x2)-y2)2 + … + (f(xk)-yk)2

我们的目标就是让这个目标函数值最小, 选定一组 <a, b> 的初始值, 然后求其梯度方向, 每次前进一个小步长, 再求梯度前进, 直到目标函数值不再下降, 说明我们已经走到了一个极值点附近, 终止迭代. 对一元线性回归, 梯度方向是对函数求偏导得到的向量方向.

另外需要注意的是, 梯度下降不一定能找到最优解, 可能会在某个局部最优解那陷进去就出不来了.

这部分更详细的推导请见参考资料里 “机器学习中的数学(1)” 一篇, 里面的公式和图做的很赞, 思路也比我清晰.

牛顿法

梯度下降一般遇到的问题迭代步长不好选, 选太小到极值点太慢, 搞太大又会在极值点附近时因为步长太大跳过去了.

牛顿法最大的贡献就是同时给出了梯度方向和迭代步长, 几乎是一步到位的求解. 方法同解方程一样, 对新的损失目标函数求解, 只是一次解可能还不够好, 需要多做几次迭代. 一般梯度下降可能需要上千轮的迭代, 而牛顿法几次迭代就能到极值点了.

最小二乘法

伟大的高斯同学提出并证明了最小二乘法是最好解答, 证明过程略… 直接看维基或百度百科上的原文吧 (数学不好伤不起).

应用

虽然这个方法看起来很简单粗暴, 但是很多时候变化确实就是线性的. 比如在很多论文和工业实践中, 大家认为同等质量的广告或搜索结果, 放在从上到下不同的位置上, 其点击率和位置的关系符合线性关系, 即 ctr(rank) = a*rank + b.

在六月的一次随笔杂记里, 提到了这样的问题:

如下式子里不同的阿拉伯数字只是一个符号, 实际表示的可能是其他数字 967621 = 3 797321 = 1 378581 = 4 422151 = 0 535951 = 1 335771 = 0

根据上述式子, 判断下式等于?
565441 = ?

假设每个式子最后做的都是加法, 并把字符 0~9 映射到 x0~x9, 则统计不同字母出现的次数就可以列线性返程, 可以将第一个式子表示为

x0*0 + x1*1 + x2*1 + x3*0 + x4*0 + x5*0 + x6*2 + x7*1 + x8*0 + x9*1 = 3

其他类推. 对这堆式子求解就可以得到不同数字对应的真实数值, 可以得到 565441 = 1. // 具体代码和方法下次给出

参考资料

// 本文开写于 9.24, 拖到 10.26 才马马虎虎完成, 最后那个题的解也没写, 各种错乱后续再修改或补充吧