浏览器中更快的NumPy

# 漫漫长路的终点：更快的浏览器NumPy 来源：https://notebook.link/blog/the-last-mile-faster-numpy/ 长期以来，在浏览器中运行NumPy意味着无法使用加速BLAS。矩阵乘法退化为简单的循环（可移植，但忽略了缓存和SIMD）。 这种情况刚刚发生了改变。Emscripten-forge（https://emscripten-forge.org/）上的NumPy包现在将OpenBLAS（https://www.openblas.net/）链接为WebAssembly，当`n = 1024`时，`np.matmul`的`float32`性能提升约**30.92倍**，`float64`性能提升约**14.90倍**。下一个OpenBLAS版本（已在Emscripten-forge上作为实验包提供，由QuantStack贡献了内核）将进一步提升性能，而可选的Relaxed SIMD（https://github.com/WebAssembly/relaxed-simd）构建在支持该特性的引擎上又增加了一层加速。 ## 什么是Emscripten-forge？（https://notebook.link/blog/the-last-mile-faster-numpy/#what-is-emscripten-forge） Emscripten-forge（https://emscripten-forge.org/）是一个面向WebAssembly的软件发行版。它与conda-forge（https://conda-forge.org/）合作，为WebAssembly重建了conda生态系统，使得编译器、运行时和共享库以具有统一ABI的**conda包**形式分发（不是Python轮子，也不是R包，而是任何语言生态系统都可以链接的原生库）。 将科学语言带到浏览器已经足够困难，到目前为止，这通常一次只针对一种语言实现。Pyodide（https://pyodide.org/）为Python开创了这一领域。后来，WebR（https://docs.r-wasm.org/webr/latest/）在相同前提下为R实现了这一点。每一个都是独立的、特定于语言的发行版。Emscripten-forge采用了不同的形式：一个语言无关的发行版。它建立在Pyodide和WebR项目的基础工作之上，并进一步扩展，不仅涵盖Python和R，还包括C++、OCaml、Lua、编译器工具链以及许多其他工具（并附带了一个包管理器）。 ## 缺失的编译器：wasm32上的Fortran（https://notebook.link/blog/the-last-mile-faster-numpy/#the-missing-compiler-fortran-on-wasm32） Fortran是基础科学软件包的核心：LAPACK（https://www.netlib.org/lapack/）、历史上大量的SciPy（https://scipy.org/）（尽管随着SciPy努力成为无Fortran的项目，情况正在改变），以及R的编译统计程序底层都是Fortran。因此，将其中任何一个带到WebAssembly都意味着需要一个能以`wasm32`为目标平台的Fortran编译器。 有一段时间，Pyodide通过变通方法解决了缺乏此类编译器的问题：其SciPy构建使用了f2c（https://www.netlib.org/f2c/），一个Fortran到C的转换器，因此Fortran源代码被转换为C并使用现有的Emscripten工具链进行编译。这个变通方案使SciPy能够在浏览器中运行，但对于R来说还不够，R的发行版无法绕过真正的Fortran编译器。正是这一限制启动了针对WebAssembly的Flang（https://flang.llvm.org/）（LLVM的Fortran前端）开发工作，这是将R通过Emscripten-forge（https://medium.com/jupyter-blog/r-in-the-browser-announcing-our-webassembly-distribution-9450e9539ed5）带到浏览器计划的一部分。 ## WebAssembly中的OpenBLAS（https://notebook.link/blog/the-last-mile-faster-numpy/#openblas-in-webassembly） OpenBLAS（https://www.openblas.net/）是一个广泛使用的、优化的BLAS（https://www.netlib.org/blas/）（基础线性代数子程序）实现。BLAS分为三个层次：**第1层**涵盖向量-向量操作（例如AXPY和点积）；**第2层**涵盖矩阵-向量操作（例如GEMV，如`A @ x`）；**第3层**涵盖矩阵-矩阵操作（例如GEMM，如`np.matmul`），这些操作主导了大规模稠密线性代数。OpenBLAS还提供了LAPACK（https://www.netlib.org/lapack/）（线性代数包）：用于求解线性系统、分解和特征问题的高级程序，大多数`np.linalg`调用最终都会委托给这些LAPACK入口点。 为WebAssembly构建Fortran工具链本身是一项多人协作的工作。Isabel Paredes（https://github.com/IsabelParedes）和Serge Guelton（https://github.com/serge-sans-paille）受George Stagg在WebAssembly上运行Fortran的工作（https://gws.phd/posts/fortran_wasm/）启发，编写了Flang的初始WebAssembly补丁；Serge Guelton（https://github.com/serge-sans-paille）还将其中部分工作上游合入了LLVM。Isabel Paredes（https://github.com/IsabelParedes）目前在Emscripten-forge上维护着最新的配方`flang_emscripten-wasm32`（https://github.com/emscripten-forge/recipes/tree/main/recipes/recipes/flang_emscripten-wasm32）。 但是，即使有了可用的Flang，仅靠工具链还不足以将OpenBLAS带到浏览器。最初的OpenBLAS和LAPACK构建无法立即使用：WebAssembly强制要求比原生目标更严格的调用约定，并且需要额外的补丁才能使库在Emscripten下编译、归档和链接。Ian Thomas（https://github.com/ianthomas23）为Flang和Emscripten适配了OpenBLAS；这些补丁位于Emscripten-forge的OpenBLAS配方中（https://github.com/emscripten-forge/recipes/tree/main/recipes/recipes_emscripten/openblas）（即上游OpenBLAS和本文测试的NumPy构建之间的桥梁）。 本文的剩余部分将走完这最后一段路：当NumPy最终能够调用BLAS时发生了什么变化，以及OpenBLAS 0.3.34和即将到来的0.3.35版本带来了什么。 ## 最后一公里：NumPy链接OpenBLAS（https://notebook.link/blog/the-last-mile-faster-numpy/#the-last-mile-numpy-links-openblas） 最后一步是将NumPy链接到加速的BLAS。之前的Emscripten-forge（https://emscripten-forge.org/）NumPy包（`emscripten-forge-4x`（https://prefix.dev/channels/emscripten-forge-4x/packages/numpy）上的`2.5.2`）未包含OpenBLAS，因此`np.matmul`、`@`和`np.linalg`退回到了可移植的C循环：忽略了缓存层次结构和SIMD的朴素实现。矩阵乘法的大部分时间都花在内存传输而非算术运算上。这是我们的**无BLAS基准**。 从`emscripten-forge/recipes#6310`（https://github.com/emscripten-forge/recipes/pull/6310）开始，`emscripten-forge-4x`（https://prefix.dev/channels/emscripten-forge-4x/packages/numpy）上的默认NumPy（`2.5.3`，构建号3）链接了带WASM SIMD（`TARGET=WASM128_GENERIC`）的OpenBLAS**0.3.34**。现在这是一个稳定的、主渠道发布版本。对`np.matmul`、`@`和`np.linalg.solve`的调用现在会分发到WebAssembly中的加速BLAS。 正是这种打包模式使其变得强大。Emscripten-forge实际上是**面向WebAssembly的conda-forge**：包从配方构建，发布在conda渠道上，并通过包管理器在共享ABI下安装：与`linux-64`上的工作流程相同，只是以`emscripten-wasm32`为目标。在这种模式下，OpenBLAS是一个独立的conda包。NumPy在运行时动态链接`libopenblas`，因此您可以升级OpenBLAS而无需重新构建NumPy。这有益于整个生态系统：像SciPy和scikit-learn这样的科学Python项目，以及像xtensor-blas（https://github.com/xtensor-stack/xtensor-blas）这样的非Python技术栈，都共享同一个库。相比之下，PyPI上的NumPy轮子在构建时内嵌了一个BLAS快照。因此，相同的NumPy 2.5.3包将在OpenBLAS**0.3.35**发布后自动使用它。 ## OpenBLAS 0.3.34与无BLAS对比（https://notebook.link/blog/the-last-mile-faster-numpy/#openblas-0334-versus-no-blas） 相对于同一Emscripten-forge技术栈但未使用BLAS实现的情况，`n = 1024`时的方形`np.matmul`分别达到了约**30.92倍**（`float32`）和**14.90倍**（`float64`）的性能（28.0 / 14.6 GFLOPS vs 0.90 / 0.98 GFLOPS）。完整尺寸网格和单线程linux-64参考数据见附录F（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-F-Full-benchmark-results-including-linux-64）；此比较的重点图表和表格见附录B（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-B-OpenBLAS-0334-versus-no-BLAS）。 **OpenBLAS 0.3.34相对于无BLAS NumPy的几何平均加速比（按Python API）** （几何平均基于两种数据类型和尺寸网格，分母为OpenBLAS 0.3.34耗时）。矩阵乘法受益最大。某些`np.linalg.*`API的增益较小（相对无BLAS为1.08–1.65倍），因为它们委托给LAPACK，而OpenBLAS中LAPACK的实现尚未针对WebAssembly进行优化或特化。`A @ x`和向量`np.dot`仍接近1倍：OpenBLAS 0.3.34未为此布局提供快速的列主序GEMV。OpenBLAS 0.3.35添加了该内核。 无BLAS时，`np.matmul`吞吐量随`n`增加而下降（朴素GEMM，缓存未命中）。使用OpenBLAS 0.3.34后，它随`n`增加而提高。基于GEMM实现的操作（`@`、二维`np.dot`、`np.tensordot`、`np.linalg.multi_dot`）遵循相同趋势。相对`n`的加速比及相关数据见附录B（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-B-OpenBLAS-0334-versus-no-BLAS）。 在`n = 1024`（`float32`）时，`np.linalg`的中位时间相对无BLAS减少了约**1.2–2.1倍**（`solve` 2.11倍，`cholesky` 1.88倍，`qr` 2.03倍，`inv` 1.66倍，`eigh` 1.60倍，`svd` 1.21倍；`float64`结果在几个百分点内一致）。对于矩阵-向量乘积，`x @ A`在0.3.34中已经使用了第2层内核（`float32`、`n = 1024`时为**14.86倍**），而`A @ x`仍接近无BLAS的约4 GFLOPS（**1.02倍**）。 ## OpenBLAS 0.3.35将会更快（https://notebook.link/blog/the-last-mile-faster-numpy/#openblas-0335-will-be-even-faster） 下一个OpenBLAS版本包含了已在`develop`分支上的WASM SIMD工作。下面的数字使用了由`emscripten-forge/recipes#6761`（https://github.com/emscripten-forge/recipes/pull/6761）发布的`openblas-experimental`包（`0.3.35.dev0`，`develop` @ `539bb47`（https://github.com/OpenMathLib/OpenBLAS/commit/539bb47f020d3278a05805f86c1ac63326b8a3d1），仍然是`TARGET=WASM128_GENERIC`）：可移植的**`simd128_hbf81ecf_3`**（默认）和可选的**`relaxed_simd_h13a9a80_3`**。WebAssembly Relaxed SIMD（https://github.com/WebAssembly/relaxed-simd）是一种引擎扩展，它允许略微宽松的浮点语义（尤其是融合乘加）以换取更快的向量数学运算；当构建启用时，OpenBLAS可以使用这些操作码。除非另有说明，本文中的**0.3.35**指的是可移植的SIMD128构建。 QuantStack向OpenBLAS上游贡献了WASM SIMD内核（Julien Jerphanion（https://github.com/jjerphan），Matthias Meschede（https://github.com/MMesch）），由Martin Kroeker（https://github.com/martin-frbg）审核和整合。这项工作涵盖了第3层GEMM/TRMM、第1/2层AXPY和GEMV、用于Node/Emscripten的数值CBLAS测试以及可选的Relaxed SIMD路径。可移植的SIMD128仍然是默认选项：Relaxed SIMD FMA默认关闭，除非您选择`relaxed_simd`构建。主要的`float32`GEMM步骤从0.3.34开始是8x4微内核；可选的Relaxed SIMD变体将在下一小节介绍。 与WASM OpenBLAS 0.3.34在`n = 1024`时相比，`np.matmul`提升了**1.79倍**（`float32`，28.0 → 50.0 GFLOPS）和**1.41倍**（`float64`，14.6 → 20.5 GFLOPS）。此比较的图表和数据见附录C（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-C-OpenBLAS-0335-SIMD128-versus-OpenBLAS-0334）；完整网格见附录F（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-F-Full-benchmark-results-including-linux-64）。 **OpenBLAS 0.3.35相对于OpenBLAS 0.3.34的几何平均加速比（按Python API）** （几何平均基于两种数据类型和尺寸网格，分母为OpenBLAS 0.3.35耗时）。相对WASM 0.3.34的最大相对增益不仅在GEMM上：**`A @ x`**通过GEMV内核从无BLAS的约4 GFLOPS跃升至**23.1 GFLOPS**（`float32`、`n = 1024`）（相对0.3.34为**5.52倍**；`x @ A`为**1.73倍**）。向量`np.dot`也有类似改善。方形`np.matmul`的几何平均加速比为**1.58倍**（在`n = 1024` `float32`时为**1.79倍**，得益于8x4 SGEMM）。在`n = 1024` `float32`的`np.linalg`上，中位时间相对0.3.34改善了约**1.3–1.9倍**（`solve` 1.38倍，`cholesky` 1.32倍，`qr` 1.49倍，`inv` 1.33倍，`eigh` 1.64倍，`svd` 1.90倍）。 ### 可选：在0.3.35基础上的Relaxed SIMD（https://notebook.link/blog/the-last-mile-faster-numpy/#optional-relaxed-simd-on-top-of-0335） 实现了WebAssembly Relaxed SIMD（https://github.com/WebAssembly/relaxed-simd）的引擎可以加载`relaxed_simd` OpenBLAS构建（`WASM_RELAXED_SIMD=1`），它使用来自#6001（https://github.com/OpenMathLib/OpenBLAS/pull/6001）/#6020（https://github.com/OpenMathLib/OpenBLAS/pull/6020）的FMA风格`v_muladd`。在本次测试的Chromium 153中，对于`n = 1024`的`np.matmul`，相对可移植的0.3.35又提升了**1.17倍**（`float32`）和**1.43倍**（`float64`）（50.0 → 58.7 GFLOPS 和 20.5 → 29.4 GFLOPS）。相对SIMD128的几何平均增益在矩阵乘法上最大（约1.3倍）；第2层和大部分`np.linalg`操作提升了约**1.0–1.2倍**。此比较的图表和数据见附录D（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-D-OpenBLAS-0335-Relaxed-SIMD-versus-OpenBLAS-0335-SIMD128）。 **OpenBLAS 0.3.35 Relaxed SIMD相对于可移植0.3.35的几何平均加速比（按Python API）** Chrome ≥114和Firefox ≥146已公开此特性；Safari需要JavaScriptCore标志`useWebAssemblyRelaxedSIMD`，否则模块将无法实例化。出于这个原因，可移植的`simd128`包仍然是默认选项（`emscripten-forge/recipes#6761`（https://github.com/emscripten-forge/recipes/pull/6761）降低了Relaxed SIMD变体的优先级）。当前引擎支持情况可在webassembly.org/features（https://webassembly.org/features/）跟踪。 ### 从无BLAS到OpenBLAS 0.3.35（https://notebook.link/blog/the-last-mile-faster-numpy/#from-no-blas-to-openblas-0335） 将各步骤结合（链接OpenBLAS 0.3.34，然后采用0.3.35 SIMD128内核，再可选地使用Relaxed SIMD），就是从之前无BLAS的Emscripten-forge NumPy开始的加速路径。 **OpenBLAS 0.3.34、0.3.35 SIMD128和0.3.35 Relaxed SIMD相对于无BLAS NumPy的几何平均加速比（按Python API）** （几何平均基于两种数据类型和尺寸网格，分子为无BLAS耗时）。每个API显示三个柱状条：**OpenBLAS 0.3.34**、**OpenBLAS 0.3.35（SIMD128）**和**OpenBLAS 0.3.35（Relaxed SIMD）**。矩阵乘法增益最大；一旦0.3.35的GEMV内核落地，第2层`A @ x`和向量`np.dot`最终开始提升；`np.linalg.*`位于1倍以上一个较小但清晰的区间内，受限于OpenBLAS中尚未针对WASM特化的LAPACK例程。Emscripten-forge上发布的NumPy已经提供了该路径中OpenBLAS 0.3.34的部分；可移植的0.3.35（来自`emscripten-forge/recipes#6761`（https://github.com/emscripten-forge/recipes/pull/6761）的`simd128_hbf81ecf_3`）今天已可在实验渠道使用，支持该特性的引擎还可选用Relaxed SIMD构建（`relaxed_simd_h13a9a80_3`）。在Relaxed SIMD下，`n = 1024`时的端到端`np.matmul`相对无BLAS达到了约**64.89倍**（`float32`）和**30.07倍**（`float64`）。此端到端比较的相对`n`加速比见附录E（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-E-All-changes-OpenBLAS-0335-Relaxed-SIMD-versus-no-BLAS）；包括linux-64在内的完整尺寸网格见附录F（https://notebook.link/blog/the-last-mile-faster-numpy/#Appendix-F-Full-benchmark-results-including-linux-64）。 ## 使用这些软件包（https://notebook.link/blog/the-last-mile-faster-numpy/#using-the-packages） 链接了OpenBLAS的NumPy位于主`emscripten-forge-4x`（https://prefix.dev/channels/emscripten-forge-4x）渠道（截至`emscripten-forge/r