Julia中对称矩阵的高效求和方法探究

问题背景 #

需要对n×n对称且对角线元素为0的矩阵A求和：

• 直接调用sum(A)速度最快，但未利用矩阵对称性，存在计算冗余；
• sum(tril(A, -1))速度明显更慢；
• 嵌套循环sum(A[i, j] for i = 1:n-1 for j = i+1:n)性能最差。

希望找到兼顾效率与对称性利用的求和方式。

高效解决方案 #

@AboAmmar提供的手动遍历实现性能优异，核心逻辑是仅计算一半非对角线元素的和，利用对称性翻倍后得到全部非对角线元素的总和（若对角线非0可额外累加，本题中可省略）。

实现代码 #

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  # 利用对称性，下三角元素和与上三角相等
    # 若对角线元素为0，可直接移除以下循环
    for i = 1:n
        s += A[i, i]
    end
    return s
end

性能测试对比 #

使用1000×1000的对称零对角线矩阵进行测试，代码如下：

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)

测试方案说明 #

• sum(A)：直接遍历全矩阵，无对称性利用，存在冗余计算；
• 2 * sum(tril(A))：计算下三角区域和后翻倍，但tril(A)会生成新矩阵，带来额外内存开销；
• sum_triu(A)：手动遍历上三角区域，结合编译器优化指令减少额外开销，无冗余计算。

关键优化点 #

• 无额外内存分配：直接操作原矩阵元素，避免像tril(A)那样生成新矩阵；
• SIMD指令加速：@simd标记让编译器使用单指令多数据指令并行处理循环；
• 消除边界检查：@inbounds跳过数组索引的边界检查，减少运行时损耗。

内容的提问来源于stack exchange，提问作者matnor