很早之前我就有一个“写一个机器学习的系列笔记”的想法了。之所以叫做“笔记”而不叫“教程”，是因为机器学习只是我的业余爱好之一，我并不精通此道。事实上，我只是学到最简单的皮毛知识，还不一定明白透彻。所以我把我将要写的这个系列定义为笔记，并且在开篇就写得明明白白，这样我就可以随心所欲地写而并不担心误人子弟。不过既然我把笔记放在网上了，那表明我也希望其他人能从我写的东西里得到启发。如果你发现我的理解有误或者文章有错，也请务必留言让我知道。

这个系列的名字叫做 npnet，意思是要用 numpy 去实现一些机器学习的算法。暂定的目标是实现线性模型，CNN， RNN 以及这些基本模型的简单复合（类似 torch.nn.Sequential）。这些算法实现只是为了学习算法背后的理论，并不是为了用于实践。如果你希望在实际应用中使用机器学习，那么你应该去学习现有的库，比如 PyTorch 和 TensorFlow。我是参考 PyTorch 来设计这些算法的。

此系列的算法实现环境是 Python 3.7.3，numpy 的版本是 1.17.2。

# import
import sys
import numpy as np
import matplotlib.pyplot as plt
py = 'Python ' + '.'.join(map(str, sys.version_info[:3]))
print('Jupyter notebook with kernel: {}'.format(py))
print('Numpy version: {}'.format(np.__version__))

Jupyter notebook with kernel: Python 3.7.3
Numpy version: 1.17.2

机器学习算法的三个模块

接下来我们用一个简单的例子——线性回归 (linear regression)——来看一下机器学习算法的主要组成部分。

x = np.arange(0, 5, 0.5)
z = 3 * x - 1
y = z + np.random.randn(len(x))
plt.scatter(x, y, c='r')
plt.plot(x, z)
plt.show()

上图中蓝色直线代表的是一个我们想知道的“真理”。虽然这个“真理”是一条直线，但是由于技术的限制或者人为因素，现实中我们只能观测到一些差不多在一条直线上的红色的点。线性回归的目标是从这些观测点推测出最合理的直线。将这个问题进一步抽象成一个数学问题：给定一些点

