如何在 Julia 中有效地计算二次型？

Question

我想在 Julia 中计算二次形式：x' Q y。
对于这些情况，最有效的计算方法是什么：

1. 没有假设。
2. Q 是对称的。
3. x 和 y 是相同的 (x = y)。
4. Q 和x = y 都是对称的。

我知道 Julia 有 dot()。但我想知道它是否比 BLAS 调用更快。

Answer 1

现有的答案都很好。但是，还有几点：

• 虽然默认情况下 Julia 确实会在这里简单地调用 BLAS，但正如 Oscar 所指出的那样，一些 BLAS 比其他的更快。特别是，在现代 x86 硬件上，MKL 通常会比 OpenBLAS 快一些。幸运的是，在 Julia 中，手动选择 BLAS 后端非常容易。只需在 Julia 1.7 或更高版本的 REPL 中输入 using MKL 即可通过 MKL.jl 包切换到 MKL 后端。

• 虽然默认情况下不使用它，但现在确实存在一些纯 Julia 线性代数包，它们可以匹配甚至击败传统的基于 Fortran 的 BLAS。特别是Octavian.jl：

Answer 2

Julia 的 LinearAlgebra 标准库具有 3 参数 dot 的本机实现，以及专门用于对称/厄米矩阵的版本。可以查看来源here和here。

您可以使用BenchmarkTools.@btime 或BenchmarkTools.@ballocated 确认它们没有分配（请记住使用$ 插入变量）。矩阵的对称性被利用了，但是看看源代码，我看不出x == y 是如何实现任何严重的加速的，除了可能节省一些数组查找。

编辑：比较BLAS版本和原生的执行速度，可以做

1.7.0> using BenchmarkTools, LinearAlgebra

1.7.0> X = rand(100,100); y=rand(100);

1.7.0> @btime $y' * $X * $y
  42.800 μs (1 allocation: 896 bytes)
1213.5489200642382

1.7.0> @btime dot($y, $X, $y)
  1.540 μs (0 allocations: 0 bytes)
1213.548920064238

这是原生版本的一大胜利。但是，对于更大的矩阵，情况会发生变化：

1.7.0> X = rand(10000,10000); y=rand(10000);

1.7.0> @btime $y' * $X * $y
  33.790 ms (2 allocations: 78.17 KiB)
1.2507105095988091e7

1.7.0> @btime dot($y, $X, $y)
  44.849 ms (0 allocations: 0 bytes)
1.2507105095988117e7

可能是因为 BLAS 使用线程，而 dot 不是多线程的。还有一些浮点差异。

Answer 3

如果您的矩阵是对称的，请使用 Symmetric 包装器来提高性能（然后调用不同的方法）：

julia> a = rand(10000); b = rand(10000);

julia> x = rand(10000, 10000); x = (x + x') / 2;

julia> y = Symmetric(x);

julia> @btime dot($a, $x, $b);
  47.000 ms (0 allocations: 0 bytes)

julia> @btime dot($a, $y, $b);
  27.392 ms (0 allocations: 0 bytes)

如果x 与y 相同，请参阅https://discourse.julialang.org/t/most-efficient-way-to-compute-a-quadratic-matrix-form/66606 以讨论选项（但总的来说，dot 似乎仍然很快）。

Answer 4

您可以使用Tullio.jl 编写优化循环，一次扫描完成所有这些。但我认为它不会明显击败 BLAS：

链式乘法也很慢，因为它不知道有更好的算法。

julia> # a, b, x are the same as in Bogumił's answer

julia> @btime dot($a, $x, $b);
  82.305 ms (0 allocations: 0 bytes)

julia> f(a, x, b) = @tullio r := a[i] * x[i,j] * b[j]
f (generic function with 1 method)

julia> @btime f($a, $x, $b);
  80.430 ms (1 allocation: 16 bytes)

添加 LoopVectorization.jl 可能是值得的：

julia> using LoopVectorization

julia> f3(a, x, b) = @tullio r := a[i] * x[i,j] * b[j]
f3 (generic function with 1 method)

julia> @btime f3($a, $x, $b);
  73.239 ms (1 allocation: 16 bytes)

但我不知道如何处理对称情况。

julia> @btime dot($a, $(Symmetric(x)), $b);
  42.896 ms (0 allocations: 0 bytes)

虽然可能有线性代数技巧可以用 Tullio.jl 智能地减少它。

在这类问题中，基准权衡就是一切。