NumPy 作为 Python 生态中最受欢迎的科学计算包，很多读者已经非常熟悉它了。
它为 Python 提供高效率的多维数组计算，并提供了一系列高等数学函数，我们可以快速搭建模型的整个计算流程。毫不负责任地说，NumPy 就是现代深度学习框架的「爸爸」。
尽管目前使用 NumPy 写模型已经不是主流，但这种方式依然不失为是理解底层架构和深度学习原理的好方法。最近，来自普林斯顿的一位博士后将 NumPy 实现的所有机器学习模型全部开源，超过 3 万行代码、30 多个模型，并提供了相应的论文和一些实现的测试效果。
粗略估计，该项目大约有 30 个主要机器学习模型，此外还有 15 个用于预处理和计算的小工具，全部.py 文件数量有 62 个之多。平均每个模型的代码行数在 500 行以上，在神经网络模型的 layer.py 文件中，代码行数接近 4000。
这，应该是目前用 NumPy 手写机器学习模型的「最高境界」吧。
项目地址为：
https://github.com/ddbourgin/numpy-ml
通过项目的代码目录，我们能发现，作者基本上把主流模型都实现了一遍，这个工作量简直惊为天人。
作者 David Bourgin 是一位大神，于 2018 年获得加州大学伯克利分校计算认知科学博士学位，随后在普林斯顿大学从事博士后研究。
尽管毕业不久，David 在顶级期刊与计算机会议上都发表了一些优秀论文。在 ICML 2019 中，其关于认知模型先验的研究就被接收为少有的 Oral 论文。
David Bourgin 就是用 NumPy 手写 ML 模型、手推反向传播的大神。这么多的工作量，当然还是需要很多参考资源的，David 会理解这些资源或实现，并以一种更易读的方式写出来。
他表示，从 autograd repo 学到了很多，但二者的不同之处在于，他显式地进行了所有梯度计算，以突出概念/数学的清晰性。当然，这么做的缺点也很明显，在每次需要微分一个新函数时，你都要写出它的公式……
估计 David Bourgin 在写完这个项目后，机器学习基础已经极其牢固了。
这个项目最大的特点是作者把机器学习模型都用 NumPy 手写了一遍，包括更显式的梯度计算和反向传播过程。可以说它就是一个机器学习框架了，只不过代码可读性会强很多。
David Bourgin 表示他一直在慢慢写或收集不同模型与模块的纯 NumPy 实现，它们跑起来可能没那么快，但是模型的具体过程一定足够直观。每当我们想了解模型 API 背后的实现，却又不想看复杂的框架代码，那么它可以作为快速的参考。
如下所示为项目文件，不同的文件夹即不同种类的代码集：
在每一个代码集下，作者都会提供不同实现的参考资料，例如模型的效果示例图、参考论文和参考链接等。
当然如此庞大的代码总会存在一些 Bug，作者也非常希望我们能一起完善这些实现。如果我们以前用纯 NumPy 实现过某些好玩的模型，那也可以直接提交 PR 请求。因为实现基本上都只依赖于 NumPy，那么环境配置就简单很多了，大家差不多都能跑得动。
作者在 GitHub 中提供了模型/模块的实现列表，列表结构基本就是代码文件的结构了。整体上，模型主要分为两部分，即传统机器学习模型与主流的深度学习模型。
其中浅层模型既有隐马尔可夫模型和提升方法这样的复杂模型，也包含了线性回归或最近邻等经典方法。而深度模型则主要从各种模块、层级、损失函数、最优化器等角度搭建代码架构，从而能快速构建各种神经网络。
除了模型外，整个项目还有一些辅助模块，包括一堆预处理相关的组件和有用的小工具。
该 repo 的模型或代码结构如下所示：
由于代码量庞大，这里整理了一些示例。
例如，实现点积注意力机制：
class DotProductAttention(LayerBase):
def __init__(self, scale=True, dropout_p=0, init="glorot_uniform", optimizer=None):
super().__init__(optimizer)
self.init = init
self.scale = scale
self.dropout_p = dropout_p
self.optimizer = self.optimizer
self._init_params()
def _fwd(self, Q, K, V):
scale = 1 / np.sqrt(Q.shape[-1]) if self.scale else 1
scores = Q @ K.swapaxes(-2, -1) * scale # attention scores
weights = self.softmax.forward(scores) # attention weights
Y = weights @ V
return Y, weights
def _bwd(self, dy, q, k, v, weights):
d_k = k.shape[-1]
scale = 1 / np.sqrt(d_k) if self.scale else 1
dV = weights.swapaxes(-2, -1) @ dy
dWeights = dy @ v.swapaxes(-2, -1)
dScores = self.softmax.backward(dWeights)
dQ = dScores @ k * scale
dK = dScores.swapaxes(-2, -1) @ q * scale
return dQ, dK, dV
在以上代码中，Q、K、V 三个向量输入到「_fwd」函数中，用于计算每个向量的注意力分数，并通过 softmax 的方式得到权重。而「_bwd」函数则计算 V、注意力权重、注意力分数、Q 和 K 的梯度，用于更新网络权重。
在一些实现中，作者也进行了测试，并给出了测试结果。如图为隐狄利克雷（Latent Dirichlet allocation，LDA）实现进行文本聚类的结果。左图为词语在特定主题中的分布热力图。右图则为文档在特定主题中的分布热力图。