$$\{(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\},$$

我们要找到一条最为拟合的直线

$$y = ax + b.$$

上述的抽象包含了机器学习的一个重要概念，模型 (model)。我们需要找的直线就是这个问题中的模型。模型中最为重要的东西就是模型的参数 (parameters)。比如直线由两个参数 $a$ 和 $b$ 确定，要找最合适的直线其实等价于寻找最为合适的参数$a$ 和 $b$。

如何确定模型的参数呢？在我们的例子中，我们如何根据给定的点来找到 $a$ 和 $b$ 呢？在回答这个问题之前，我们还必须需要选择一个标准 (criterion) 来判断参数 $a$ 和 $b$ 的优劣。这是机器学习中第二个重要的概念。以我们的线性回归为例，对于某个观测点 $(x_i, y_i)$，模型给出的值是 $\tilde{y}_i = ax_i + b$。我们希望 $\tilde{y}_i$ 越接近 $y_i$ 越好，换句话说我们希望两者差的绝对值 $|\tilde{y}_i - y_i|$ 越小越好。我们希望所有观测点都能接近模型给出的预测值，因此我们取所有的观测值与预测值之间差的绝对值的平均值来给出一个标准来判断 $a$ 和 $b$ 的好坏：

$$L_1(a, b) = \frac{1}{n} \sum_{i = 1}^n |\tilde{y}_i - y_i| = \frac{1}{n} \sum_{i = 1}^n |ax_i +b - y_i|.$$

$L_1(a, b)$ 的值越小表示 $a$ 和 $b$ 越好。我们希望找到 $a$ 和 $b$ 使得 $L_1(a, b)$ 的值达到最小。

正像参数 $a$ 和 $b$ 有优劣之分，标准也有好坏之分。在微积分里我们学过使用导数去求函数的最值。$L_1(a, b)$ 是一个判断的标准，但是其中包含的绝对值函数 ($y = |x|$) 在 $x = 0$ 处并不可导，所以我们不能简单地使用微积分的知识去找它的最小值。这种情况下，我们不认为 $L_1$ 是一个好的标准。我们知道 $|\tilde{y}_i - y_i|$小等价于 $\frac{1}{2}(\tilde{y}_i - y_i)^2$小，而 $y = \frac{1}{2}x^2$ 是一个处处可导的函数，因此我们有一个更好的标准：

$$L_2(a, b) = \frac{1}{2n} \sum_{i = 1}^n (\tilde{y}_i - y_i)^2 = \frac{1}{2n} \sum_{i = 1}^n (ax_i +b - y_i)^2.$$

$L_2$ 是机器学习中常用的标准，它的名字为均方误差 (mean square error)。公式中出现的 $\frac{1}{2}$ 只是为了方便：$y = \frac{1}{2}x^2$ 的导数是 $y = x$。从这两个标准我们可以大概看出一些门道来。在寻找标准的时候我们一般是先衡量一个观测点与预测值的局部错误，而后取所有局部错误的平均值作为一个衡量模型（或者参数）的标准。这些衡量模型好坏的标准也常被称为损失函数 (loss function)。

现在我们可以根据设定的标准来找到最好的模型了。这个过程是在优化 (optimize) 使得标准给出的错误值最小。优化是机器学习中第三个最重要的组成部分。这个组成部分也经常被称作学习 (learning)，或者训练 (training)。我们可以使用微积分的方法来优化 $L_2(a, b)$。

$$\begin{align*} &\frac{\partial L_2}{\partial a} = \frac{1}{n} \sum_{i = 1}^n (ax_i+b-y_i)x_i = 0, \\ &\frac{\partial L_2}{\partial b} = \frac{1}{n} \sum_{i = 1}^n (ax_i+b-y_i) = 0. \end{align*}$$

由此我们解出

$$\begin{align*} &a = \frac{\left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n y_i\right) - n\left(\sum_{i=1}^n x_iy_i\right)}{\left(\sum_{i=1}^n x_i\right)^2 - n \left(\sum_{i=1}^n x_i^2\right)},\\ &b = \frac{\left(\sum_{i=1}^n x_i\right)\left(\sum_{i=1}^n x_iy_i\right) - \left(\sum_{i=1}^n y_i\right)\left(\sum_{i=1}^n x_i^2\right)}{\left(\sum_{i=1}^n x_i\right)^2 - n \left(\sum_{i=1}^n x_i^2\right)}. \end{align*}$$

我们使用上述公式求出 $a$ 和 $b$ 的值：

n = len(x)
sx = np.sum(x)
sy = np.sum(y)
sxy = np.sum(x * y)
sx2 = np.sum(x ** 2)

a = (sx * sy - n * sxy) / (sx ** 2 - n * sx2)
b = (sx * sxy - sy * sx2) / (sx ** 2 - n * sx2)

print(f'a: {a:.5f}, b: {b:.5f}')

a: 2.90060, b: -1.38144

我们来直观地看一下我们找到的模型：

plt.scatter(x, y, c='r')
plt.plot(x, z, 'b')
plt.plot(x, a * x + b, 'g')
plt.show()

虽然我们找到的模型（绿色的直线）并不是“真理”（蓝色的直线），但是在观测的范围内，我们的模型跟真理非常靠近。值得注意的是，对于两条不同且不平行的直线，只要 $x$ 值足够大，它们的 $y$ 值会相差很远。所以我们在使用机器学习的时候要特别注意模型的适用范围。

回到优化这个问题。如果每个模型的参数我们都能通过微积分的方法给出每个参数关于观测点的公式，那么这是极好的。在残酷复杂的现实中，这种幻想是天真的。只要模型稍微复杂些，求解参数公式就变得不可能了。那么对于复杂的模型，我们有什么办法根据标准去优化它的参数呢？答案还在微积分中：梯度下降法 (gradient descent) 以及链式法则 (chain rule)。

假设 $f(x)$ 是一个多元函数（$x$ 是一个向量）。它的梯度 $\nabla f$ 能指出 $f$ 在每个点增加最快的方向。相对地，$-\nabla f$ 则给出 $f$ 在每个点减小最快的方向。比如在 $x = x_0$ 这个点上，$x$ 要往 $-\nabla f(x_0)$ 这个方向走能使 $f$ 减小最多。梯度下降法就是根据这个想法给出的优化方法：

