从Julia中的短向量快速随机选择

Question

我有一个简单的函数出现在我的Julia代码中的几个地方,并且在循环中运行了数百万次.该函数基本上就是这样rand([1,-1,im,-im]),它选择了四个可能给定值中的一个.我注意到这个函数在我的庞大循环中需要花费大量时间,因此,我尝试以稍快的方式编写它:

function qpsk()
  temp1 = ifelse(rand(Bool), 1+0im, -1+0im)
  temp2 = ifelse(rand(Bool), 1+0im,  0+1im)
  temp1*temp2
end

然后,它通常像这样调用:

sig = complex(zeros(N))
for i = 1:N
  sig[i] = qpsk()
end

现在,有没有办法进一步优化这个功能,或者使用另一种更快的方法？感谢您的帮助.

---

对当前答案的评论:

@DanGetz(22行??)的答案并没有解决问题,因为目前,Julia并不像使用显式循环那样擅长向量.另外,我简单的qpsk2(s)下面1行比Dan的原始答案中的那些"神秘"22行代码快2倍(创建了一个矢量,这增加了更多的时间).

但问题仍然存在,为什么他们没有实现类似qpsk1下面的东西？为什么我的原始qpsk分支比qpsk4(s)下面的直接快3倍？

如果有经验的人喜欢参加,我在下面添加了更多版本来指导讨论.

qpsk1(s) = s[1+(rand(Int8)&3)]          # Blazingly fast
qpsk2(s) = s[1+rand(Bool)+2rand(Bool)]  # Very fast
qpsk3(s) = s[rand(1:4,1)]               # Compiler issue here?
qpsk4(s) = s[rand(1:4)]                 # Why slow?
qpsk5(s) = rand([s])                    # Ridiculously slow!!
function test_orig(n)                   # Test qpsk(), very fast(branching!), why?
  for i = 1:n
    qpsk()
  end
end

using StaticArrays
function test(func, n)                  # Test all qpsk1 --> qpsk5
  s = SVector(1,-1,im,-im)
  for i=1:n
    func(s)
  end
end

@time test(qpsk1,10^8)  0.554994 seconds (5 allocations: 176 bytes)
@time test(qpsk2,10^8)  0.755286 seconds (5 allocations: 176 bytes)
@time test(qpsk3,10^8) 13.431529 seconds (400 M allocations: 26.822 GiB, 20.68% gc time)
@time test(qpsk4,10^8)  2.520085 seconds (5 allocations: 176 bytes)
@time test(qpsk5,10^8) 10.881852 seconds (200 M allocations: 20.862 GiB, 19.76% gc time)
@time test_orig(10^8)   0.771778 seconds (5 allocations: 176 bytes)
@time nqpsk2(10^8);     1.402830 seconds (9 allocations: 1.490 GiB, 6.39% gc time)

Answer 1

答案摘要

[(-1)^b1*im^b2 for (b1,b2) in zip(rand!(BitVector(N)),rand!(BitVector(N)))]

更快地生成长度为N的向量.

回答

计算随机位是大部分工作,因此从使用RandomNumbers.jl的注释中探索Chris的想法值得一试.另外,我们可以使用@ rickhg12hs的想法从生成的每个随机数中提取更多位.无论如何,一起生成一个值块对于更好的优化至关重要.

例如,下面的代码(nqpsk1使用qpsk从提问作为基准.nqpsk2是建议的改善):

function qpsk()
  temp1 = ifelse(rand(Bool), 1+0im, -1+0im)
  temp2 = ifelse(rand(Bool), 1+0im,  0+1im)
  temp1*temp2
end

nqpsk1(n::Int) = [qpsk() for i=1:n]

nqpsk2(n::Int) = begin
    res = zeros(Int,2*n)
    blocks = n >>> 4                 # use blocks of 16 values
    btail = n & 0x000000000000000f   # in case n is not a multiple of 16
    pos = 1
    @inbounds for i=1:blocks
        bits = rand(UInt32)          # get random bits for a whole block
        for j=1:16
            b1 = Bool(bits & 1)
            bits >>>= 1
            b2 = Bool(bits & 1)
            bits >>>= 1
            res[pos+b1] = (-1)^b2
            pos += 2
        end
    end
    @inbounds for i=1:btail
        res[pos+rand(Bool)] = (-1)^rand(Bool)
        pos += 2
    end
    return reinterpret(Complex{Int64},res)
end

在我的设置上实现了> 4倍的改进(Julia 0.7):

julia> using BenchmarkTools

julia> @btime nqpsk1(320);
  8.791 ?s (323 allocations: 15.19 KiB)

julia> @btime nqpsk2(320);
  1.056 ?s (3 allocations: 5.20 KiB)

更新

只有适度的速度折扣(和一些分配),但更好看的代码:

function nqpsk3(n::Int)
    res = zeros(Int,2n)
    rv1 = rand!(BitVector(n))
    rv2 = rand!(BitVector(n))
    @inbounds for (b1,b2,i) in zip(rv1,rv2,1:2:2n)
        res[i+b1] = (-1)^b2
    end
    return reinterpret(Complex{Int},res)
end

基准:

julia> @btime nqpsk3(320);
  1.780 ?s (11 allocations: 5.83 KiB)

附录

并且单 - (包裹) - 线版本也可以(2.48μs):

nqpsk4(n) = [(1+0im,-1+0im,0+im,0-im)[2b1+b2+1] for
  (b1,b2) in zip(rand!(BitVector(n)),rand!(BitVector(n)))]

最后,真正的单行版本(1.96μs):

nqpsk5(n) = [(-1)^b1*im^b2 for (b1,b2) in zip(rand!(BitVector(n)),rand!(BitVector(n)))]