对 APL 等数组语言的有原则性重新思考

# 对 APL 类数组语言的 principled 重新思考 **Dercuano 来源**：<https://dercuano.github.io/notes/principled-apl.html> Kragen Javier Sitaker，2015-05-16（更新于 2019-09-30）（阅读时间约 31 分钟） Matlab、R 和 APL 等数组处理语言的力量和便利性从何而来？我们能否在不支付这些语言所伴随的高可读性成本和易错性成本的情况下获得这种能力？我觉得我终于对这里发生的事情有了一个解释，而且我们不仅能从中获得更好的错误检查，还能获得一种比传统数组处理语言更具表达力（以“更少的 token 拥有更大的力量”为义）、更易读且更高效的编程语言。虽然我认为本笔记本身是自包含的，但该想法的初步发展见于 2007 或 2008 年左右的《索引集合推断或针对带有索引族编程的域推断》（<https://dercuano.github.io/notes/matrix-multiply.html>）。 ## 基本思想 计算机程序中的变量拥有一个变化的值。有时代码运行时，变量具有一个值；其他时候，它则有另一个值。（在许多语言中，即使在同一运行期间它也可以改变，但这与我在本文探讨的内容无关。）这些值是一些任意输入集的函数，这些输入通常隐含在其上下文中。在图像处理中，亮度值可能随 X、Y、颜色通道和帧号变化；在统计处理中，经过时间的测量值可能随试验次数变化；在渲染中，像素颜色可能随着色算法、相机位置以及输入场景中所有物体的状态而变化。 很难判断它们依赖于什么，事实上，假设变量相对于给定输入是恒定的，而实际上它应当变化，这是一个常见的错误来源。（你可以认为这是缓存困难性的根源。）我们可以将这些输入视为多维空间中的维度，该空间中的每个点代表一个可能的宇宙——或许是一个非常小的、只包含少数变量的可能宇宙，但无论如何它是一个宇宙。数组语言允许我们将此空间具体化为运行时变量，从而编写作用于整个子空间的程序，而不是像 C 或 Java 等 ALGOL 系语言那样采用逐点和一维的方法。 ## 光线追踪中的一些示例 考虑这段来自[我的第一个光线追踪器](http://canonical.org/~kragen/sw/aspmisc/my-very-first-raytracer)的 C 代码： ```c static color pixel_color(world here, int ww, int hh, int xx, int yy) { vec pv = { (double)xx/ww - 0.5, (double)yy/hh - 0.5, 1 }; ray rr = { {0}, normalize(pv) }; return trace(here, rr, 1.0); } ``` 这段代码被多次调用，其中 `world`、`ww` 和 `hh` 变量相同，但 `xx` 和 `yy` 的值不同；但仅凭观察代码无法看出这一点。同样，当调用 `trace` 时，它使用相同的 `here` 值、不同的 `rr` 值以及恒定的重要性值 1.0 进行多次调用。但必须将其写为函数（或至少是一个循环），以允许 `xx` 和 `yy` 的这种灵活性。 从某种意义上说，在这段代码中，`xx` 不代表单个整数，而是代表大量可能的值，甚至可能是无限系列的值；它不是一个标量，而是关于我们正在查看哪个像素的函数。当我们写 `xx` 除以 `ww` 时，我们并不是在写单个浮点数除以单个整数，而是在写人类历史中 `xx` 将转换成的所有浮点数除以人类历史中 `ww` 将包含的所有可能图像宽度。或者，如果我们限制视角到单次执行，那条除法指令最终将被用来将所有水平图像像素坐标除以图像宽度——事实上冗余地执行多次，每行一次；因此这是一条隐含表示向量操作的机器指令。 但这是一种相当不寻常的 C 和机器代码的解释学。C 代码强制特定的求值顺序：在第一个 `trace()` 调用完成之前，它无法开始求值第二个调用，并且在第一个 `pixel_color` 调用完成之前，也无法求值第二个调用。但这可能不是找出像素应该是什么颜色的最高效方式。你会注意到，C 代码在允许这些变量变化方面也有相当多的语法开销；它们必须声明为参数。 如果，相反，我们使用一种编程语言，其中这些变量在源文本中且在运行时内存里明确地代表了整个可能性向量——一种更加原则化的 APL，会怎样？也许我们可以这样写： ```apl pixel_color = trace(here, ray(vec(0,0,0), normalize(vec(xx/ww-.5, yy/hh-.5, 1))), 1.0); ``` （由于无需命名临时结构体，我们也获得了一些简洁性。） 具有类似 APL 求值策略的语言可以推断出 `xx` 和 `yy` 独立变化，而 `ww` 和 `hh` 根本不变，从而为我们正在标准化的 `vec`s 生成一个 $\rho$`xx` × $\rho$`yy` 的可能性空间，其中 $\rho$ 是给出矩阵形状的 APL 算子。（关于如何做到这一点的更多细节将在下一节中给出。） 我认为这是 APL 的相容性和广播规则的最终哲学证明；例如，如果你正在对 640 列和 480 行的像素进行光线追踪，那么对所有那些像素恒定的值只是一个标量；随列变化但不随行变化的值将是一个 640 元素向量；随行变化但不随列变化的值将是一个 480 元素向量；随像素变化的值必然是一个 640×480 矩阵。因此，将 `xx` 除以 `ww` 或将 `yy` 除以 `hh` 是有意义的，但对它们二者一起做任何事情都需要外积操作（在 APL 中是显式的），或者将那些变化轴之一减少为某种虚拟变量，比如求和索引之类。 但当然 APL 无法判断你代表场景中三个球体 X 坐标的 3 元素向量并没有真正与你代表相机 X、Y 和 Z 坐标的三元素向量兼容维度。数组维度类型是一种弱类型，类似于当前流行的将所有内容都键入为字符串的做法。这就是为什么外积在 APL 中必然是显式的，而我认为一种更原则化的数组处理语言可以推断出其中大部分。 这是来自同一程序的另一个例子： ```c static vec scale(vec vv, sc c) { vec rv = { vv.x*c, vv.y*c, vv.z*c }; return rv; } static ray reflect(ray rr, vec start, vec normal) { // 将射线方向投影到法线上 vec proj = scale(normal, dot(rr.dir, normal)); // 减去两倍该值以将射线移动到表面另一侧 vec reflected_dir = sub(rr.dir, scale(proj, 2)); ray rv = { add(start, scale(reflected_dir, 0.001)), reflected_dir }; return rv; } ``` （那里的 0.001 调整因子是为了防止反射射线因舍入误差而从内部再次击中同一表面。） 这里的 `scale` 函数显然只存在于因为 C 不是数组处理语言。或者说真的是这样吗？如果我们尝试编写一个处理多个法线的 `reflect`，使用三个单独的数组来存储法线的 X、Y 和 Z 分量并不完全疯狂。仅取 `reflect` 的第一行，并将其翻译为看似 Fortran 风格的 C 代码： ```c static void scale(sc vx[], sc vy[], sc vz[], sc c[], sc vox[], sc voy[], sc voz[], int n) { for (int ii = 0; ii < n; ii++) { vox[ii] = vx[ii] * c; voy[ii] = vy[ii] * c; voz[ii] = vz[ii] * c; } } static void proj(sc vx[], sc vy[], sc vz[], sc dx, sc dy, sc dz, sc vox[], sc voy[], sc voz[], int n) { sc dots[n]; for (int ii = 0; ii < n; ii++) { dots[ii] = vx[ii] * dx + vy[ii] * dy + vz[ii] * dz; } scale(vx, vy, vz, dots, vox, voy, voz, n); } ``` 你可能不同意这并不完全疯狂，但希望你能同意这是一种 Fortran 程序员认为不疯狂的方法。此外你可以看到，对于这类事情，APL 相对于 Fortran 的巨大好处在于，你至少有一些希望改变主意关于你想要多个射线、法线还是两者都有，因为在 Fortran 中使循环和索引隐式意味着你不必修改代码以索引 `dx` 数组而不索引 `vx` 数组。（此外，对于某些尺寸组合和 CPU 型号，此版本更可预测的内存访问可能实际上使其更快，尽管其计算与访问内存位置的比例较小。显然，如果你的循环是隐式的而非隐式操作受到解释开销的影响，如 Numpy 或正常 APL 实现中，数组方法将大大快得多。） 然而，三维空间中的坐标确实是合理地被看作从 0 到 2 或从 1 到 3 编号的东西，而不是同一事物的不相关属性。那么你可能希望这里的法线数组表示为 $n \times 3$ 数组而不是 $n$ 的三个数组或具有三个字段的结构体数组。然后像 `scale` 这样的函数完全消失，但你需要某种方法来指定在计算点积时要对哪个维度求和。APL `\+/` 默认操作最后维度，并可以通过数字索引指定不同维度，如 `+\[1]`。这在我看来似乎是临时拼凑且不可读的，就像 APL 的许多部分一样。但如果你已经命名了你的变化轴，其中一个区分了 XYZ，那么你可以非常合理地说 `XYZ.sum()` 或 `+\[XYZ]`，这很清晰，将 XYZ 变化转化为虚拟变量；如果将其应用于具有多于一个 XYZ 区分（由显式外积运算符引入）或没有 XYZ 区分的某种聚合，你会得到一个错误。然后你可以写下 ```apl proj = normal * XYZ.sum(raydir * normal) ``` 代替上面那 20 行废话，并且进一步保持抽象，无论你有单个法线和许多射线方向、单个射线方向和许多法线、每个对应于不同射线方向的许多法线，甚至是 APL 永远无法隐含处理的许多射线方向和许多法线，而我们想要其笛卡尔积。然后也许你可以这样写整个函数： ```apl proj = normal * XYZ.sum(raydir * normal) reflected_dir = raydir - 2 * proj reflect = ray(start + reflected_dir * 0.001, reflected_dir) ``` 当涉及到在不同维度上隐式广播操作时，C 语言等同于数组语言一样简洁——除了它为了给你变量所需要的数据键入和抽象开销之外。但由于 C 值只有隐含地、以晦涩的解释学方式成为向量或可能性宇宙，因此难以编写类似于上面的 `XYZ.sum` 函数的内容；相反，我们被迫编写显式循环，或者像这里的情况一样，编写文本重复性很高的代码。 ## 更加严谨：带有隐式参数的功能语义 好吧，“严谨”和“语义”可能是夸大其词。但我至少要在这里概述一种严谨的语义。假设，我们不认为诸如 `normal` 这样的变量持有一些有限大小的标量、向量或矩阵，而是考虑它们持有可计算函数，但其域不一定已知或有限。这不是一个大飞跃：我们可以将向量 `[6 832 4]` 视为三个整数域上的函数 $f$：当作为 $f(0)$ 调用时返回 6，当作为 $f(1)$ 调用时返回 832，或者当作为 $f(2)$ 调用时返回 4。 在这种解释下，我们将通常的算术运算符提升到操作单参数函数的通常方式：例如，`*` 是我们通常在 $\lambda$-演算中写的函数 $\lambda f.\lambda g.\lambda p.(f p)*(g p)$，或者在 Python 中写作 `lambda f, g: lambda p: f(p)*g(p)`；我们认为常数如 `2` 表示常函数，比如 `K 2` 或 `lambda p: 2`。 但是这个神秘的 `p` 参数是什么？它是前面提到的多维可能性空间中的上下文或点，通常是隐含的，所以它需要包括我们之前讨论的所有变化轴；为了得到传统的 APL 语义，你可能想要它类似于数字堆栈，但现在字典比堆栈、数组或列表更时尚，所以我们认为它类似于按符号索引的 Lisp 关联列表，每个符号表示一些变化轴。 现在我们可以解释这一行： ```apl proj = normal * XYZ.sum(raydir * normal) ``` 意味着（用 Python 表示）： ```python def proj(p): return normal(p) * sum(raydir(q) * normal(q) for q in XYZ.possibilities_augmenting(p)) ``` 这里 `possibilities_augmenting` 是 `XYZ` 维度的一个方法，它为你提供点 `p` 的版本，其中所有可能的 `XYZ` 值替换进去。因此，第一次调用 `normal` 可能返回某个特定法线的 `x`、`y` 或 `z` 分量；但所有这些都将乘以相同的点积。 当然，这并不是建议必须以此方式计算，那样会极其浪费；它旨在规范输入和输出之间的关系。 这提出了之前提到的 `vec` 函数的实现，在 C 程序中是一种 `struct` 类型： ```python def vec(x, y, z): values = {XYZ.x: x, XYZ.y: y, XYZ.z: z} return lambda p: values[p[XYZ]](p.without(XYZ)) ``` 也就是说，由 `vec` 生成的函数消耗输入的 `XYZ` 维度，并分发到作为其 X、Y 和 Z 分量的三个传入函数中的任何一个。因此，来自 `pixel_color` 函数的这个表达式： ```apl vec(xx/ww-.5, yy/hh-.5, 1) ``` 当用 `z` 调用时，将分派到表示 `1` 的常函数；当用 `y` 调用时，将分派到表示 `yy/hh-.5` 的函数，其自身将分派到 `yy`（在此程序中随像素变化），到 `hh`（在程序运行期间根本不变化），以及另一个返回 0.5 的常函数。 另一个有用的高阶函数是“重命名”、“别名”、“轴旋转”或“重塑”函数，它将一个轴转变为另一个轴： ```python def rename(a, b, f): return lambda p: f(p.without(a).with(b, p[a])) ``` 从空间考虑，这防止 `f` 在轴 `a` 上变化，将 `p` 沿 `a` 的运动旋转为沿新轴 `b` 的运动。从关系角度考虑，它将关系 `f` 的输入列 `b` 重命名为 `a`。这是执行显式外积所需的运算符；如果 `f` 和 `g` 都是轴 `b` 上的向量，则 `rename(a,b,f)+g` 给出它们的外积和，`f` 的值沿新轴 `a` 变化，`g` 的值沿轴 `b` 如前变化。（此 `rename` 函数也在某种意义上提供了一般轴转置。） 此“上下文”或“点”对象可能携带大量大多数函数懒得访问的上下文属性，并且可以传递给它们被调用的地方而不显式提及它们。（如果我们考虑 N 元关系而不是函数，你可以将此“点”视为查询-示例部分记录。但这与上面描述的功能语义不太一致。） 理论上，如果所有组合通过提升运算符、`rename`、轴构造函数如 `vec` 等组成的组件函数