Clojure 数据科学（三）

最新推荐文章于 2025-07-20 12:39:36 发布

绝不原创的飞龙

最新推荐文章于 2025-07-20 12:39:36 发布

阅读量914

点赞数 10

CC 4.0 BY-SA版权

分类专栏：默认分类文章标签：默认分类

License CC BY-NC-SA 4.0 / 自豪地采用谷歌翻译

本文链接：https://blog.youkuaiyun.com/wizardforcel/article/details/149481297

默认分类专栏收录该内容

5217 篇文章

订阅专栏

原文：annas-archive.org/md5/bd3fb3fcef82cc10763b55da27a85129

译者：飞龙

协议：CC BY-NC-SA 4.0

第五章大数据

	“更多即不同。”
	–Philip Warren Anderson

在前几章中，我们使用回归技术将模型拟合到数据中。例如，在第三章，相关性中，我们构建了一个线性模型，利用普通最小二乘法和正态方程，通过运动员的身高和对数体重拟合了一条直线。在第四章，分类中，我们使用 Incanter 的优化命名空间，最小化逻辑代价函数，并构建了一个泰坦尼克号乘客的分类器。在本章中，我们将以适合更大数据量的方式应用类似的分析。

我们将使用一个相对简单的数据集，只有 100,000 条记录。这不是大数据（数据大小为 100 MB，足以在一台机器的内存中舒适地存储），但它足够大，可以展示大规模数据处理的常用技术。本章将以 Hadoop（分布式计算的流行框架）为案例，重点讨论如何通过并行化将算法扩展到非常大的数据量。我们将介绍 Clojure 提供的两个与 Hadoop 一起使用的库——Tesser和Parkour。

然而，在深入讨论 Hadoop 和分布式数据处理之前，我们将首先看看一些相同的原则，这些原则使得 Hadoop 在大规模上有效，并且这些原则也可以应用于单机的数据处理，借助现代计算机的并行能力。

下载代码和数据

本章使用的数据来自美国国税局（IRS）按邮政编码划分的个人收入数据。数据包含按州、邮政编码和收入等级分类的收入和税务项目。

文件大小为 100 MB，可以从www.irs.gov/pub/irs-soi/12zpallagi.csv下载到示例代码的数据目录中。由于该文件包含美国国税局的所得税统计（SoI），我们已将文件重命名为soi.csv以供示例使用。

注意

本章的示例代码可以从 Packt Publishing 的网站或github.com/clojuredatascience/ch5-big-data获取。

和往常一样，提供了一个脚本，可以为你下载并重命名数据。可以在项目目录内的命令行中运行：

script/download-data.sh

如果你运行此操作，文件将自动下载并重命名。

检查数据

下载数据后，查看文件第一行的列标题。访问文件第一行的一种方法是将文件加载到内存中，按换行符分割，并取第一个结果。Clojure 核心库中的函数slurp将返回整个文件作为字符串：

