高效地产生反对称矩阵

综合讨论区
1

我的项目需要反复从一个给定矢量构建一个反对称矩阵，例如可以通过长度为 L*(L-1)/2 的数组生成一个 L*L 的反对称矩阵，最容易想到的是下面的代码：

function test1(d::Int64, D::Int64, v::Vector{Float64})
    L = d * D
    A = zeros(eltype(v),(L, L))
    count = 1
    for j = 1:L
        for i=1:L
               if i > j
                   A[i,j] = v[count]
                   A[j,i] = -v[count]
                   count += 1
               end
          end
    end
    return A
end

但考虑到 Julia 以列为主的存储形式，for 循环中的 A[j,i] = -v[count] 跨越了不同的列，因而似乎这不是优化的实现？

请问各位有什么优化的想法吗？ 谢谢各位！

2

加上 @inbounds 可以轻微提升(7.7->6.2μs)，如果 L 不变 A 能重复利用，可以在外边创建，减少开销 (6.2->5.8)

const A = zeros(Float64, (100, 100)) # L = 100
function test1(v::AbstractVector, L::Int)
    # A = zeros(Float64, (L, L))
    count = 1
    @inbounds for i = 1:L, j = 1:L
        if i > j
            A[i, j] = v[count]
            A[j, i] = -v[count]
            count += 1
        end
    end
    return A
end

using BenchmarkTools
v = rand(4950)
@btime test1($v, 100);

关于内存顺序，我分别注释了 A[i, j] 和 A[j, i] 但计算开销区别很小，其他因素可能更重要，比如 count += 1 使这段代码不能并行进行。

1 个赞
3

有点班门弄斧…

function test2(d::Int64, D::Int64, v::Vector{Float64})
  L = d * D
  A = zeros(eltype(v),(L, L))
  
  for i = 2:L
    for j = 1:i-1
      A[i,j] = -v[(i-2)*(i-1)/2+j]
    end
  end
  return A-transpose(A)
end

100维的反对称矩阵，大概快5us(15.658-> 10.524)

image
image1078×110 15.9 KB

而且应该更容易写并行吧…

2 个赞
4

这段应该写反了， i 不变 j 先走 是行优先遍历（除号 / 可以用 ÷ 或 >> 避免浮点数）。

等号右边的 i, j 反过来就得到上三角的索引了，可以避免加法运算。

A = Array{eltype(v)}(undef, L, L)
@inbounds for i in 1:L, j in 1:L
    A[j, i] = if i > j
        v[((i-1)*(i-2))>>1+j]
    elseif i < j
        -v[((j-1)*(j-2))>>1+i]
    else
        0
    end
end

题主如果对构建规则的要求宽松，可以考虑用方阵 v 代替向量再跟负转置求和。

ADD:

@LDA111222 :joy: 抱歉没表述好，我是想说用索引来并行的点子不错，然后补充了可能的改进。。。

5

谢谢，@inbounds 以及 const A 的声明的确有帮助，后者可以把开销降为 0

6

谢谢！A-transpose(A) 中的 transpose(A) 会把开销翻倍；我在 L = 1000 下测试似乎速度反而变慢？

8

啊 我只试了小规模的例子…
抱歉了我不太懂（只是提醒一下可以绕开count +=1 这样并行更容易实现吧）
我去看下transpose 的源码（逃

1 个赞
9

没事，谢谢，我也只是来学习一下。