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

Question

我想计算一个二次形式：朱莉娅中的x'Q y。
对于这种情况，最有效的计算方法是什么：

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

我知道朱莉娅有点（）。但是我想知道它是否比BLAS呼叫更快。

Answer 1

Julia的LinearAlgebra stdlib有3-参数点的本机实现，也是一个专门用于对称/厄米特矩阵的版本。您可以在这里和这里查看源代码。

您可以使用BenchmarkTools确认它们没有进行分配@b时间或基准工具@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 2

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

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 3

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

>

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

虽然默认情况下不使用它，但现在确实存在两个纯Julia线性代数包，可以匹配甚至击败传统的基于Fortran的BLAS。尤其是屋大维。jl：