(defn ex-5-1 []
  (-> (slurp "data/soi.csv")
      (str/split #"\n")
      (first)))

该文件大约 100 MB 大小。当加载到内存并转换为对象表示时，数据在内存中将占用更多空间。当我们只对第一行感兴趣时，这样的做法尤其浪费。

幸运的是，如果我们利用 Clojure 的懒序列，我们不必将整个文件加载到内存中。我们可以返回文件的引用，然后逐行读取：

(defn ex-5-2 []
  (-> (io/reader "data/soi.csv")
      (line-seq)
      (first)))

在前面的代码中，我们使用 clojure.java.io/reader 返回文件的引用。同时，我们使用 clojure.core 函数 line-seq 返回文件的懒序列。通过这种方式，我们可以读取比可用内存更大的文件。

上一个函数的结果如下：

"STATEFIPS,STATE,zipcode,AGI_STUB,N1,MARS1,MARS2,MARS4,PREP,N2,NUMDEP,A00100,N00200,A00200,N00300,A00300,N00600,A00600,N00650,A00650,N00900,A00900,SCHF,N01000,A01000,N01400,A01400,N01700,A01700,N02300,A02300,N02500,A02500,N03300,A03300,N00101,A00101,N04470,A04470,N18425,A18425,N18450,A18450,N18500,A18500,N18300,A18300,N19300,A19300,N19700,A19700,N04800,A04800,N07100,A07100,N07220,A07220,N07180,A07180,N07260,A07260,N59660,A59660,N59720,A59720,N11070,A11070,N09600,A09600,N06500,A06500,N10300,A10300,N11901,A11901,N11902,A11902"

文件中有 77 个字段，因此我们不会全部列出。前四个字段是：

STATEFIPS：这是联邦信息处理系统（FIPS）代码。
STATE：这是州的两字母代码。
zipcode：这是 5 位数的邮政编码。
AGI_STUB：这是调整后总收入的一部分，按以下方式分箱：
1. $1 至 $25,000
2. $25,000 至 $50,000
3. $50,000 至 $75,000
4. $75,000 至 $100,000
5. $100,000 至 $200,000
6. $200,000 或更多

我们感兴趣的其他字段如下：

N1：提交的报税表数量
MARS2：提交的联合报税表数量
NUMDEP：受抚养人数量
N00200：包含薪水和工资的报税表数量
N02300：包含失业补偿的报税表数量

如果你感兴趣，可以在 www.irs.gov/pub/irs-soi/12zpdoc.doc 查阅完整的列定义列表。

计数记录

我们的文件当然很宽，但它高吗？我们想要确定文件中的总行数。创建了懒序列后，这只是计算序列长度的问题：

(defn ex-5-3 []
  (-> (io/reader "data/soi.csv")
      (line-seq)
      (count)))

上述示例返回 166,905 行，包括标题行，因此我们知道文件中实际上有 166,904 行。

count 函数是计算序列中元素数量的最简单方法。对于向量（以及实现计数接口的其他类型），这是最有效的方法，因为集合已经知道包含多少个元素，因此无需重新计算。然而，对于懒序列，唯一确定序列中包含多少元素的方法是从头到尾逐步遍历它。

Clojure 中的 count 实现是用 Java 编写的，但等效的 Clojure 版本是对序列进行 reduce 操作，如下所示：

(defn ex-5-4 []
  (->> (io/reader "data/soi.csv")
       (line-seq)
       (reduce (fn [i x]
                 (inc i)) 0)))

我们传递给reduce的前一个函数接受一个计数器i和来自序列的下一个元素x。对于每个x，我们简单地增加计数器i。reduce函数接受一个初始值零，代表“无”的概念。如果没有要合并的行，则会返回零。

从版本 1.5 开始，Clojure 提供了reducers库（clojure.org/reducers），它提供了一种通过牺牲内存效率来换取速度的替代方式来执行合并操作。

reducers 库

我们之前实现的count操作是一个顺序算法。每一行按顺序处理，直到序列耗尽。但是，这个操作并没有要求它必须按这种方式进行。

我们可以将行数分成两个序列（理想情况下长度大致相等），然后独立地对每个序列进行合并操作。当我们完成时，只需将每个序列的总行数相加，即可得到文件中的总行数：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_100.jpg

如果每个Reduce在自己的处理单元上运行，那么这两个计数操作将会并行执行。其他条件相同的情况下，算法将运行得更快，是原来的两倍。这就是clojure.core.reducers库的目标之一——通过利用多核处理器，让算法能够在单台机器上实现并行化。

使用 Reducers 的并行折叠

由 reducers 库实现的reduce的并行版本称为fold。为了使用fold，我们必须提供一个合并函数，它将接收我们减少过的序列结果（部分行计数）并返回最终结果。由于我们的行计数是数字，所以合并函数就是+。

注意

Reducers 是 Clojure 标准库的一部分，无需作为外部依赖添加。

调整后的示例，使用clojure.core.reducers作为r，如下所示：

(defn ex-5-5 []
  (->> (io/reader "data/soi.csv")
       (line-seq)
       (r/fold + (fn [i x]
                   (inc i)))))

合并函数+已作为第一个参数传递给fold，我们未更改的reduce函数则作为第二个参数传入。我们不再需要传递初始值零——fold会通过调用没有参数的合并函数来获取初始值。我们之前的示例之所以有效，是因为+在没有参数的情况下已经返回零：

(defn ex-5-6 []
  (+))

;; 0

因此，要参与折叠，合并函数必须有两个实现：一个没有参数，返回标识值，另一个有两个参数，用于合并这两个参数。当然，不同的折叠操作将需要不同的合并函数和标识值。例如，乘法的标识值是1。

我们可以将用标识值初始化计算的过程，迭代地在xs序列上执行合并，并将所有合并结果作为一棵树的输出值：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_110.jpg

当然，可能有超过两个的归约需要合并。fold 的默认实现会将输入集合拆分成 512 个元素的块。我们的 166,000 元素的序列因此会生成 325 次归约以供合并。由于树形表示图占用的空间较大，我们将很快用尽页面的可用空间，因此让我们改为以更简洁的方式来可视化这个过程——作为一个两步的归约和合并过程。

第一步在集合中的所有块上执行并行归约。第二步则对中间结果执行串行归约，以得到最终结果：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_120.jpg

上述表示展示了如何对多个 xs 序列进行归约，这里以圆圈表示，最终输出以方块表示。这些方块将串行合并，以产生最终结果，用星形表示。

使用 iota 加载大文件

在懒序列上调用 fold 需要 Clojure 将序列加载到内存中，然后将序列分块以进行并行执行。对于每行计算较小的情况，协调开销会超过并行化带来的好处。我们可以通过使用一个名为 iota 的库来稍微改善这种情况 (github.com/thebusby/iota)。

注意

iota 库通过使用内存映射文件将文件直接加载到适合使用 reducers 进行折叠的数据结构中，可以处理比可用内存更大的文件。

将 iota 代替我们的 line-seq 函数后，我们的行数统计变成了：

(defn ex-5-7 []
  (->> (iota/seq "data/soi.csv")
       (r/fold + (fn [i x]
                   (inc i)))))

到目前为止，我们只处理了未格式化的行序列，但如果我们要做的不仅仅是计算行数，我们希望将其解析为更有用的数据结构。这是 Clojure 的 reducers 可以帮助我们提高代码效率的另一个领域。

创建一个 reducers 处理管道

我们已经知道文件是以逗号分隔的，因此让我们首先创建一个函数，将每行转换为一个字段向量。除了前两个字段，其他字段都包含数值数据，因此我们在此过程中将其解析为 double 类型：

(defn parse-double [x]
  (Double/parseDouble x))

(defn parse-line [line]
  (let [[text-fields double-fields] (->> (str/split line #",")
                                         (split-at 2))]
    (concat text-fields
            (map parse-double double-fields))))

我们使用 map 的 reducers 版本，将我们的 parse-line 函数依次应用到文件中的每一行：

(defn ex-5-8 []
   (->> (iota/seq "data/soi.csv")
        (r/drop 1)
        (r/map parse-line)
        (r/take 1)
        (into [])))

;; [("01" "AL" 0.0 1.0 889920.0 490850.0 ...)]

最终的 into 函数调用将 reducers 的内部表示（一个可归约的集合）转换为一个 Clojure 向量。上面的示例应该返回一个包含 77 个字段的序列，表示文件中第一行（头部后）的内容。

目前我们只是去掉了列名，但如果我们能利用这些列名返回每条记录的映射表示，将列名与字段值关联起来，那就太好了。映射的键将是列名，值将是解析后的字段。clojure.core 函数 zipmap 将从两个序列（一个用于键，一个用于值）创建一个映射：

(defn parse-columns [line]
  (->> (str/split line #",")
       (map keyword)))

(defn ex-5-9 []
  (let [data (iota/seq "data/soi.csv")
        column-names (parse-columns (first data))]
    (->> (r/drop 1 data)
         (r/map parse-line)
         (r/map (fn [fields]
                  (zipmap column-names fields)))
         (r/take 1)
         (into []))))

此函数返回每行的映射表示，这是一种更加用户友好的数据结构：

[{:N2 1505430.0, :A19300 181519.0, :MARS4 256900.0 ...}]

Clojure 的 reducers 的一个很棒的功能是，在上述计算中， r/map, r/drop 和 r/take 的调用被组合成一个减少，将在数据上进行单次遍历。随着操作数量的增加，这变得尤为重要。

假设我们想要过滤掉零邮政编码。我们可以扩展 reducer 管道如下：

(defn ex-5-10 []
  (let [data (iota/seq "data/soi.csv")
        column-names (parse-columns (first data))]
    (->> (r/drop 1 data)
         (r/map parse-line)
         (r/map (fn [fields]
                  (zipmap column-names fields)))
         (r/remove (fn [record]
                     (zero? (:zipcode record))))
         (r/take 1)
         (into []))))

r/remove 步骤现在也与 r/map, r/drop 和 r/take 调用一起运行。随着数据量的增加，避免不必要地多次迭代数据变得越来越重要。使用 Clojure 的 reducers 确保我们的计算编译成单次遍历。

使用 reducers 的柯里化减少

为了使过程更加清晰，我们可以为之前的每个步骤创建一个 柯里化 版本。解析行，从字段创建记录并过滤掉零邮政编码。函数的柯里化版本是一个等待集合的减少：

(def line-formatter
  (r/map parse-line))

(defn record-formatter [column-names]
  (r/map (fn [fields]
           (zipmap column-names fields))))

(def remove-zero-zip
  (r/remove (fn [record]
              (zero? (:zipcode record)))))

在每种情况下，我们调用 reducers 的函数之一，但没有提供集合。响应是函数的柯里化版本，稍后可以应用于集合。这些柯里化函数可以使用 comp 组合成单个 parse-file 函数：

(defn load-data [file]
  (let [data (iota/seq file)
        col-names  (parse-columns (first data))
        parse-file (comp remove-zero-zip
                         (record-formatter col-names)
                         line-formatter)]
    (parse-file (rest data))))

只有当 parse-file 函数与序列一起调用时，管道才会实际执行。

使用 reducers 进行统计折叠

数据解析完成后，是时候进行一些描述性统计了。假设我们想知道 IRS 按邮政编码提交的平均退回数量（列 N1）。有一种方法可以做到这一点——这是本书中多次尝试过的方法——通过将值加起来并除以计数来实现。我们的第一次尝试可能如下所示：

(defn ex-5-11 []
  (let [data (load-data "data/soi.csv")
        xs (into [] (r/map :N1 data))]
    (/ (reduce + xs)
       (count xs))))

;; 853.37

虽然这样做可以工作，但速度相对较慢。我们需要对数据进行三次迭代：一次创建 xs，一次计算总和，一次计算数量。数据集越大，我们支付的时间成本越大。理想情况下，我们希望能够在单次数据遍历中计算均值，就像我们之前的 parse-file 函数一样。如果能并行执行更好。

结合性质

在继续之前，花点时间反思以下代码为什么不能达到我们想要的效果：

(defn mean
  ([] 0)
  ([x y] (/ (+ x y) 2)))

我们的 mean 函数是两个参数的函数。没有参数时，它返回零，是 mean 计算的单位元。有两个参数时，它返回它们的平均值：

(defn ex-5-12 []
  (->> (load-data "data/soi.csv")
       (r/map :N1)
       (r/fold mean)))

;; 930.54

前面的示例对 N1 数据使用我们的 mean 函数进行了折叠，并产生了与之前不同的结果。如果我们能够扩展出前三个 xs 的计算，我们可能会看到类似以下代码的内容：

(mean (mean (mean 0 a) b) c)

这是一个不好的主意，因为 mean 函数不是关联的。对于关联函数，以下条件成立：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_01.jpg

加法是结合律的，但乘法和除法不是。所以mean函数也不是结合的。将mean函数与以下简单加法进行对比：

(+ 1 (+ 2 3))

这与以下结果相同：

(+ (+ 1 2) 3)

+的参数如何分区并不重要。结合性是用于对数据集进行归约的函数的一个重要特性，因为根据定义，先前计算的结果会作为下一个计算的输入。

将mean函数转化为结合函数的最简单方法是分别计算和与计数。由于和与计数是结合的，它们可以并行地计算数据。mean函数可以通过将两者相除来简单计算。

使用 fold 计算均值

我们将使用两个自定义函数mean-combiner和mean-reducer来创建一个 fold。这需要定义三个实体：

fold 的身份值
reduce 函数
combine 函数

我们在前一节中发现了结合律的好处，因此我们希望只使用结合操作来更新中间的mean，并分别计算和与计数。表示这两个值的一种方式是一个包含：count和：sum两个键的映射。表示零的值就是和为零，计数为零，或者如下所示的映射：{:count 0 :sum 0}。

combine 函数mean-combiner在没有参数时会提供初始值。两个参数的 combiner 需要将每个参数的:count和:sum相加。我们可以通过使用+来合并这些映射：

(defn mean-combiner
  ([] {:count 0 :sum 0})
  ([a b] (merge-with + a b)))

mean-reducer函数需要接受一个累积值（无论是身份值还是先前计算结果），并将新的x纳入其中。我们只需通过递增:count并将x加到累积的:sum中来实现：

(defn mean-reducer [acc x]
  (-> acc
      (update-in [:count] inc)
      (update-in [:sum] + x)))

前述两个函数已经足够完全定义我们的mean fold：

(defn ex-5-13 []
  (->> (load-data "data/soi.csv")
       (r/map :N1)
       (r/fold mean-combiner
               mean-reducer)))

;; {:count 166598, :sum 1.4216975E8}

结果给出了我们计算N1均值所需的所有信息，而这个计算只需一次遍历数据。计算的最后一步可以通过以下mean-post-combiner函数来执行：

(defn mean-post-combiner [{:keys [count sum]}]
  (if (zero? count) 0 (/ sum count)))

(defn ex-5-14 []
  (->> (load-data "data/soi.csv")
       (r/map :N1)
       (r/fold mean-combiner
               mean-reducer)
       (mean-post-combiner)))

;; 853.37

幸运的是，结果与我们之前计算的均值一致。

使用 fold 计算方差

接下来，我们希望计算N1值的方差。记住，方差是均值的平方差：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_02.jpg

为了将其实现为 fold，我们可能会编写如下代码：

(defn ex-5-15 []
   (let [data (->> (load-data "data/soi.csv")
                   (r/map :N1))
         mean-x (->> data
                     (r/fold mean-combiner
                             mean-reducer)
                     (mean-post-combine))
         sq-diff (fn [x] (i/pow (- x mean-x) 2))]
     (->> data
          (r/map sq-diff)
          (r/fold mean-combiner
                  mean-reducer)
          (mean-post-combine))))

;; 3144836.86

首先，我们使用刚才构建的 fold 计算该系列的mean值。然后，我们定义一个关于x和sq-diff的函数，计算x与mean值之间的平方差。我们将其映射到平方差上，并再次调用mean fold，最终得到方差结果。

因此，我们对数据进行了两次完整的遍历，首先计算均值，然后计算每个 x 与 mean 值的差异。看起来计算方差必然是一个顺序算法：似乎无法进一步减少步骤，只能通过单次遍历数据来计算方差。

事实上，方差计算可以用一个折叠表达。为此，我们需要追踪三项内容：计数、（当前）均值和平方差之和：

(defn variance-combiner
  ([] {:count 0 :mean 0 :sum-of-squares 0})
  ([a b]
   (let [count (+ (:count a) (:count b))]
     {:count count
      :mean (/ (+ (* (:count a) (:mean a))
                  (* (:count b) (:mean b)))
               count)
      :sum-of-squares (+ (:sum-of-squares a)
                         (:sum-of-squares b)
                         (/ (* (- (:mean b)
                                  (:mean a))
                               (- (:mean b)
                                  (:mean a))
                               (:count a)
                               (:count b))
                            count))})))

我们的合并函数在前面的代码中展示。单位值是一个映射，所有三项的值都设为零。零元合并器返回这个值。

二元合并器需要合并两个传入值的计数、均值和平方和。合并计数很简单——我们只需要将它们相加。均值则稍微复杂一些：我们需要计算两个均值的加权平均。如果一个均值是基于较少的记录计算的，那么它在合并均值时应该占较少的份额：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_03.jpg

合并平方和是最复杂的计算。合并平方和时，我们还需要添加一个因子，以考虑到 a 和 b 的平方和可能是基于不同的均值计算得出的：

(defn variance-reducer [{:keys [count mean sum-of-squares]} x]
  (let [count' (inc count)
        mean'  (+ mean (/ (- x mean) count'))]
    {:count count'
     :mean mean'
     :sum-of-squares (+ sum-of-squares
                        (* (- x mean') (- x mean)))}))

缩减器要简单得多，并解释了方差折叠如何通过一次数据扫描计算方差。对于每个新记录，mean 值是根据之前的 mean 和当前的 count 重新计算的。然后，我们将平方和增加为均值前后差异的乘积，考虑到这条新记录。

最终结果是一个包含 count、mean 和总 sum-of-squares 的映射。由于方差只是 sum-of-squares 除以 count，因此我们的 variance-post-combiner 函数是相对简单的：

(defn variance-post-combiner [{:keys [count mean sum-of-squares]}]
   (if (zero? count) 0 (/ sum-of-squares count)))

将这三个函数结合起来，得到如下结果：

(defn ex-5-16 []
  (->> (load-data "data/soi.csv")
       (r/map :N1)
       (r/fold variance-combiner
               variance-reducer)
       (variance-post-combiner)))

;; 3144836.86

由于标准差只是方差的平方根，因此我们只需要稍微修改过的 variance-post-combiner 函数，也可以计算标准差。

使用 Tesser 进行数学折叠

我们现在应该理解如何使用折叠来计算简单算法的并行实现。希望我们也能理解找到高效解决方案所需的巧思，这样能够在数据上进行最少次数的迭代。

幸运的是，Clojure 库 Tesser (github.com/aphyr/tesser) 包含了常见数学折叠的实现，包括均值、标准差和协方差。为了了解如何使用 Tesser，我们考虑一下 IRS 数据集中两个字段的协方差：工资和薪金 A00200，以及失业补偿 A02300。

使用 Tesser 计算协方差

我们在第三章，相关性中遇到了协方差，它是用来衡量两组数据如何一起变化的指标。公式如下所示：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_04.jpg

tesser.math中包含了协方差 fold。在下面的代码中，我们将tesser.math作为m，tesser.core作为t来引用：

(defn ex-5-17 []
  (let [data (into [] (load-data "data/soi.csv"))]
    (->> (m/covariance :A02300 :A00200)
         (t/tesser (t/chunk 512 data )))))

;; 3.496E7

m/covariance函数期望接收两个参数：一个返回x值的函数，另一个返回y值的函数。由于关键字作为提取相应值的函数，我们只需将关键字作为参数传递。

Tesser 的工作方式类似于 Clojure 的 reducers，但有一些小差异。Clojure 的fold负责将我们的数据分割为子序列进行并行执行。而在 Tesser 中，我们必须显式地将数据分成小块。由于这是我们将要反复做的事情，让我们创建一个名为chunks的辅助函数：

(defn chunks [coll]
  (->> (into [] coll)
       (t/chunk 1024)))

在本章的大部分内容中，我们将使用chunks函数将输入数据拆分为1024条记录的小组。

交换律

Clojure 的 reducers 与 Tesser 的 folds 之间的另一个区别是，Tesser 不保证输入顺序会被保持。除了是结合律的，如我们之前讨论的，Tesser 的函数还必须满足交换律。交换律函数是指无论其参数的顺序如何变化，结果都相同的函数：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_05.jpg

加法和乘法是交换律的，但减法和除法不是。交换律是分布式数据处理函数的一个有用属性，因为它减少了子任务之间协调的需求。当 Tesser 执行合并函数时，它可以在任何先返回结果的 reducer 函数上进行，而不需要等到第一个完成。

让我们将load-data函数改写为prepare-data函数，它将返回一个交换律的 Tesser fold。它执行与我们之前基于 reducer 的函数相同的步骤（解析文本文件中的一行，格式化记录为 map 并删除零 ZIP 代码），但它不再假设列头是文件中的第一行——first是一个明确要求有序数据的概念：

(def column-names
  [:STATEFIPS :STATE :zipcode :AGI_STUB :N1 :MARS1 :MARS2 ...])

(defn prepare-data []
  (->> (t/remove #(.startsWith % "STATEFIPS"))
       (t/map parse-line)
       (t/map (partial format-record column-names))
       (t/remove  #(zero? (:zipcode %)))))

现在，所有的准备工作都已在 Tesser 中完成，我们可以直接将iota/seq的结果作为输入。这在本章后面我们将代码分布式运行在 Hadoop 时尤其有用：

(defn ex-5-18 []
  (let [data (iota/seq "data/soi.csv")]
    (->> (prepare-data)
         (m/covariance :A02300 :A00200)
         (t/tesser (chunks data)))))

;; 3.496E7

在第三章，相关性中，我们看到在简单线性回归中，当只有一个特征和一个响应变量时，相关系数是协方差与标准差乘积的比值：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_06.jpg

Tesser 还包括计算一对属性相关性的函数，作为一个 fold 来实现：

(defn ex-5-19 []
  (let [data (iota/seq "data/soi.csv")]
    (->> (prepare-data)
         (m/correlation :A02300 :A00200)
         (t/tesser (chunks data)))))

;; 0.353

这两个变量之间存在适度的正相关。让我们构建一个线性模型，用工资和薪金（A00200）来预测失业补偿（A02300）的值。

使用 Tesser 进行简单线性回归

Tesser 当前没有提供线性回归的 fold，但它提供了我们实现线性回归所需的工具。在第三章相关性中，我们看到如何通过输入的方差、协方差和均值来计算简单线性回归模型的系数——斜率和截距：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_07.jpg https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_08.jpg

斜率 b 是协方差除以 X 的方差。截距是确保回归线穿过两个序列均值的值。因此，理想情况下，我们能够在数据的单个 fold 中计算出这四个变量。Tesser 提供了两个 fold 组合器，t/fuse 和 t/facet，用于将更基本的 folds 组合成更复杂的 folds。

在我们有一个输入记录并且需要并行运行多个计算的情况下，我们应该使用 t/fuse。例如，在以下示例中，我们将均值和标准差的 fold 融合为一个单独的 fold，同时计算这两个值：

(defn ex-5-20 []
  (let [data (iota/seq "data/soi.csv")]
    (->> (prepare-data)
         (t/map :A00200)
         (t/fuse {:A00200-mean (m/mean)
                  :A00200-sd   (m/standard-deviation)})
         (t/tesser (chunks data)))))

;; {:A00200-sd 89965.99846545042, :A00200-mean 37290.58880658831}

在这里，我们需要对映射中的所有字段进行相同的计算；因此，我们应该使用 t/facet：

(defn ex-5-21 []
  (let [data (iota/seq "data/soi.csv")]
    (->> (prepare-data)
         (t/map #(select-keys % [:A00200 :A02300]))
         (t/facet)
         (m/mean)
         (t/tesser (chunks data)))))

;; {:A02300 419.67862159209596, :A00200 37290.58880658831}

在前面的代码中，我们只选择了记录中的两个值（A00200 和 A02300），并同时计算了它们的 mean 值。回到执行简单线性回归的挑战——我们有四个数值需要计算，因此我们将它们 fuse（结合）起来：

(defn calculate-coefficients [{:keys [covariance variance-x
                                      mean-x mean-y]}]
  (let [slope     (/ covariance variance-x)
        intercept (- mean-y (* mean-x slope))]
    [intercept slope]))

(defn ex-5-22 []
  (let [data (iota/seq "data/soi.csv")
        fx :A00200
        fy :A02300]
    (->> (prepare-data)
         (t/fuse {:covariance (m/covariance fx fy)
                  :variance-x (m/variance (t/map fx))
                  :mean-x (m/mean (t/map fx))
                  :mean-y (m/mean (t/map fx))})
         (t/post-combine calculate-coefficients)
         (t/tesser (chunks data)))))

;; [37129.529236553506 0.0043190406799462925]

fuse 非常简洁地将我们想要执行的计算结合在一起。此外，它允许我们指定一个 post-combine 步骤，作为 fuse 的一部分。我们可以直接将其作为 fold 的一个整体来指定，而不是将结果交给另一个函数来完成输出。post-combine 步骤接收四个结果，并从中计算出斜率和截距，返回两个系数作为一个向量。

计算相关性矩阵

我们只比较了两个特征，看看它们之间的相关性，但 Tesser 使得查看多个目标特征之间的相互关联变得非常简单。我们将目标特征提供为一个映射，其中特征名称对应于输入记录中某个函数，这个函数返回所需的特征。例如，在第三章相关性中，我们本来会取身高的对数。在这里，我们将简单地提取每个特征，并为它们提供易于理解的名称：

(defn ex-5-23 []
  (let [data (iota/seq "data/soi.csv")
        attributes {:unemployment-compensation :A02300
                    :salary-amount             :A00200
                    :gross-income              :AGI_STUB
                    :joint-submissions         :MARS2
                    :dependents                :NUMDEP}]
    (->> (prepare-data)
         (m/correlation-matrix attributes)
         (t/tesser (chunks data)))))

Tesser 将计算每对特征之间的相关性，并将结果以映射的形式返回。映射以包含每对特征名称的元组（两个元素的向量）为键，相关值为它们之间的相关性。

如果你现在运行前面的示例，你会发现某些变量之间有较高的相关性。例如，:dependents和:unemployment-compensation之间的相关性为0.821。让我们建立一个线性回归模型，将所有这些变量作为输入。

梯度下降的多元回归

当我们在第三章，相关性中运行多元线性回归时，我们使用了正规方程和矩阵来快速得出多元线性回归模型的系数。正规方程如下重复：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_09.jpg

正规方程使用矩阵代数来非常快速且高效地得出最小二乘估计。当所有数据都能装入内存时，这是一个非常方便且简洁的方程。然而，当数据超出单台机器可用的内存时，计算变得笨重。原因在于矩阵求逆。计算https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_10.jpg并不是可以在数据上折叠一次完成的——输出矩阵中的每个单元格都依赖于输入矩阵中的许多其他单元格。这些复杂的关系要求矩阵必须以非顺序的方式处理。

解决线性回归问题以及许多其他相关机器学习问题的另一种方法是梯度下降技术。梯度下降将问题重新构造为一个迭代算法的解决方案——这个算法并不是在一个计算量极大的步骤中计算答案，而是通过一系列更小的步骤逐渐接近正确答案。

在上一章中，我们遇到了梯度下降，当时我们使用了 Incanter 的minimize函数来计算出使我们逻辑回归分类器成本最低的参数。随着数据量的增加，Incanter 不再是执行梯度下降的可行方案。在下一节中，我们将看到如何使用 Tesser 自行运行梯度下降。

梯度下降更新规则

梯度下降通过反复应用一个函数来工作，这个函数将参数朝着其最优值的方向移动。为了应用这个函数，我们需要知道当前参数下的成本函数的梯度。

计算梯度公式涉及微积分内容，超出了本书的范围。幸运的是，最终的公式并不难理解：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_11.jpg

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_12.jpg是我们代价函数J(β)相对于索引j的参数的偏导数或梯度。因此，我们可以看到，代价函数相对于索引j的参数的梯度等于我们预测值与真实值y之间的差乘以索引j的特征值x。

因为我们希望沿着梯度下降，所以我们希望从当前的参数值中减去梯度的某个比例。因此，在每一步梯度下降中，我们会执行以下更新：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_13.jpg

这里，:=是赋值操作符，α是一个叫做学习率的因子。学习率控制我们希望在每次迭代中根据梯度对参数进行调整的大小。如果我们的预测ŷ几乎与实际值y相符，那么对参数的调整就不大。相反，较大的误差会导致对参数进行更大的调整。这个规则被称为Widrow-Hoff 学习规则或Delta 规则。

梯度下降学习率

正如我们所见，梯度下降是一个迭代算法。学习率，通常用α表示，决定了梯度下降收敛到最终答案的速度。如果学习率太小，收敛过程会非常缓慢。如果学习率过大，梯度下降就无法找到接近最优值的结果，甚至可能会从正确答案发散出去：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_130.jpg

在前面的图表中，一个较小的学习率会导致算法在多次迭代中收敛得非常慢。虽然算法最终能够达到最小值，但它需要经过比理想情况更多的步骤，因此可能会花费相当长的时间。相反，在下图中，我们可以看到一个过大学习率的效果。参数估计在每次迭代之间变化如此剧烈，以至于它们实际上超出了最优值，并且从最小值开始发散：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_140.jpg

梯度下降算法要求我们在数据集上反复迭代。在选择正确的α值后，每次迭代应该会逐步提供更接近理想参数的估计值。当迭代之间的变化非常小，或者达到了预定的迭代次数时，我们可以选择终止算法。

特征缩放

随着更多特征被添加到线性模型中，适当缩放特征变得非常重要。如果特征的尺度差异非常大，梯度下降的表现会很差，因为无法为所有特征选择一个合适的学习率。

我们可以执行的简单缩放是从每个值中减去均值并除以标准差。这将使得值趋于零均值，通常在-3和3之间波动：

(defn feature-scales [features]
  (->> (prepare-data)
       (t/map #(select-keys % features))
       (t/facet)
       (t/fuse {:mean (m/mean)
                :sd   (m/standard-deviation)})))

上述代码中的feature-factors函数使用t/facet来计算所有输入特征的mean值和标准差：

(defn ex-5-24 []
  (let [data (iota/seq "data/soi.csv")
        features [:A02300 :A00200 :AGI_STUB :NUMDEP :MARS2]]
    (->> (feature-scales features)
         (t/tesser (chunks data)))))

;; {:MARS2 {:sd 533.4496892658647, :mean 317.0412009748016}...}

如果你运行上述示例，你会看到feature-scales函数返回的不同均值和标准差。由于我们的特征缩放和输入记录表示为 maps，我们可以使用 Clojure 的merge-with函数一次性对所有特征进行缩放：

(defn scale-features [factors]
  (let [f (fn [x {:keys [mean sd]}]
            (/ (- x mean) sd))]
    (fn [x]
      (merge-with f x factors))))

同样，我们可以使用unscale-features进行至关重要的反转操作：

(defn unscale-features [factors]
  (let [f (fn [x {:keys [mean sd]}]
            (+ (* x sd) mean))]
    (fn [x]
      (merge-with f x factors))))

让我们对特征进行缩放并查看第一个特征。Tesser 不允许我们在没有 reducer 的情况下执行折叠操作，因此我们将暂时恢复使用 Clojure 的 reducers：

(defn ex-5-25 []
  (let [data     (iota/seq "data/soi.csv")
        features [:A02300 :A00200 :AGI_STUB :NUMDEP :MARS2]
        factors (->> (feature-scales features)
                     (t/tesser (chunks data)))]
    (->> (load-data "data/soi.csv")
         (r/map #(select-keys % features ))
         (r/map (scale-features factors))
         (into [])
         (first))))

;; {:MARS2 -0.14837567114357617, :NUMDEP 0.30617757526890155,
;;  :AGI_STUB -0.714280814223704, :A00200 -0.5894942801950217,
;;  :A02300 0.031741856083514465}

这一步骤将帮助梯度下降在我们的数据上进行优化。

特征提取

尽管我们在本章中使用 maps 来表示输入数据，但在运行梯度下降时，将特征表示为矩阵会更加方便。让我们编写一个函数，将输入数据转换为一个包含xs和y的 map。y轴将是标量响应值，xs将是缩放后的特征值矩阵。

如前几章所示，我们将偏置项添加到返回的特征矩阵中：

(defn feature-matrix [record features]
  (let [xs (map #(% record) features)]
    (i/matrix (cons 1 xs))))

(defn extract-features [fy features]
  (fn [record]
    {:y  (fy record)
     :xs (feature-matrix record features)}))

我们的feature-matrix函数仅接受一个记录输入和要转换为矩阵的特征。我们在extract-features中调用此函数，该函数返回一个可以应用于每个输入记录的函数：

(defn ex-5-26 []
  (let [data     (iota/seq "data/soi.csv")
        features [:A02300 :A00200 :AGI_STUB :NUMDEP :MARS2]
        factors (->> (feature-scales features)
                     (t/tesser (chunks data)))]
    (->> (load-data "data/soi.csv")
         (r/map (scale-features factors))
         (r/map (extract-features :A02300 features))
         (into [])
         (first))))

;; {:y 433.0, :xs  A 5x1 matrix
;;  -------------
;;  1.00e+00
;; -5.89e-01
;; -7.14e-01
;;  3.06e-01
;; -1.48e-01
;; }

上述示例展示了数据被转换为适合进行梯度下降的格式：一个包含y响应变量的 map 和包含偏置项的值矩阵。

创建自定义 Tesser 折叠

梯度下降的每次迭代都会根据成本函数确定的值调整系数。成本函数是通过对数据集中每个参数的误差求和计算得出的，因此，拥有一个对矩阵元素进行逐项求和的折叠函数将非常有用。

而 Clojure 通过 reducer、combiner 以及从 combiner 获得的身份值来表示折叠，Tesser 的折叠则通过六个协作函数来表达。Tesser 的m/mean折叠的实现如下：

{:reducer-identity  (constantly [0 0])
 :reducer           (fn reducer [[s c] x]
                     [(+ s x) (inc c)])
 :post-reducer      identity
 :combiner-identity (constantly [0 0])
 :combiner          (fn combiner [x y] (map + x y))
 :post-combiner     (fn post-combiner [x]
                      (double (/ (first x)
                                 (max 1 (last x)))))}

Tesser 选择将reducer身份与combiner函数分开表示，并且还包含另外三个函数：combiner-identity、post-reducer和post-combiner函数。Tesser 的mean折叠将一对数字（计数和累积和）表示为两个数字的向量，但在其他方面，它与我们自己的折叠类似。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_150.jpg

我们已经看到如何使用post-combiner函数，与我们在本章前面提到的mean-post-combiner和variance-post-combiner函数配合使用。

创建一个矩阵求和折叠

要创建一个自定义的matrix-sum折叠，我们需要一个单位值。我们在第三章《相关性》中遇到过单位矩阵，但这是矩阵乘法的单位矩阵，而不是加法的单位矩阵。如果+的单位值是零（因为将零加到一个数字上不会改变它），那么矩阵加法的单位矩阵就是一个全零矩阵。

我们必须确保矩阵的大小与我们想要相加的矩阵相同。所以，让我们使用矩阵的行数和列数来参数化我们的matrix-sum折叠。我们无法提前知道需要多大，因为单位函数会在折叠中的任何操作之前被调用：

(defn matrix-sum [nrows ncols]
  (let [zeros-matrix (i/matrix 0 nrows ncols)]
    {:reducer-identity (constantly zeros-matrix)
     :reducer i/plus
     :combiner-identity (constantly zeros-matrix)
     :combiner i/plus}))

上面的示例是完整的matrix-sum折叠定义。我们没有提供post-combiner和post-reducer函数；因为如果省略这些，默认它们是单位函数，这正是我们想要的。我们可以使用新的折叠来计算输入中所有特征的总和，如下所示：

(defn ex-5-27 []
   (let [columns [:A02300 :A00200 :AGI_STUB :NUMDEP :MARS2]
         data    (iota/seq "data/soi.csv")]
     (->> (prepare-data)
          (t/map (extract-features :A02300 columns))
          (t/map :xs)
          (t/fold (matrix-sum (inc (count columns)) 1))
          (t/tesser (chunks data)))))

;; A 6x1 matrix
;; -------------
;; 1.67e+05
;; 6.99e+07
;; 6.21e+09
;; ...
;; 5.83e+05
;; 9.69e+07
;; 5.28e+07

计算矩阵的总和让我们更接近执行梯度下降的目标。让我们使用新的折叠计算总模型误差，前提是我们有一些初始系数。

计算总模型误差

让我们再看一看梯度下降的 delta 规则：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_14.jpg

对于每个参数* j ，我们根据总体预测误差 ŷ - y *的某个比例来调整该参数，并乘以特征。因此，较大的特征比较小的特征承担更多的成本，并相应地调整更大的量。为了在代码中实现这一点，我们需要计算：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_15.jpg

这是所有输入记录中特征与预测误差乘积的总和。正如我们之前所做的那样，我们的预测值* y 将使用以下公式为每个输入记录 x *计算：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_16.jpg

系数β在所有输入记录中都是相同的，因此让我们创建一个calculate-error函数。给定转置的系数β^T，我们返回一个函数来计算！计算总模型误差。由于* x 是一个矩阵， ŷ - y *是一个标量，结果将是一个矩阵：

(defn calculate-error [coefs-t]
  (fn [{:keys [y xs]}]
    (let [y-hat (first (i/mmult coefs-t xs))
          error (- y-hat y)]
      (i/mult xs error))))

为了计算整个数据集的误差总和，我们可以简单地在calculate-error步骤后连接我们之前定义的matrix-sum函数：

(defn ex-5-28 []
  (let [columns [:A02300 :A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount  (inc (count columns))
        coefs   (vec (replicate fcount 0))
        data    (iota/seq "data/soi.csv")]
    (->> (prepare-data)
         (t/map (extract-features :A02300 columns))
         (t/map (calculate-error (i/trans coefs)))
         (t/fold (matrix-sum fcount 1))
         (t/tesser (chunks data)))))

;; A 6x1 matrix
;;  -------------
;; -6.99e+07
;; -2.30e+11
;; -8.43e+12
;;  ...
;; -1.59e+08
;; -2.37e+11
;; -8.10e+10

请注意，所有特征的梯度都是负的。这意味着，为了下降梯度并生成更好的模型系数估计，必须增加参数。

创建一个矩阵均值折叠

前面的代码中定义的更新规则实际上要求将代价的平均值分配给每一个特征。这意味着我们需要计算sum和count。我们不想对数据执行两次单独的遍历。因此，正如我们之前做的那样，我们将这两个折叠合并成一个：

(defn ex-5-29 []
  (let [columns [:A02300 :A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount  (inc (count columns))
        coefs   (vec (replicate fcount 0))
        data    (iota/seq "data/soi.csv")]
    (->> (prepare-data)
         (t/map (extract-features :A02300 columns))
         (t/map (calculate-error (i/trans coefs)))
         (t/fuse {:sum   (t/fold (matrix-sum fcount 1))
                  :count (t/count)})
         (t/post-combine (fn [{:keys [sum count]}]
                           (i/div sum count)))
         (t/tesser (chunks data)))))

fuse函数将返回：sum和：count的映射，因此我们将在结果上调用post-combine。post-combine函数指定了一个在折叠结束时运行的函数，该函数简单地将总和除以计数。

另外，我们可以创建另一个自定义折叠，返回一个矩阵序列的均值，而不是总和。它与之前定义的matrix-sum折叠有很多相似之处，但与我们在本章前面计算的mean折叠一样，我们还需要跟踪处理过的记录数：

(defn matrix-mean [nrows ncols]
  (let [zeros-matrix (i/matrix 0 nrows ncols)]
    {:reducer-identity  (constantly [zeros-matrix 0])
     :reducer           (fn [[sum counter] x]
                          [(i/plus sum x) (inc counter)])
     :combiner-identity (constantly [zeros-matrix 0])
     :combiner          (fn [[sum-a count-a] [sum-b count-b]]
                          [(i/plus sum-a sum-b)
                           (+ count-a count-b)])
     :post-combiner     (fn [[sum count]]
                          (i/div sum count))}))

减少器的恒等式是一个包含[zeros-matrix 0]的向量。每次减少都会将值添加到矩阵总和中，并将计数器加一。每个合并步骤都会将两个矩阵以及它们的计数相加，以得出所有分区的总和和总计数。最后，在post-combiner步骤中，均值被计算为sum与count的比值。

尽管自定义折叠的代码比我们合并的sum和count解决方案要更长，但我们现在有了一种计算矩阵均值的一般方法。这使得示例更加简洁易读，而且我们可以像这样在误差计算代码中使用它：

(defn ex-5-30 []
  (let [features [:A02300 :A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount   (inc (count features))
        coefs    (vec (replicate fcount 0))
        data     (iota/seq "data/soi.csv")]
    (->> (prepare-data)
         (t/map (extract-features :A02300 features))
         (t/map (calculate-error (i/trans coefs)))
         (t/fold (matrix-mean fcount 1))
         (t/tesser (chunks data)))))

;;  A 5x1 matrix
;;  -------------
;;  4.20e+01
;;  3.89e+01
;;  -3.02e+01
;;  9.02e+01
;;  6.62e+01

创建自定义折叠的小额额外努力，使得调用代码的意图变得更容易理解。

应用梯度下降的单步

计算代价的目标是确定调整每个系数的幅度。一旦我们计算出平均代价，正如我们之前所做的那样，我们需要更新对系数β的估计。这些步骤合起来代表梯度下降的单次迭代：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_14.jpg

我们可以在post-combiner步骤中返回更新后的系数，该步骤利用平均代价、alpha 的值和之前的系数。我们来创建一个实用函数update-coefficients，它将接收系数和 alpha，并返回一个函数，该函数将在给定总模型代价的情况下计算新的系数：

(defn update-coefficients [coefs alpha]
  (fn [cost]
    (->> (i/mult cost alpha)
         (i/minus coefs))))

在前面的函数到位后，我们就有了将批量梯度下降更新规则打包的所有必要工具：

(defn gradient-descent-fold [{:keys [fy features factors
                                     coefs alpha]}]
  (let [zeros-matrix (i/matrix 0 (count features) 1)]
    (->> (prepare-data)
         (t/map (scale-features factors))
         (t/map (extract-features fy features))
         (t/map (calculate-error (i/trans coefs)))
         (t/fold (matrix-mean (inc (count features)) 1))
         (t/post-combine (update-coefficients coefs alpha)))))

(defn ex-5-31 []
  (let [features [:A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount   (inc (count features))
        coefs    (vec (replicate fcount 0))
        data     (chunks (iota/seq "data/soi.csv"))
        factors  (->> (feature-scales features)
                      (t/tesser data))
        options {:fy :A02300 :features features
                 :factors factors :coefs coefs :alpha 0.1}]
    (->> (gradient-descent-fold options)
         (t/tesser data))))

;; A 6x1 matrix
;; -------------
;; -4.20e+02
;; -1.38e+06
;; -5.06e+07
;; -9.53e+02
;; -1.42e+06
;; -4.86e+05

结果矩阵表示梯度下降第一次迭代后的系数值。

运行迭代梯度下降

梯度下降是一个迭代算法，通常需要运行多次才能收敛。对于一个大型数据集，这可能非常耗时。

为了节省时间，我们在数据目录中包含了一个名为soi-sample.csv的soi.csv随机样本。较小的文件大小使得我们能够在合理的时间尺度内运行迭代梯度下降。以下代码执行 100 次梯度下降迭代，并绘制每次迭代之间的参数值在xy-plot上的变化：

(defn descend [options data]
  (fn [coefs]
    (->> (gradient-descent-fold (assoc options :coefs coefs))
         (t/tesser data))))

(defn ex-5-32 []
  (let [features [:A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount   (inc (count features))
        coefs    (vec (replicate fcount 0))
        data     (chunks (iota/seq "data/soi-sample.csv"))
        factors  (->> (feature-scales features)
                      (t/tesser data))
        options  {:fy :A02300 :features features
                  :factors factors :coefs coefs :alpha 0.1}
        iterations 100
        xs (range iterations)
        ys (->> (iterate (descend options data) coefs)
                (take iterations))]
    (-> (c/xy-plot xs (map first ys)
                   :x-label "Iterations"
                   :y-label "Coefficient")
        (c/add-lines xs (map second ys))
        (c/add-lines xs (map #(nth % 2) ys))
        (c/add-lines xs (map #(nth % 3) ys))
        (c/add-lines xs (map #(nth % 4) ys))
        (i/view))))

如果你运行这个示例，你应该会看到一个类似于下图的图表：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_160.jpg

在上面的图表中，你可以看到在 100 次迭代过程中，参数是如何趋向相对稳定的值的。

使用 Hadoop 扩展梯度下降

批量梯度下降每次迭代运行所需的时间由数据大小和计算机的处理器数量决定。尽管多个数据块是并行处理的，但数据集很大，且处理器是有限的。通过并行计算，我们已经实现了速度提升，但如果我们将数据集的大小翻倍，运行时间也会翻倍。

Hadoop 是过去十年中出现的多个系统之一，旨在并行化超出单台机器能力的工作。Hadoop 并不是将代码分发到多个处理器上执行，而是负责在多个服务器上运行计算。实际上，Hadoop 集群可以，也确实有很多，包含成千上万的服务器。

Hadoop 由两个主要子系统组成——Hadoop 分布式文件系统（HDFS）——和作业处理系统MapReduce。HDFS 将文件存储为块。一个文件可能由许多块组成，且这些块通常会在多个服务器之间进行复制。通过这种方式，Hadoop 能够存储任何单一服务器无法处理的庞大数据量，并且通过复制确保数据在硬件故障时能够可靠存储。正如其名所示，MapReduce 编程模型围绕 map 和 reduce 步骤构建。每个作业至少包含一个 map 步骤，并可以选择性地指定一个 reduce 步骤。一个完整的作业可能由多个 map 和 reduce 步骤串联而成。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_170.jpg

在 reduce 步骤是可选的这一点上，Hadoop 相比 Tesser 在分布式计算方面有一个稍微更灵活的方法。在本章以及未来章节中，我们将进一步探讨 Hadoop 提供的更多功能。Tesser 确实使我们能够将折叠转换为 Hadoop 作业，接下来我们就来做这个。

使用 Tesser 和 Parkour 在 Hadoop 上运行梯度下降

Tesser 的 Hadoop 功能可以在tesser.hadoop命名空间中找到，我们将其包含为h。Hadoop 命名空间中的主要公共 API 函数是h/fold。

fold 函数期望接收至少四个参数，分别代表 Hadoop 作业的配置、我们想要处理的输入文件、Hadoop 用于存储临时文件的工作目录，以及我们希望执行的 fold，它被引用为 Clojure var。任何额外的参数将在执行时作为参数传递给 fold。

使用 var 来表示我们的 fold 的原因是，启动 fold 的函数调用可能发生在与实际执行它的计算机完全不同的计算机上。在分布式环境中，var 和参数必须完全指定函数的行为。我们通常不能依赖其他可变的局部状态（例如，一个原子变量的值，或闭包内的变量值）来提供任何额外的上下文。

Parkour 分布式源和接收器

我们希望 Hadoop 作业处理的数据可能也存在于多台机器上，并且在 HDFS 上分块分布。Tesser 使用一个名为 Parkour 的库（github.com/damballa/parkour/）来处理访问可能分布式的数据源。我们将在本章和下一章更详细地学习 Parkour，但目前，我们只会使用 parkour.io.text 命名空间来引用输入和输出文本文件。

虽然 Hadoop 设计为在多个服务器上运行和分布式执行，但它也可以在 本地模式 下运行。本地模式适用于测试，并使我们能够像操作 HDFS 一样与本地文件系统交互。我们将从 Parkour 中使用的另一个命名空间是 parkour.conf。这将允许我们创建一个默认的 Hadoop 配置，并在本地模式下运行它：

(defn ex-5-33 []
  (->> (text/dseq "data/soi.csv")
       (r/take 2)
       (into [])))

在前面的示例中，我们使用 Parkour 的 text/dseq 函数创建了一个 IRS 输入数据的表示。返回值实现了 Clojure 的 reducers 协议，因此我们可以对结果使用 r/take。

使用 Hadoop 运行功能规模的 `fold`

Hadoop 在执行任务时需要一个位置来写入其临时文件，如果我们尝试覆盖现有目录，它会发出警告。由于在接下来的几个示例中我们将执行多个作业，让我们创建一个小的工具函数，返回一个带有随机生成名称的新文件。

(defn rand-file [path]
  (io/file path (str (long (rand 0x100000000)))))

(defn ex-5-34 []
  (let [conf     (conf/ig)
        input    (text/dseq "data/soi.csv")
        workdir  (rand-file "tmp")
        features [:A00200 :AGI_STUB :NUMDEP :MARS2]]
    (h/fold conf input workdir #'feature-scales features)))

Parkour 提供了一个默认的 Hadoop 配置对象，使用简写（conf/ig）。它将返回一个空的配置。默认值已经足够，我们不需要提供任何自定义配置。

注意

我们的所有 Hadoop 作业将把临时文件写入项目 tmp 目录中的一个随机目录。如果你担心磁盘空间，可以稍后删除这个文件夹。

如果你现在运行前面的示例，你应该会得到类似以下的输出：

;; {:MARS2 317.0412009748016, :NUMDEP 581.8504423822615,
;; :AGI_STUB 3.499939975269811, :A00200 37290.58880658831}

尽管返回值与我们之前获得的值相同，但我们现在在幕后使用 Hadoop 来处理我们的数据。尽管如此，注意 Tesser 将从我们的折叠操作返回一个单一的 Clojure 数据结构作为响应。

使用 Hadoop 运行梯度下降

由于tesser.hadoop折叠返回的 Clojure 数据结构与tesser.core折叠类似，定义一个利用我们缩放后的特征的梯度下降函数非常简单：

(defn hadoop-gradient-descent [conf input-file workdir]
  (let [features [:A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount  (inc (count features))
        coefs   (vec (replicate fcount 0))
        input   (text/dseq input-file)
        options {:column-names column-names
                 :features features
                 :coefs coefs
                 :fy :A02300
                 :alpha 1e-3}
        factors (h/fold conf input (rand-file workdir)
                        #'feature-scales
                        features)
        descend (fn [coefs]
                  (h/fold conf input (rand-file workdir)
                          #'gradient-descent-fold
                          (merge options {:coefs coefs
                                          :factors factors})))]
    (take 5 (iterate descend coefs))))

上述代码定义了一个hadoop-gradient-descent函数，该函数迭代执行descend函数5次。每次descend迭代都会基于gradient-descent-fold函数计算改进后的系数。最终返回值是梯度下降5次后的系数向量。

我们在以下示例中运行整个 IRS 数据的作业：

(defn ex-5-35 []
  (let [workdir  "tmp"
        out-file (rand-file workdir)]
    (hadoop-gradient-descent (conf/ig) "data/soi.csv" workdir)))

在经过几次迭代后，你应该看到类似以下的输出：

;; ([0 0 0 0 0]
;; (20.9839310796048 46.87214911003046 -7.363493937722712
;;  101.46736841329326 55.67860863427868)
;; (40.918665605227744 56.55169901254631 -13.771345753228694
;;  162.1908841131747 81.23969785586247)
;; (59.85666340457121 50.559130068258995 -19.463888245285332
;;  202.32407094149158 92.77424653758085)
;; (77.8477613139478 38.67088624825574 -24.585818946408523
;;  231.42399118694212 97.75201693843269))

我们已经看到如何使用本地的分布式技术计算梯度下降。现在，让我们看看如何在我们自己的集群上运行这个过程。

为 Hadoop 集群准备我们的代码

Hadoop 的 Java API 定义了Tool和相关的ToolRunner类，用于帮助在 Hadoop 集群上执行任务。Tool类是 Hadoop 为通用命令行应用程序定义的名称，它与 Hadoop 框架进行交互。通过创建我们自己的工具，我们就创建了一个可以提交到 Hadoop 集群的命令行应用程序。

由于这是一个 Java 框架，Hadoop 期望与我们代码的类表示进行交互。因此，定义我们工具的命名空间需要包含:gen-class声明，指示 Clojure 编译器从我们的命名空间创建一个类：

(ns cljds.ch5.hadoop
  (:gen-class)
  ...)

默认情况下，:gen-class将期望命名空间定义一个名为-main的主函数。这将是 Hadoop 用我们的参数调用的函数，因此我们可以简单地将调用委托给一个实际执行我们任务的函数：

(defn -main [& args]
  (tool/run hadoop-gradient-descent args))

Parkour 提供了一个 Clojure 接口，用于与 Hadoop 的许多类进行交互。在这种情况下，parkour.tool/run包含了我们在 Hadoop 上运行分布式梯度下降函数所需的一切。通过上述示例，我们需要指示 Clojure 编译器提前（AOT）编译我们的命名空间，并指定我们希望将项目的主类定义为哪个类。我们可以通过将:aot和:main声明添加到project.clj函数中来实现：

{:main cljds.ch5.hadoop
 :aot [cljds.ch5.hadoop]}

在示例代码中，我们将这些设置为：uberjar配置的一部分，因为在将作业发送到集群之前，我们的最后一步是将其打包为 uberjar 文件。

构建 uberjar

JAR 文件包含可执行的 Java 代码。一个 uberjar 文件包含可执行的 Java 代码，以及运行所需的所有依赖项。uberjar 提供了一种方便的方式来打包代码，以便在分布式环境中运行，因为作业可以在机器之间传递，同时携带它的依赖项。虽然这会导致较大的作业负载，但它避免了确保所有集群中的机器预先安装特定作业依赖项的需要。要使用Leiningen创建 uberjar 文件，在项目目录下执行以下命令行：

lein uberjar

完成此操作后，目标目录中将创建两个文件。一个文件ch5-0.1.0.jar包含项目的编译代码。这与使用lein jar生成的文件相同。此外，uberjar 会生成ch5-0.1.0-standalone.jar文件。这个文件包含了项目代码的 AOT 编译版本以及项目的依赖项。生成的文件虽然较大，但它包含了 Hadoop 作业运行所需的所有内容。

将 uberjar 提交到 Hadoop

一旦我们创建了一个 uberjar 文件，就可以将其提交给 Hadoop。拥有一个本地的 Hadoop 安装并不是跟随本章示例的前提条件，我们也不会在此描述安装步骤。

注意

Hadoop 安装指南的链接已提供在本书的 wiki 页面：wiki.clojuredatascience.com。

然而，如果你已经在本地模式下安装并配置了 Hadoop，现在就可以在命令行上运行示例作业。由于主类指定的工具也接受两个参数——工作目录和输入文件——因此这些参数也需要提供：

hadoop jar target/ch5-0.1.0-standalone.jar data/soi.csv tmp

如果命令成功运行，你可能会看到 Hadoop 进程输出的日志信息。经过一段时间，你应该能看到作业输出的最终系数。

尽管目前执行需要更多时间，但我们的 Hadoop 作业有一个优势，那就是它可以分布在一个可以随着数据量不断扩展的集群上。

随机梯度下降

我们刚才看到的梯度下降计算方法通常被称为批量梯度下降，因为每次对系数的更新都是在对所有数据进行单次批量迭代的过程中完成的。对于大量数据来说，每次迭代可能会非常耗时，等待收敛可能需要很长时间。

梯度下降的另一种方法叫做随机梯度下降或SGD。在这种方法中，系数的估算会随着输入数据的处理不断更新。随机梯度下降的更新方法如下：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_14.jpg

实际上，这与批量梯度下降是完全相同的。应用的不同之处纯粹在于表达式 https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_17.jpg 是在 小批次——即数据的随机子集——上计算的。小批次的大小应该足够大，以便代表输入记录的公平样本——对于我们的数据，一个合理的小批次大小可能是 250。

随机梯度下降通过将整个数据集分成小批次并逐个处理它们来获得最佳估计。由于每个小批次的输出就是我们希望用于下一个小批次的系数（以便逐步改进估计），因此该算法本质上是顺序的。

随机梯度下降相较于批量梯度下降的优势在于，它可以在对数据集进行一次迭代后就获得良好的估计。对于非常大的数据集，甚至可能不需要处理所有小批次，就能达到良好的收敛效果。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_180.jpg

我们可以通过利用组合器串行应用的事实，使用 Tesser 实现 SGD，并将每个块视为一个小批次，从中可以计算出系数。这意味着我们的归约步骤就是恒等函数——我们不需要进行任何归约操作。

相反，让我们利用这个机会来学习如何在 Parkour 中构建一个 Hadoop 作业。在深入了解 Parkour 之前，让我们先看看如何使用我们已经掌握的知识实现随机梯度下降：

(defn stochastic-gradient-descent [options data]
  (let [batches (->> (into [] data)
                     (shuffle)
                     (partition 250))
        descend (fn [coefs batch]
                  (->> (gradient-descent-fold
                        (assoc options :coefs coefs))
                       (t/tesser (chunks batch))))]
    (reductions descend (:coefs options) batches)))

上述代码将输入集合分成较小的 250 元素组。对每个小批次运行梯度下降并更新系数。梯度下降的下一次迭代将在下一个批次上使用新的系数，并且对于适当的 alpha 值，生成改进的推荐结果。

以下代码将在数百个批次中绘制输出：

(defn ex-5-36 []
  (let [features [:A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount   (inc (count features))
        coefs    (vec (replicate fcount 0))
        data     (chunks (iota/seq "data/soi.csv"))
        factors  (->> (feature-scales features)
                      (t/tesser data))
        options  {:fy :A02300 :features features
                  :factors factors :coefs coefs :alpha 1e-3}
        ys       (stochastic-gradient-descent options data)
        xs       (range (count ys))]
    (-> (c/xy-plot xs (map first ys)
                   :x-label "Iterations"
                   :y-label "Coefficient")
        (c/add-lines xs (map #(nth % 1) ys))
        (c/add-lines xs (map #(nth % 2) ys))
        (c/add-lines xs (map #(nth % 3) ys))
        (c/add-lines xs (map #(nth % 4) ys))
        (i/view))))

我们提供的学习率比批量梯度下降的值小 100 倍。这有助于确保包含离群点的小批次不会使参数偏离其最优值。由于每个小批次固有的方差，随机梯度下降的输出将不会完全收敛到最优参数，而是会围绕最小值波动。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_190.jpg

上面的图片展示了随机梯度下降更为随机的效果；特别是小批次之间方差对参数估计的影响。尽管学习率大大降低，我们仍能看到与包含离群点的数据批次相对应的尖峰。

使用 Parkour 的随机梯度下降

在本章的其余部分，我们将直接使用 Parkour 构建 Hadoop 作业。Parkour 比 Tesser 更加暴露 Hadoop 的底层能力，这既是优点也是缺点。尽管 Tesser 使得在 Hadoop 中编写 fold 操作并应用于大型数据集变得非常容易，但 Parkour 需要我们更多地理解 Hadoop 的计算模型。

尽管 Hadoop 的 MapReduce 方法体现了我们在本章中遇到的许多原则，但它与 Tesser 的抽象在几个关键方面有所不同：

Hadoop 假定待处理的数据是键/值对
Hadoop 不要求在 map 之后进行 reduce 阶段
Tesser 在整个输入序列上进行折叠，Hadoop 在各组数据上进行归约
Hadoop 的值组由分区器定义
Tesser 的 combine 阶段发生在 reduce 之后，Hadoop 的 combine 阶段发生在 reduce 之前

最后一项特别令人遗憾。我们为 Clojure reducer 和 Tesser 学到的术语，在 Hadoop 中被反转了：在 Hadoop 中，combiner 会在数据被发送到 reducer 之前聚合 mapper 的输出。

我们可以在以下图示中看到广泛的流程，其中 mapper 的输出被合并成中间表示并在发送到 reducers 之前进行排序。每个 reducer 对整个数据的子集进行归约。combine 步骤是可选的，事实上，在我们的随机梯度下降作业中我们不需要一个：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_200.jpg

无论是否有 combine 步骤，数据都会在发送到 reducers 之前按组进行排序，分组策略由分区器定义。默认的分区方案是按键/值对的键进行分区（不同的键在前面的图示中由不同的灰度表示）。事实上，也可以使用任何自定义的分区方案。

如你所见，Parkour 和 Hadoop 都不假定输出是单一结果。由于 Hadoop 的 reduce 操作默认通过分组键定义的组进行，reduce 的结果可以是一个包含多个记录的数据集。在前面的图示中，我们为数据中的每个键展示了三个不同结果的情况。

定义 mapper

我们将定义的 Hadoop 作业的第一个组件是 mapper。mapper 的角色通常是接受一块输入记录并以某种方式转换它们。可以指定一个没有 reducer 的 Hadoop 作业；在这种情况下，mapper 的输出就是整个作业的输出。

Parkour 允许我们将映射器的操作定义为 Clojure 函数。该函数的唯一要求是，它需要将输入数据（无论是来自源文件，还是来自前一个 MapReduce 步骤的输出）作为最后一个参数传入。如果需要，还可以提供额外的参数，只要输入数据是最后一个参数：

(defn parse-m
  {::mr/source-as :vals
   ::mr/sink-as   :vals}
  [fy features factors lines]
  (->> (skip-header lines)
       (r/map parse-line)
       (r/map (partial format-record column-names))
       (r/map (scale-features factors))
       (r/map (extract-features fy features))
       (into [])
       (shuffle)
       (partition 250)))

前面的例子中，map 函数中的 parse-m（按照惯例，Parkour 映射器的后缀是 -m）负责将输入的单行数据转换成特征表示。我们重新使用了本章前面定义的许多函数：parse-line、format-record、scale-features 和 extract-features。Parkour 将以可简化集合的形式向映射器函数提供输入数据，因此我们将函数通过 r/map 进行链式调用。

随机梯度下降需要按小批量处理数据，因此我们的映射器负责将数据划分为 250 行的小组。在调用 partition 之前，我们进行 shuffle，以确保数据的顺序是随机的。

Parkour 造型函数

我们还向 parse-m 函数提供了元数据，形式为 {::mr/source-as :vals ::mr/sink-as :vals} 的映射。这两个命名空间关键字分别引用 parkour.mapreduce/source-as 和 parkour.mapreduce/sink-as，它们是指示 Parkour 在将数据提供给函数之前应该如何构造数据以及它可以期望返回的数据的结构。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_05_210.jpg

Parkour 映射器的有效选项有 :keyvals、:keys 和 :vals。前面的图示展示了三个键/值对的短序列的效果。通过请求将数据源作为 :vals，我们得到一个仅包含键/值对中值部分的序列。

定义一个归约器

在 Parkour 中定义一个归约器与定义映射器相同。最后一个参数必须是输入（现在是来自先前映射步骤的输入），但可以提供额外的参数。我们的随机梯度下降的 Parkour 归约器如下所示：

(defn sum-r
  {::mr/source-as :vals
   ::mr/sink-as   :vals}
  [fcount alpha batches]
  (let [initial-coefs (vec (replicate fcount 0))
        descend-batch (fn [coefs batch]
                        (->> (t/map (calculate-error
                                     (i/trans coefs)))
                             (t/fold (matrix-mean fcount 1))
                             (t/post-combine
                              (update-coefficients coefs alpha))
                             (t/tesser (chunks batch))))]
    (r/reduce descend-batch initial-coefs batches)))

我们的输入数据与之前一样作为可简化集合提供，因此我们使用 Clojure 的 reducers 库来进行迭代。我们使用 r/reduce 而不是 r/fold，因为我们不想对数据进行并行的归约处理。实际上，使用 Hadoop 的原因是我们可以独立控制映射阶段和归约阶段的并行度。现在，我们已经定义了映射和归约步骤，可以通过使用 parkour.graph 命名空间中的函数将它们组合成一个作业。

使用 Parkour 图定义 Hadoop 作业

graph 命名空间是 Parkour 用来定义 Hadoop 作业的主要 API。每个作业至少需要一个输入、一个映射器和一个输出，我们可以通过 Clojure 的 -> 宏将这些规格链式连接。首先定义一个非常简单的作业，它从我们的映射器获取输出并立即写入磁盘：

(defn hadoop-extract-features [conf workdir input output]
  (let [fy       :A02300
        features [:A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount   (inc (count features))
        input   (text/dseq input)
        factors (h/fold conf input (rand-file workdir)
                        #'feature-scales
                        features)
        conf (conf/ig)]
    (-> (pg/input input)
        (pg/map #'parse-m fy features factors)
        (pg/output (text/dsink output))
        (pg/execute conf "extract-features-job"))))

(defn ex-5-37 []
  (let [workdir  "tmp"
        out-file (rand-file workdir)]
    (hadoop-extract-features (conf/ig) "tmp"
                             "data/soi.csv" out-file)
    (str out-file)))

;; "tmp/1935333306"

前面的示例返回的响应应该是项目 tmp 目录中的一个子目录，Hadoop 会将文件写入该目录。如果你进入该目录，你应该会看到几个文件。在我的电脑上，我看到了四个文件——_SUCCESS、part-m-00000、part-m-00001 和 part-m-00002。_SUCCESS 文件的存在表示我们的作业已成功完成。part-m-xxxxx 文件是输入文件的块。

文件数量为三，说明 Hadoop 创建了三个映射器来处理我们的输入数据。如果我们以分布式模式运行，可能会并行创建这些文件。如果你打开其中一个文件，你应该看到一长串 clojure.lang.LazySeq@657d118e。由于我们写入的是文本文件，所以这是我们映射器输出的文本表示。

使用 Parkour 图链式连接映射器和 reducers

我们真正想做的是将映射和减少步骤链式连接，使它们一个接一个地执行。为此，我们必须插入一个中间步骤，分区器，并告诉分区器如何序列化我们的 clojure.lang.LazySeqs。

后者可以通过借用 Tesser 来实现，Tesser 实现了使用 Fressian 对任意 Clojure 数据结构进行序列化和反序列化的功能。在下一章中，我们将更深入地了解 Parkour 如何提供支持，帮助我们为分区数据创建良定义的架构，但现在仅仅让分区器通过编码后的数据就足够了。

注意

Fressian 是一种可扩展的二进制数据格式。你可以通过 github.com/clojure/data.fressian 的文档了解更多信息。

我们的键将编码为 FressianWritable，而我们的键没有明确指定（我们像处理 vals 一样处理我们的映射数据）。Hadoop 对 nil 的表示是 NullWritable 类型。我们通过以下方式将两者导入到我们的命名空间中：

(:import [org.apache.hadoop.io NullWritable]
         [tesser.hadoop_support FressianWritable])

在导入完成后，我们可以完全指定我们的作业：

(defn hadoop-sgd [conf workdir input-file output]
  (let [kv-classes [NullWritable FressianWritable]
        fy       :A02300
        features [:A00200 :AGI_STUB :NUMDEP :MARS2]
        fcount   (inc (count features))
        input   (text/dseq input-file)
        factors (h/fold conf input (rand-file workdir)
                        #'feature-scales
                        features)
        conf (conf/assoc! conf "mapred.reduce.tasks" 1)]
    (-> (pg/input input)
        (pg/map #'parse-m fy features factors)
        (pg/partition kv-classes)
        (pg/reduce #'sum-r fcount 1e-8)
        (pg/output (text/dsink output))
        (pg/execute conf "sgd-job"))))

我们需要确保只有一个 reducer 来处理我们的迷你批次（尽管 SGD 有一些变种允许我们平均多个随机梯度下降的结果，但我们希望得到一组接近最优的系数）。我们将使用 Parkour 的 conf 命名空间，通过 assoc! mapred.reduce.tasks 设置为 1。

在映射和减少步骤之间，我们指定分区器并传递在函数顶部定义的 kv-classes 函数。最后的示例仅仅运行这个作业：

(defn ex-5-38 []
  (let [workdir  "tmp"
        out-file (rand-file workdir)]
    (hadoop-sgd (conf/ig) "tmp" "data/soi.csv" out-file)
    (str out-file)))

;; "tmp/4046267961"

如果你进入作业返回的目录，你应该现在看到一个只包含两个文件的目录——_SUCCESS 和 part-r-00000。每个文件是一个 reducer 的输出，因此在只有一个 reducer 的情况下，我们最终得到了一个 part-r-xxxxx 文件。这个文件内部包含了使用随机梯度下降计算的线性模型系数。

总结

在本章中，我们学习了一些分布式数据处理的基本技术，并了解了本地用于数据处理的函数——映射（map）和归约（reduce），它们是处理即使是非常大量数据的强大方法。我们了解了 Hadoop 如何通过在数据的较小子集上运行函数，来扩展超越任何单一服务器的能力，这些子集的输出最终会组合起来生成结果。一旦你理解了其中的权衡，这种“分而治之”的数据处理方法，就是一种简单且非常通用的大规模数据分析方式。

我们看到了简单折叠操作在使用 Clojure 的 reducers 和 Tesser 时处理数据的强大功能和局限性。我们还开始探索 Parkour 如何暴露 Hadoop 更深层的能力。

在下一章，我们将看到如何使用 Hadoop 和 Parkour 来解决一个特定的机器学习挑战——对大量文本文档进行聚类。

第六章聚类

	具有共同特质的事物总会很快地寻找到它们自己的同类。
	–马库斯·奥勒留

在前几章中，我们涵盖了多个学习算法：线性回归和逻辑回归，C4.5，朴素贝叶斯和随机森林。在每种情况下，我们需要通过提供特征和期望输出来训练算法。例如，在线性回归中，期望输出是奥运游泳选手的体重，而对于其他算法，我们提供了一个类别：乘客是否生还。这些是监督学习算法的示例：我们告诉算法期望的输出，它将尝试学习生成相似的模型。

还有一类学习算法称为无监督学习。无监督算法能够在没有参考答案集的情况下操作数据。我们甚至可能不知道数据内部的结构；算法将尝试自行确定这种结构。

聚类是无监督学习算法的一个例子。聚类分析的结果是输入数据的分组，这些数据在某种方式上更相似于彼此。该技术是通用的：任何具有概念上相似性或距离的实体都可以进行聚类。例如，我们可以通过共享粉丝的相似性来聚类社交媒体账户组，或者通过测量受访者问卷答案的相似性来聚类市场调研的结果。

聚类的一个常见应用是识别具有相似主题的文档。这为我们提供了一个理想的机会来讨论文本处理，本章将介绍一系列处理文本的特定技术。

下载数据

本章使用Reuters-21578数据集：这是一组在 1987 年通过路透社新闻线发布的文章。它是最广泛用于测试文本分类和分类的数据集之一。文章文本和 Reuters-21578 集合中的注释版权归 Reuters Ltd.所有。Reuters Ltd.和 Carnegie Group, Inc.已同意仅为研究目的免费分发这些数据。

注意

您可以从 Packt Publishing 的网站或github.com/clojuredatascience/ch6-clustering下载本章的示例代码。

如往常一样，在示例代码中有一个脚本用于下载并解压文件到数据目录。您可以在项目目录中使用以下命令运行它：

script/download-data.sh

另外，在写作时，Reuters 数据集可以从kdd.ics.uci.edu/databases/reuters21578/reuters21578.tar.gz下载。本章的其余部分将假设文件已被下载并安装到项目的data目录中。

提取数据

运行前面的脚本后，文章将被解压到目录data/reuters-sgml中。每个.sgm文件包含大约 1,000 篇短文章，这些文章使用标准通用标记语言（SGML）被包装在 XML 样式的标签中。我们不需要自己编写解析器，而是可以利用已经在 Lucene 文本索引器中写好的解析器。

(:import [org.apache.lucene.benchmark.utils ExtractReuters])

(defn sgml->txt [in-path out-path]
  (let [in-file  (clojure.java.io/file in-path)
        out-file (clojure.java.io/file out-path)]
    (.extract (ExtractReuters. in-file out-file))))

在这里，我们利用 Clojure 的 Java 互操作性，简单地调用 Lucene 的ExtractReuters类中的提取方法。每篇文章都被提取为一个独立的文本文件。

可以通过执行以下代码运行：

lein extract-reuters

在项目目录内的命令行中运行。输出将是一个新目录data/reuters-text，其中包含超过 20,000 个独立的文本文件。每个文件包含一篇单独的路透社新闻文章。

如果磁盘空间不足，现在可以删除reuters-sgml和reuters21578.tar.gz文件：reuters-text目录中的内容是本章中我们唯一会使用的文件。现在让我们看看其中的几个文件。

检查数据

1987 年是“黑色星期一”那一年。10 月 19 日，全球股市暴跌，道琼斯工业平均指数下降了 508 点，降至 1738.74。像reut2-020.sgm-962.txt中包含的文章描述了这一事件：

19-OCT-1987 16:14:37.57

WALL STREET SUFFERS WORST EVER SELLOFF

Wall Street tumbled to its worst point loss ever and the worst percentage decline since the First World War as a frenzy of stock selling stunned even the most bearish market participants. "Everyone is in awe and the word 'crash' is on everyone's mind," one trader said.     The Dow Jones industrial average fell 508 points to 1738, a level it has not been at since the Autumn of 1986\.     Volume soared to 603 mln shares, almost doubling the previous record of 338 mln traded just last Friday. Reuter &#3;

这篇文章的结构代表了语料库中大多数文章的结构。第一行是时间戳，表示文章的发布时间，后面是一个空行。文章有一个标题，通常——但并非总是——是大写字母，然后是另一个空行。最后是文章正文文本。与处理此类半结构化文本时常见的情况一样，文章中存在多个空格、奇怪的字符和缩写。

其他文章只是标题，例如在reut2-020.sgm-761.txt中：

20-OCT-1987 17:09:34.49
REAGAN SAYS HE SEES NO RECESSION

这些文件是我们将进行聚类分析的对象。

聚类文本

聚类是找到彼此相似的对象组的过程。目标是，簇内的对象应该比簇外的对象更加相似。与分类一样，聚类并不是一个特定的算法，而是解决一般问题的算法类别。

尽管有多种聚类算法，但它们都在某种程度上依赖于距离度量。为了让算法判断两个对象是否属于同一簇或不同簇，它必须能够确定它们之间的距离（或者，如果你愿意，可以认为是相似度）的定量度量。这就需要一个数字化的距离度量：距离越小，两个对象之间的相似度就越大。

由于聚类是一种可以应用于多种数据类型的通用技术，因此有大量可能的距离度量方法。尽管如此，大多数数据都可以通过一些常见的抽象表示：集合、空间中的点或向量。对于这些，每种情况都有一种常用的度量方法。

单词集合与 Jaccard 指数

如果你的数据可以表示为事物的集合，则可以使用 Jaccard 指数，也称为Jaccard 相似度。它在概念上是最简单的度量之一：它是集合交集除以集合并集，或者是集合中共同元素的数量占总独特元素数量的比例：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_01.jpg

许多事物可以表示为集合。社交网络上的账户可以表示为朋友或关注者的集合，顾客可以表示为购买或查看过的商品集合。对于我们的文本文件，集合表示可以简单地是使用的唯一单词的集合。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_100.jpg

Jaccard 指数在 Clojure 中非常简单计算：

(:require [clojure.set :as set])

(defn jaccard-similarity [a b]
  (let [a (set a)
        b (set b)]
    (/ (count (set/intersection a b))
       (count (set/union a b)))))

(defn ex-6-1 []
  (let [a [1 2 3]
        b [2 3 4]]
    (jaccard a b)))

;; => 1/2

它的优点是，即使集合的基数不同，距离度量依然有意义。在前面的图示中，A“比”B要“大”，但交集除以并集仍然能公平地反映它们的相似度。要将 Jaccard 指数应用于文本文件，我们需要将它们转化为单词集合。这就是词元化的过程。

词元化 Reuters 文件

词元化是将一串文本拆分成更小的单元以便进行分析的技术名称。常见的方法是将文本字符串拆分成单个单词。一个明显的分隔符是空格，因此 "tokens like these" 会变成 ["tokens" "like" "these"]。

(defn tokenize [s]
  (str/split s #"\W+"))

这既方便又简单，但不幸的是，语言是微妙的，很少有简单的规则可以普遍适用。例如，我们的词元化器将撇号视为空格：

(tokenize "doesn't handle apostrophes")
;; ["doesn" "t" "handle" "apostrophes"]

连字符也被视为空格：

(tokenize "good-looking user-generated content")
;; ["good" "looking" "user" "generated" "content"]

而删除它们则会改变句子的含义。然而，并非所有的连字符都应该被保留：

(tokenize "New York-based")
;; ["New" "York" "based"]

"New"、"York" 和 "based" 正确地表示了短语的主体，但最好将 "New York" 归为一个词，因为它代表了一个特定的名称，应该完整保留。另一方面，York-based 单独作为词元没有意义。

简而言之，文本是杂乱的，从自由文本中可靠地解析出意义是一个极其丰富和活跃的研究领域。特别是，对于从文本中提取名称（例如，“纽约”），我们需要考虑术语使用的上下文。通过其语法功能标注句子中词汇的技术被称为词性标注器。

注意

如需了解更多关于高级分词和词性标注的信息，请参见clojure-opennlp库：github.com/dakrone/clojure-opennlp。

在本章中，我们有足够多的文档可以使用，因此我们将继续使用我们的简单分词器。我们会发现，尽管它存在一些缺点，它仍然能够足够好地从文档中提取意义。

让我们编写一个函数，根据文档的文件名返回文档的标记：

(defn tokenize-reuters [content]
  (-> (str/replace content  #"^.*\n\n" "")
      (str/lower-case)
      (tokenize)))

(defn reuters-terms [file]
  (-> (io/resource file)
      (slurp)
      (tokenize-reuters)))

我们正在去除文件顶部的时间戳，并在分词之前将文本转换为小写。在下一部分，我们将看到如何衡量分词后的文档相似性。

应用 Jaccard 指数到文档上

在对输入文档进行分词之后，我们可以将结果的标记序列简单地传递给我们之前定义的jaccard-similarity函数。让我们比较一下来自路透社语料库的几篇文档相似性：

(defn ex-6-2 []
  (let [a (set (reuters-terms "reut2-020.sgm-761.txt"))
        b (set (reuters-terms "reut2-007.sgm-750.txt"))
        s (jaccard a b)]
    (println "A:" a)
    (println "B:" b)
    (println "Similarity:" s)))

A: #{recession says reagan sees no he}
B: #{bill transit says highway reagan and will veto he}
Similarity: 1/4

Jaccard 指数输出一个介于零和一之间的数字，因此它认为这两篇文档在标题中的词汇相似度为 25%。注意我们丢失了标题中词汇的顺序。没有我们稍后会讲到的其他技巧，Jaccard 指数只关注两个集合中共同的项目。我们丢失的另一个方面是一个词在文档中出现的次数。重复出现同一个词的文档，可能会在某种意义上视该词为更重要。例如，reut2-020.sgm-932.txt的标题如下：

19-OCT-1987 16:41:40.58
NYSE CHAIRMAN JOHN PHELAN SAYS NYSE WILL OPEN TOMORROW ON TIME

NYSE 在标题中出现了两次。我们可以推测这个标题特别关注纽约证券交易所，可能比只提到一次 NYSE 的标题更为重要。

词袋模型和欧几里得距离

一种可能优于词集方法的改进是词袋模型方法。这种方法保留了文档中各个词汇的词频。通过距离度量可以将词频纳入考虑，从而可能得到更准确的结果。

最常见的距离概念之一是欧几里得距离度量。在几何学中，欧几里得度量是我们计算空间中两点之间距离的方式。在二维空间中，欧几里得距离由毕达哥拉斯公式给出：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_02.jpg

这表示两个点之间的差异是它们之间直线距离的长度。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_110.jpg

这可以扩展到三维：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_03.jpg

并且推广到 n 维度：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_04.jpg

其中 A[i] 和 B[i] 是在维度 i 上 A 或 B 的值。因此，距离度量就是两个文档之间的整体相似性，考虑了每个单词出现的频率。

(defn euclidean-distance [a b]
  (->> (map (comp i/sq -) a b)
       (apply +)
       (i/sqrt)))

由于每个单词现在代表空间中的一个维度，我们需要确保在计算欧几里得距离时，我们是在比较每个文档中同一维度上的大小。否则，我们可能会字面意义上比较“苹果”和“橙子”。

将文本表示为向量

与 Jaccard 指数不同，欧几里得距离依赖于将单词一致地排序为各个维度。词频（term frequency）表示文档在一个大型多维空间中的位置，我们需要确保在比较值时，比较的是正确维度中的值。我们将文档表示为术语频率向量。

想象一下，文档中可能出现的所有单词都被赋予一个唯一的编号。例如，单词“apple”可以被赋值为 53，单词“orange”可以被赋值为 21,597。如果所有数字都是唯一的，它们可以对应于单词在术语向量中出现的索引。

这些向量的维度可能非常大。最大维度数是向量的基数。对应于单词的索引位置的元素值通常是该单词在文档中出现的次数。这被称为术语频率（tf）加权。

为了能够比较文本向量，重要的是相同的单词始终出现在向量中的相同索引位置。这意味着我们必须为每个创建的向量使用相同的单词/索引映射。这个单词/索引映射就是我们的词典。

创建词典

为了创建一个有效的词典，我们需要确保两个单词的索引不会冲突。做到这一点的一种方法是使用一个单调递增的计数器，每添加一个单词到词典中，计数器就增加一次。单词被添加时的计数值将成为该单词的索引。为了线程安全地同时将单词添加到词典并递增计数器，我们可以使用原子操作：

(def dictionary
  (atom {:count 0
         :words {}}))

(defn add-term-to-dict [dict word]
  (if (contains? (:terms dict) word)
    dict
    (-> dict
        (update-in [:terms] assoc word (get dict :count))
        (update-in [:count] inc))))

(defn add-term-to-dict! [dict term]
  (doto dict
    (swap! add-term-to-dict term)))

为了对原子进行更新，我们必须在swap!函数中执行我们的代码。

(add-term-to-dict! dictionary "love")

;; #<Atom@261d1f0a: {:count 1, :terms {"love" 0}}>

添加另一个单词将导致计数增加：

(add-term-to-dict! dictionary "music")

;; #<Atom@261d1f0a: {:count 2, :terms {"music" 1, "love" 0}}>

并且重复添加相同单词不会产生任何效果：

(add-term-to-dict! dictionary "love")

;; #<Atom@261d1f0a: {:count 2, :terms {"music" 1, "love" 0}}>

在原子操作中执行此更新可以确保即使在多个线程同时更新词典时，每个单词也能获得自己的索引。

(defn build-dictionary! [dict terms]
  (reduce add-term-to-dict! dict terms))

构建整个词典就像在提供的词典原子上使用add-term-to-dict!函数并对一组术语进行简化操作一样简单。

创建术语频率向量

为了计算欧几里得距离，首先让我们从字典和文档中创建一个向量。这将使我们能够轻松地比较文档之间的术语频率，因为它们将占据向量的相同索引。

(defn term-id [dict term]
  (get-in @dict [:terms term]))

(defn term-frequencies [dict terms]
  (->> (map #(term-id dict %) terms)
       (remove nil?)
       (frequencies)))

(defn map->vector [dictionary id-counts]
  (let [zeros (vec (replicate (:count @dictionary) 0))]
    (-> (reduce #(apply assoc! %1 %2) (transient zeros) id-counts)
        (persistent!))))

(defn tf-vector [dict document]
  (map->vector dict (term-frequencies dict document)))

term-frequencies 函数为文档中的每个术语创建一个术语 ID 到频率计数的映射。map->vector 函数简单地接受此映射，并在由术语 ID 给出的向量索引位置关联频率计数。由于可能有许多术语，且向量可能非常长，因此我们使用 Clojure 的 transient! 和 persistent! 函数暂时创建一个可变向量以提高效率。

让我们打印 reut2-020.sgm-742.txt 的文档、字典和生成的向量：

(defn ex-6-3 []
  (let [doc  (reuters-terms "reut2-020.sgm-742.txt")
        dict (build-dictionary! dictionary doc)]
    (println "Document:" doc)
    (println "Dictionary:" dict)
    (println "Vector:" (tf-vector dict doc))))

输出结果如下所示（格式已调整以提高可读性）：

;; Document: [nyse s phelan says nyse will continue program
;;            trading curb until volume slows]
;; Dictionary: #<Atom@bb156ec: {:count 12, :terms {s 1, curb 8,
;;             phelan 2, says 3, trading 7, nyse 0, until 9,
;;             continue 5, volume 10, will 4, slows 11,
;;             program 6}}>
;; Vector: [2 1 1 1 1 1 1 1 1 1 1 1]

输入中有 12 个术语，字典中有 12 个术语，并返回了一个包含 12 个元素的向量。

(defn print-distance [doc-a doc-b measure]
  (let [a-terms (reuters-terms doc-a)
        b-terms (reuters-terms doc-b)
        dict (-> dictionary
                 (build-dictionary! a-terms)
                 (build-dictionary! b-terms))
        a (tf-vector dict a-terms)
        b (tf-vector dict b-terms)]
    (println "A:" a)
    (println "B:" b)
    (println "Distance:" (measure a b))))

(defn ex-6-4 []
  (print-distance "reut2-020.sgm-742.txt"
                  "reut2-020.sgm-932.txt"
                  euclidean-distance))

;; A: [2 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0]
;; B: [2 0 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1]
;; Distance: 3.7416573867739413

与 Jaccard 指数类似，欧几里得距离不能低于零。然而，不同于 Jaccard 指数，欧几里得距离的值可以无限增长。

向量空间模型与余弦距离

向量空间模型可以视为词集（set-of-words）和词袋（bag-of-words）模型的推广。与词袋模型类似，向量空间模型将每个文档表示为一个向量，每个元素表示一个术语。每个索引位置的值是该词的重要性度量，可能是也可能不是术语频率。

如果你的数据从概念上表示一个向量（也就是说，一个特定方向上的大小），那么余弦距离可能是最合适的选择。余弦距离度量通过计算两个元素的向量表示之间夹角的余弦值来确定它们的相似度。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_120.jpg

如果两个向量指向相同的方向，那么它们之间的夹角为零，而零的余弦值为 1。余弦相似度可以通过以下方式定义：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_05.jpg

这是一个比我们之前讨论的方程更为复杂的方程。它依赖于计算两个向量的点积及其大小。

(defn cosine-similarity [a b]
  (let [dot-product (->> (map * a b)
                         (apply +))
        magnitude (fn [d]
                    (->> (map i/sq d)
                         (apply +)
                         (i/sqrt)))]
    (/ dot-product (* (magnitude a) (magnitude b)))))

余弦相似度的例子如下所示：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_130.jpg

余弦相似度通常作为高维空间中的相似度度量使用，其中每个向量包含很多零，因为它计算起来非常高效：只需要考虑非零维度。由于大多数文本文档只使用了所有单词中的一小部分（因此在大多数维度上是零），余弦度量通常用于文本聚类。

在向量空间模型中，我们需要一个一致的策略来衡量每个词项的重要性。在词集合模型中，所有词项被平等对待。这相当于将该点的向量值设置为 1。在词袋模型中，计算了词项的频率。我们目前将继续使用词频，但很快我们将看到如何使用一种更复杂的重要性度量，称为词频-逆文档频率（TF-IDF）。

(defn ex-6-5 []
  (print-distance "reut2-020.sgm-742.txt"
                  "reut2-020.sgm-932.txt"
                  cosine-similarity))

;; A: [2 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0]
;; B: [2 0 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1]
;; Distance: 0.5012804118276031

余弦值越接近1，这两个实体越相似。为了将余弦相似度转化为距离度量，我们可以简单地将余弦相似度从1中减去。

尽管前面提到的所有度量方法对相同输入产生不同的度量值，但它们都满足一个约束，即* A 和 B 之间的距离应该与 B 和 A *之间的差异相同。通常，相同的底层数据可以被转换为表示集合（Jaccard）、空间中的一个点（欧几里得）或一个向量（余弦）。有时，唯一的方法是尝试并查看结果如何，才能知道哪个是正确的。

出现在一个文档中的唯一单词数量通常比处理中的文档集合中任何文档中出现的唯一单词数量要少。因此，这些高维文档向量通常是非常稀疏的。

去除停用词

许多头条新闻之间的相似性是由一些经常出现的词汇所产生的，这些词汇对内容的意义贡献不大。例如，“a”，“says”和“and”。我们应该过滤掉这些词，以避免产生虚假的相似性。

考虑以下两个习语：

"音乐是爱的食粮"
"战争是历史的动力"

我们可以使用以下 Clojure 代码来计算它们之间的余弦相似度：

(defn ex-6-6 []
  (let [a (tokenize "music is the food of love")
        b (tokenize "war is the locomotive of history")]
    (add-documents-to-dictionary! dictionary [a b])
    (cosine-similarity (tf-vector dictionary a)
                       (tf-vector dictionary b))))

;; 0.5

尽管这两个文档之间唯一的共同单词是is、the和of，它们的相似度为0.5。理想情况下，我们希望移除这些词。

词干提取

现在让我们考虑一个替代表达：

"音乐是爱的食粮"
"你很有音乐天赋，真好"

让我们也来比较它们的余弦相似度：

(defn ex-6-7 []
  (let [a (tokenize "music is the food of love")
        b (tokenize "it's lovely that you're musical")]
    (add-documents-to-dictionary! dictionary [a b])
    (cosine-similarity (tf-vector dictionary a)
                       (tf-vector dictionary b))))

;; 0.0

尽管这两个句子都提到音乐和积极情感，但这两个短语的余弦相似度为零：两个短语之间没有共同的单词。这是有道理的，但并没有表达我们通常想要的行为，即捕捉“概念”之间的相似性，而不是精确使用的词汇。

解决这个问题的一种方法是词干提取，它将单词简化为其词根。具有共同意义的单词更可能提取到相同的词根。Clojure 库的词干提取器（github.com/mattdw/stemmers）可以为我们完成这项工作，幸运的是，它们也会去除停用词。

(defn ex-6-8 []
  (let [a (stemmer/stems "music is the food of love")
        b (stemmer/stems "it's lovely that you're musical")]
    (add-documents-to-dictionary! dictionary [a b])
    (cosine-similarity (tf-vector dictionary a)
                       (tf-vector dictionary b))))

;; 0.8164965809277259

好得多了。经过词干提取和去除停用词后，短语之间的相似度从 0.0 降低到 0.82。这是一个很好的结果，因为尽管句子使用了不同的词汇，它们表达的情感是相关的。

使用 k-means 和 Incanter 进行聚类

最终，在对输入文档进行标记化、词干提取和向量化处理，并且选择了不同的距离度量方式后，我们可以开始对数据进行聚类。我们将首先探讨的聚类算法是 k-means 聚类。

k-means 是一个迭代算法，步骤如下：

随机选择 k 个聚类中心点。
将每个数据点分配到与其最近的聚类中心点所属的聚类中。
调整每个聚类中心点的位置，使其位于其分配数据点的均值位置。
重复进行，直到收敛或达到最大迭代次数。

该过程在以下图示中展示，适用于 k=3 个聚类：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_140.jpg

在上图中，我们可以看到，在第一次迭代时，初始的聚类中心点并没有很好地代表数据的结构。尽管这些点显然被分成了三组，但初始的聚类中心点（以十字标记表示）都分布在图形的上方区域。这些点的颜色表示它们与最近的聚类中心的关系。随着迭代的进行，我们可以看到聚类中心点如何逐步靠近每个点组的“自然”位置，即组中心。

在定义主要的 k-means 函数之前，先定义几个实用的辅助函数是有用的：一个用于计算聚类中心点的函数，另一个则是将数据分组到各自聚类的函数。

(defn centroid [xs]
  (let [m (i/trans (i/matrix xs))]
    (if (> (i/ncol m) 1)
      (i/matrix (map s/mean m))
     m)))

(defn ex-6-9 []
  (let [m (i/matrix [[1 2 3]
                     [2 2 5]])]
    (centroid m)))

;; A 3x1 matrix
;;  -------------
;; 1.50e+00
;; 2.00e+00
;; 4.00e+00

centroid 函数简单地计算输入矩阵每列的均值。

(defn clusters [cluster-ids data]
  (->> (map vector cluster-ids data)
       (conj-into {})
       (vals)
       (map i/matrix))) 

(defn ex-6-10 []
  (let [m (i/matrix [[1 2 3]
                     [4 5 6]
                     [7 8 9]])]
    (clusters [0 1 0] m)))

;; A 1x3 matrix
;; -------------
;; 4.00e+00  5.00e+00  6.00e+00 
;;  A 2x3 matrix
;; -------------
;; 7.00e+00  8.00e+00  9.00e+00 
;; 1.00e+00  2.00e+00  3.00e+00

clusters 函数根据提供的聚类 ID 将较大的矩阵拆分为一系列较小的矩阵。聚类 ID 被提供为与聚类点长度相同的元素序列，列出该序列中每个点对应的聚类 ID。共享相同聚类 ID 的项目将被分到一起。通过这两个函数，我们得到了完整的 k-means 函数：

(defn k-means [data k]
  (loop [centroids (s/sample data :size k)
         previous-cluster-ids nil]
    (let [cluster-id (fn [x]
                       (let [distance  #(s/euclidean-distance x %)
                             distances (map distance centroids)]
                         (->> (apply min distances)
                              (index-of distances))))
          cluster-ids (map cluster-id data)]
      (if (not= cluster-ids previous-cluster-ids)
        (recur (map centroid (clusters cluster-ids data))
               cluster-ids)
        clusters))))

我们首先通过抽样输入数据随机选择 k 个聚类中心点。接着，使用循环/递归不断更新聚类中心点，直到 previous-cluster-ids 和 cluster-ids 相同。此时，所有文档都没有移动到其他聚类，因此聚类过程已收敛。

对路透社文档进行聚类

现在让我们使用 k-means 函数对路透社文档进行聚类。我们先让算法简单一些，选取一些较大的文档样本。较大的文档更容易让算法识别它们之间的相似性。我们将最低字符数设定为 500 字符。这意味着我们的输入文档至少有一个标题和几句正文内容。

(defn ex-6-11 []
  (let [documents (fs/glob "data/reuters-text/*.txt")
        doc-count 100
        k 5
        tokenized (->> (map slurp documents)
                       (remove too-short?)
                       (take doc-count)
                       (map stem-reuters))]
    (add-documents-to-dictionary! dictionary tokenized)
    (-> (map #(tf-vector dictionary %) tokenized)
        (k-means k))))

我们使用fs库（github.com/Raynes/fs）通过调用fs/glob，并使用匹配所有文本文件的模式，来创建文件列表进行聚类。我们删除那些太短的文件，标记前 100 个，并将它们添加到字典中。我们为输入创建tf向量，然后对它们调用k-means。

如果你运行前面的示例，你将得到一组聚类文档向量，但这些并不太有用。让我们创建一个summary函数，利用字典报告每个聚类中最常见的术语。

(defn cluster-summary [dict clusters top-term-count]
  (for [cluster clusters]
    (let [sum-terms (if (= (i/nrow cluster) 1)
                      cluster
                      (->> (i/trans cluster)
                           (map i/sum)
                           (i/trans)))
          popular-term-ids (->> (map-indexed vector sum-terms)
                                (sort-by second >)
                                (take top-term-count)
                                (map first))
          top-terms (map #(id->term dict %) popular-term-ids)]
      (println "N:" (i/nrow cluster))
      (println "Terms:" top-terms))))

(defn ex-6-12 []
  (cluster-summary dictionary (ex-6-11) 5))

k-均值算法本质上是一个随机算法，对质心的初始位置敏感。我得到了以下输出，但你的结果几乎肯定会不同：

;; N: 2
;; Terms: (rocket launch delta satellit first off weather space)
;;  N: 4
;; Terms: (said will for system 000 bank debt from bond farm)
;; N: 12
;; Terms: (said reuter for iranian it iraq had new on major)
;; N: 62
;; Terms: (said pct dlr for year mln from reuter with will)
;; N: 20
;; Terms: (said for year it with but dlr mln bank week)

不幸的是，我们似乎没有得到很好的结果。第一个聚类包含两篇关于火箭和太空的文章，第三个聚类似乎由关于伊朗的文章组成。大多数文章中最常见的词汇是“said”。

使用 TF-IDF 进行更好的聚类

词汇 频率-逆文档频率（TF-IDF）是一种通用方法，用于在文档向量中加权术语，以便在整个数据集中流行的术语不会像那些较不常见的术语那样被高估。这捕捉了直观的信念——也是我们之前观察到的——即“said”这样的词汇并不是构建聚类的强有力基础。

Zipf 定律

Zipf 定律指出，任何词汇的频率与其在频率表中的排名成反比。因此，最频繁的词汇出现的频率大约是第二常见词汇的两倍，第三常见词汇的三倍，依此类推。让我们看看这一规律是否适用于我们的路透社语料库：

(defn ex-6-13 []
  (let [documents (fs/glob "data/reuters-text/*.txt")
        doc-count 1000
        top-terms 25
        term-frequencies (->> (map slurp documents)
                              (remove too-short?)
                              (take doc-count)
                              (mapcat tokenize-reuters)
                              (frequencies)
                              (vals)
                              (sort >)
                              (take top-terms))]
    (-> (c/xy-plot (range (inc top-terms)) term-frequencies
                   :x-label "Terms"
                   :y-label "Term Frequency")
        (i/view))))

使用前述代码，我们可以计算出前 1,000 篇路透社文档中前 25 个最流行词汇的频率图。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_150.jpg

在前 1,000 篇文档中，最流行的词汇出现了将近 10,000 次。第 25^(th)个最流行的词汇总共出现了大约 1,000 次。实际上，数据显示，词汇在路透社语料库中的出现频率比频率表中的位置所暗示的要高。这很可能是由于路透社语料库的公告性质，导致相同的短词被反复使用。

计算 TF-IDF 权重

计算 TF-IDF 只需要对我们已经创建的代码进行两个修改。首先，我们必须追踪给定术语出现在哪些文档中。其次，在构建文档向量时，我们必须适当地加权该术语。

既然我们已经创建了一个术语字典，我们不妨将每个术语的文档频率存储在其中。

(defn inc-df! [dictionary term-id]
  (doto dictionary
    (swap! update-in [:df term-id] (fnil inc 0))))

(defn build-df-dictionary! [dictionary document]
  (let [terms    (distinct document)
        dict     (build-dictionary! dictionary document)
        term-ids (map #(term-id dictionary %) document)]
    (doseq [term-id term-ids]
      (inc-df! dictionary term-id))
    dict))

build-df-dictionary 函数之前接受一个字典和一个术语序列。我们从不同的术语中构建字典，并查找每个术语的 term-id。最后，我们遍历术语 ID 并为每个 ID 增加 :df。

如果一个文档包含词语 w[1]、…、w[n]，那么词语 w[i] 的逆文档频率定义为：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_06.jpg

即，词语出现在文档中的倒数。如果一个词语在一组文档中频繁出现，它的 DF 值较大，而 IDF 值较小。在文档数量很大的情况下，通常通过乘以一个常数（通常是文档总数 N）来对 IDF 值进行归一化，因此 IDF 公式看起来像这样：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_07.jpg

词语 w[i] 的 TF-IDF 权重 W[i] 由词频和逆文档频率的乘积给出：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_08.jpg

然而，前述公式中的 IDF 值仍然不理想，因为对于大型语料库，IDF 项的范围通常远大于 TF，并且可能会压倒 TF 的效果。为了减少这个问题并平衡 TF 和 IDF 项的权重，通常的做法是使用 IDF 值的对数：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_09.jpg

因此，词语 w[i] 的 TF-IDF 权重 w[i] 变为：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_10.jpg

这是经典的 TF-IDF 权重：常见词语的权重较小，而不常见的词语的权重较大。确定文档主题的重要词语通常具有较高的 TF 和适中的 IDF，因此二者的乘积成为一个较大的值，从而在结果向量中赋予这些词语更多的权重。

(defn document-frequencies [dict terms]
  (->> (map (partial term-id dict) terms)
       (select-keys (:df @dict))))

(defn tfidf-vector [dict doc-count terms]
  (let [tf (term-frequencies dict terms)
        df (document-frequencies dict (distinct terms))
        idf   (fn [df] (i/log (/ doc-count df)))
        tfidf (fn [tf df] (* tf (idf df)))]
    (map->vector dict (merge-with tfidf tf df))))

前述代码计算了先前定义的 term-frequencies 和从字典中提取的 document-frequencies 的 TF-IDF 值。

使用 TF-IDF 的 k-means 聚类

在进行前述调整后，我们可以计算路透社文档的 TF-IDF 向量。以下示例是基于新的 tfidf-vector 函数修改的 ex-6-12：

(defn ex-6-14 []
  (let [documents (fs/glob "data/reuters-text/*.txt")
        doc-count 100
        k 5
        tokenized (->> (map slurp documents)
                       (remove too-short?)
                       (take doc-count)
                       (map stem-reuters))]
    (reduce build-df-dictionary! dictionary tokenized)
    (-> (map #(tfidf-vector dictionary doc-count %) tokenized)
        (k-means k)
        (cluster-summary dictionary 10))))

前述代码与之前的示例非常相似，但我们已经替换了新的 build-df-dictionary 和 tfidf-vector 函数。如果你运行该示例，应该会看到比之前稍微改进的输出：

N: 5
Terms: (unquot unadjust year-on-year novemb grew sundai labour m-3 ahead 120)
N: 15
Terms: (rumor venezuela azpurua pai fca keat ongpin boren gdp moder)
N: 16
Terms: (regan drug lng soviet bureau deleg gao dean fdic algerian)
N: 46
Terms: (form complet huski nrc rocket north underwrit card oat circuit)
N: 18
Terms: (freez cocoa dec brown bean sept seixa telex argentin brown-forman)

尽管由于词语已被词干化，顶级词语可能难以解释，但这些词语代表了每个聚类中最为不寻常的常见词。注意，“said”不再是所有聚类中评分最高的词语。

更好的 n-gram 聚类

从前面列出的词汇表来看，应该清楚地意识到，通过将文档简化为无序的词汇序列，我们丧失了多少信息。如果没有句子的上下文，我们很难准确把握每个聚类的含义。

然而，向量空间模型本身并没有什么内在的限制，阻止我们保持输入词语的顺序。我们可以简单地创建一个新的术语来表示多个词的组合。这个组合词，可能表示多个连续的输入词，被称为 n-gram。

一个 n-gram 的例子可能是“new york”（纽约）或“stock market”（股市）。事实上，因为它们包含两个词，所以这些被称为 bigrams（二元组）。n-grams 可以是任意长度的。n-gram 越长，携带的上下文信息越多，但它的出现也就越为罕见。

n-grams 与 **瓦片法（shingling）**的概念密切相关。当我们对 n-grams 进行瓦片化处理时，我们是在创建重叠的词组序列。瓦片化（shingling）一词来源于这些词语像屋顶瓦片一样重叠的方式。

(defn n-grams [n words]
  (->> (partition n 1 words)
       (map (partial str/join " ")))) 

(defn ex-6-15 []
  (let [terms (reuters-terms "reut2-020.sgm-761.txt")]
    (n-grams 2 terms)))

;; ("reagan says" "says he" "he sees" "sees no" "no recession")

目前，使用 2-grams 就能让我们（例如）区分数据集中“coconut”一词的不同用法：“coconut oil”（椰子油）、“coconut planters”（椰子种植者）、“coconut plantations”（椰子种植园）、“coconut farmers”（椰子农民）、“coconut association”（椰子协会）、“coconut authority”（椰子管理机构）、“coconut products”（椰子产品）、“coconut exports”（椰子出口）、“coconut industry”（椰子产业）和相当吸引人的“coconut chief”（椰子首领）。这些词对定义了不同的概念——有时是微妙的不同——我们可以在不同的文档中捕捉并进行比较。

我们可以通过结合不同长度的 n-gram 的结果，来实现 n-gram 和瓦片化的双重优势：

(defn multi-grams [n words]
  (->> (range 1 (inc n))
       (mapcat #(n-grams % words))))

(defn ex-6-16 []
  (let [terms (reuters-terms "reut2-020.sgm-761.txt")]
    (multi-grams 4 terms)))

;; ("reagan" "says" "he" "sees" "no" "recession" "reagan says" 
;; "says he" "he sees" "sees no" "no recession" "reagan says he" 
;; "says he sees" "he sees no" "sees no recession" "reagan says he 
;; sees" "says he sees no" "he sees no recession")

虽然词干提取（stemming）和停用词移除（stop word removal）起到了缩减字典规模的作用，而使用 TF-IDF 则提高了每个词在文档中的权重效用，但生成 n-grams 的效果是大大增加了我们需要处理的词汇数量。

特征的爆炸性增长将立刻使我们在 Incanter 中实现 k-means 的效果超出负荷。幸运的是，有一个叫做 Mahout 的机器学习库，专门设计用于在海量数据上运行像 k-means 这样的算法。

使用 Mahout 进行大规模聚类

Mahout (mahout.apache.org/) 是一个机器学习库，旨在分布式计算环境中使用。该库的 0.9 版本支持 Hadoop，并且是我们将在此使用的版本。

注意

在本文撰写时，Mahout 0.10 刚刚发布，并且同样支持 Spark。Spark 是一个替代性的分布式计算框架，我们将在下一章介绍。

我们在上一章看到，Hadoop 的一个抽象概念是序列文件：Java 键和值的二进制表示。许多 Mahout 的算法期望在序列文件上操作，我们需要创建一个序列文件作为 Mahout 的 k-means 算法的输入。Mahout 的 k-means 算法也期望将其输入作为向量，表示为 Mahout 的向量类型之一。

尽管 Mahout 包含提取向量的类和实用程序，我们将借此机会演示如何将 Parkour 和 Mahout 一起使用。这样不仅能让我们对创建的向量有更细粒度的控制，还可以展示更多 Parkour 在指定 Hadoop 作业时的能力。

将文本文件转换为序列文件

然而，我们不会定义自定义作业来将文本文档转换为序列文件表示：Mahout 已经定义了一个有用的 SequenceFilesFromDirectory 类，用来转换文本文件目录。我们将使用这个工具来创建一个代表整个 reuters-txt 目录内容的单一文件。

尽管序列文件可能物理上分布在不同的块中（例如在 HDFS 上），但它在逻辑上是一个文件，表示所有输入文档作为键/值对。键是文件名，值是文件的文本内容。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_160.jpg

以下代码将处理转换：

(:import [org.apache.mahout.text
           SequenceFilesFromDirectory])

(defn text->sequencefile [in-path out-path]
  (SequenceFilesFromDirectory/main
   (into-array String (vector "-i" in-path
                              "-o" out-path
                              "-xm" "sequential"
                              "-ow"))))

(defn ex-6-17 []
  (text->sequencefile "data/reuters-text"
                      "data/reuters-sequencefile"))

SequenceFilesFromDirectory 是 Mahout 的一个实用类，是一套可以在命令行调用的类之一。

提示

由于运行前面的示例是后续示例的先决条件，它也可以在命令行上运行：

lein create-sequencefile

我们直接调用 main 函数，传递我们通常在命令行上传递的参数，作为字符串数组。

使用 Parkour 创建 Mahout 向量

现在我们已经有了 Reuters 语料库的序列文件表示，我们需要将每个文档（现在表示为一个单一的键/值对）转换为向量。我们之前已经看到如何使用共享字典（被建模为 Clojure 原子）来实现这一点。原子确保即使在多线程环境中，每个不同的术语也会得到自己的 ID。

我们将使用 Parkour 和 Hadoop 来生成向量，但这带来了一个挑战。由于 MapReduce 编程的特点是映射器并行工作且不共享状态，我们如何为每个单词分配一个唯一的 ID 呢？Hadoop 并未提供类似于 Clojure 原子的机制来在集群中的节点之间共享可变状态，实际上，最小化共享状态是扩展分布式应用程序的关键。

因此，创建一个共享的唯一 ID 集合对我们的 Parkour 作业来说是一个有趣的挑战：让我们看看如何以分布式方式为我们的字典生成唯一的 ID。

创建分布式唯一 ID

在我们查看 Hadoop 特定的解决方案之前，值得注意的是，创建一个集群范围内唯一标识符的一个简单方法是创建一个通用唯一标识符，或 UUID。

(defn uuid []
  (str (java.util.UUID/randomUUID)))

这会创建一个形式为 3a65c7db-6f41-4087-a2ec-8fe763b5f185 的长字节字符串，这几乎可以保证不会与世界上任何地方生成的其他 UUID 冲突。

虽然这对于生成唯一 ID 有效，但可能的 ID 数量是天文数字，而 Mahout 的稀疏向量表示需要以整数的形式初始化向量的基数。使用 uuid 生成的 ID 太大了。此外，这并不能帮助我们协调 ID 的创建：集群中的每台机器都会生成不同的 UUID 来表示相同的术语。

解决这个问题的一种方法是使用术语本身来生成唯一 ID。如果我们使用一致性哈希函数从每个输入术语创建一个整数，集群中的所有机器都会生成相同的 ID。由于良好的哈希函数可能会为唯一的输入术语产生唯一的输出，这个技巧可能会很好地工作。虽然会有一些哈希冲突（即两个词哈希到相同的 ID），但这应该只是整体中的一个小比例。

注意

哈希特征本身以创建唯一 ID 的方法通常被称为“哈希技巧”。尽管它通常用于文本向量化，但它可以应用于任何涉及大量特征的问题。

然而，生成跨整个集群唯一的区分 ID 的挑战，给了我们一个讨论 Hadoop 中 Parkour 所揭示的一个有用功能的机会：分布式缓存。

使用 Hadoop 分布式唯一 ID

如果我们要计算唯一的集群范围 ID，考虑一下我们的 Parkour 映射器和化简器可能是什么样子。映射器很简单：我们希望计算每个遇到的术语的文档频率，因此以下映射器简单地为每个唯一术语返回一个向量：向量的第一个元素（键）是术语本身，第二个元素（值）是 1。

(defn document-count-m
  {::mr/source-as :vals}
  [documents]
  (->> documents
       (r/mapcat (comp distinct stemmer/stems))
       (r/map #(vector % 1))))

化简器的任务是将这些术语的键/值对（术语与文档计数）进行归约，以确保每个唯一术语都有一个唯一的 ID。做这件事的一个简单方法是确保集群上只有一个化简器。由于所有术语都会传递给这个单一的进程，化简器可以简单地保持一个内部计数器，并像我们之前用 Clojure 原子做的那样为每个术语分配一个 ID。然而，这并没有利用 Hadoop 的分布式能力。

我们尚未介绍 Parkour 的一个特点，即每个映射器（mapper）和归约器（reducer）内部可以访问的运行时上下文。Parkour 将 parkour.mapreduce/*context* 动态变量绑定到执行我们映射器和归约器的 Hadoop 任务的任务上下文中。任务上下文包含以下属性（其中之一）：

属性	类型	描述
`mapred.job.id`	字符串	作业的 ID
`mapred.task.id`	int	任务尝试的 ID
`mapred.task.partition`	int	任务在作业中的 ID

其中最后一个属性 mapred.task.partition，是 Hadoop 分配的任务编号，保证是一个单调递增的唯一整数，且在整个集群中唯一。这个数字就是我们任务的全局偏移量。在每个任务内部，我们还可以保持一个本地偏移量，并在处理每个单词时输出全局和本地偏移量。全局偏移量和本地偏移量一起为集群中的每个术语提供了一个唯一的标识符。

以下图表展示了在三个独立的映射器上处理的八个术语的过程：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_170.jpg

每个映射器只知道它自己的分区号和术语的本地偏移量。然而，这两个数字就是计算唯一全局 ID 所需的全部信息。前面提到的 计算偏移量 盒子确定了每个任务分区的全局偏移量是什么。分区 1 的全局偏移量为 0。分区 2 的全局偏移量为 3，因为分区 1 处理了 3 个单词。分区 3 的偏移量为 5，因为分区 1 和 2 处理了总共 5 个单词，依此类推。

要使上述方法工作，我们需要知道三件事：映射器的全局偏移量、术语的本地偏移量，以及每个映射器处理的术语总数。这三组数字可以用来为每个术语定义一个全局唯一的集群 ID。生成这三组数字的归约器定义如下。它引入了一些新概念，稍后我们将详细讨论。

(defn unique-index-r
  {::mr/source-as :keyvalgroups,
   ::mr/sink-as dux/named-keyvals}
  [coll]
  (let [global-offset (conf/get-long mr/*context*
                                     "mapred.task.partition" -1)]
    (tr/mapcat-state
     (fn [local-offset [word doc-counts]]
       [(inc local-offset)
        (if (identical? ::finished word)
          [[:counts [global-offset local-offset]]]
          [[:data [word [[global-offset local-offset]
                         (apply + doc-counts)]]]])])
     0 (r/mapcat identity [coll [[::finished nil]]]))))

归约器执行的第一步是获取 global-offset，即此归约器对应的任务分区。我们使用 mapcat-state，这是在 transduce 库中定义的一个函数（github.com/brandonbloom/transduce），来构建一系列元组，格式为 [[:data ["apple" [1 4]] [:data ["orange" [1 5]] ...]，其中数字向量 [1 4] 分别表示全局和本地偏移量。最后，当我们到达此归约任务的末尾时，我们会将一个元组以 [:counts [1 5]] 的格式添加到序列中。这代表了该特定归约器分区 1 的最终本地计数 5。因此，单个归约器计算了我们计算所有术语 ID 所需的三项元素。

提供给::mr/source-as的关键字是我们之前没有遇到过的。在上一章中，我们看到了如何通过：keyvals、:keys和:vals等塑形选项，让 Parkour 知道我们希望如何提供数据，以及我们将返回的数据结构。对于聚合器，Parkour 描述了一个更全面的塑形函数集，考虑到输入可能是分组的。下图展示了可用的选项：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_175.jpg

提供给::mr/sink-as的选项我们也没有遇到过。parkour.io.dux命名空间提供了输出去多路复用的选项。实际上，这意味着通过将 sink 指定为dux/named-keyvals，单个聚合器可以写入多个不同的输出。换句话说，我们在数据管道中引入了一个分支：部分数据写入一个分支，其余数据写入另一个分支。

设置了dux/named-keyvals的 sink 规范后，我们元组的第一个元素将被解释为写入的目标；元组的第二个元素将被视为要写入的键值对。因此，我们可以将:data（本地和全局偏移量）写入一个目标，将:counts（每个映射器处理的术语数量）写入另一个目标。

接下来展示的是使用我们定义的映射器和聚合器的作业。与上一章中指定的 Parkour 作业类似，我们将输入、映射、分区、归约和输出步骤串联起来。

(defn df-j [dseq]
  (-> (pg/input dseq)
      (pg/map #'document-count-m)
      (pg/partition (mra/shuffle [:string :long]))
      (pg/reduce #'unique-index-r)
      (pg/output :data (mra/dsink [:string index-value])
                 :counts (mra/dsink [:long :long]))))

前面的代码和我们之前看到的作业规范之间有两个主要区别。首先，我们的输出指定了两个命名的 sink：每个聚合器的输出各一个。其次，我们使用parkour.io.avro命名空间作为mra来为我们的数据指定模式，使用(mra/dsink [:string long-pair])。

在上一章中，我们使用了 Tesser 的FressianWritable将任意 Clojure 数据结构序列化到磁盘。这是因为FressianWritable的内容不需要 Hadoop 解析：该值是完全不透明的。使用 Parkour 时，我们可以定义自定义的键/值对类型。由于 Hadoop 需要将键和值作为独立实体进行解析（用于读取、分区和写入序列文件），Parkour 允许我们使用parkour.io.avro命名空间定义“元组模式”，该模式明确地定义了键和值的类型。long-pair是一个自定义模式，用于将本地和全局偏移量存储在单个元组中。

(def long-pair (avro/tuple-schema [:long :long]))

由于模式是可组合的，我们可以在定义输出模式时引用long-pair模式：(mra/dsink [:string long-pair])。

注意

Parkour 使用Acbracad库来通过 Avro 序列化 Clojure 数据结构。有关序列化选项的更多信息，请参考 Abracad 文档，网址：github.com/damballa/abracad。

让我们来看看 Parkour 暴露的 Hadoop 的另一个功能，它使得我们的术语 ID 任务比其他方式更高效：分布式缓存。

使用分布式缓存共享数据

如我们在上一节中讨论的那样，如果我们知道每个映射器中每个单词的本地偏移量，并且我们知道每个映射器处理了多少记录，那么我们就能计算出每个单词的唯一连续 ID。

几页前展示的图示包含了两个中央框，每个框上标有计算偏移量和全局 ID。这些框直接对应于我们接下来要展示的函数：

(defn global-id [offsets [global-offset local-offset]]
  (+ local-offset (get offsets global-offset)))

(defn calculate-offsets [dseq]
  (->> (into [] dseq)
       (sort-by first)
       (reductions (fn [[_ t] [i n]]
                     [(inc i) (+ t n)])
                   [0 0])
       (into {})))

一旦我们计算出了用于生成唯一 ID 的偏移量映射，我们希望这些映射能够作为共享资源对所有的映射任务和归约任务可用。既然我们已经以分布式的方式生成了偏移量，我们也希望以分布式的方式进行使用。

分布式缓存是 Hadoop 让任务能够访问公共数据的一种方式。这比通过可能昂贵的数据连接共享少量数据（足够小可以存储在内存中的数据）要高效得多。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_180.jpg

在从分布式缓存读取数据之前，我们需要向它写入一些数据。这可以通过 Parkour 的parkour.io.dval命名空间来实现：

(defn unique-word-ids [conf df-data df-counts]
  (let [offsets-dval (-> (calculate-offsets df-counts)
                         (dval/edn-dval))]
    (-> (pg/input df-data)
        (pg/map #'word-id-m offsets-dval)
        (pg/output (mra/dsink [word-id]))
        (pg/fexecute conf `word-id)
        (->> (r/map parse-idf)
             (into {}))
        (dval/edn-dval))))

在这里，我们使用dval/edn-dval函数将两组数据写入分布式缓存。第一组数据是刚刚定义的calculate-offsets函数的结果，它将传递给word-id-m映射器使用。写入分布式缓存的第二组数据是它们的输出。我们将看到如何在word-id-m函数中生成这些数据，如下所示：

(defn word-id-m
  {::mr/sink-as :keys}
  [offsets-dval coll]
  (let [offsets @offsets-dval]
    (r/map
     (fn [[word word-offset]]
       [word (global-id offsets word-offset)])
     coll)))

dval/edn-dval返回的值实现了IDRef接口。这意味着我们可以使用 Clojure 的deref函数（或者@反引用宏）来获取它所包装的值，就像我们操作 Clojure 的原子值一样。第一次对分布式值进行反引用时，会从分布式缓存中下载数据到本地映射器缓存中。一旦数据在本地可用，Parkour 会负责重新构建我们以 EDN 格式写入的数据结构（偏移量的映射）。

从输入文档构建 Mahout 向量

在前面的章节中，我们绕了一些路介绍了几个新的 Parkour 和 Hadoop 概念，但现在我们终于可以为 Mahout 构建文本向量，并为每个术语使用唯一的 ID。为了简洁起见，部分代码被省略，但整个任务可以在cljds.ch6.vectorizer示例代码命名空间中查看。

如前所述，Mahout 的k-means 实现要求我们使用其向量类之一提供输入的向量表示。由于我们的字典很大，并且大多数文档只使用其中少数术语，因此我们将使用稀疏向量表示。以下代码利用了一个dictionary分布式值，为每个输入文档创建一个org.apache.mahout.math.RandomAccessSparseVector：

(defn create-sparse-tfidf-vector [dictionary [id doc]]
  (let [vector (RandomAccessSparseVector. (count dictionary))]
    (doseq [[term tf] (-> doc stemmer/stems frequencies)]
      (let [term-info (get dictionary term)
            id  (:id term-info)
            idf (:idf term-info)]
        (.setQuick vector id (* tf idf))))
    [id vector]))

(defn create-tfidf-vectors-m [dictionary coll]
  (let [dictionary @dictionary]
    (r/map #(create-sparse-tfidf-vector dictionary %) coll)))

最后，我们利用create-tfidf-vectors-m函数，它将我们所学的内容汇聚成一个单一的 Hadoop 作业：

(defn tfidf [conf dseq dictionary-path vector-path]
  (let [doc-count (->> dseq (into []) count)
        [df-data df-counts] (pg/execute (df-j dseq) conf df)
        dictionary-dval (make-dictionary conf df-data
                                         df-counts doc-count)]
    (write-dictionary dictionary-path dictionary-dval)
    (-> (pg/input dseq)
        (pg/map #'create-tfidf-vectors-m dictionary-dval)
        (pg/output (seqf/dsink [Text VectorWritable] vector-path))
        (pg/fexecute conf `vectorize))))

这个任务处理字典的创建，将字典写入分布式缓存，然后使用我们刚定义的映射器，将每个输入文档转换为 Mahout 向量。为了确保与 Mahout 的序列文件兼容，我们将最终输出的键/值类设置为Text和VectorWritable，其中键是文档的原始文件名，值是文档内容的 Mahout 向量表示。

我们可以通过运行以下命令来调用此作业：

(defn ex-6-18 []
  (let [input-path  "data/reuters-sequencefile" 
        output-path "data/reuters-vectors"]
    (vectorizer/tfidf-job (conf/ig) input-path output-path)))

该作业将字典写入dictionary-path（我们稍后还需要它），并将向量写入vector-path。

提示

由于运行前面的例子是后续示例的先决条件，它也可以通过命令行访问：

lein create-vectors

接下来，我们将学习如何使用这些向量来实际执行 Mahout 的聚类。

使用 Mahout 运行 k-means 聚类

现在，我们已经有了一个适合 Mahout 使用的向量序列文件，接下来就该在整个数据集上实际运行k-means 聚类了。与我们本地的 Incanter 版本不同，Mahout 在处理完整的 Reuters 语料库时不会遇到任何问题。

与SequenceFilesFromDirectory类一样，我们已经为 Mahout 的另一个命令行程序KMeansDriver创建了一个封装器。Clojure 变量名使得我们更容易理解每个命令行参数的作用。

(defn run-kmeans [in-path clusters-path out-path k]
  (let [distance-measure  "org.apache.mahout.common.distance.CosineDistanceMeasure"
        max-iterations    100
        convergence-delta 0.001]
    (KMeansDriver/main
     (->> (vector "-i"  in-path
                  "-c"  clusters-path
                  "-o"  out-path
                  "-dm" distance-measure
                  "-x"  max-iterations
                  "-k"  k
                  "-cd" convergence-delta
                  "-ow"
                  "-cl")
          (map str)
          (into-array String)))))

我们提供字符串org.apache.mahout.common.distance.CosineDistanceMeasure，以指示驱动程序我们希望使用 Mahout 的余弦距离度量实现。Mahout 还包括EuclideanDistanceMeasure和TanimotoDistanceMeasure（类似于 Jaccard 距离，是 Jaccard 指数的补集，但将作用于向量而非集合）。还有几种其他的距离度量可以选择；请参考 Mahout 文档以了解所有可用选项。

有了前面的run-kmeans函数后，我们只需要告诉 Mahout 在哪里访问我们的文件。与上一章一样，我们假设 Hadoop 在本地模式下运行，所有文件路径都相对于项目根目录：

(defn ex-6-19 []
  (run-kmeans "data/reuters-vectors/vectors"
              "data/kmeans-clusters/clusters"
              "data/kmeans-clusters"
              10))

这个例子可能会运行一段时间，因为 Mahout 需要对我们的大数据集进行迭代。

查看 k-means 聚类结果

完成后，我们希望能看到每个聚类的聚类总结，就像我们在 Incanter 实现中所做的那样。幸运的是，Mahout 定义了一个ClusterDumper类，正是用来做这个事情的。我们需要提供聚类的位置，当然，还需要提供字典的位置。提供字典意味着输出将返回每个聚类的顶部术语。

(defn run-cluster-dump [in-path dict-path points-dir out-path]
  (let [distance-measure
        "org.apache.mahout.common.distance.CosineDistanceMeasure"]
    (ClusterDumper/main
     (->> (vector "-i" in-path
                  "-o" out-path
                  "-d" dict-path
                  "--pointsDir" points-dir
                  "-dm" distance-measure
                  "-dt" "sequencefile"
                  "-b" "100"
                  "-n" "20"
                  "-sp" "0"
                  "--evaluate")
          (map str)
          (into-array String)))))

接下来，我们定义实际调用run-cluster-dump函数的代码：

(defn path-for [path]
  (-> (fs/glob path)
      (first)
      (.getAbsolutePath)))

(defn ex-6-20 []
  (run-cluster-dump
   (path-for "data/kmeans-clusters/clusters-*-final")
   "data/reuters-vectors/dictionary/part-r-00000"
   "data/kmeans-clusters/clusteredPoints"
   "data/kmeans-clusterdump"))

我们再次使用me.raynes.fs库来确定最终聚类所在的目录。Mahout 会在包含最终聚类的目录名后附加-final，但我们事先并不知道哪个目录会是这个最终目录。fs/glob函数会查找与clusters-*-final模式匹配的目录，并将*替换为实际目录名称中包含的迭代编号。

解释聚类输出

如果你在任何文本编辑器中打开之前示例创建的文件data/kmeans-clusterdump，你将看到表示 Mahout 聚类的顶部术语的输出。该文件可能很大，但下面提供了一个摘录：

:VL-11417{n=312 c=0.01:0.039, 0.02:0.030, 0.07:0.047, 0.1:0.037, 0.10:0.078, 0.11:0.152, 0.12:0.069,
  Top Terms:
    tonnes              =>   2.357810452962533
    department          =>   1.873890568048526
    wheat               =>  1.7797807546762319
    87                  =>  1.6685682321206117
    u.s                 =>   1.634764205186795
    mln                 =>  1.5050923755535712
    agriculture         =>  1.4595903158187866
    ccc                 =>  1.4314624499051998
    usda                =>  1.4069041441648433
    dlrs                =>  1.2770121846443567

第一行包含关于聚类的信息：ID（在本例中为VL-11417），后面跟着包含聚类大小和聚类质心位置的大括号。由于文本已转换为权重和数字 ID，单独解读质心是不可行的。不过，质心描述下方的顶部术语暗示了聚类的内容；它们是聚类汇聚的核心术语。

VL-12535{n=514 c=[0:0.292, 0.25:0.015, 0.5:0.012, 00:0.137, 00.46:0.018, 00.50:0.036, 00.91:0.018, 0
  Top Terms:
    president           =>   3.330068911559851
    reagan              =>   2.485271333256584
    chief               =>  2.1148699971952327
    senate              =>   1.876725117983985
    officer             =>  1.8531712558019022
    executive           =>  1.7373591731030653
    bill                =>  1.6326750159727461
    chairman            =>  1.6280977206471365
    said                =>  1.6279512813119108
    house               =>  1.5771017798189988

前面提到的两个聚类暗示了数据集中存在的两个明确主题，尽管由于k-means 算法的随机性，你的聚类可能会有所不同。

根据初始质心和算法运行的迭代次数，你可能会看到某些聚类在某些方面看起来“更好”或“更差”。这将基于对聚类词汇如何组合在一起的直觉反应。但通常，仅通过查看顶部术语并不能清楚地判断聚类的效果如何。无论如何，直觉并不是判断无监督学习算法质量的可靠方法。我们理想的情况是有一个定量指标来衡量聚类的效果。

聚类评估指标

在我们上一节查看的文件底部，你会看到一些统计信息，表明数据的聚类效果如何：

Inter-Cluster Density: 0.6135607681542804
Intra-Cluster Density: 0.6957348405534836

这两个数字可以看作是我们在[第二章，推理和第三章，相关性中看到的组内方差和组间方差的等价物。理想情况下，我们希望聚类内的方差较低（或密度较高），而聚类间的密度较低。

聚类间密度

聚类间密度是聚类质心之间的平均距离。好的聚类通常不会有太靠近的质心。如果它们太近，那就意味着聚类正在创建具有相似特征的组，并可能在区分聚类成员方面存在难以支持的情况。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_190.jpg

因此，理想情况下，我们的聚类应产生具有大聚类间距离的聚类。

聚类内部密度

相反，聚类内密度是衡量聚类紧凑性的指标。理想情况下，聚类会识别出彼此相似的项目组。紧凑的聚类表示聚类中的所有项目彼此高度相似。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_200.jpg

因此，最佳的聚类结果会产生紧凑、独特的聚类，具有高聚类内密度和低聚类间密度。

然而，数据到底应该有多少个聚类并不总是清楚。考虑以下示例，它展示了同一数据集以不同聚类数分组的情况。很难有足够的信心判断理想的聚类数是多少。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_210.jpg

尽管前面的示意图是人为构造的，但它展示了聚类数据时的一个普遍问题。通常没有一个明确的“最佳”聚类数。最有效的聚类将很大程度上取决于数据的最终用途。

然而，我们可以通过确定某些质量评分如何随着聚类数量的变化来推断可能的较优 k 值。质量评分可以是像聚类间密度或聚类内密度这样的统计数据。当聚类数接近理想值时，我们期望该质量评分的值会提高。相反，当聚类数偏离理想值时，我们期望质量评分会下降。因此，为了合理地估算数据集中多少个聚类是合理的，我们应该为不同的 k 值多次运行算法。

使用 Parkour 计算均方根误差

最常见的聚类质量度量之一是平方误差和（SSE）。对于每个点，误差是测量到最近聚类质心的距离。因此，总的聚类 SSE 是聚类点到其相应质心的所有聚类的和：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_11.jpg

其中 µ[i] 是聚类 S[i] 中点的质心，k 是聚类的总数，n 是点的总数。

因此，在 Clojure 中计算 RMSE 时，我们需要能够将聚类中的每个点与其对应的聚类质心关联起来。Mahout 将聚类质心和聚类点保存在两个不同的文件中，因此在下一节中我们将把它们合并。

加载聚类点和质心

给定一个父目录（例如 data/reuters-kmeans/kmeans-10），以下函数将使用 Parkour 的 seqf/dseq 函数从序列文件加载键/值对，并将点加载到按聚类 ID 索引的向量映射中。在这种情况下，键是聚类 ID（整数），值是 TF-IDF 向量。

(defn load-cluster-points [dir]
  (->> (points-path dir)
       (seqf/dseq)
       (r/reduce
        (fn [accum [k v]]
          (update-in accum [k] conj v)) {})))

上述函数的输出是一个按聚类 ID 索引的映射，其值是聚类点的序列。同样，以下函数将把每个聚类转换为一个按聚类 ID 索引的映射，其值是包含 :id 和 :centroid 键的映射。

(defn load-cluster-centroids [dir]
  (let [to-tuple (fn [^Cluster kluster]
                   (let [id (.getId kluster)]
                     [id  {:id id
                           :centroid (.getCenter kluster)}]))]
    (->> (centroids-path dir)
         (seqf/dseq)
         (r/map (comp to-tuple last))
         (into {}))))

拥有两个按聚类 ID 索引的映射意味着将聚类的点与聚类中心结合起来，只需对映射调用 merge-with 并提供一个自定义的合并函数即可。在以下代码中，我们将聚类的点合并到包含聚类 :id 和 :centroid 的映射中。

(defn assoc-points [cluster points]
  (assoc cluster :points points))

(defn load-clusters [dir]
  (->> (load-cluster-points dir)
       (merge-with assoc-points
                   (load-cluster-centroids dir))
       (vals)))

最终输出是一个单一的映射，按聚类 ID 索引，其中每个值都是一个包含 :id、:centroid 和 :points 的映射。我们将在下一节中使用这个映射来计算聚类的 RMSE。

计算聚类 RMSE

为了计算 RMSE，我们需要能够确定每个点与其关联的聚类中心之间的距离。由于我们使用了 Mahout 的 CosineDistanceMeasure 来执行初始聚类，因此我们也应该使用余弦距离来评估聚类。事实上，我们可以直接利用 Mahout 的实现。

(def measure
  (CosineDistanceMeasure.))

(defn distance [^DistanceMeasure measure a b]
  (.distance measure a b))

(defn centroid-distances [cluster]
  (let [centroid (:centroid cluster)]
    (->> (:points cluster)
         (map #(distance measure centroid %)))))

(defn squared-errors [cluster]
  (->> (centroid-distances cluster)
       (map i/sq)))

(defn root-mean-square-error [clusters]
  (->> (mapcat squared-errors clusters)
       (s/mean)
       (i/sqrt)))

如果将 RMSE 与聚类数绘制成图，你会发现随着聚类数的增加，RMSE 会逐渐下降。单一聚类将具有最高的 RMSE 错误（原始数据集与均值的方差），而最低的 RMSE 将是每个点都在自己的聚类中的退化情况（RMSE 为零）。显然，这两个极端都无法很好地解释数据的结构。然而，RMSE 不是线性下降的。当聚类数从 1 增加时，它会急剧下降，但一旦超过“自然”聚类数后，下降的速度就会变慢。

因此，判断理想聚类数的一种方法是绘制 RMSE 随聚类数变化的图表。这就是所谓的 肘部法。

使用肘部法确定最佳的 k 值

为了使用肘部法确定 k 的值，我们需要多次重新运行 k-均值聚类。以下代码实现了对 2 到 21 之间的所有 k 值进行聚类。

(defn ex-6-21 []
  (doseq [k (range 2 21)
          :let [dir (str "data/kmeans-clusters-" k)]]
    (println dir)
    (run-kmeans "data/reuters-vectors/vectors"
                (str dir "/clusters")
                dir k)))

这可能需要一些时间，所以不妨去泡杯热饮：println 语句会记录每次聚类的运行情况，让你知道进展如何。在我的笔记本电脑上，整个过程大约需要 15 分钟。

完成后，你应该能够运行示例，生成每个聚类值的 RMSE 散点图：

(defn ex-6-22 []
  (let [ks (range 2 21)
        ys (for [k ks
                 :let [dir (str "data/kmeans-clusters-" k)
                       clusters (load-clusters dir)]]
             (root-mean-square-error clusters))]
    (-> (c/scatter-plot ks ys
                        :x-label "k"
                        :y-label "RMSE")
        (i/view))))

这应该返回一个类似于以下的图表：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_220.jpg

上面的散点图展示了 RMSE 与簇数之间的关系。应该能看出，当k超过大约 13 个簇时，RMSE 的变化速率放慢，进一步增加簇数会带来递减的收益。因此，上图表明对于我们的路透社数据，约 13 个簇是一个不错的选择。

肘部法则提供了一种直观的方式来确定理想的簇数，但在实践中有时难以应用。这是因为我们必须解释每个k对应的 RMSE 曲线的形状。如果k很小，或者 RMSE 包含较多噪声，可能不容易看出肘部在哪里，甚至是否存在肘部。

注意

由于聚类是一种无监督学习算法，我们在此假设簇的内部结构是验证聚类质量的唯一手段。如果已知真实的簇标签，则可以使用外部验证度量（例如熵）来验证模型的成功，这在第四章，分类中我们曾遇到过。

其他聚类评估方案旨在提供更清晰的方式来确定精确的簇数。我们将讨论的两种方法是邓恩指数和戴维斯-鲍尔丁指数。它们都是内部评估方案，意味着它们只关注聚类数据的结构。每种方法都旨在以不同的方式识别产生最紧凑、最分离的簇的聚类。

使用邓恩指数确定最优 k

邓恩指数提供了另一种选择最优k的方式。邓恩指数不考虑聚类数据中的平均误差，而是考虑两种“最坏情况”的比率：两个簇中心之间的最小距离，除以最大簇直径。因此，较高的邓恩指数表示更好的聚类，因为通常我们希望簇间距离较大，而簇内距离较小。

对于k个簇，我们可以通过以下方式表达邓恩指数：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_12.jpg

其中，δ(C[i],C[j])是两个簇C[i]和C[j]之间的距离，而https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_13.jpg表示最大簇的大小（或散布）。

计算一个簇的散布有几种可能的方法。我们可以计算簇内最远两个点之间的距离，或者计算簇内所有数据点之间成对距离的平均值，或者计算每个数据点到簇中心的距离的平均值。在以下代码中，我们通过计算从簇中心到各点的中位距离来确定大小。

(defn cluster-size [cluster]
  (-> cluster
      centroid-distances
      s/median))

(defn dunn-index [clusters]
  (let [min-separation (->> (combinations clusters 2)
                            (map #(apply separation %))
                            (apply min))
        max-cluster-size (->> (map cluster-size clusters)
                              (apply max))]
    (/ min-separation max-cluster-size)))

上面的代码使用了 clojure.math.combinatorics 中的 combinations 函数（github.com/clojure/math.combinatorics/）来生成所有聚类的懒序列对。

(defn ex-6-23 []
  (let [ks (range 2 21)
        ys (for [k ks
                 :let [dir (str "data/kmeans-clusters-" k)
                       clusters (load-clusters dir)]]
             (dunn-index clusters))]
    (-> (c/scatter-plot ks ys
                        :x-label "k"
                        :y-label "Dunn Index")
        (i/view))))

我们在上面的代码中使用了 dunn-index 函数来为 k=2 到 k=20 的聚类生成散点图：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_230.jpg

更高的 Dunn 指数表示更好的聚类效果。因此，最佳聚类似乎是 k=2，接下来是 k=6，然后是 k=12 和 k=13，它们紧随其后。我们来尝试一种替代的聚类评估方案，并比较结果。

使用 Davies-Bouldin 指数确定最优的 k 值

Davies-Bouldin 指数是一种替代的评估方案，用于衡量聚类中所有值的大小和分离度的平均比率。对于每个聚类，找到一个替代聚类，使得聚类大小之和与聚类间距离的比率最大化。Davies-Bouldin 指数被定义为所有聚类的此值的平均值：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_14.jpg https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_15.jpg

其中 δ(C[i],C[j]) 是两个聚类中心 C[i] 和 C[j]* 之间的距离，S[i] 和 S[j]* 是散布度。我们可以使用以下代码计算 Davies-Bouldin 指数：

(defn scatter [cluster]
  (-> (centroid-distances cluster)
      (s/mean)))

(defn assoc-scatter [cluster]
  (assoc cluster :scatter (scatter cluster)))

(defn separation [a b]
  (distance measure (:centroid a) (:centroid b)))

(defn davies-bouldin-ratio [a b]
  (/ (+ (:scatter a)
        (:scatter b))
     (separation a b)))

(defn max-davies-bouldin-ratio [[cluster & clusters]]
  (->> (map #(davies-bouldin-ratio cluster %) clusters)
       (apply max)))

(defn rotations [xs]
  (take (count xs)
        (partition (count xs) 1 (cycle xs))))

(defn davies-bouldin-index [clusters]
  (let [ds (->> (map assoc-scatter clusters)
                (rotations)
                (map max-davies-bouldin-ratio))]
    (s/mean ds)))

现在，让我们在散点图中绘制 k=2 到 k=20 的 Davies-Bouldin 指数：

(defn ex-6-24 []
  (let [ks (range 2 21)
        ys (for [k ks
                 :let [dir (str "data/kmeans-clusters-" k)
                       clusters (load-clusters dir)]]
             (davies-bouldin-index clusters))]
    (-> (c/scatter-plot ks ys
                        :x-label "k"
                        :y-label "Davies-Bouldin Index")
        (i/view))))

这将生成以下图表：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_240.jpg

与 Dunn 指数不同，Davies-Bouldin 指数对于良好的聚类方案是最小化的，因为一般来说，我们寻找的是紧凑的聚类，并且聚类间的距离较大。前面的图表表明，k=2 是理想的聚类大小，其次是 k=13。

k-均值算法的缺点

k-均值算法是最流行的聚类算法之一，因为它相对易于实现，并且可以很好地扩展到非常大的数据集。尽管它很受欢迎，但也有一些缺点。

k-均值算法是随机的，并不能保证找到全局最优的聚类解。事实上，该算法对离群点和噪声数据非常敏感：最终聚类的质量可能高度依赖于初始聚类中心的位置。换句话说，k-均值算法通常会发现局部最优解，而非全局最优解。

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_250.jpg

上图说明了如何根据不良的初始聚类质心，k均值可能会收敛到局部最小值。如果初始聚类质心的位置合适，非最优聚类仍然可能发生，因为 k 均值倾向于选择具有相似大小和密度的聚类。在聚类大小和密度不大致相等时，k 均值可能无法收敛到最自然的聚类：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_260.jpg

此外，k 均值算法强烈倾向于选择“球形”的聚类。形状较为复杂的聚类往往不能被 k 均值算法很好地识别。

在下一章，我们将看到各种降维技术如何帮助解决这些问题。但在那之前，我们先来培养一种直觉，理解另一种定义距离的方式：作为衡量某个元素距离一“组”物体的远近。

马哈拉诺比斯距离度量

在本章开始时，我们通过展示 Jaccard、欧几里得和余弦距离如何与数据表示相关，看到了一些距离度量在特定数据下可能比其他度量更合适。选择距离度量和聚类算法时需要考虑的另一个因素是数据的内部结构。请考虑以下散点图：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_270.jpg

很“显然”，箭头所指的点与其他点不同。我们可以清楚地看到它远离其他点的分布，因此代表了一个异常点。然而，如果我们计算所有点到均值（“质心”）的欧几里得距离，这个点将被其他距离相等或更远的点所掩盖：

(defn ex-6-25 []
  (let [data (dataset-with-outlier)
        centroid  (i/matrix [[0 0]])
        distances (map #(s/euclidean-distance centroid %) data)]
    (-> (c/bar-chart (range 202) distances
                     :x-label "Points"
                     :y-label "Euclidean Distance") 
        (i/view))))

上述代码生成了以下图表：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_280.jpg

马哈拉诺比斯距离在计算距离时考虑了变量之间的协方差。在二维空间中，我们可以将欧几里得距离想象成一个从质心发出的圆：圆上所有的点与质心的距离相等。马哈拉诺比斯距离将这个圆拉伸并扭曲，以纠正不同变量的尺度差异，并考虑它们之间的相关性。我们可以在以下示例中看到这种效果：

(defn ex-6-26 []
  (let [data (dataset-with-outlier)
        distances    (map first (s/mahalanobis-distance data))]
    (-> (c/bar-chart (range 202) distances
                     :x-label "Points"
                     :y-label "Mahalanobis Distance")
        (i/view))))

上述代码使用了 incanter.stats 提供的函数来绘制相同数据点的马哈拉诺比斯距离。结果显示在下图中：

https://github.com/OpenDocCN/freelearn-ds-pt5-zh/raw/master/docs/clj-ds/img/7180OS_06_290.jpg

该图清晰地标识出一个点与其他点相比，距离要远得多。这与我们对这个点应该被认为比其他点更远的看法一致。

维度灾难

然而，马哈拉诺比斯距离度量无法克服一个事实，这就是所谓的维度灾难。随着数据集中的维度数量增加，每个点都趋向于和其他点一样远。我们可以通过下面的代码简单地演示这一点：

(defn ex-6-27 []
  (let [distances (for [d (range 2 100)
                        :let [data (->> (dataset-of-dimension d)
                                        (s/mahalanobis-distance)
                                        (map first))]]
                    [(apply min data) (apply max data)])]
    (-> (c/xy-plot (range 2 101) (map first distances)
                   :x-label "Number of Dimensions"
                   :y-label "Distance Between Points"
                   :series-label "Minimum Distance"
                   :legend true)
        (c/add-lines (range 2 101) (map second distances)
                     :series-label "Maximum Distance")
        (i/view))))