从零开始的数据科学第二版（五）

原文：zh.annas-archive.org/md5/48ab308fc34189a6d7d26b91b72a6df9

译者：飞龙

协议：CC BY-NC-SA 4.0

第十九章：深度学习

略知一二是危险的；要么深入探索，要么不要触碰那泊里亚之泉。

亚历山大·蒲柏

深度学习最初是指“深度”神经网络的应用（即具有多个隐藏层的网络），尽管实际上这个术语现在包含了各种各样的神经网络架构（包括我们在第十八章中开发的“简单”神经网络）。

在本章中，我们将在之前的工作基础上继续，并查看更广泛的神经网络。为此，我们将介绍一些允许我们以更一般方式思考神经网络的抽象概念。

张量

之前，我们区分了向量（一维数组）和矩阵（二维数组）。当我们开始处理更复杂的神经网络时，我们还需要使用更高维度的数组。

在许多神经网络库中，n维数组被称为张量，这也是我们将它们称为的方式。（有一些学术严谨的数学原因不将n维数组称为张量；如果你是这样的学究，我们注意到你的反对。）

如果我要写一本关于深度学习的整本书，我会实现一个功能齐全的Tensor类，重载 Python 的算术运算符，并能处理各种其他操作。这样的实现将需要一个完整的章节。在这里，我们将简单处理，并说一个Tensor只是一个list。在某种意义上是正确的——我们所有的向量、矩阵和更高维度的模拟 都 是列表。但在另一方面则不正确——大多数 Python 的list不是我们所说的n维数组。

注意

理想情况下，你想做这样的事情：

# A Tensor is either a float, or a List of Tensors
Tensor = Union[float, List[Tensor]]

然而，Python 不允许你定义这样的递归类型。即使它允许，那个定义仍然是错误的，因为它允许像这样的坏“张量”：

[[1.0, 2.0],
 [3.0]]

其行具有不同的大小，这使得它不是一个n维数组。

所以，正如我所说的，我们将简单地作弊：

Tensor = list

而且，我们将编写一个辅助函数来找到张量的形状：

from typing import List

def shape(tensor: Tensor) -> List[int]:
    sizes: List[int] = []
    while isinstance(tensor, list):
        sizes.append(len(tensor))
        tensor = tensor[0]
    return sizes

assert shape([1, 2, 3]) == [3]
assert shape([[1, 2], [3, 4], [5, 6]]) == [3, 2]

因为张量可以具有任意数量的维度，我们通常需要递归地处理它们。在一维情况下我们会做一件事，在更高维情况下我们会递归处理：

def is_1d(tensor: Tensor) -> bool:
    """
 If tensor[0] is a list, it's a higher-order tensor.
 Otherwise, tensor is 1-dimensional (that is, a vector).
 """
    return not isinstance(tensor[0], list)

assert is_1d([1, 2, 3])
assert not is_1d([[1, 2], [3, 4]])

我们可以利用这一点编写一个递归的tensor_sum函数：

def tensor_sum(tensor: Tensor) -> float:
    """Sums up all the values in the tensor"""
    if is_1d(tensor):
        return sum(tensor)  # just a list of floats, use Python sum
    else:
        return sum(tensor_sum(tensor_i)      # Call tensor_sum on each row
                   for tensor_i in tensor)   # and sum up those results.

assert tensor_sum([1, 2, 3]) == 6
assert tensor_sum([[1, 2], [3, 4]]) == 10

如果你不习惯递归思维，你应该思考直到理解为止，因为我们将在本章中始终使用相同的逻辑。但是，我们将创建几个辅助函数，这样我们就不必在每个地方重写这个逻辑。首先，将一个函数逐元素地应用于单个张量：

from typing import Callable

def tensor_apply(f: Callable[[float], float], tensor: Tensor) -> Tensor:
    """Applies f elementwise"""
    if is_1d(tensor):
        return [f(x) for x in tensor]
    else:
        return [tensor_apply(f, tensor_i) for tensor_i in tensor]

assert tensor_apply(lambda x: x + 1, [1, 2, 3]) == [2, 3, 4]
assert tensor_apply(lambda x: 2 * x, [[1, 2], [3, 4]]) == [[2, 4], [6, 8]]

我们可以用它来编写一个函数，创建一个与给定张量形状相同的零张量：

def zeros_like(tensor: Tensor) -> Tensor:
    return tensor_apply(lambda _: 0.0, tensor)

assert zeros_like([1, 2, 3]) == [0, 0, 0]
assert zeros_like([[1, 2], [3, 4]]) == [[0, 0], [0, 0]]

我们还需要将函数应用于两个张量对应的元素（它们最好具有完全相同的形状，尽管我们不会检查这一点）：

def tensor_combine(f: Callable[[float, float], float],
                   t1: Tensor,
                   t2: Tensor) -> Tensor:
    """Applies f to corresponding elements of t1 and t2"""
    if is_1d(t1):
        return [f(x, y) for x, y in zip(t1, t2)]
    else:
        return [tensor_combine(f, t1_i, t2_i)
                for t1_i, t2_i in zip(t1, t2)]

import operator
assert tensor_combine(operator.add, [1, 2, 3], [4, 5, 6]) == [5, 7, 9]
assert tensor_combine(operator.mul, [1, 2, 3], [4, 5, 6]) == [4, 10, 18]

层的抽象

在上一章中，我们构建了一个简单的神经网络，允许我们堆叠两层神经元，每层计算 sigmoid(dot(weights, inputs))。

虽然这可能是对实际神经元功能的理想化表示，但实际上，我们希望允许更多种类的操作。也许我们希望神经元记住它们以前的输入。也许我们想使用不同的激活函数而不是 sigmoid。而且通常情况下，我们希望使用超过两层的网络。（我们的 feed_forward 函数实际上可以处理任意数量的层，但我们的梯度计算不能。）

在本章中，我们将构建用于实现各种神经网络的机制。我们的基本抽象将是 Layer，它知道如何将某个函数应用到其输入上，并知道如何反向传播梯度。

我们在第十八章构建的神经网络可以理解为一个“线性”层，后跟一个“sigmoid”层，然后是另一个线性层和另一个 sigmoid 层。在这些术语中没有区分它们，但这样做将允许我们尝试更通用的结构：

from typing import Iterable, Tuple

class Layer:
    """
 Our neural networks will be composed of Layers, each of which
 knows how to do some computation on its inputs in the "forward"
 direction and propagate gradients in the "backward" direction.
 """
    def forward(self, input):
        """
 Note the lack of types. We're not going to be prescriptive
 about what kinds of inputs layers can take and what kinds
 of outputs they can return.
 """
        raise NotImplementedError

    def backward(self, gradient):
        """
 Similarly, we're not going to be prescriptive about what the
 gradient looks like. It's up to you the user to make sure
 that you're doing things sensibly.
 """
        raise NotImplementedError

    def params(self) -> Iterable[Tensor]:
        """
 Returns the parameters of this layer. The default implementation
 returns nothing, so that if you have a layer with no parameters
 you don't have to implement this.
 """
        return ()

    def grads(self) -> Iterable[Tensor]:
        """
 Returns the gradients, in the same order as params().
 """
        return ()

forward 和 backward 方法需要在我们的具体子类中实现。一旦我们建立了神经网络，我们将希望使用梯度下降来训练它，这意味着我们希望使用其梯度更新网络中的每个参数。因此，我们要求每一层能够告诉我们它的参数和梯度。

一些层（例如，将 sigmoid 应用于每个输入的层）没有需要更新的参数，因此我们提供了一个处理这种情况的默认实现。

让我们看看那层：

from scratch.neural_networks import sigmoid

class Sigmoid(Layer):
    def forward(self, input: Tensor) -> Tensor:
        """
 Apply sigmoid to each element of the input tensor,
 and save the results to use in backpropagation.
 """
        self.sigmoids = tensor_apply(sigmoid, input)
        return self.sigmoids

    def backward(self, gradient: Tensor) -> Tensor:
        return tensor_combine(lambda sig, grad: sig * (1 - sig) * grad,
                              self.sigmoids,
                              gradient)

这里有几件事情需要注意。一是在前向传播期间，我们保存了计算出的 sigmoid 值，以便稍后在反向传播中使用。我们的层通常需要执行这种操作。

其次，您可能想知道 sig * (1 - sig) * grad 是从哪里来的。这只是微积分中的链式法则，对应于我们先前神经网络中的 output * (1 - output) * (output - target) 项。

最后，您可以看到我们如何使用 tensor_apply 和 tensor_combine 函数。我们的大多数层将类似地使用这些函数。

线性层

我们还需要从第十八章复制神经网络所需的“线性”层，该层代表神经元中的 dot(weights, inputs) 部分。

这个层将有参数，我们希望用随机值初始化它们。

结果表明，初始参数值可能极大地影响网络的训练速度（有时甚至影响网络是否能训练）。如果权重过大，它们可能会在激活函数梯度接近零的范围内产生大的输出。具有零梯度部分的网络必然无法通过梯度下降学习任何内容。

因此，我们将实现三种不同的方案来随机生成我们的权重张量。第一种是从[0, 1]上的随机均匀分布中选择每个值，即使用random.random()。第二种（默认）是从标准正态分布中随机选择每个值。第三种是使用Xavier 初始化，其中每个权重都从均值为 0、方差为 2 / (num_inputs + num_outputs)的正态分布中随机抽取。事实证明，这通常对神经网络权重效果很好。我们将使用random_uniform函数和random_normal函数来实现这些：

import random

from scratch.probability import inverse_normal_cdf

def random_uniform(*dims: int) -> Tensor:
    if len(dims) == 1:
        return [random.random() for _ in range(dims[0])]
    else:
        return [random_uniform(*dims[1:]) for _ in range(dims[0])]

def random_normal(*dims: int,
                  mean: float = 0.0,
                  variance: float = 1.0) -> Tensor:
    if len(dims) == 1:
        return [mean + variance * inverse_normal_cdf(random.random())
                for _ in range(dims[0])]
    else:
        return [random_normal(*dims[1:], mean=mean, variance=variance)
                for _ in range(dims[0])]

assert shape(random_uniform(2, 3, 4)) == [2, 3, 4]
assert shape(random_normal(5, 6, mean=10)) == [5, 6]

然后将它们全部包装在一个random_tensor函数中：

def random_tensor(*dims: int, init: str = 'normal') -> Tensor:
    if init == 'normal':
        return random_normal(*dims)
    elif init == 'uniform':
        return random_uniform(*dims)
    elif init == 'xavier':
        variance = len(dims) / sum(dims)
        return random_normal(*dims, variance=variance)
    else:
        raise ValueError(f"unknown init: {init}")

现在我们可以定义我们的线性层了。我们需要用输入的维度来初始化它（这告诉我们每个神经元需要多少个权重），输出的维度（这告诉我们应该有多少个神经元），以及我们想要的初始化方案：

from scratch.linear_algebra import dot

class Linear(Layer):
    def __init__(self,
                 input_dim: int,
                 output_dim: int,
                 init: str = 'xavier') -> None:
        """
 A layer of output_dim neurons, each with input_dim weights
 (and a bias).
 """
        self.input_dim = input_dim
        self.output_dim = output_dim

        # self.w[o] is the weights for the oth neuron
        self.w = random_tensor(output_dim, input_dim, init=init)

        # self.b[o] is the bias term for the oth neuron
        self.b = random_tensor(output_dim, init=init)

注意

如果你想知道初始化方案有多重要，这一章中一些网络如果使用与我使用的不同初始化方法，我根本无法训练它们。

forward方法很容易实现。我们将得到每个神经元的一个输出，将其放入一个向量中。每个神经元的输出只是其权重与输入的dot积，加上偏置：

    def forward(self, input: Tensor) -> Tensor:
        # Save the input to use in the backward pass.
        self.input = input

        # Return the vector of neuron outputs.
        return [dot(input, self.w[o]) + self.b[o]
                for o in range(self.output_dim)]

backward方法更复杂一些，但如果你懂微积分，它并不难：

    def backward(self, gradient: Tensor) -> Tensor:
        # Each b[o] gets added to output[o], which means
        # the gradient of b is the same as the output gradient.
        self.b_grad = gradient

        # Each w[o][i] multiplies input[i] and gets added to output[o].
        # So its gradient is input[i] * gradient[o].
        self.w_grad = [[self.input[i] * gradient[o]
                        for i in range(self.input_dim)]
                       for o in range(self.output_dim)]

        # Each input[i] multiplies every w[o][i] and gets added to every
        # output[o]. So its gradient is the sum of w[o][i] * gradient[o]
        # across all the outputs.
        return [sum(self.w[o][i] * gradient[o] for o in range(self.output_dim))
                for i in range(self.input_dim)]

注意

在一个“真正”的张量库中，这些（以及许多其他）操作将被表示为矩阵或张量乘法，这些库被设计成能够非常快速地执行。我们的库非常慢。

最后，我们确实需要实现params和grads。我们有两个参数和两个相应的梯度：

    def params(self) -> Iterable[Tensor]:
        return [self.w, self.b]

    def grads(self) -> Iterable[Tensor]:
        return [self.w_grad, self.b_grad]

神经网络作为层序列

我们希望将神经网络视为层序列，因此让我们想出一种将多个层组合成一个的方法。得到的神经网络本身就是一个层，它以明显的方式实现了Layer方法：

from typing import List

class Sequential(Layer):
    """
 A layer consisting of a sequence of other layers.
 It's up to you to make sure that the output of each layer
 makes sense as the input to the next layer.
 """
    def __init__(self, layers: List[Layer]) -> None:
        self.layers = layers

    def forward(self, input):
        """Just forward the input through the layers in order."""
        for layer in self.layers:
            input = layer.forward(input)
        return input

    def backward(self, gradient):
        """Just backpropagate the gradient through the layers in reverse."""
        for layer in reversed(self.layers):
            gradient = layer.backward(gradient)
        return gradient

    def params(self) -> Iterable[Tensor]:
        """Just return the params from each layer."""
        return (param for layer in self.layers for param in layer.params())

    def grads(self) -> Iterable[Tensor]:
        """Just return the grads from each layer."""
        return (grad for layer in self.layers for grad in layer.grads())

因此，我们可以将用于 XOR 的神经网络表示为：

xor_net = Sequential([
    Linear(input_dim=2, output_dim=2),
    Sigmoid(),
    Linear(input_dim=2, output_dim=1),
    Sigmoid()
])

但我们仍然需要一些更多的机制来训练它。

损失和优化

之前，我们为我们的模型编写了单独的损失函数和梯度函数。在这里，我们将想要尝试不同的损失函数，所以（像往常一样）我们将引入一个新的Loss抽象，它封装了损失计算和梯度计算：

class Loss:
    def loss(self, predicted: Tensor, actual: Tensor) -> float:
        """How good are our predictions? (Larger numbers are worse.)"""
        raise NotImplementedError

    def gradient(self, predicted: Tensor, actual: Tensor) -> Tensor:
        """How does the loss change as the predictions change?"""
        raise NotImplementedError

我们已经多次使用过损失函数，即平方误差的和，所以实现它应该很容易。唯一的技巧是我们需要使用tensor_combine：

class SSE(Loss):
    """Loss function that computes the sum of the squared errors."""
    def loss(self, predicted: Tensor, actual: Tensor) -> float:
        # Compute the tensor of squared differences
        squared_errors = tensor_combine(
            lambda predicted, actual: (predicted - actual) ** 2,
            predicted,
            actual)

        # And just add them up
        return tensor_sum(squared_errors)

    def gradient(self, predicted: Tensor, actual: Tensor) -> Tensor:
        return tensor_combine(
            lambda predicted, actual: 2 * (predicted - actual),
            predicted,
            actual)

（稍后我们将看一个不同的损失函数。）

最后需要弄清楚的是梯度下降。在整本书中，我们通过一个训练循环手动进行所有的梯度下降，类似于以下操作：

theta = gradient_step(theta, grad, -learning_rate)

这对我们来说不太适用，原因有几个。首先是我们的神经网络将有很多参数，我们需要更新所有这些参数。其次是我们希望能够使用更聪明的梯度下降的变体，而不想每次都重新编写它们。

因此，我们将引入一个（你猜对了）Optimizer 抽象，梯度下降将是一个具体的实例：

class Optimizer:
    """
 An optimizer updates the weights of a layer (in place) using information
 known by either the layer or the optimizer (or by both).
 """
    def step(self, layer: Layer) -> None:
        raise NotImplementedError

之后，使用 tensor_combine 再次实现梯度下降就很容易了：

class GradientDescent(Optimizer):
    def __init__(self, learning_rate: float = 0.1) -> None:
        self.lr = learning_rate

    def step(self, layer: Layer) -> None:
        for param, grad in zip(layer.params(), layer.grads()):
            # Update param using a gradient step
            param[:] = tensor_combine(
                lambda param, grad: param - grad * self.lr,
                param,
                grad)

可能令人惊讶的唯一一件事是“切片赋值”，这反映了重新分配列表不会改变其原始值的事实。也就是说，如果你只是做了 param = tensor_combine(. . .)，你会重新定义局部变量 param，但你不会影响存储在层中的原始参数张量。然而，如果你分配给切片 [:]，它实际上会改变列表内的值。

这里有一个简单的示例来演示：

tensor = [[1, 2], [3, 4]]

for row in tensor:
    row = [0, 0]
assert tensor == [[1, 2], [3, 4]], "assignment doesn't update a list"

for row in tensor:
    row[:] = [0, 0]
assert tensor == [[0, 0], [0, 0]], "but slice assignment does"

如果你在 Python 方面有些经验不足，这种行为可能会让你感到惊讶，所以要沉思一下，并尝试自己的例子，直到它变得清晰为止。

为了展示这种抽象的价值，让我们实现另一个使用 动量 的优化器。其思想是我们不希望对每一个新的梯度过于反应，因此我们保持已看到的梯度的运行平均，每次新的梯度更新它，并朝平均梯度的方向迈出一步：

class Momentum(Optimizer):
    def __init__(self,
                 learning_rate: float,
                 momentum: float = 0.9) -> None:
        self.lr = learning_rate
        self.mo = momentum
        self.updates: List[Tensor] = []  # running average

    def step(self, layer: Layer) -> None:
        # If we have no previous updates, start with all zeros
        if not self.updates:
            self.updates = [zeros_like(grad) for grad in layer.grads()]

        for update, param, grad in zip(self.updates,
                                       layer.params(),
                                       layer.grads()):
            # Apply momentum
            update[:] = tensor_combine(
                lambda u, g: self.mo * u + (1 - self.mo) * g,
                update,
                grad)

            # Then take a gradient step
            param[:] = tensor_combine(
                lambda p, u: p - self.lr * u,
                param,
                update)

因为我们使用了 Optimizer 抽象，我们可以轻松地在不同的优化器之间切换。

示例：重新思考 XOR

让我们看看使用我们的新框架来训练一个能计算 XOR 的网络有多容易。我们首先重新创建训练数据：

# training data
xs = [[0., 0], [0., 1], [1., 0], [1., 1]]
ys = [[0.], [1.], [1.], [0.]]

然后我们定义网络，尽管现在我们可以省略最后的 sigmoid 层：

random.seed(0)

net = Sequential([
    Linear(input_dim=2, output_dim=2),
    Sigmoid(),
    Linear(input_dim=2, output_dim=1)
])

现在我们可以编写一个简单的训练循环，除了现在我们可以使用 Optimizer 和 Loss 的抽象。这使得我们可以轻松尝试不同的优化方法：

import tqdm

optimizer = GradientDescent(learning_rate=0.1)
loss = SSE()

with tqdm.trange(3000) as t:
    for epoch in t:
        epoch_loss = 0.0

        for x, y in zip(xs, ys):
            predicted = net.forward(x)
            epoch_loss += loss.loss(predicted, y)
            gradient = loss.gradient(predicted, y)
            net.backward(gradient)

            optimizer.step(net)

        t.set_description(f"xor loss {epoch_loss:.3f}")

这应该会快速训练，并且你应该会看到损失下降。现在我们可以检查权重了：

for param in net.params():
    print(param)

对于我的网络，我发现大致上：

hidden1 = -2.6 * x1 + -2.7 * x2 + 0.2  # NOR
hidden2 =  2.1 * x1 +  2.1 * x2 - 3.4  # AND
output =  -3.1 * h1 + -2.6 * h2 + 1.8  # NOR

因此，如果 hidden1 激活，那么没有一个输入是 1。如果 hidden2 激活，那么两个输入都是 1。如果 output 激活，那么既不是 hidden 输出是 1，也不是两个输入都是 1。确实，这正是 XOR 的逻辑。

注意，这个网络学习了不同于我们在 第十八章 中训练的网络的特征，但它仍然能够执行相同的操作。

其他激活函数

sigmoid 函数因为几个原因已经不再流行。其中一个原因是 sigmoid(0) 等于 1/2，这意味着输入总和为 0 的神经元具有正的输出。另一个原因是对于非常大和非常小的输入，它的梯度非常接近 0，这意味着它的梯度可能会“饱和”，它的权重可能会被困住。

一个流行的替代方案是tanh（“双曲正切”），这是一个不同的 S 型函数，范围从–1 到 1，并且如果其输入为 0，则输出为 0。tanh(x)的导数就是1 - tanh(x) ** 2，这样写起来很容易：

import math

def tanh(x: float) -> float:
    # If x is very large or very small, tanh is (essentially) 1 or -1.
    # We check for this because, e.g., math.exp(1000) raises an error.
    if x < -100:  return -1
    elif x > 100: return 1

    em2x = math.exp(-2 * x)
    return (1 - em2x) / (1 + em2x)

class Tanh(Layer):
    def forward(self, input: Tensor) -> Tensor:
        # Save tanh output to use in backward pass.
        self.tanh = tensor_apply(tanh, input)
        return self.tanh

    def backward(self, gradient: Tensor) -> Tensor:
        return tensor_combine(
            lambda tanh, grad: (1 - tanh ** 2) * grad,
            self.tanh,
            gradient)

在更大的网络中，另一个流行的替代方案是Relu，对于负输入为 0，对于正输入为恒等：

class Relu(Layer):
    def forward(self, input: Tensor) -> Tensor:
        self.input = input
        return tensor_apply(lambda x: max(x, 0), input)

    def backward(self, gradient: Tensor) -> Tensor:
        return tensor_combine(lambda x, grad: grad if x > 0 else 0,
                              self.input,
                              gradient)

还有许多其他的函数。我鼓励你在你的网络中尝试它们。

示例：FizzBuzz 再探讨

现在，我们可以使用我们的“深度学习”框架来重新生成我们在“示例：Fizz Buzz”中的解决方案。让我们设置数据：

from scratch.neural_networks import binary_encode, fizz_buzz_encode, argmax

xs = [binary_encode(n) for n in range(101, 1024)]
ys = [fizz_buzz_encode(n) for n in range(101, 1024)]

并创建网络：

NUM_HIDDEN = 25

random.seed(0)

net = Sequential([
    Linear(input_dim=10, output_dim=NUM_HIDDEN, init='uniform'),
    Tanh(),
    Linear(input_dim=NUM_HIDDEN, output_dim=4, init='uniform'),
    Sigmoid()
])

当我们训练时，让我们也跟踪一下训练集上的准确率：

def fizzbuzz_accuracy(low: int, hi: int, net: Layer) -> float:
    num_correct = 0
    for n in range(low, hi):
        x = binary_encode(n)
        predicted = argmax(net.forward(x))
        actual = argmax(fizz_buzz_encode(n))
        if predicted == actual:
            num_correct += 1

    return num_correct / (hi - low)

optimizer = Momentum(learning_rate=0.1, momentum=0.9)
loss = SSE()

with tqdm.trange(1000) as t:
    for epoch in t:
        epoch_loss = 0.0

        for x, y in zip(xs, ys):
            predicted = net.forward(x)
            epoch_loss += loss.loss(predicted, y)
            gradient = loss.gradient(predicted, y)
            net.backward(gradient)

            optimizer.step(net)

        accuracy = fizzbuzz_accuracy(101, 1024, net)
        t.set_description(f"fb loss: {epoch_loss:.2f} acc: {accuracy:.2f}")

# Now check results on the test set
print("test results", fizzbuzz_accuracy(1, 101, net))

经过 1000 次训练迭代后，模型在测试集上达到了 90%的准确率；如果你继续训练，它应该能表现得更好。（我认为只用 25 个隐藏单元不可能训练到 100%的准确率，但如果增加到 50 个隐藏单元，这是可能的。）

Softmax 函数和交叉熵

我们在前面章节使用的神经网络以Sigmoid层结束，这意味着其输出是介于 0 和 1 之间的向量。特别是，它可以输出一个完全是 0 的向量，或者一个完全是 1 的向量。然而，在分类问题中，我们希望输出一个 1 代表正确类别，0 代表所有不正确的类别。通常我们的预测不会那么完美，但我们至少希望预测一个实际的类别概率分布。

例如，如果我们有两个类别，而我们的模型输出[0, 0]，很难理解这意味着什么。它不认为输出属于任何类别？

但是如果我们的模型输出[0.4, 0.6]，我们可以将其解释为预测我们的输入属于第一类的概率为 0.4，属于第二类的概率为 0.6。

为了实现这一点，我们通常放弃最后的Sigmoid层，而是使用softmax函数，将实数向量转换为概率向量。我们对向量中的每个数计算exp(x)，得到一组正数。然后，我们将这些正数除以它们的和，得到一组加起来为 1 的正数——即概率向量。

如果我们试图计算，比如exp(1000)，我们会得到一个 Python 错误，所以在计算exp之前，我们要减去最大值。这样做结果是相同的概率；在 Python 中这样计算更安全：

def softmax(tensor: Tensor) -> Tensor:
    """Softmax along the last dimension"""
    if is_1d(tensor):
        # Subtract largest value for numerical stability.
        largest = max(tensor)
        exps = [math.exp(x - largest) for x in tensor]

        sum_of_exps = sum(exps)                 # This is the total "weight."
        return [exp_i / sum_of_exps             # Probability is the fraction
                for exp_i in exps]              # of the total weight.
    else:
        return [softmax(tensor_i) for tensor_i in tensor]

一旦我们的网络产生了概率，我们通常使用另一种称为交叉熵（有时称为“负对数似然”）的损失函数。

你可能记得，在“最大似然估计”中，我们通过引用最小二乘在线性回归中的使用来证明（在某些假设下）最小二乘系数最大化了观察数据的似然。

在这里我们可以做类似的事情：如果我们的网络输出是概率，交叉熵损失表示观察数据的负对数似然，这意味着最小化该损失等同于最大化（因而最大化）训练数据的似然。

通常情况下，我们不会将softmax函数作为神经网络本身的一部分。这是因为事实证明，如果softmax是你的损失函数的一部分，但不是网络本身的一部分，那么损失关于网络输出的梯度计算非常容易。

class SoftmaxCrossEntropy(Loss):
    """
 This is the negative-log-likelihood of the observed values, given the
 neural net model. So if we choose weights to minimize it, our model will
 be maximizing the likelihood of the observed data.
 """
    def loss(self, predicted: Tensor, actual: Tensor) -> float:
        # Apply softmax to get probabilities
        probabilities = softmax(predicted)

        # This will be log p_i for the actual class i and 0 for the other
        # classes. We add a tiny amount to p to avoid taking log(0).
        likelihoods = tensor_combine(lambda p, act: math.log(p + 1e-30) * act,
                                     probabilities,
                                     actual)

        # And then we just sum up the negatives.
        return -tensor_sum(likelihoods)

    def gradient(self, predicted: Tensor, actual: Tensor) -> Tensor:
        probabilities = softmax(predicted)

        # Isn't this a pleasant equation?
        return tensor_combine(lambda p, actual: p - actual,
                              probabilities,
                              actual)

如果现在我使用SoftmaxCrossEntropy损失训练相同的 Fizz Buzz 网络，我发现它通常训练速度快得多（即在更少的周期内）。这可能是因为找到使softmax到给定分布的权重比找到使sigmoid到给定分布的权重更容易。

换句话说，如果我需要预测类别 0（一个向量，第一个位置为 1，其余位置为 0），在linear + sigmoid的情况下，我需要第一个输出为一个较大的正数，其余的输出为较大的负数。然而，在softmax的情况下，我只需要第一个输出比其余的输出大。显然，第二种情况发生的方式更多，这表明找到使其成为可能的权重应该更容易：

random.seed(0)

net = Sequential([
    Linear(input_dim=10, output_dim=NUM_HIDDEN, init='uniform'),
    Tanh(),
    Linear(input_dim=NUM_HIDDEN, output_dim=4, init='uniform')
    # No final sigmoid layer now
])

optimizer = Momentum(learning_rate=0.1, momentum=0.9)
loss = SoftmaxCrossEntropy()

with tqdm.trange(100) as t:
    for epoch in t:
        epoch_loss = 0.0

        for x, y in zip(xs, ys):
            predicted = net.forward(x)
            epoch_loss += loss.loss(predicted, y)
            gradient = loss.gradient(predicted, y)
            net.backward(gradient)

            optimizer.step(net)

        accuracy = fizzbuzz_accuracy(101, 1024, net)
        t.set_description(f"fb loss: {epoch_loss:.3f} acc: {accuracy:.2f}")

# Again check results on the test set
print("test results", fizzbuzz_accuracy(1, 101, net))

Dropout

像大多数机器学习模型一样，神经网络容易过拟合其训练数据。我们之前已经看到了一些缓解这种情况的方法；例如，在“正则化”中，我们对大的权重进行了惩罚，这有助于防止过拟合。

常见的神经网络正则化方法之一是使用dropout。在训练时，我们随机关闭每个神经元（即将其输出替换为 0），关闭的概率固定。这意味着网络不能学会依赖任何单个神经元，这似乎有助于防止过拟合。

在评估时，我们不希望关闭任何神经元，因此Dropout层需要知道它是在训练还是不训练。此外，在训练时，Dropout层仅传递其输入的一些随机部分。为了在评估期间使其输出可比较，我们将使用相同的比例缩减输出（均匀地）：

class Dropout(Layer):
    def __init__(self, p: float) -> None:
        self.p = p
        self.train = True

    def forward(self, input: Tensor) -> Tensor:
        if self.train:
            # Create a mask of 0s and 1s shaped like the input
            # using the specified probability.
            self.mask = tensor_apply(
                lambda _: 0 if random.random() < self.p else 1,
                input)
            # Multiply by the mask to dropout inputs.
            return tensor_combine(operator.mul, input, self.mask)
        else:
            # During evaluation just scale down the outputs uniformly.
            return tensor_apply(lambda x: x * (1 - self.p), input)

    def backward(self, gradient: Tensor) -> Tensor:
        if self.train:
            # Only propagate the gradients where mask == 1.
            return tensor_combine(operator.mul, gradient, self.mask)
        else:
            raise RuntimeError("don't call backward when not in train mode")

我们将使用这个来帮助防止我们的深度学习模型过拟合。

示例：MNIST

MNIST 是一个手写数字数据集，每个人都用它来学习深度学习。

这种数据以一种有些棘手的二进制格式提供，因此我们将安装mnist库来处理它。（是的，这部分从技术上讲并非“从头开始”。）

python -m pip install mnist

然后我们可以加载数据：

import mnist

# This will download the data; change this to where you want it.
# (Yes, it's a 0-argument function, that's what the library expects.)
# (Yes, I'm assigning a lambda to a variable, like I said never to do.)
mnist.temporary_dir = lambda: '/tmp'

# Each of these functions first downloads the data and returns a numpy array.
# We call .tolist() because our "tensors" are just lists.
train_images = mnist.train_images().tolist()
train_labels = mnist.train_labels().tolist()

assert shape(train_images) == [60000, 28, 28]
assert shape(train_labels) == [60000]

让我们绘制前 100 张训练图像，看看它们的样子（图 19-1）：

import matplotlib.pyplot as plt

fig, ax = plt.subplots(10, 10)

for i in range(10):
    for j in range(10):
        # Plot each image in black and white and hide the axes.
        ax[i][j].imshow(train_images[10 * i + j], cmap='Greys')
        ax[i][j].xaxis.set_visible(False)
        ax[i][j].yaxis.set_visible(False)

plt.show()

图 19-1. MNIST 图像

您可以看到它们确实看起来像手写数字。

注意

我第一次尝试显示图像时，结果是黄色数字在黑色背景上。我既不聪明也不够细心，不知道我需要添加 cmap=*Greys* 才能获得黑白图像；我在 Stack Overflow 上搜索找到了解决方案。作为一名数据科学家，你将变得非常擅长这种工作流程。

我们还需要加载测试图像：

test_images = mnist.test_images().tolist()
test_labels = mnist.test_labels().tolist()

assert shape(test_images) == [10000, 28, 28]
assert shape(test_labels) == [10000]

每个图像是 28 × 28 像素，但我们的线性层只能处理一维输入，因此我们将它们展平（同时除以 256 以使它们在 0 到 1 之间）。此外，如果我们的输入平均值为 0，我们的神经网络将更好地训练，因此我们将减去平均值：

# Compute the average pixel value
avg = tensor_sum(train_images) / 60000 / 28 / 28

# Recenter, rescale, and flatten
train_images = [[(pixel - avg) / 256 for row in image for pixel in row]
                for image in train_images]
test_images = [[(pixel - avg) / 256 for row in image for pixel in row]
               for image in test_images]

assert shape(train_images) == [60000, 784], "images should be flattened"
assert shape(test_images) == [10000, 784], "images should be flattened"

# After centering, average pixel should be very close to 0
assert -0.0001 < tensor_sum(train_images) < 0.0001

我们还希望对目标进行独热编码，因为我们有 10 个输出。首先让我们编写一个 one_hot_encode 函数：

def one_hot_encode(i: int, num_labels: int = 10) -> List[float]:
    return [1.0 if j == i else 0.0 for j in range(num_labels)]

assert one_hot_encode(3) == [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
assert one_hot_encode(2, num_labels=5) == [0, 0, 1, 0, 0]

然后将其应用到我们的数据中：

train_labels = [one_hot_encode(label) for label in train_labels]
test_labels = [one_hot_encode(label) for label in test_labels]

assert shape(train_labels) == [60000, 10]
assert shape(test_labels) == [10000, 10]

我们抽象化的一个优势是，我们可以用多种模型来使用相同的训练/评估循环。因此，让我们首先编写它。我们会传入我们的模型、数据、损失函数和（如果我们在训练）一个优化器。

它将遍历我们的数据，跟踪性能，并且（如果我们传入了优化器）更新我们的参数：

import tqdm

def loop(model: Layer,
         images: List[Tensor],
         labels: List[Tensor],
         loss: Loss,
         optimizer: Optimizer = None) -> None:
    correct = 0         # Track number of correct predictions.
    total_loss = 0.0    # Track total loss.

    with tqdm.trange(len(images)) as t:
        for i in t:
            predicted = model.forward(images[i])             # Predict.
            if argmax(predicted) == argmax(labels[i]):       # Check for
                correct += 1                                 # correctness.
            total_loss += loss.loss(predicted, labels[i])    # Compute loss.

            # If we're training, backpropagate gradient and update weights.
            if optimizer is not None:
                gradient = loss.gradient(predicted, labels[i])
                model.backward(gradient)
                optimizer.step(model)

            # And update our metrics in the progress bar.
            avg_loss = total_loss / (i + 1)
            acc = correct / (i + 1)
            t.set_description(f"mnist loss: {avg_loss:.3f} acc: {acc:.3f}")

作为基准，我们可以使用我们的深度学习库训练（多类别）逻辑回归模型，这只是一个单一线性层，后跟一个 softmax。该模型（实质上）只是寻找 10 个线性函数，如果输入代表，比如说，一个 5，那么第 5 个线性函数将产生最大的输出。

通过我们的 60,000 个训练样本的一次遍历应该足以学习模型：

random.seed(0)

# Logistic regression is just a linear layer followed by softmax
model = Linear(784, 10)
loss = SoftmaxCrossEntropy()

# This optimizer seems to work
optimizer = Momentum(learning_rate=0.01, momentum=0.99)

# Train on the training data
loop(model, train_images, train_labels, loss, optimizer)

# Test on the test data (no optimizer means just evaluate)
loop(model, test_images, test_labels, loss)

这个准确率约为 89%。让我们看看是否可以通过深度神经网络做得更好。我们将使用两个隐藏层，第一个有 30 个神经元，第二个有 10 个神经元。我们将使用我们的 Tanh 激活函数：

random.seed(0)

# Name them so we can turn train on and off
dropout1 = Dropout(0.1)
dropout2 = Dropout(0.1)

model = Sequential([
    Linear(784, 30),  # Hidden layer 1: size 30
    dropout1,
    Tanh(),
    Linear(30, 10),   # Hidden layer 2: size 10
    dropout2,
    Tanh(),
    Linear(10, 10)    # Output layer: size 10
])

我们可以只使用相同的训练循环！

optimizer = Momentum(learning_rate=0.01, momentum=0.99)
loss = SoftmaxCrossEntropy()

# Enable dropout and train (takes > 20 minutes on my laptop!)
dropout1.train = dropout2.train = True
loop(model, train_images, train_labels, loss, optimizer)

# Disable dropout and evaluate
dropout1.train = dropout2.train = False
loop(model, test_images, test_labels, loss)

我们的深度模型在测试集上的准确率超过了 92%，这比简单的逻辑回归模型有了显著提升。

MNIST 网站 描述了多种优于这些模型的模型。其中许多模型可以使用我们迄今为止开发的机制实现，但在我们的列表作为张量的框架中训练时间将非常长。一些最佳模型涉及到 卷积 层，这很重要，但不幸的是对于数据科学入门书籍来说有些超出范围。

保存和加载模型

这些模型需要很长时间来训练，因此如果我们能保存它们以免每次都重新训练就太好了。幸运的是，我们可以使用 json 模块轻松地将模型权重序列化到文件中。

对于保存，我们可以使用 Layer.params 收集权重，将它们放入列表中，并使用 json.dump 将该列表保存到文件中：

import json

def save_weights(model: Layer, filename: str) -> None:
    weights = list(model.params())
    with open(filename, 'w') as f:
        json.dump(weights, f)

将权重加载回来只需要多做一点工作。我们只需使用json.load从文件中获取权重列表，然后使用切片赋值将权重设置到我们的模型中。

(具体来说，这意味着我们必须自己实例化模型，然后加载权重。另一种方法是保存模型架构的某种表示，并使用它来实例化模型。这并不是一个糟糕的想法，但这将需要大量的代码和对所有我们的Layer进行更改，所以我们将坚持使用更简单的方法。)

在加载权重之前，我们希望检查它们与我们加载到其中的模型参数具有相同的形状。（这是为了防止例如将保存的深度网络的权重加载到浅网络中，或类似的问题。）

def load_weights(model: Layer, filename: str) -> None:
    with open(filename) as f:
        weights = json.load(f)

    # Check for consistency
    assert all(shape(param) == shape(weight)
               for param, weight in zip(model.params(), weights))

    # Then load using slice assignment
    for param, weight in zip(model.params(), weights):
        param[:] = weight

注意

JSON 将数据存储为文本，这使其成为极其低效的表示形式。在实际应用中，您可能会使用pickle序列化库，它将事物序列化为更有效的二进制格式。在这里，我决定保持简单和人类可读性。

您可以从该书的 GitHub 存储库下载我们训练的各种网络的权重。

深入探索

现在深度学习非常火热，在本章中我们只是简单介绍了一下。关于几乎任何您想了解的深度学习方面，都有许多好书和博客文章（以及很多很多糟糕的博客文章）。

• 深度学习 这本经典教材，作者是 Ian Goodfellow、Yoshua Bengio 和 Aaron Courville（MIT Press），可以在线免费获取。这本书非常好，但涉及到相当多的数学。

• Francois Chollet 的Python 深度学习（Manning）是了解 Keras 库的绝佳入门书籍，我们的深度学习库就是基于这种模式设计的。

• 我自己大多使用PyTorch进行深度学习。它的网站有大量的文档和教程。

第二十章：聚类

当我们有这样的聚类时

使我们变得高尚而不是疯狂

罗伯特·赫里克

本书中的大多数算法都属于监督学习算法，即它们从一组带标签的数据开始，并将其用作对新的未标记数据进行预测的基础。然而，聚类是无监督学习的一个例子，我们在其中使用完全未标记的数据（或者我们的数据具有标签但我们忽略它们）。

思路

每当您查看某些数据源时，很可能数据会以某种方式形成聚类。显示百万富翁住在哪里的数据集可能在比佛利山庄和曼哈顿等地形成聚类。显示人们每周工作多少小时的数据集可能在周围形成一个大约为 40 的聚类（如果它来自于一项法律规定每周至少工作 20 小时的州，那么它可能还有另一个大约在 19 左右的聚类）。登记选民的人口统计数据集可能形成各种各样的聚类（例如，“足球妈妈”，“无聊的退休者”，“失业的千禧一代”），这些聚类被民意调查员和政治顾问认为是相关的。

与我们已经研究过的一些问题不同，通常没有“正确”的聚类。另一种聚类方案可能会将一些“失业的千禧一代”与“研究生”分组，而将其他一些与“父母的地下室居民”分组。这两种方案都不一定更正确——而是每一种都可能更优，关于自身的“聚类有多好？”度量标准而言。

此外，这些簇不会自行标记。您需要通过查看每个簇下面的数据来完成。

模型

对我们来说，每个input将是d维空间中的一个向量，通常我们将其表示为数字列表。我们的目标将是识别相似输入的簇，并（有时）为每个簇找到一个代表值。

例如，每个输入可以是表示博客文章标题的数值向量，此时目标可能是找到相似帖子的簇，也许是为了理解我们的用户正在博客的内容。或者想象一下，我们有一张包含数千个(红色，绿色，蓝色)颜色的图片，并且我们需要丝网印刷它的 10 种颜色版本。聚类可以帮助我们选择最小化总“颜色误差”的 10 种颜色。

最简单的聚类方法之一是k-means，在其中聚类数k是预先选择的，然后目标是以最小化每个点到其分配的簇的平均值的距离的总平方和的方式将输入分区到集合S 1 , ... , S k。

有很多方法可以将n个点分配到k个簇中，这意味着找到最佳聚类是一个非常困难的问题。我们将采用一个通常能找到良好聚类的迭代算法：

1. 以一组k-means 开始，这些点位于d维空间中。

2. 将每个点分配给最接近它的均值。

3. 如果没有点的分配发生变化，请停止并保留这些簇。

4. 如果某些点的分配发生了变化，请重新计算均值并返回到步骤 2。

使用来自第四章的vector_mean函数，创建一个执行此操作的类非常简单。

首先，我们将创建一个辅助函数，用于衡量两个向量在多少坐标上不同。我们将使用这个函数来跟踪我们的训练进度：

from scratch.linear_algebra import Vector

def num_differences(v1: Vector, v2: Vector) -> int:
    assert len(v1) == len(v2)
    return len([x1 for x1, x2 in zip(v1, v2) if x1 != x2])

assert num_differences([1, 2, 3], [2, 1, 3]) == 2
assert num_differences([1, 2], [1, 2]) == 0

我们还需要一个函数，它根据一些向量及其对簇的分配计算簇的均值。有可能某个簇没有分配到任何点。我们不能对空集合取平均值，因此在这种情况下，我们将随机选择一个点作为该簇的“均值”：

from typing import List
from scratch.linear_algebra import vector_mean

def cluster_means(k: int,
                  inputs: List[Vector],
                  assignments: List[int]) -> List[Vector]:
    # clusters[i] contains the inputs whose assignment is i
    clusters = [[] for i in range(k)]
    for input, assignment in zip(inputs, assignments):
        clusters[assignment].append(input)

    # if a cluster is empty, just use a random point
    return [vector_mean(cluster) if cluster else random.choice(inputs)
            for cluster in clusters]

现在我们可以编写我们的聚类器代码了。像往常一样，我们将使用tqdm来跟踪我们的进度，但在这里，我们不知道需要多少次迭代，因此我们使用itertools.count，它创建一个无限迭代器，当完成时我们会从中return出来：

import itertools
import random
import tqdm
from scratch.linear_algebra import squared_distance

class KMeans:
    def __init__(self, k: int) -> None:
        self.k = k                      # number of clusters
        self.means = None

    def classify(self, input: Vector) -> int:
        """return the index of the cluster closest to the input"""
        return min(range(self.k),
                   key=lambda i: squared_distance(input, self.means[i]))

    def train(self, inputs: List[Vector]) -> None:
        # Start with random assignments
        assignments = [random.randrange(self.k) for _ in inputs]

        with tqdm.tqdm(itertools.count()) as t:
            for _ in t:
                # Compute means and find new assignments
                self.means = cluster_means(self.k, inputs, assignments)
                new_assignments = [self.classify(input) for input in inputs]

                # Check how many assignments changed and if we're done
                num_changed = num_differences(assignments, new_assignments)
                if num_changed == 0:
                    return

                # Otherwise keep the new assignments, and compute new means
                assignments = new_assignments
                self.means = cluster_means(self.k, inputs, assignments)
                t.set_description(f"changed: {num_changed} / {len(inputs)}")

让我们看看这是如何工作的。

例子：见面会

为了庆祝 DataSciencester 的增长，您的用户奖励副总裁希望为您的本地用户组织几次面对面见面会，提供啤酒、披萨和 DataSciencester T 恤。您知道所有本地用户的位置（见图 20-1），她希望您选择方便每个人参加的见面地点。

图 20-1. 您的本地用户位置

根据你的视角不同，你可能看到两个或三个簇。（从视觉上很容易，因为数据只有二维。如果是更多维度，只凭眼睛判断会更难。）

想象一下，她预算足够举办三次见面会。你去电脑上尝试这个：

random.seed(12)                   # so you get the same results as me
clusterer = KMeans(k=3)
clusterer.train(inputs)
means = sorted(clusterer.means)   # sort for the unit test

assert len(means) == 3

# Check that the means are close to what we expect
assert squared_distance(means[0], [-44, 5]) < 1
assert squared_distance(means[1], [-16, -10]) < 1
assert squared_distance(means[2], [18, 20]) < 1

你会找到三个以[-44, 5]，[-16, -10]和[18, 20]为中心的簇，并且寻找靠近这些位置的见面场所（见图 20-2）。

图 20-2. 用户位置分为三个簇

您将结果展示给副总裁，她通知您现在只有预算足够举办两次见面会。

“没问题”，你说：

random.seed(0)
clusterer = KMeans(k=2)
clusterer.train(inputs)
means = sorted(clusterer.means)

assert len(means) == 2
assert squared_distance(means[0], [-26, -5]) < 1
assert squared_distance(means[1], [18, 20]) < 1

如图 20-3 所示，一个见面会仍应接近[18, 20]，但另一个现在应该接近[-26, -5]。

图 20-3. 用户位置分为两个簇

选择k

在上一个示例中，k的选择受到我们控制之外的因素的驱动。一般情况下，这不会发生。有各种选择k的方法。其中一个相对容易理解的方法涉及绘制作为k的平方误差和（每个点与其簇的平均值之间的平方误差）的函数，并查看图表“弯曲”的位置：

from matplotlib import pyplot as plt

def squared_clustering_errors(inputs: List[Vector], k: int) -> float:
    """finds the total squared error from k-means clustering the inputs"""
    clusterer = KMeans(k)
    clusterer.train(inputs)
    means = clusterer.means
    assignments = [clusterer.classify(input) for input in inputs]

    return sum(squared_distance(input, means[cluster])
               for input, cluster in zip(inputs, assignments))

这可以应用到我们之前的例子中：

# now plot from 1 up to len(inputs) clusters

ks = range(1, len(inputs) + 1)
errors = [squared_clustering_errors(inputs, k) for k in ks]

plt.plot(ks, errors)
plt.xticks(ks)
plt.xlabel("k")
plt.ylabel("total squared error")
plt.title("Total Error vs. # of Clusters")
plt.show()

查看图 20-4，这种方法与我们最初的直觉观察相符，认为三是“正确”的聚类数目。

图 20-4. 选择一个k

示例：聚类颜色

Swag 的副总裁设计了引人注目的 DataSciencester 贴纸，他希望你在聚会上分发。不幸的是，你的贴纸打印机每张最多只能打印五种颜色。由于艺术副总裁正在休假，Swag 的副总裁问你是否有办法修改他的设计，使其只包含五种颜色。

计算机图像可以表示为像素的二维数组，其中每个像素本身是一个三维向量(red, green, blue)，表示其颜色。

创建图片的五种颜色版本，因此，涉及：

1. 选择五种颜色。

2. 为每个像素分配其中的一种颜色。

原来这是k-means 聚类的一个很好的任务，它可以将像素在红-绿-蓝色彩空间中分成五个聚类。然后，如果我们将每个聚类中的像素重新着色为平均颜色，我们就完成了。

首先，我们需要一种方法将图像加载到 Python 中。我们可以通过 matplotlib 来实现这一点，前提是我们首先安装 pillow 库：

python -m pip install pillow

然后我们只需使用matplotlib.image.imread：

image_path = r"girl_with_book.jpg"    # wherever your image is
import matplotlib.image as mpimg
img = mpimg.imread(image_path) / 256  # rescale to between 0 and 1

在幕后，img是一个 NumPy 数组，但是对于我们的目的，我们可以将其视为列表的列表的列表。

img[i][j]是第i行第j列的像素，每个像素是一个列表[red, green, blue]，数字介于 0 和 1 之间，表示该像素的颜色：

top_row = img[0]
top_left_pixel = top_row[0]
red, green, blue = top_left_pixel

特别是，我们可以获得所有像素的扁平化列表，如下所示：

# .tolist() converts a NumPy array to a Python list
pixels = [pixel.tolist() for row in img for pixel in row]

然后将它们提供给我们的聚类器：

clusterer = KMeans(5)
clusterer.train(pixels)   # this might take a while

完成后，我们只需构造一个新的具有相同格式的图像：

def recolor(pixel: Vector) -> Vector:
    cluster = clusterer.classify(pixel)        # index of the closest cluster
    return clusterer.means[cluster]            # mean of the closest cluster

new_img = [[recolor(pixel) for pixel in row]   # recolor this row of pixels
           for row in img]                     # for each row in the image

并显示它，使用plt.imshow：

plt.imshow(new_img)
plt.axis('off')
plt.show()

在黑白书籍中展示彩色结果很困难，但图 20-5 显示了将全彩色图片转换为灰度版本以及使用此过程减少至五种颜色的输出。

图 20-5. 原始图片及其 5-means 去色化结果

自底向上的分层聚类

聚类的另一种方法是从底部向上“增长”聚类。我们可以这样做：

1. 使每个输入成为自己的一个簇。

2. 只要还有多个剩余的聚类，就找到最接近的两个聚类并将它们合并。

最后，我们将拥有一个包含所有输入的巨大聚类。如果我们跟踪合并顺序，我们可以通过取消合并来重新创建任意数量的聚类。例如，如果我们想要三个聚类，我们可以撤销最后两个合并。

我们将使用聚类的一个非常简单的表示。我们的值将存储在leaf簇中，并将其表示为NamedTuple：

from typing import NamedTuple, Union

class Leaf(NamedTuple):
    value: Vector

leaf1 = Leaf([10,  20])
leaf2 = Leaf([30, -15])

我们将使用这些来增长merged聚类，我们也将其表示为NamedTuple：

class Merged(NamedTuple):
    children: tuple
    order: int

merged = Merged((leaf1, leaf2), order=1)

Cluster = Union[Leaf, Merged]

注意

这是另一种情况，Python 的类型注解让我们感到失望。你想用Tuple[Cluster, Cluster]作为Merged.children的类型提示，但mypy不允许这样的递归类型。

我们稍后会讨论合并顺序，但首先让我们创建一个递归返回所有值的帮助函数，这些值包含在（可能已合并的）簇中：

def get_values(cluster: Cluster) -> List[Vector]:
    if isinstance(cluster, Leaf):
        return [cluster.value]
    else:
        return [value
                for child in cluster.children
                for value in get_values(child)]

assert get_values(merged) == [[10, 20], [30, -15]]

为了合并最接近的簇，我们需要一些关于簇之间距离的概念。我们将使用两个簇中元素之间的最小距离，这将合并最接近接触的两个簇（但有时会产生不太紧密的链式簇）。如果我们想要紧凑的球状簇，我们可能会改用最大距离，因为它会合并适合最小球中的两个簇。这两种选择都很常见，同样常见的是平均距离：

from typing import Callable
from scratch.linear_algebra import distance

def cluster_distance(cluster1: Cluster,
                     cluster2: Cluster,
                     distance_agg: Callable = min) -> float:
    """
 compute all the pairwise distances between cluster1 and cluster2
 and apply the aggregation function _distance_agg_ to the resulting list
 """
    return distance_agg([distance(v1, v2)
                         for v1 in get_values(cluster1)
                         for v2 in get_values(cluster2)])

我们将使用合并顺序插槽来跟踪我们执行合并的顺序。较小的数字将表示较晚的合并。这意味着当我们想要拆分簇时，我们会从最低的合并顺序到最高的顺序进行。由于Leaf簇从未合并过，我们将给它们分配无穷大，即最大可能的值。由于它们没有.order属性，因此我们将创建一个辅助函数：

def get_merge_order(cluster: Cluster) -> float:
    if isinstance(cluster, Leaf):
        return float('inf')  # was never merged
    else:
        return cluster.order

类似地，由于Leaf簇没有子节点，因此我们将为此创建并添加一个辅助函数：

from typing import Tuple

def get_children(cluster: Cluster):
    if isinstance(cluster, Leaf):
        raise TypeError("Leaf has no children")
    else:
        return cluster.children

现在我们准备创建聚类算法：

def bottom_up_cluster(inputs: List[Vector],
                      distance_agg: Callable = min) -> Cluster:
    # Start with all leaves
    clusters: List[Cluster] = [Leaf(input) for input in inputs]

    def pair_distance(pair: Tuple[Cluster, Cluster]) -> float:
        return cluster_distance(pair[0], pair[1], distance_agg)

    # as long as we have more than one cluster left...
    while len(clusters) > 1:
        # find the two closest clusters
        c1, c2 = min(((cluster1, cluster2)
                      for i, cluster1 in enumerate(clusters)
                      for cluster2 in clusters[:i]),
                      key=pair_distance)

        # remove them from the list of clusters
        clusters = [c for c in clusters if c != c1 and c != c2]

        # merge them, using merge_order = # of clusters left
        merged_cluster = Merged((c1, c2), order=len(clusters))

        # and add their merge
        clusters.append(merged_cluster)

    # when there's only one cluster left, return it
    return clusters[0]

它的使用非常简单：

base_cluster = bottom_up_cluster(inputs)

这将生成一个如下所示的聚类：

  0  1  2  3  4  5  6  7  8  9  10 11 12 13 14 15 16 17 18
──┬──┬─────┬────────────────────────────────┬───────────┬─ [19, 28]
  │  │     │                                │           └─ [21, 27]
  │  │     │                                └─ [20, 23]
  │  │     └─ [26, 13]
  │  └────────────────────────────────────────────┬─ [11, 15]
  │                                               └─ [13, 13]
  └─────┬─────┬──┬───────────┬─────┬─ [-49, 0]
        │     │  │           │     └─ [-46, 5]
        │     │  │           └─ [-41, 8]
        │     │  └─ [-49, 15]
        │     └─ [-34, 1]
        └───────────┬──┬──┬─────┬─ [-22, -16]
                    │  │  │     └─ [-19, -11]
                    │  │  └─ [-25, -9]
                    │  └─────────────────┬─────┬─────┬─ [-11, -6]
                    │                    │     │     └─ [-12, -8]
                    │                    │     └─ [-14, 5]
                    │                    └─ [-18, -3]
                    └─────────────────┬─ [-13, -19]
                                      └─ [-9, -16]

顶部的数字表示“合并顺序”。因为我们有 20 个输入，所以需要 19 次合并才能得到这一个簇。第一次合并通过组合叶子[19, 28]和[21, 27]创建了簇 18。最后一次合并创建了簇 0。

如果你只想要两个簇，你可以在第一个分叉点“0”处分割，创建一个包含六个点的簇和另一个包含剩余点的簇。对于三个簇，你将继续到第二个分叉点“1”，它指示要将第一个簇拆分为包含([19, 28], [21, 27], [20, 23], [26, 13])的簇和包含([11, 15], [13, 13])的簇。依此类推。

尽管如此，我们通常不想看到像这样的难看文本表示。相反，让我们编写一个函数，通过执行适当数量的拆分来生成任意数量的簇：

def generate_clusters(base_cluster: Cluster,
                      num_clusters: int) -> List[Cluster]:
    # start with a list with just the base cluster
    clusters = [base_cluster]

    # as long as we don't have enough clusters yet...
    while len(clusters) < num_clusters:
        # choose the last-merged of our clusters
        next_cluster = min(clusters, key=get_merge_order)
        # remove it from the list
        clusters = [c for c in clusters if c != next_cluster]

        # and add its children to the list (i.e., unmerge it)
        clusters.extend(get_children(next_cluster))

    # once we have enough clusters...
    return clusters

例如，如果我们想生成三个簇，我们只需执行以下操作：

three_clusters = [get_values(cluster)
                  for cluster in generate_clusters(base_cluster, 3)]

可以轻松绘制的部分：

for i, cluster, marker, color in zip([1, 2, 3],
                                     three_clusters,
                                     ['D','o','*'],
                                     ['r','g','b']):
    xs, ys = zip(*cluster)  # magic unzipping trick
    plt.scatter(xs, ys, color=color, marker=marker)

    # put a number at the mean of the cluster
    x, y = vector_mean(cluster)
    plt.plot(x, y, marker='$' + str(i) + '$', color='black')

plt.title("User Locations -- 3 Bottom-Up Clusters, Min")
plt.xlabel("blocks east of city center")
plt.ylabel("blocks north of city center")
plt.show()

这与k-means 产生了非常不同的结果，如图 20-6 所示。

图 20-6. 使用最小距离生成的三个自下而上的簇

如前所述，这是因为在cluster_distance中使用min倾向于生成链式簇。如果我们改用max（生成紧密簇），它将与 3-means 结果相同（见图 20-7）。

注意

先前的bottom_up_clustering实现相对简单，但效率惊人地低下。特别是，它在每一步重新计算每对输入之间的距离。更高效的实现可能会预先计算每对输入之间的距离，然后在cluster_distance内执行查找。真正高效的实现可能还会记住前一步骤的cluster_distance。

图 20-7. 使用最大距离的三个自底向上聚类

进一步探索

• scikit-learn 有一个完整的模块，sklearn.cluster，其中包含几种聚类算法，包括KMeans和Ward层次聚类算法（其合并集群的标准与我们的不同）。

• SciPy 提供了两种聚类模型：scipy.cluster.vq，实现k-均值；以及scipy.cluster.hierarchy，提供多种层次聚类算法。

第二十一章：自然语言处理

他们已经在语言的盛宴中大快朵颐，窃取了残羹剩饭。

威廉·莎士比亚

自然语言处理（NLP）指的是涉及语言的计算技术。这是一个广泛的领域，但我们将看几种简单和复杂的技术。

词云

在 第 1 章，我们计算了用户兴趣的单词计数。一个可视化单词和计数的方法是 词云，它以比例大小艺术化地呈现单词。

通常，数据科学家们对词云并不看重，主要是因为单词的排列除了“这里是我能放置一个词的空间”之外没有其他意义。

如果你不得不创建一个词云，考虑一下是否可以让坐标轴传达某种信息。例如，想象一下，对于某些数据科学相关的流行术语，你有两个在 0 到 100 之间的数字——第一个表示它在职位发布中出现的频率，第二个表示它在简历中出现的频率：

data = [ ("big data", 100, 15), ("Hadoop", 95, 25), ("Python", 75, 50),
         ("R", 50, 40), ("machine learning", 80, 20), ("statistics", 20, 60),
         ("data science", 60, 70), ("analytics", 90, 3),
         ("team player", 85, 85), ("dynamic", 2, 90), ("synergies", 70, 0),
         ("actionable insights", 40, 30), ("think out of the box", 45, 10),
         ("self-starter", 30, 50), ("customer focus", 65, 15),
         ("thought leadership", 35, 35)]

词云的方法就是在页面上以酷炫的字体排列这些词（图 21-1）。

图 21-1 术语词云

这看起来很整齐，但实际上并没有告诉我们什么。一个更有趣的方法可能是将它们散布开来，使得水平位置表示发布的流行度，垂直位置表示简历的流行度，这将产生一个传达几个见解的可视化效果（图 21-2）：

from matplotlib import pyplot as plt

def text_size(total: int) -> float:
    """equals 8 if total is 0, 28 if total is 200"""
    return 8 + total / 200 * 20

for word, job_popularity, resume_popularity in data:
    plt.text(job_popularity, resume_popularity, word,
             ha='center', va='center',
             size=text_size(job_popularity + resume_popularity))
plt.xlabel("Popularity on Job Postings")
plt.ylabel("Popularity on Resumes")
plt.axis([0, 100, 0, 100])
plt.xticks([])
plt.yticks([])
plt.show()

图 21-2 更有意义（虽然不够吸引人）的词云

n-Gram 语言模型

DataSciencester 搜索引擎市场副总裁希望创建成千上万个关于数据科学的网页，以便您的网站在与数据科学相关的搜索结果中排名更高。（你试图向她解释搜索引擎算法已经足够聪明，这实际上不会起作用，但她拒绝听取。）

当然，她不想写成千上万个网页，也不想支付一大群“内容战略师”来完成。相反，她问你是否可以以某种程序化的方式生成这些网页。为此，我们需要某种语言建模的方法。

一种方法是从一组文档语料库开始，并学习语言的统计模型。在我们的案例中，我们将从迈克·劳凯德斯的文章《什么是数据科学？》开始。

如同 第 9 章，我们将使用 Requests 和 Beautiful Soup 库来获取数据。这里有几个值得注意的问题。

首先，文本中的撇号实际上是 Unicode 字符 u"\u2019"。我们将创建一个辅助函数来将其替换为正常的撇号：

def fix_unicode(text: str) -> str:
    return text.replace(u"\u2019", "'")

第二个问题是，一旦我们获取了网页的文本，我们将希望将其拆分为一系列的单词和句号（以便我们可以知道句子的结束位置）。我们可以使用re.findall来实现这一点：

import re
from bs4 import BeautifulSoup
import requests

url = "https://www.oreilly.com/ideas/what-is-data-science"
html = requests.get(url).text
soup = BeautifulSoup(html, 'html5lib')

content = soup.find("div", "article-body")   # find article-body div
regex = r"[\w']+|[\.]"                       # matches a word or a period

document = []

for paragraph in content("p"):
    words = re.findall(regex, fix_unicode(paragraph.text))
    document.extend(words)

我们当然可以（而且很可能应该）进一步清理这些数据。文档中仍然有一些多余的文本（例如，第一个单词是Section），我们已经在中间句点上分割了（例如，在Web 2.0中），还有一些标题和列表散布在其中。话虽如此，我们将按照文档的样子进行处理。

现在我们将文本作为一系列单词，我们可以按以下方式对语言进行建模：给定一些起始词（比如book），我们查看源文档中跟随它的所有单词。我们随机选择其中一个作为下一个单词，并重复这个过程，直到我们遇到一个句号，这表示句子的结束。我们称之为bigram 模型，因为它完全由原始数据中 bigram（单词对）的频率决定。

起始词呢？我们可以随机从跟在句号后面的单词中选择一个。首先，让我们预先计算可能的单词转换。记住zip在其输入的任何一个完成时停止，因此zip(document, document[1:])给出了document中连续元素的精确配对：

from collections import defaultdict

transitions = defaultdict(list)
for prev, current in zip(document, document[1:]):
    transitions[prev].append(current)

现在我们准备好生成句子了：

def generate_using_bigrams() -> str:
    current = "."   # this means the next word will start a sentence
    result = []
    while True:
        next_word_candidates = transitions[current]    # bigrams (current, _)
        current = random.choice(next_word_candidates)  # choose one at random
        result.append(current)                         # append it to results
        if current == ".": return " ".join(result)     # if "." we're done

它生成的句子是胡言乱语，但如果你试图听起来像数据科学，它们就是你可能会放在你的网站上的那种胡言乱语。例如：

如果你知道你想要对数据进行排序的数据源网页友好的人在热门话题上作为 Hadoop 中的数据，那么数据科学需要一本书来演示为什么可视化是数据的可视化是但我们在 Python 语言中使用许多商业磁盘驱动器上的大量相关性，并创建更可管理的形式进行连接，然后使用它来解决数据的问题。

Bigram 模型

通过查看三元组，连续三个单词的三元组，我们可以使句子不那么胡言乱语。 （更一般地说，你可以查看由n个连续单词组成的n-gram，但对于我们来说，三个就足够了。）现在的转换将取决于前两个单词：

trigram_transitions = defaultdict(list)
starts = []

for prev, current, next in zip(document, document[1:], document[2:]):

    if prev == ".":              # if the previous "word" was a period
        starts.append(current)   # then this is a start word

    trigram_transitions[(prev, current)].append(next)

现在注意，我们现在必须单独跟踪起始词。我们几乎可以以相同的方式生成句子：

def generate_using_trigrams() -> str:
    current = random.choice(starts)   # choose a random starting word
    prev = "."                        # and precede it with a '.'
    result = [current]
    while True:
        next_word_candidates = trigram_transitions[(prev, current)]
        next_word = random.choice(next_word_candidates)

        prev, current = current, next_word
        result.append(current)

        if current == ".":
            return " ".join(result)

这样产生的句子更好，比如：

事后看来 MapReduce 看起来像是一场流行病，如果是这样，那么这是否给我们提供了新的见解，即经济如何运作这不是一个我们几年前甚至可以问的问题已经被工具化。

三元模型

当然，它们听起来更好，因为在每一步生成过程中的选择更少，而在许多步骤中只有一个选择。这意味着我们经常生成句子（或至少是长短语），这些句子在原始数据中原封不动地出现过。拥有更多的数据会有所帮助；如果我们收集了关于数据科学的多篇文章中的n-gram，它也会更有效。

语法

一种不同的语言建模方法是使用语法，即生成可接受句子的规则。在小学时，你可能学过词性及其如何组合。例如，如果你有一个非常糟糕的英语老师，你可能会说一个句子必然由名词后跟一个动词。如果你有名词和动词的列表，你可以根据这个规则生成句子。

我们将定义一个稍微复杂的语法：

from typing import List, Dict

# Type alias to refer to grammars later
Grammar = Dict[str, List[str]]

grammar = {
    "_S"  : ["_NP _VP"],
    "_NP" : ["_N",
             "_A _NP _P _A _N"],
    "_VP" : ["_V",
             "_V _NP"],
    "_N"  : ["data science", "Python", "regression"],
    "_A"  : ["big", "linear", "logistic"],
    "_P"  : ["about", "near"],
    "_V"  : ["learns", "trains", "tests", "is"]
}

我编制了以下约定：以下划线开头的名称指的是需要进一步扩展的规则，而其他名称是不需要进一步处理的终端。

因此，例如，"_S"是“句子”规则，它生成一个"_NP"（“名词短语”）规则后跟一个"_VP"（“动词短语”）规则。

动词短语规则可以生成"_V"（“动词”）规则，或者动词规则后跟名词短语规则。

注意"_NP"规则包含在其一个产生中。语法可以是递归的，这使得像这样的有限语法可以生成无限多个不同的句子。

我们如何从这个语法生成句子？我们将从包含句子规则["_S"]的列表开始。然后，我们将通过用其产生之一随机替换每个规则来重复扩展每个规则。当我们有一个完全由终端组成的列表时，我们停止。

例如，这样的进展可能看起来像：

['_S']
['_NP','_VP']
['_N','_VP']
['Python','_VP']
['Python','_V','_NP']
['Python','trains','_NP']
['Python','trains','_A','_NP','_P','_A','_N']
['Python','trains','logistic','_NP','_P','_A','_N']
['Python','trains','logistic','_N','_P','_A','_N']
['Python','trains','logistic','data science','_P','_A','_N']
['Python','trains','logistic','data science','about','_A', '_N']
['Python','trains','logistic','data science','about','logistic','_N']
['Python','trains','logistic','data science','about','logistic','Python']

我们如何实现这一点？嗯，首先，我们将创建一个简单的辅助函数来识别终端：

def is_terminal(token: str) -> bool:
    return token[0] != "_"

接下来，我们需要编写一个函数，将标记列表转换为句子。我们将寻找第一个非终结符标记。如果找不到一个，那意味着我们有一个完成的句子，我们就完成了。

如果我们找到一个非终结符，然后我们随机选择它的一个产生式。如果该产生式是一个终端（即一个词），我们只需用它替换标记。否则，它是一系列以空格分隔的非终结符标记，我们需要split然后在当前标记中插入。无论哪种方式，我们都会在新的标记集上重复这个过程。

将所有这些放在一起，我们得到：

def expand(grammar: Grammar, tokens: List[str]) -> List[str]:
    for i, token in enumerate(tokens):
        # If this is a terminal token, skip it.
        if is_terminal(token): continue

        # Otherwise, it's a nonterminal token,
        # so we need to choose a replacement at random.
        replacement = random.choice(grammar[token])

        if is_terminal(replacement):
            tokens[i] = replacement
        else:
            # Replacement could be, e.g., "_NP _VP", so we need to
            # split it on spaces and splice it in.
            tokens = tokens[:i] + replacement.split() + tokens[(i+1):]

        # Now call expand on the new list of tokens.
        return expand(grammar, tokens)

    # If we get here, we had all terminals and are done.
    return tokens

现在我们可以开始生成句子了：

def generate_sentence(grammar: Grammar) -> List[str]:
    return expand(grammar, ["_S"])

尝试改变语法——添加更多单词，添加更多规则，添加你自己的词性——直到你准备生成公司所需的多个网页为止。

当语法反向使用时，语法实际上更有趣。给定一个句子，我们可以使用语法解析句子。这然后允许我们识别主语和动词，并帮助我们理解句子的意义。

使用数据科学生成文本是一个很棒的技巧；使用它来理解文本更加神奇。（参见“进一步探索”可以用于此目的的库。）

一个旁注：吉布斯采样

从一些分布生成样本很容易。我们可以得到均匀随机变量：

random.random()

和正常的随机变量一样：

inverse_normal_cdf(random.random())

但是一些分布很难进行抽样。吉布斯抽样 是一种从多维分布生成样本的技术，当我们只知道一些条件分布时使用。

例如，想象掷两个骰子。让 x 是第一个骰子的值，y 是两个骰子的和，想象你想生成大量 (x, y) 对。在这种情况下，直接生成样本是很容易的：

from typing import Tuple
import random

def roll_a_die() -> int:
    return random.choice([1, 2, 3, 4, 5, 6])

def direct_sample() -> Tuple[int, int]:
    d1 = roll_a_die()
    d2 = roll_a_die()
    return d1, d1 + d2

但是假设你只知道条件分布。知道 x 的值时，y 的分布很简单——如果你知道 x 的值，y 同样可能是 x + 1、x + 2、x + 3、x + 4、x + 5 或 x + 6：

def random_y_given_x(x: int) -> int:
    """equally likely to be x + 1, x + 2, ... , x + 6"""
    return x + roll_a_die()

另一个方向更为复杂。例如，如果你知道 y 是 2，则必然 x 是 1（因为使两个骰子的和为 2 的唯一方法是它们都是 1）。如果你知道 y 是 3，则 x 等可能是 1 或 2。同样，如果 y 是 11，则 x 必须是 5 或 6：

def random_x_given_y(y: int) -> int:
    if y <= 7:
        # if the total is 7 or less, the first die is equally likely to be
        # 1, 2, ..., (total - 1)
        return random.randrange(1, y)
    else:
        # if the total is 7 or more, the first die is equally likely to be
        # (total - 6), (total - 5), ..., 6
        return random.randrange(y - 6, 7)

吉布斯抽样的工作原理是我们从任意（有效的）x 和 y 的值开始，然后反复替换，用在y 条件下随机选择的值替换 x，并在 x 条件下随机选择的值替换 y。经过若干次迭代，x 和 y 的结果值将代表无条件联合分布的一个样本：

def gibbs_sample(num_iters: int = 100) -> Tuple[int, int]:
    x, y = 1, 2 # doesn't really matter
    for _ in range(num_iters):
        x = random_x_given_y(y)
        y = random_y_given_x(x)
    return x, y

你可以检查这是否给出与直接样本相似的结果：

def compare_distributions(num_samples: int = 1000) -> Dict[int, List[int]]:
    counts = defaultdict(lambda: [0, 0])
    for _ in range(num_samples):
        counts[gibbs_sample()][0] += 1
        counts[direct_sample()][1] += 1
    return counts

我们将在下一节中使用这种技术。

主题建模

当我们在第一章中构建“您可能认识的数据科学家”推荐系统时，我们简单地查找人们声明的兴趣的完全匹配。

更复杂的方法是尝试理解用户兴趣的主题。一种称为潜在狄利克雷分配（LDA）的技术通常用于识别一组文档中的常见主题。我们将其应用于由每个用户兴趣组成的文档。

LDA 与我们在第十三章中构建的朴素贝叶斯分类器有一些相似之处，因为它假设文档的概率模型。对于我们的目的，该模型假设以下内容，我们将略过更复杂的数学细节：

• 有一些固定数量 K 的主题。

• 有一个随机变量为每个主题分配与之相关联的单词概率分布。你应该将这个分布看作是给定主题 k 下看到单词 w 的概率。

• 还有一个随机变量为每个文档分配一个主题的概率分布。你应该将这个分布看作是文档 d 中主题的混合。

• 文档中的每个单词是通过首先随机选择一个主题（从文档的主题分布中）然后随机选择一个单词（从主题的单词分布中）生成的。

特别是，我们有一个documents的集合，每个文档都是一个单词的list。并且我们有一个相应的document_topics集合，它为每个文档中的每个单词分配一个主题（这里是 0 到K-1 之间的数字）。

因此，第四个文档中的第五个单词是：

documents[3][4]

选择该单词的主题是：

document_topics[3][4]

这非常明确地定义了每个文档在主题上的分布，并且隐含地定义了每个主题在单词上的分布。

通过比较主题 1 生成该单词的次数与主题 1 生成任何单词的次数，我们可以估计主题 1 生成某个单词的可能性。（类似地，在第十三章中建立垃圾邮件过滤器时，我们比较了每个单词在垃圾邮件中出现的次数与垃圾邮件中出现的总字数。）

虽然这些主题只是数字，但我们可以通过查看它们赋予最高权重的单词来为它们命名描述性名称。我们只需以某种方式生成document_topics。这就是吉布斯抽样发挥作用的地方。

我们首先随机地为每个文档中的每个单词分配一个主题。现在我们逐个单词地遍历每个文档。对于该单词和文档，我们为每个主题构造依赖于该文档中主题的（当前）分布和该主题中单词的（当前）分布的权重。然后我们使用这些权重来对该单词抽样一个新的主题。如果我们多次迭代这个过程，我们最终会得到从主题-单词分布和文档-主题分布的联合样本。

起步，我们需要一个函数根据任意一组权重随机选择一个索引：

def sample_from(weights: List[float]) -> int:
    """returns i with probability weights[i] / sum(weights)"""
    total = sum(weights)
    rnd = total * random.random()      # uniform between 0 and total
    for i, w in enumerate(weights):
        rnd -= w                       # return the smallest i such that
        if rnd <= 0: return i          # weights[0] + ... + weights[i] >= rnd

例如，如果给定权重 [1, 1, 3]，那么它将返回 0 的概率为五分之一，返回 1 的概率为五分之一，返回 2 的概率为三分之五。让我们编写一个测试：

from collections import Counter

# Draw 1000 times and count
draws = Counter(sample_from([0.1, 0.1, 0.8]) for _ in range(1000))
assert 10 < draws[0] < 190   # should be ~10%, this is a really loose test
assert 10 < draws[1] < 190   # should be ~10%, this is a really loose test
assert 650 < draws[2] < 950  # should be ~80%, this is a really loose test
assert draws[0] + draws[1] + draws[2] == 1000

我们的文档是我们用户的兴趣，看起来像：

documents = [
    ["Hadoop", "Big Data", "HBase", "Java", "Spark", "Storm", "Cassandra"],
    ["NoSQL", "MongoDB", "Cassandra", "HBase", "Postgres"],
    ["Python", "scikit-learn", "scipy", "numpy", "statsmodels", "pandas"],
    ["R", "Python", "statistics", "regression", "probability"],
    ["machine learning", "regression", "decision trees", "libsvm"],
    ["Python", "R", "Java", "C++", "Haskell", "programming languages"],
    ["statistics", "probability", "mathematics", "theory"],
    ["machine learning", "scikit-learn", "Mahout", "neural networks"],
    ["neural networks", "deep learning", "Big Data", "artificial intelligence"],
    ["Hadoop", "Java", "MapReduce", "Big Data"],
    ["statistics", "R", "statsmodels"],
    ["C++", "deep learning", "artificial intelligence", "probability"],
    ["pandas", "R", "Python"],
    ["databases", "HBase", "Postgres", "MySQL", "MongoDB"],
    ["libsvm", "regression", "support vector machines"]
]

我们将尝试找到：

K = 4

主题。为了计算抽样权重，我们需要跟踪几个计数。让我们首先为它们创建数据结构。

• 每个主题分配给每个文档的次数是：

  # a list of Counters, one for each document
  document_topic_counts = [Counter() for _ in documents]
  

• 每个单词分配到每个主题的次数是：

  # a list of Counters, one for each topic
  topic_word_counts = [Counter() for _ in range(K)]
  

• 分配给每个主题的总字数是：

  # a list of numbers, one for each topic
  topic_counts = [0 for _ in range(K)]
  

• 每个文档包含的总字数是：

  # a list of numbers, one for each document
  document_lengths = [len(document) for document in documents]
  

• 不同单词的数量是：

  distinct_words = set(word for document in documents for word in document)
  W = len(distinct_words)
  

• 以及文档的数量：

  D = len(documents)
  

一旦我们填充了这些数据，我们可以如下找到例如与主题 1 相关联的documents[3]中的单词数：

document_topic_counts[3][1]

并且我们可以找到nlp与主题 2 相关联的次数如下：

topic_word_counts[2]["nlp"]

现在我们准备定义我们的条件概率函数。就像第十三章中那样，每个函数都有一个平滑项，确保每个主题在任何文档中被选择的概率都不为零，并且每个单词在任何主题中被选择的概率也不为零：

def p_topic_given_document(topic: int, d: int, alpha: float = 0.1) -> float:
    """
 The fraction of words in document 'd'
 that are assigned to 'topic' (plus some smoothing)
 """
    return ((document_topic_counts[d][topic] + alpha) /
            (document_lengths[d] + K * alpha))

def p_word_given_topic(word: str, topic: int, beta: float = 0.1) -> float:
    """
 The fraction of words assigned to 'topic'
 that equal 'word' (plus some smoothing)
 """
    return ((topic_word_counts[topic][word] + beta) /
            (topic_counts[topic] + W * beta))

我们将使用这些函数来创建更新主题的权重：

def topic_weight(d: int, word: str, k: int) -> float:
    """
 Given a document and a word in that document,
 return the weight for the kth topic
 """
    return p_word_given_topic(word, k) * p_topic_given_document(k, d)

def choose_new_topic(d: int, word: str) -> int:
    return sample_from([topic_weight(d, word, k)
                        for k in range(K)])

有坚实的数学原因解释为什么 topic_weight 被定义为它的方式，但是它们的细节会让我们走得太远。希望至少直观地理解，鉴于一个词和它的文档，选择任何主题的可能性取决于该主题对文档的可能性以及该词对该主题的可能性。

这就是我们需要的所有机制。我们从将每个单词分配给一个随机主题开始，并相应地填充我们的计数器：

random.seed(0)
document_topics = [[random.randrange(K) for word in document]
                   for document in documents]

for d in range(D):
    for word, topic in zip(documents[d], document_topics[d]):
        document_topic_counts[d][topic] += 1
        topic_word_counts[topic][word] += 1
        topic_counts[topic] += 1

我们的目标是获取主题-词分布和文档-主题分布的联合样本。我们使用一种基于吉布斯抽样的形式来完成这一过程，该过程使用之前定义的条件概率：

import tqdm

for iter in tqdm.trange(1000):
    for d in range(D):
        for i, (word, topic) in enumerate(zip(documents[d],
                                              document_topics[d])):

            # remove this word / topic from the counts
            # so that it doesn't influence the weights
            document_topic_counts[d][topic] -= 1
            topic_word_counts[topic][word] -= 1
            topic_counts[topic] -= 1
            document_lengths[d] -= 1

            # choose a new topic based on the weights
            new_topic = choose_new_topic(d, word)
            document_topics[d][i] = new_topic

            # and now add it back to the counts
            document_topic_counts[d][new_topic] += 1
            topic_word_counts[new_topic][word] += 1
            topic_counts[new_topic] += 1
            document_lengths[d] += 1

主题是什么？它们只是数字 0、1、2 和 3。如果我们想要它们的名称，我们必须自己添加。让我们看看每个主题的五个权重最高的词汇（表 21-1）：

for k, word_counts in enumerate(topic_word_counts):
    for word, count in word_counts.most_common():
        if count > 0:
            print(k, word, count)

表 21-1. 每个主题的最常见词汇

主题 0	主题 1	主题 2	主题 3
Java	R	HBase	回归分析
大数据	统计	Postgres	libsvm
Hadoop	Python	MongoDB	scikit-learn
深度学习	概率	卡桑德拉	机器学习
人工智能	熊猫	NoSQL	神经网络

基于这些，我可能会分配主题名称：

topic_names = ["Big Data and programming languages",
               "Python and statistics",
               "databases",
               "machine learning"]

在这一点上，我们可以看到模型如何将主题分配给每个用户的兴趣：

for document, topic_counts in zip(documents, document_topic_counts):
    print(document)
    for topic, count in topic_counts.most_common():
        if count > 0:
            print(topic_names[topic], count)
    print()

这给出了：

['Hadoop', 'Big Data', 'HBase', 'Java', 'Spark', 'Storm', 'Cassandra']
Big Data and programming languages 4 databases 3
['NoSQL', 'MongoDB', 'Cassandra', 'HBase', 'Postgres']
databases 5
['Python', 'scikit-learn', 'scipy', 'numpy', 'statsmodels', 'pandas']
Python and statistics 5 machine learning 1

等等。考虑到我们在一些主题名称中需要使用的“和”，可能我们应该使用更多的主题，尽管最可能我们没有足够的数据成功学习它们。

单词向量

最近自然语言处理的许多进展涉及深度学习。在本章的其余部分中，我们将使用我们在第十九章中开发的机制来看一些这样的进展。

一个重要的创新涉及将单词表示为低维向量。这些向量可以进行比较、相加、输入到机器学习模型中，或者任何你想做的事情。它们通常具有良好的特性；例如，相似的单词倾向于有相似的向量。也就是说，通常单词 big 的向量与单词 large 的向量非常接近，因此一个操作单词向量的模型可以（在某种程度上）免费处理类似的词语。

经常向量也会展示出令人愉悦的算术特性。例如，在某些模型中，如果你取 king 的向量，减去 man 的向量，再加上 woman 的向量，你将得到一个非常接近 queen 向量的向量。思考这对于单词向量实际上“学到”了什么，可能会很有趣，尽管我们在这里不会花时间讨论这一点。

对于一个庞大的词汇表来说，设计这样的向量是一项困难的任务，所以通常我们会从文本语料库中 学习 它们。有几种不同的方案，但在高层次上，任务通常看起来像这样：

1. 获取一堆文本。

2. 创建一个数据集，目标是预测给定附近单词的单词（或者，预测给定单词的附近单词）。

3. 训练一个神经网络在这个任务上表现良好。

4. 将训练好的神经网络的内部状态作为单词向量。

特别是，由于任务是根据附近的单词预测单词，出现在类似上下文中的单词（因此具有类似的附近单词）应该具有类似的内部状态，因此也应该具有相似的单词向量。

我们将使用 余弦相似度（cosine similarity）来衡量“相似性”，它是一个介于-1 和 1 之间的数值，用于衡量两个向量指向相同方向的程度：

from scratch.linear_algebra import dot, Vector
import math

def cosine_similarity(v1: Vector, v2: Vector) -> float:
    return dot(v1, v2) / math.sqrt(dot(v1, v1) * dot(v2, v2))

assert cosine_similarity([1., 1, 1], [2., 2, 2]) == 1, "same direction"
assert cosine_similarity([-1., -1], [2., 2]) == -1,    "opposite direction"
assert cosine_similarity([1., 0], [0., 1]) == 0,       "orthogonal"

让我们学习一些词向量，看看它是如何工作的。

首先，我们需要一个玩具数据集。通常使用的单词向量通常是通过在数百万甚至数十亿个单词上训练而来的。由于我们的玩具库无法处理那么多数据，我们将创建一个具有某些结构的人工数据集：

colors = ["red", "green", "blue", "yellow", "black", ""]
nouns = ["bed", "car", "boat", "cat"]
verbs = ["is", "was", "seems"]
adverbs = ["very", "quite", "extremely", ""]
adjectives = ["slow", "fast", "soft", "hard"]

def make_sentence() -> str:
    return " ".join([
        "The",
        random.choice(colors),
        random.choice(nouns),
        random.choice(verbs),
        random.choice(adverbs),
        random.choice(adjectives),
        "."
    ])

NUM_SENTENCES = 50

random.seed(0)
sentences = [make_sentence() for _ in range(NUM_SENTENCES)]

这将生成许多具有类似结构但不同单词的句子；例如，“绿色的船似乎相当慢。” 在这种设置下，颜色将主要出现在“相似”的上下文中，名词也是如此，依此类推。因此，如果我们成功地分配了单词向量，颜色应该会得到相似的向量，依此类推。

注意

在实际使用中，您可能会有数百万个句子的语料库，在这种情况下，您将从句子中获得“足够的”上下文。在这里，我们只有 50 个句子，我们必须使它们有些人为的。

如前所述，我们将需要对我们的单词进行一位有效编码，这意味着我们需要将它们转换为 ID。我们将引入一个Vocabulary类来跟踪这个映射：

from scratch.deep_learning import Tensor

class Vocabulary:
    def __init__(self, words: List[str] = None) -> None:
        self.w2i: Dict[str, int] = {}  # mapping word -> word_id
        self.i2w: Dict[int, str] = {}  # mapping word_id -> word

        for word in (words or []):     # If words were provided,
            self.add(word)             # add them.

    @property
    def size(self) -> int:
        """how many words are in the vocabulary"""
        return len(self.w2i)

    def add(self, word: str) -> None:
        if word not in self.w2i:        # If the word is new to us:
            word_id = len(self.w2i)     # Find the next id.
            self.w2i[word] = word_id    # Add to the word -> word_id map.
            self.i2w[word_id] = word    # Add to the word_id -> word map.

    def get_id(self, word: str) -> int:
        """return the id of the word (or None)"""
        return self.w2i.get(word)

    def get_word(self, word_id: int) -> str:
        """return the word with the given id (or None)"""
        return self.i2w.get(word_id)

    def one_hot_encode(self, word: str) -> Tensor:
        word_id = self.get_id(word)
        assert word_id is not None, f"unknown word {word}"

        return [1.0 if i == word_id else 0.0 for i in range(self.size)]

这些都是我们可以手动完成的事情，但将它们放在一个类中很方便。我们可能应该测试它：

vocab = Vocabulary(["a", "b", "c"])
assert vocab.size == 3,              "there are 3 words in the vocab"
assert vocab.get_id("b") == 1,       "b should have word_id 1"
assert vocab.one_hot_encode("b") == [0, 1, 0]
assert vocab.get_id("z") is None,    "z is not in the vocab"
assert vocab.get_word(2) == "c",     "word_id 2 should be c"
vocab.add("z")
assert vocab.size == 4,              "now there are 4 words in the vocab"
assert vocab.get_id("z") == 3,       "now z should have id 3"
assert vocab.one_hot_encode("z") == [0, 0, 0, 1]

我们还应该编写简单的辅助函数来保存和加载词汇表，就像我们为深度学习模型所做的那样：

import json

def save_vocab(vocab: Vocabulary, filename: str) -> None:
    with open(filename, 'w') as f:
        json.dump(vocab.w2i, f)       # Only need to save w2i

def load_vocab(filename: str) -> Vocabulary:
    vocab = Vocabulary()
    with open(filename) as f:
        # Load w2i and generate i2w from it
        vocab.w2i = json.load(f)
        vocab.i2w = {id: word for word, id in vocab.w2i.items()}
    return vocab

我们将使用一个称为 skip-gram 的词向量模型，它以单词作为输入，并生成可能性，表明哪些单词可能会在附近出现。我们将向其提供训练对 (单词, 附近单词) 并尝试最小化 SoftmaxCrossEntropy 损失。

注意

另一个常见的模型，连续词袋（continuous bag-of-words，CBOW），将附近的单词作为输入，并尝试预测原始单词。

让我们设计我们的神经网络。在其核心将是一个 嵌入（embedding）层，它以单词 ID 作为输入并返回一个单词向量。在内部，我们可以简单地使用查找表来实现这一点。

然后，我们将单词向量传递给一个具有与我们词汇表中单词数量相同的输出的 Linear 层。与以前一样，我们将使用 softmax 将这些输出转换为附近单词的概率。当我们使用梯度下降训练模型时，我们将更新查找表中的向量。训练完成后，该查找表为我们提供了单词向量。

让我们创建那个嵌入层。实际上，我们可能希望嵌入除单词之外的其他内容，因此我们将构建一个更通用的Embedding层。（稍后我们将编写一个TextEmbedding子类，专门用于词向量。）

在其构造函数中，我们将提供我们嵌入向量的数量和维度，因此它可以创建嵌入（最初将是标准的随机正态分布）：

from typing import Iterable
from scratch.deep_learning import Layer, Tensor, random_tensor, zeros_like

class Embedding(Layer):
    def __init__(self, num_embeddings: int, embedding_dim: int) -> None:
        self.num_embeddings = num_embeddings
        self.embedding_dim = embedding_dim

        # One vector of size embedding_dim for each desired embedding
        self.embeddings = random_tensor(num_embeddings, embedding_dim)
        self.grad = zeros_like(self.embeddings)

        # Save last input id
        self.last_input_id = None

在我们的情况下，我们一次只嵌入一个词。然而，在其他模型中，我们可能希望嵌入一个词序列并返回一个词向量序列。（例如，如果我们想要训练前面描述的 CBOW 模型。）因此，另一种设计将采用单词 ID 序列。我们将坚持一次只处理一个，以简化事务。

    def forward(self, input_id: int) -> Tensor:
        """Just select the embedding vector corresponding to the input id"""
        self.input_id = input_id    # remember for use in backpropagation

        return self.embeddings[input_id]

对于反向传播，我们将获得对应于所选嵌入向量的梯度，并且我们需要为self.embeddings构建相应的梯度，对于除所选之外的每个嵌入，其梯度都为零：

    def backward(self, gradient: Tensor) -> None:
        # Zero out the gradient corresponding to the last input.
        # This is way cheaper than creating a new all-zero tensor each time.
        if self.last_input_id is not None:
            zero_row = [0 for _ in range(self.embedding_dim)]
            self.grad[self.last_input_id] = zero_row

        self.last_input_id = self.input_id
        self.grad[self.input_id] = gradient

因为我们有参数和梯度，我们需要重写那些方法：

    def params(self) -> Iterable[Tensor]:
        return [self.embeddings]

    def grads(self) -> Iterable[Tensor]:
        return [self.grad]

如前所述，我们需要一个专门用于词向量的子类。在这种情况下，我们的嵌入数量由我们的词汇决定，所以让我们直接传入它：

class TextEmbedding(Embedding):
    def __init__(self, vocab: Vocabulary, embedding_dim: int) -> None:
        # Call the superclass constructor
        super().__init__(vocab.size, embedding_dim)

        # And hang onto the vocab
        self.vocab = vocab

所有其他内置方法都可以原样工作，但我们将添加一些特定于文本处理的方法。例如，我们希望能够检索给定单词的向量。（这不是Layer接口的一部分，但我们始终可以根据需要向特定层添加额外的方法。）

    def __getitem__(self, word: str) -> Tensor:
        word_id = self.vocab.get_id(word)
        if word_id is not None:
            return self.embeddings[word_id]
        else:
            return None

这个 dunder 方法将允许我们使用索引检索单词向量：

word_vector = embedding["black"]

我们还希望嵌入层告诉我们给定单词的最接近的单词：

    def closest(self, word: str, n: int = 5) -> List[Tuple[float, str]]:
        """Returns the n closest words based on cosine similarity"""
        vector = self[word]

        # Compute pairs (similarity, other_word), and sort most similar first
        scores = [(cosine_similarity(vector, self.embeddings[i]), other_word)
                  for other_word, i in self.vocab.w2i.items()]
        scores.sort(reverse=True)

        return scores[:n]

我们的嵌入层只输出向量，我们可以将其馈送到Linear层中。

现在我们准备好组装我们的训练数据。对于每个输入单词，我们将选择其左边的两个单词和右边的两个单词作为目标单词。

让我们从将句子转换为小写并拆分为单词开始：

import re

# This is not a great regex, but it works on our data.
tokenized_sentences = [re.findall("[a-z]+|[.]", sentence.lower())
                       for sentence in sentences]

在此之后，我们可以构建一个词汇表：

# Create a vocabulary (that is, a mapping word -> word_id) based on our text.
vocab = Vocabulary(word
                   for sentence_words in tokenized_sentences
                   for word in sentence_words)

现在我们可以创建训练数据：

from scratch.deep_learning import Tensor, one_hot_encode

inputs: List[int] = []
targets: List[Tensor] = []

for sentence in tokenized_sentences:
    for i, word in enumerate(sentence):          # For each word
        for j in [i - 2, i - 1, i + 1, i + 2]:   # take the nearby locations
            if 0 <= j < len(sentence):           # that aren't out of bounds
                nearby_word = sentence[j]        # and get those words.

                # Add an input that's the original word_id
                inputs.append(vocab.get_id(word))

                # Add a target that's the one-hot-encoded nearby word
                targets.append(vocab.one_hot_encode(nearby_word))

利用我们建立的机制，现在很容易创建我们的模型：

from scratch.deep_learning import Sequential, Linear

random.seed(0)
EMBEDDING_DIM = 5  # seems like a good size

# Define the embedding layer separately, so we can reference it.
embedding = TextEmbedding(vocab=vocab, embedding_dim=EMBEDDING_DIM)

model = Sequential([
    # Given a word (as a vector of word_ids), look up its embedding.
    embedding,
    # And use a linear layer to compute scores for "nearby words."
    Linear(input_dim=EMBEDDING_DIM, output_dim=vocab.size)
])

使用来自第十九章的工具，训练我们的模型非常容易：

from scratch.deep_learning import SoftmaxCrossEntropy, Momentum, GradientDescent

loss = SoftmaxCrossEntropy()
optimizer = GradientDescent(learning_rate=0.01)

for epoch in range(100):
    epoch_loss = 0.0
    for input, target in zip(inputs, targets):
        predicted = model.forward(input)
        epoch_loss += loss.loss(predicted, target)
        gradient = loss.gradient(predicted, target)
        model.backward(gradient)
        optimizer.step(model)
    print(epoch, epoch_loss)            # Print the loss
    print(embedding.closest("black"))   # and also a few nearest words
    print(embedding.closest("slow"))    # so we can see what's being
    print(embedding.closest("car"))     # learned.

当你看着这个训练过程时，你会看到颜色变得越来越接近，形容词变得越来越接近，名词也变得越来越接近。

模型训练完成后，探索最相似的单词是件有趣的事情：

pairs = [(cosine_similarity(embedding[w1], embedding[w2]), w1, w2)
         for w1 in vocab.w2i
         for w2 in vocab.w2i
         if w1 < w2]
pairs.sort(reverse=True)
print(pairs[:5])

这对我来说结果如下：

[(0.9980283554864815, 'boat', 'car'),
 (0.9975147744587706, 'bed', 'cat'),
 (0.9953153441218054, 'seems', 'was'),
 (0.9927107440377975, 'extremely', 'quite'),
 (0.9836183658415987, 'bed', 'car')]

（显然bed和cat并不真正相似，但在我们的训练句子中它们似乎相似，并且模型捕捉到了这一点。）

我们还可以提取前两个主成分并将它们绘制出来：

from scratch.working_with_data import pca, transform
import matplotlib.pyplot as plt

# Extract the first two principal components and transform the word vectors
components = pca(embedding.embeddings, 2)
transformed = transform(embedding.embeddings, components)

# Scatter the points (and make them white so they're "invisible")
fig, ax = plt.subplots()
ax.scatter(*zip(*transformed), marker='.', color='w')

# Add annotations for each word at its transformed location
for word, idx in vocab.w2i.items():
    ax.annotate(word, transformed[idx])

# And hide the axes
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)