$$x_0 = x_0 - \alpha \nabla f(x_0).$$

我们从一个随机的点 $x_0$ 出发，重复使用上述的公式迭代 $x_0$ 的值来找到一个使 $f$ 最小的 $x_0$。公式中的 $\alpha$ 是一个大于0的数，代表往 $-\nabla f(x_0)$ 方向走出的距离。其术语是学习率 (learning rate)。

在我们线性回归的例子中，我们想优化的是 $L_2$。它的梯度是 $\nabla L_2 = (\frac{\partial L_2}{\partial a}, \frac{\partial L_2}{\partial b})$。我们之前已经算过这两个偏导数了：

$$\begin{align*} &\frac{\partial L_2}{\partial a} = \frac{1}{n} \sum_{i = 1}^n (ax_i+b-y_i)x_i = \frac{1}{n}\left(\sum_{i=1}^n x_i^2\right)a + \frac{1}{n}\left(\sum_{i=1}^n x_i\right)b - \frac{1}{n}\left(\sum_{i=1}^n x_iy_i\right), \\ &\frac{\partial L_2}{\partial b} = \frac{1}{n} \sum_{i = 1}^n (ax_i+b-y_i) = \frac{1}{n}\left(\sum_{i=1}^n x_i\right)a + b - \frac{1}{n}\left(\sum_{i=1}^n y_i\right). \end{align*}$$

梯度下降法给出的迭代公式是

$$\begin{align*} & a_0 = a_0 - \alpha \frac{\partial L_2}{\partial a}(a_0, b_0), \\ & b_0 = b_0 - \alpha \frac{\partial L_2}{\partial b}(a_0, b_0). \end{align*}$$

我们随机化 $a_0$ 和 $b_0$，使用 $\alpha = 0.1$ 迭代100次。然后我们比较一下用梯度下降法得到的 $a_0$ 和 $b_0$，以及用公式得到的 $a$ 和 $b$。

lr = 0.1
a0, b0 = np.random.rand(2)
for _ in range(100):
    ga = (sx2 * a0 + sx * b0 - sxy) / n
    gb = (sx * a0 - sy) / n + b0
    a0 -= lr * ga
    b0 -= lr * gb
print(f'a0: {a0:.5f}, b0: {b0:.5f}')
print(f'a: {a:.5f}, b: {b:.5f}')

a0: 2.83866, b0: -1.19253
a: 2.90060, b: -1.38144

我们从这个具体的例子看到梯度下降法是求解模型参数的可行性。若是选定了损失函数 $L$，对于每个参数 $a$，我们只需知道 $L$ 相对于 $a$ 的偏导数 $\frac{\partial L}{\partial a}$ 就可以使用梯度下降法了。对于一个复杂的模型，我们需要使用链式法则来求出所有参数的偏导数。链式法则是关于复合函数 (composition function) 的求导方法：

$$\frac{\partial f(g(x))}{\partial x} = \frac{\partial f}{\partial x}(g(x)) \cdot \frac{\partial g}{\partial x}(x)$$

用一个简略的图，一个复合函数是如下图

$$\xrightarrow{x} g(x) \xrightarrow{u} f(x) \xrightarrow{v}$$

以 $x$ 为输入，$v$ 为输出的函数。其中 $u = g(x), v = f(u)$。那么链式法则是说

$$\frac{\partial v}{\partial x} = \frac{\partial f}{\partial x}(u) \cdot \frac{\partial g}{\partial x}(x) $$

机器学习中的复杂模型可以抽象地看成有多个简单模型复合而成的。一个简单的模型可以抽象地看作一个函数 $f(x; a)$。其中 $x$ 代表模型的输入，而分号后边的变量 $a$ 则代表模型的参数。我们来看一个有两个简单模型复合而成的模型：

$$\xrightarrow{x} f(x; a) \xrightarrow{u} g(x; b) \xrightarrow{v} L(x) \xrightarrow{w}$$

