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

Question

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

没有假设。
Q是对称的。
x和y是一样的(x = y).
Q 和 x = y.

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

Answer 1

如果您的矩阵是对称的，请使用 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 2

现有答案都很好。但是，补充几点：

虽然默认情况下 Julia 确实会在此处简单地调用 BLAS，但正如 Oscar 指出的那样，一些 BLAS 比其他的更快。特别是，在现代 x86 硬件上，MKL 通常会比 OpenBLAS 快一些。幸运的是，在 Julia 中手动选择 BLAS 后端非常容易。只需在 Julia 1.7 或更高版本的 REPL 中键入 using MKL 即可通过 MKL.jl 包将您切换到 MKL 后端。
虽然默认情况下不使用它，但现在确实存在几个纯 Julia 线性代数包，它们可以匹配甚至击败传统的基于 Fortran 的 BLAS。特别是 Octavian.jl:

Answer 3

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

Chained multiplication is also very slow, because it doesn't know there's a better algorithm.

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 智能地减少它。

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

Answer 4

Julia 的 LinearAlgebra stdlib 具有 3 参数 dot 的本机实现，还有一个专门用于 symmetric/hermitian 矩阵的版本。可以查看来源here and 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 不是多线程的。还有一些浮点数差异。

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

How can I efficiently calculate a quadratic form in Julia?

performance

julia

quadratic-programming