plt.show()

这表明相似的单词确实聚集在一起（参见图 21-3）：

图 21-3. 单词向量

如果你感兴趣，训练 CBOW 词向量并不难。你需要做一些工作。首先，你需要修改Embedding层，使其接受一个ID 列表作为输入，并输出一个嵌入向量列表。然后你需要创建一个新的层（Sum？），它接受一个向量列表并返回它们的总和。

每个单词表示一个训练示例，其中输入是周围单词的单词 ID，目标是单词本身的独热编码。

修改后的Embedding层将周围的单词转换为向量列表，新的Sum层将向量列表合并为单个向量，然后Linear层可以生成分数，这些分数可以经过softmax处理，得到表示“在这个上下文中最可能的单词”的分布。

我发现 CBOW 模型比跳字模型更难训练，但我鼓励你去尝试一下。

循环神经网络

我们在前一节中开发的单词向量通常用作神经网络的输入。做这件事的一个挑战是句子的长度不同：你可以将一个三个单词的句子想象为一个[3, embedding_dim]张量，而一个十个单词的句子想象为一个[10, embedding_dim]张量。为了，比如，将它们传递给Linear层，我们首先需要处理第一个可变长度维度。

一个选择是使用Sum层（或者一个接受平均值的变体）；然而，句子中单词的顺序通常对其含义很重要。以一个常见的例子来说，“狗咬人”和“人咬狗”是两个非常不同的故事！

