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

Question

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

performancejuliaquadratic-programming

23

我想计算一个二次形式：x' Q y，使用Julia编程语言。
以下是在不同情况下最有效的计算方式：

没有假设。
Q 是对称的。
x 和 y 相同 (x = y)。
同时满足 Q 对称和 x = y。

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

- Eric Johnson

如果你想更多了解Julia中的二次形式，请访问：https://forem.julialang.org/vinodv/exploring-quadratic-forms-with-julia-k6e - undefined

4个回答

13

Julia的LinearAlgebra标准库有原生实现的三参数dot，还有一个专门针对对称/共轭矩阵的版本。您可以在这里查看源代码和这里。

您可以使用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没有实现多线程。还存在一些浮点数差异。

- DNF

1

你能否检查一下使用@tturbo而不是@simd的@dot版本？我很想看看LoopVectorization在这里的表现如何。 - Oscar Smith

1

仅仅将 @simd 替换为 @turbo 可以略微提高速度，而 @tturbo 与 @turbo 看起来相同。这些宏放置在内循环中。我认为将外循环线程化会更有益。 - DNF

4

去掉 iszero 检查（以允许 @tturbo 在外部循环中运行）后，我得到了一个与 BLAS 在大规模情况下一样快的结果，并且比小规模情况下的 dot 快 3 倍。 - Oscar Smith

@OscarSmith，能否分享你使用的准确代码？ - Eric Johnson

10

如果您的矩阵是对称的，请使用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仍然比较快）。

- Bogumił Kamiński

dot 产生一个标量。 - Bogumił Kamiński

2

@Eric Johnson 不，与 y'*X*y 不同，dot(y,X,y) 似乎根本不分配内存。您可以使用 BenchmarkTools.@btime 进行检查。 - DNF

@DNF，你知道有没有BLAS函数可以处理上述问题吗？我认为使用BLAS意味着它会执行“矩阵x向量”然后进行内积运算。 - Eric Johnson

@OscarSmith通常比我更了解代码库，所以我有些犹豫，但是@less dot(a, X, b)显示了一种本地实现，而@btime显示零分配。我会再深入研究一下，然后发表答案。 - DNF

1

我可能应该听从自己的建议，在说实现方式之前查看源代码。对于提供错误信息感到抱歉 :) - Oscar Smith

显示剩余6条评论

7

你可以使用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智能地缩短时间，但在这种问题中，权衡基准测试非常重要。

- phipsgabler

这真的很棒。或许支持对称矩阵应该成为 Tulio.jl 添加的一个特性。 - Eric Johnson

1

如果我理解正确的话，Tullio总是将向量化转换为嵌套循环，因此这并没有太多意义。你可以将Symmetric要执行的优化实现表示为一个einsum表达式，然后再使用Tullio。 - phipsgabler

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- cbk · Accepted Answer

现有的答案都很好。但是还有一些额外的要点：

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

- 虽然默认情况下它没有被使用，但现在已经存在一些纯Julia线性代数包，可以匹敌甚至击败传统的基于Fortran的BLAS。特别是，Octavian.jl: