Julia中对称矩阵的快速求和

Question

我想求一个维度为n乘以n的矩阵A中的所有元素的和。矩阵是对称的，对角线上有0。我发现最快的方法就是sum(A)。然而，这似乎是浪费，因为它没有使用我只需要计算矩阵的下三角的事实。然而，sum(tril(A, -1))明显慢，sum(A[i, j] for i = 1:n-1 for j = i+1:n)更慢。有没有更有效的方法来求和矩阵？

编辑：@AboAmmar的解决方案表现良好。以下是要比较的代码(对角线分别求和，如果对角线上只有零，则可以删除(：

using BenchmarkTools
using LinearAlgebra
function sum_triu(A)
m, n = size(A)
@assert m == n
s = zero(eltype(A))
for j = 2:n
@simd for i = 1:j-1
s += @inbounds A[i,j]
end
end
s *= 2
for i = 1:n
s += A[i, i]
end
return s
end
N = 1000
A = Symmetric(rand(0:9,N,N))
A -= diagm(diag(A))
@btime sum(A)
@btime 2 * sum(tril(A))
@btime sum_triu(A)

Answer 1

对于n=1000矩阵，这比sum快2.7倍。确保在循环之前添加一个@simd，并使用@inbounds。此外，使用正确的循环顺序进行快速内存访问。

function sum_triu(A)
m, n = size(A)
@assert m == n
s = zero(eltype(A))
for j = 1:n
@simd for i = 1:j 
s += @inbounds A[i,j]
end
end
return 2 * s
end

在我的电脑上运行的示例：

sum_triu(A) = 499268.7328022966
sum(A) = 499268.73280229873
93.000 μs (0 allocations: 0 bytes)
249.900 μs (0 allocations: 0 bytes)

Answer 2

怎么样

2 * sum(LowerTriangular(A))

help?> LA.LowerTriangular
LowerTriangular(A::AbstractMatrix)
Construct a LowerTriangular view of the matrix A.

tril创建一个新的矩阵，用于分配内存。由于LowerTriangular是现有矩阵中的view，因此没有内存分配。