julia - Julia 中对称矩阵的快速求和

Question

我想对维度为 n 乘以 n 的矩阵 A 中的所有元素求和。该矩阵是对称的并且对角线上有 0。我发现最快的方法就是求和(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)