处理这个问题的另一种方法是使用循环神经网络（RNNs），它们具有它们在输入之间保持的隐藏状态。在最简单的情况下，每个输入与当前隐藏状态结合以产生输出，然后将其用作新的隐藏状态。这允许这些网络在某种意义上“记住”它们看到的输入，并建立到依赖于所有输入及其顺序的最终输出。

我们将创建一个非常简单的 RNN 层，它将接受单个输入（例如句子中的一个单词或一个单词中的一个字符），并在调用之间保持其隐藏状态。

回想一下，我们的Linear层有一些权重，w，和一个偏置，b。它接受一个向量input并使用逻辑生成不同的向量作为output：

output[o] = dot(w[o], input) + b[o]

这里我们将要整合我们的隐藏状态，因此我们将有两组权重——一组用于应用于input，另一组用于前一个hidden状态：

output[o] = dot(w[o], input) + dot(u[o], hidden) + b[o]

接下来，我们将使用output向量作为新的hidden值。这并不是一个巨大的改变，但它将使我们的网络能够做出奇妙的事情。

from scratch.deep_learning import tensor_apply, tanh

class SimpleRnn(Layer):
    """Just about the simplest possible recurrent layer."""
    def __init__(self, input_dim: int, hidden_dim: int) -> None:
        self.input_dim = input_dim
        self.hidden_dim = hidden_dim

        self.w = random_tensor(hidden_dim, input_dim, init='xavier')
        self.u = random_tensor(hidden_dim, hidden_dim, init='xavier')
        self.b = random_tensor(hidden_dim)

        self.reset_hidden_state()

    def reset_hidden_state(self) -> None:
        self.hidden = [0 for _ in range(self.hidden_dim)]

