julia - 矩阵列上的内存高效 sortperm

Question

我有一个大矩阵，我想申请 sortperm到该矩阵的每一列。天真的做法是

order = sortperm(X[:,j])

这制作了一个副本。这似乎是一种耻辱，所以我想我会尝试使用 SubArray :

order = sortperm(sub(X,1:n,j))

但这甚至更慢。为了笑我试过了

order = sortperm(1:n,by=i->X[i,j])

但这当然是可怕的。执行此操作的最快方法是什么？

这是一些基准代码:

getperm1(X,n,j) = sortperm(X[:,j])
getperm2(X,n,j) = sortperm(sub(X,1:n,j))
getperm3(X,n) = mapslices(sortperm, X, 1)
n = 1000000
X = rand(n, 10)
for f in [getperm1, getperm2]
    println(f)
    for it in 1:5
        gc()
        @time f(X,n,5)
    end
end
for f in [getperm3]
    println(f)
    for it in 1:5
        gc()
        @time getperm3(X,n)
    end
end

结果:

getperm1
elapsed time: 0.258576164 seconds (23247944 bytes allocated)
elapsed time: 0.141448346 seconds (16000208 bytes allocated)
elapsed time: 0.137306078 seconds (16000208 bytes allocated)
elapsed time: 0.137385171 seconds (16000208 bytes allocated)
elapsed time: 0.139137529 seconds (16000208 bytes allocated)
getperm2
elapsed time: 0.433251141 seconds (11832620 bytes allocated)
elapsed time: 0.33970986 seconds (8000624 bytes allocated)
elapsed time: 0.339840795 seconds (8000624 bytes allocated)
elapsed time: 0.342436716 seconds (8000624 bytes allocated)
elapsed time: 0.342867431 seconds (8000624 bytes allocated)
getperm3
elapsed time: 1.766020534 seconds (257397404 bytes allocated, 1.55% gc time)
elapsed time: 1.43763525 seconds (240007488 bytes allocated, 1.85% gc time)
elapsed time: 1.41373546 seconds (240007488 bytes allocated, 1.82% gc time)
elapsed time: 1.42215519 seconds (240007488 bytes allocated, 1.83% gc time)
elapsed time: 1.419174037 seconds (240007488 bytes allocated, 1.83% gc time)

凡mapslices版本是 getperm1 的 10 倍版本，如您所料。

值得指出的是，至少在我的机器上，copy+sortperm 选项并不比相同长度的向量上的 sortperm 慢多少，但不需要内存分配，所以最好避免它。

Answer 1

在一些非常特殊的情况下(比如连续查看 Array)，您可以使用指针魔法击败 SubArray 的性能:

function colview(X::Matrix,j::Int)
    n = size(X,1)
    offset = 1+n*(j-1) # The linear start position
    checkbounds(X, offset+n-1)
    pointer_to_array(pointer(X, offset), (n,))
end

getperm4(X,n,j) = sortperm(colview(X,j))

函数colview将返回一个成熟的Array与原始 X 共享其数据.请注意，这是一个 可怕的想法 因为返回的数组正在引用 Julia 仅通过 X 跟踪的数据.这意味着如果 X在列“ View ”数据访问将因段错误而崩溃之前超出范围。

结果:

getperm1
elapsed time: 0.317923176 seconds (15 MB allocated)
elapsed time: 0.252215996 seconds (15 MB allocated)
elapsed time: 0.215124686 seconds (15 MB allocated)
elapsed time: 0.210062109 seconds (15 MB allocated)
elapsed time: 0.213339974 seconds (15 MB allocated)
getperm2
elapsed time: 0.509172302 seconds (7 MB allocated)
elapsed time: 0.509961218 seconds (7 MB allocated)
elapsed time: 0.506399583 seconds (7 MB allocated)
elapsed time: 0.512562736 seconds (7 MB allocated)
elapsed time: 0.506199265 seconds (7 MB allocated)
getperm4
elapsed time: 0.225968056 seconds (7 MB allocated)
elapsed time: 0.220587707 seconds (7 MB allocated)
elapsed time: 0.219854355 seconds (7 MB allocated)
elapsed time: 0.226289377 seconds (7 MB allocated)
elapsed time: 0.220391515 seconds (7 MB allocated)

我没有研究为什么 SubArray 的性能更差，但它可能只是来自每次内存访问时的额外指针取消引用。非常值得注意的是，分配实际上花费的时间很少 - getperm1 的时间变化更大，但它仍然偶尔会胜过 getperm4!我认为这是由于 Array 中的一些额外的指针数学引起的使用共享数据的内部实现。还有一些疯狂的缓存行为...... getperm1 在重复运行时变得明显更快。