如果我们想用梯度下降法，我们必须要知道 $\frac{\partial w}{\partial a}$ $\frac{\partial w}{\partial b}$。根据链式法则，我们得到

$$\begin{align*} &\frac{\partial w}{\partial b} = \frac{\partial L}{\partial x}(v) \cdot \frac{\partial g}{\partial b}(u; b), \\ &\frac{\partial w}{\partial a} = \frac{\partial L}{\partial x}(v) \cdot \frac{\partial g}{\partial x}(u; b) \cdot \frac{\partial f}{\partial a}(x; a). \end{align*}$$

观察上述求偏导数的公式，我们发现在求某个模型参数的偏导数的时候我们需要后面的导数以及模型本身对参数的偏导数。比如 $\frac{\partial w}{\partial b}$ 需要后面的导数 $\frac{\partial L}{\partial x}(v)$ 以及模型本身对参数的偏导数 $\frac{\partial g}{\partial b}(u; b)$；$\frac{\partial w}{\partial a}$ 需要后面的导数 $\frac{\partial L}{\partial x}(v) \cdot \frac{\partial g}{\partial x}(u; b)$ 以及模型本身对参数的偏导数 $\frac{\partial f}{\partial a}(x; a)$。这样的观察让我们能简单地依次算出所有参数的偏导数。这个算法叫做反向传播 (backpropagation)。

总结一下，机器学习有三个最重要的模块：模型、标准以及优化。我们的接下来的算法实现会将这三个模块用三个 Python 类 (class) 表示。

代码复用

最后我们讲一下代码复用的事情。这个系列的博客是使用 Jupyter Notebook 写的，因此我们需要一种可以导入 Jupyter Notebook 的方法来复用之前的代码。我参考官方的文档整合出 import_npnet.py 的文件，其中的代码如下：

%%writefile import_npnet.py
import nbformat
import types
from IPython import get_ipython
from IPython.core.interactiveshell import InteractiveShell
from IPython.utils import io as IPythonIO


def import_npnet(n):
    shell = InteractiveShell.instance()
    mod = types.ModuleType('npnet')
    mod.__dict__['get_ipython'] = get_ipython
    mod.__dict__['npnet'] = mod
    save_user_ns = shell.user_ns
    shell.user_ns = mod.__dict__

    try:
        with IPythonIO.capture_output() as captured:
            for i in range(n):
                nb_path = f'npnet{i}.ipynb'
                nb = nbformat.read(nb_path, 4)
                for cell in nb.cells:
                    if cell.cell_type == 'code':
                        code = shell.input_transformer_manager.transform_cell(
                            cell.source)

                        # import the designated cells only,
                        # those marked with '# import' in the 1st line
                        line = list(code.split('\n'))[0].strip()
                        if line and line[0] == '#' and 'import' in line:
                            exec(code, mod.__dict__)
    finally:
        shell.user_ns = save_user_ns

    return mod

Overwriting import_npnet.py

这个文件中只有一个函数 import_npnet，它接受一个整数 n 作为输入，代表要导入这个系列之前 n 篇笔记的代码。每一篇笔记都是一个 Jupyter Notebook ，名字格式是 npnet{i}.ipynb，其中 {i} 是笔记的编号。比如本笔记保存在 npnet0.ipynb 。对于一篇笔记，如果某个代码 Cell 的第一行是 Python 注释，而且包含 import 这个词，那么这个代码 Cell 将会被导入。这种设定是非常有必要的：以后我们要测试我们的算法，而机器学习的算法一般运行的时间都挺长的，我们希望 import_npnet 能跳过这些耗时测试代码。这个函数会返回由这些代码生成的一个模块。一般来说，n 的值为当前笔记的编号。比如说当前是编号为0的笔记，我可以使用一下代码达到 import 的效果。

%run import_npnet.py
npnet = import_npnet(0)
print(npnet)
print(dir(npnet))


['__doc__', '__loader__', '__name__', '__package__', '__spec__', 'get_ipython', 'npnet']

当然啦，因为这是第0篇笔记，所以以上命令并不会有什么实际效果。

python numpy npnet

[npnet] 序言

机器学习算法的三个模块

代码复用