你可以看到，我们将隐藏状态初始为一个 0 向量，并提供一个函数，供使用网络的人调用以重置隐藏状态。

在这种设置下，forward函数相对直接（至少，如果你记得并理解我们的Linear层是如何工作的话）：

    def forward(self, input: Tensor) -> Tensor:
        self.input = input              # Save both input and previous
        self.prev_hidden = self.hidden  # hidden state to use in backprop.

        a = [(dot(self.w[h], input) +           # weights @ input
              dot(self.u[h], self.hidden) +     # weights @ hidden
              self.b[h])                        # bias
             for h in range(self.hidden_dim)]

        self.hidden = tensor_apply(tanh, a)  # Apply tanh activation
        return self.hidden                   # and return the result.

backward传递类似于我们Linear层中的传递，只是需要计算额外的u权重的梯度：

    def backward(self, gradient: Tensor):
        # Backpropagate through the tanh
        a_grad = [gradient[h] * (1 - self.hidden[h] ** 2)
                  for h in range(self.hidden_dim)]

        # b has the same gradient as a
        self.b_grad = a_grad

        # Each w[h][i] is multiplied by input[i] and added to a[h],
        # so each w_grad[h][i] = a_grad[h] * input[i]
        self.w_grad = [[a_grad[h] * self.input[i]
                        for i in range(self.input_dim)]
                       for h in range(self.hidden_dim)]

        # Each u[h][h2] is multiplied by hidden[h2] and added to a[h],
        # so each u_grad[h][h2] = a_grad[h] * prev_hidden[h2]
        self.u_grad = [[a_grad[h] * self.prev_hidden[h2]
                        for h2 in range(self.hidden_dim)]
                       for h in range(self.hidden_dim)]

        # Each input[i] is multiplied by every w[h][i] and added to a[h],
        # so each input_grad[i] = sum(a_grad[h] * w[h][i] for h in ...)
        return [sum(a_grad[h] * self.w[h][i] for h in range(self.hidden_dim))
                for i in range(self.input_dim)]

最后，我们需要重写params和grads方法：

    def params(self) -> Iterable[Tensor]:
        return [self.w, self.u, self.b]

    def grads(self) -> Iterable[Tensor]:
        return [self.w_grad, self.u_grad, self.b_grad]

警告

这个“简单”的 RNN 实在太简单了，你可能不应该在实践中使用它。

我们的SimpleRnn有几个不理想的特性。其中一个是每次调用它时，它的整个隐藏状态都用来更新输入。另一个是每次调用它时，整个隐藏状态都会被覆盖。这两点使得训练变得困难；特别是，它使模型难以学习长期依赖性。

因此，几乎没有人使用这种简单的 RNN。相反，他们使用更复杂的变体，如 LSTM（“长短期记忆”）或 GRU（“门控循环单元”），这些变体有更多的参数，并使用参数化的“门”来允许每个时间步只更新一部分状态（并且只使用一部分状态）。

这些变体并没有什么特别的困难；然而，它们涉及更多的代码，我认为阅读起来并不会相应更具有教育性。本章的代码可以在GitHub找到，其中包括了 LSTM 的实现。我鼓励你去看看，但这有点乏味，所以我们在这里不再详细讨论。

我们实现的另一个怪癖是，它每次只处理一个“步骤”，并且需要我们手动重置隐藏状态。一个更实用的 RNN 实现可以接受输入序列，将其隐藏状态在每个序列开始时设为 0，并生成输出序列。我们的实现肯定可以修改成这种方式；同样地，这将需要更多的代码和复杂性，而对理解的帮助不大。

示例：使用字符级别的循环神经网络

新任品牌副总裁并不是亲自想出DataSciencester这个名字的，因此他怀疑，换一个更好的名字可能会更有利于公司的成功。他请你使用数据科学来提出替换的候选名。

RNN 的一个“可爱”的应用包括使用字符（而不是单词）作为它们的输入，训练它们学习某个数据集中微妙的语言模式，然后使用它们生成该数据集的虚构实例。

例如，你可以训练一个 RNN 来学习另类乐队的名称，使用训练好的模型来生成新的假另类乐队的名称，然后手动选择最有趣的名称并分享在 Twitter 上。太有趣了！

见过这个技巧很多次后，你不再认为它很聪明，但你决定试试看。

经过一番调查，你发现创业加速器 Y Combinator 发布了其最成功的 100（实际上是 101）家初创企业的列表，这看起来是一个很好的起点。检查页面后，你发现公司名称都位于<b class="h4">标签内，这意味着你可以轻松使用你的网络爬虫技能来获取它们：

from bs4 import BeautifulSoup
import requests

url = "https://www.ycombinator.com/topcompanies/"
soup = BeautifulSoup(requests.get(url).text, 'html5lib')

# We get the companies twice, so use a set comprehension to deduplicate.
companies = list({b.text
                  for b in soup("b")
                  if "h4" in b.get("class", ())})
assert len(companies) == 101

和往常一样，页面可能会发生变化（或消失），这种情况下这段代码就不起作用了。如果是这样，你可以使用你新学到的数据科学技能来修复它，或者直接从书的 GitHub 站点获取列表。

那么我们的计划是什么呢？我们将训练一个模型来预测名称的下一个字符，给定当前字符和表示到目前为止所有字符的隐藏状态。

和往常一样，我们将预测字符的概率分布，并训练我们的模型以最小化SoftmaxCrossEntropy损失。

一旦我们的模型训练好了，我们可以使用它生成一些概率，根据这些概率随机抽取一个字符，然后将该字符作为下一个输入。这将允许我们使用学到的权重生成公司名称。

要开始，我们应该从名称中构建一个Vocabulary：

vocab = Vocabulary([c for company in companies for c in company])

此外，我们将使用特殊的标记来表示公司名称的开始和结束。这允许模型学习哪些字符应该开始一个公司名称，以及何时一个公司名称结束。

我们将只使用正则表达式字符来表示开始和结束，这些字符（幸运的是）不会出现在我们的公司名称列表中：

START = "^"
STOP = "$"

# We need to add them to the vocabulary too.
vocab.add(START)
vocab.add(STOP)

对于我们的模型，我们将对每个字符进行独热编码，通过两个SimpleRnn传递它们，然后使用Linear层生成每个可能的下一个字符的分数：

HIDDEN_DIM = 32  # You should experiment with different sizes!

rnn1 =  SimpleRnn(input_dim=vocab.size, hidden_dim=HIDDEN_DIM)
rnn2 =  SimpleRnn(input_dim=HIDDEN_DIM, hidden_dim=HIDDEN_DIM)
linear = Linear(input_dim=HIDDEN_DIM, output_dim=vocab.size)

model = Sequential([
    rnn1,
    rnn2,
    linear
])

假设我们已经训练好了这个模型。让我们编写一个函数，使用来自“主题建模”的sample_from函数来生成新的公司名称：

from scratch.deep_learning import softmax

def generate(seed: str = START, max_len: int = 50) -> str:
    rnn1.reset_hidden_state()  # Reset both hidden states
    rnn2.reset_hidden_state()
    output = [seed]            # Start the output with the specified seed

    # Keep going until we produce the STOP character or reach the max length
    while output[-1] != STOP and len(output) < max_len:
        # Use the last character as the input
        input = vocab.one_hot_encode(output[-1])

        # Generate scores using the model
        predicted = model.forward(input)

        # Convert them to probabilities and draw a random char_id
        probabilities = softmax(predicted)
        next_char_id = sample_from(probabilities)

        # Add the corresponding char to our output
        output.append(vocab.get_word(next_char_id))

    # Get rid of START and END characters and return the word
    return ''.join(output[1:-1])

终于，我们准备好训练我们的字符级 RNN。这会花费一些时间！

loss = SoftmaxCrossEntropy()
optimizer = Momentum(learning_rate=0.01, momentum=0.9)

for epoch in range(300):
    random.shuffle(companies)  # Train in a different order each epoch.
    epoch_loss = 0             # Track the loss.
    for company in tqdm.tqdm(companies):
        rnn1.reset_hidden_state()  # Reset both hidden states.
        rnn2.reset_hidden_state()
        company = START + company + STOP   # Add START and STOP characters.

        # The rest is just our usual training loop, except that the inputs
        # and target are the one-hot-encoded previous and next characters.
        for prev, next in zip(company, company[1:]):
            input = vocab.one_hot_encode(prev)
            target = vocab.one_hot_encode(next)
            predicted = model.forward(input)
            epoch_loss += loss.loss(predicted, target)
            gradient = loss.gradient(predicted, target)
            model.backward(gradient)
            optimizer.step(model)

    # Each epoch, print the loss and also generate a name.
    print(epoch, epoch_loss, generate())

    # Turn down the learning rate for the last 100 epochs.
    # There's no principled reason for this, but it seems to work.
    if epoch == 200:
        optimizer.lr *= 0.1

训练后，模型生成了一些实际的名称（这并不奇怪，因为模型具有相当的容量，但训练数据并不多），以及与训练名称略有不同的名称（Scripe, Loinbare, Pozium），看起来确实创意十足的名称（Benuus, Cletpo, Equite, Vivest），以及类似单词但是有点垃圾的名称（SFitreasy, Sint ocanelp, GliyOx, Doorboronelhav）。

不幸的是，像大多数字符级 RNN 输出一样，这些名称只是略有创意，品牌副总裁最终无法使用它们。

如果我将隐藏维度提升到 64，我将从列表中获得更多名称的原样；如果我将其降至 8，我将得到大多数垃圾名称。所有这些模型尺寸的词汇表和最终权重都可以在书的 GitHub 站点上找到，并且你可以使用load_weights和load_vocab来自己使用它们。

正如前面提到的，本章的 GitHub 代码还包含了 LSTM 的实现，你可以自由地将其替换为我们公司名称模型中的 SimpleRnn。

进一步探索

• NLTK 是一个流行的 Python 自然语言处理工具库。它有自己的整本 书籍，可以在线阅读。

• gensim 是一个用于主题建模的 Python 库，比我们从头开始的模型更可靠。

• spaCy 是一个用于“Python 中的工业级自然语言处理”库，也非常受欢迎。

• Andrej Karpathy 有一篇著名的博文，“递归神经网络的非理性有效性”，非常值得一读。

• 我的日常工作涉及构建 AllenNLP，一个用于进行自然语言处理研究的 Python 库。（至少在本书付印时是这样。）该库超出了本书的范围，但你可能会觉得它很有趣，而且还有一个很酷的交互式演示展示了许多最先进的 NLP 模型。