在 Julia 中从大型数组中对列集进行子采样的最快方法答案

【问题标题】：Fastest way to subsample sets of columns from large arrays in Julia在 Julia 中从大型数组中对列集进行子采样的最快方法
【发布时间】：2016-11-24 08:27:00
【问题描述】：

我有兴趣对大型数据矩阵的不同随机采样子集执行统计，现在我的代码中的一个瓶颈是实际的子采样。这对我来说似乎很奇怪，因为对子采样数据有相当多的 O(N^2) 距离计算。我根本没有修改子样本，只是看着它们。

using Distributions, Distances

function test_subsetting(X; batch_size=500, nloops=100)
    nfeatures, nsamples = size(X)
    ref_samples = X[:,1:10]
    batch_inds = zeros(batch_size)
    batch = zeros(nfeatures,batch_size)
    d_matrix = zeros(batch_size,batch_size)
    for i = 1:nloops
        batch_inds = sort(sample(1:nsamples, batch_size, replace = false))
        batch = X[:,batch_inds]
        d_matrix = pairwise(SqEuclidean(), batch, ref_samples)
    end
end

当我在 5000 个特征乘 50000 个样本矩阵上对其进行测试时：

X_test = randn(5000,50000);

我发现我在_unsafe_getindex 中花费了大约一半的时间在multidimensional.jl 中，另一半在计算距离。

有没有更有效的方法来解决这个问题？

【问题讨论】：

如果您首先完全打乱矩阵 X，然后从中绘制后续批次的 n 列（如果需要，一旦到达末尾就重新打乱它），那么这可以加快速度。你也可以考虑并行化事情。这种采样非常适合一次在多个工作人员上运行此功能。
并行化是一个有趣的想法，但是我在每个批次上运行的统计信息会影响下一轮下一批的统计信息（抱歉，这在我的示例中不清楚），所以我不要认为它是微不足道的并行化。理想情况下，我想将其扩展到内存无法容纳的东西，所以我不重新洗牌是要走的路。
我想我确实需要创建batch，或者可能是一个替代品，它实际上并不为新数组分配新内存，而是将数据保存在X 中并伪造它-- 我收集了sub 和view 做这样的事情，但我无法让他们使用distances.jl 的pairwise 函数。
重新并行化 - 至少您可以让一个进程生成样本，另一个进程分析它们。
不错，这是个好主意。

标签： arrays subset julia

【解决方案1】：

这在 julia 0.5 上适用于我：

julia> using Distances, Distributions

julia> X = randn(500,1000);

julia> S = sample(1:1000,500,replace=false);

julia> M = view(X, :, S);

julia> S2 = sample(1:1000,500,replace=false);

julia> R = view(X, :, S2);

julia> pairwise(SqEuclidean(), M, R)
500×500 Array{Float64,2}:
  994.67 ...
...

在 julia 0.5 上的view 在 julia 0.4 上称为 slice（或 sub，在这里它们将是相同的）。不要与ArrayViews.view 混淆，后者的功能类似，但实现完全不同。

从理论上讲，您似乎应该能够将view 替换为slice，但似乎在julia 0.4 上缺少At_mul_B! 方法。所以你可能会被困在制作副本中。

【讨论】：

我想是时候升级到 0.5 了。感谢您阐明不同的 view 实现。

【解决方案2】：

如果您乐于“预生成”一个不错的随机索引大矩阵，然后在运行时引用该矩阵，则可以完全缩短“采样”时间。您甚至可以在每次使用之前对该矩阵进行一些“行和列”改组，而成本最低。

另外，你为什么需要排序？这肯定会破坏随机抽样的意义，以及引入不必要的计算吗？

【讨论】：

我在对样本进行排序时认为，如果列按顺序排列，数组访问可能会稍快一些。子样本的顺序无关紧要，因为对它们执行的统计在排列下是不变的。使用大的预先生成的主列表是一个有趣的想法，但瓶颈不是生成样本索引，而是使用这些索引来访问数据矩阵中条目的子集。跨度>
batch = X[:,batch_inds] 行似乎过慢了——我想知道是否有人想过使用 view 或 sub 或类似的东西。
@Rory 然后摆脱sort。排序是一项非常昂贵的操作，您这样做的目的是实现 micro 优化！此外，它基于一个误解，实际上它不会做任何事情。 Julia 是列主序的，大致意思是单个列中的所有元素（即该列的行元素）都连续存储在内存中。这意味着在相邻列中的同一行中的两个元素实际上在内存地址方面彼此相距很远。所以你的排序“优化”在这里几乎没有什么区别。
@Rory 也是，我发现 batch = X[:,batch_inds] 不太可能是罪魁祸首，因为如果有什么 julia 被优化要做的就是访问矩阵元素......如果我不得不猜测，我会说batch_inds 可能更可能在需要时才实际评估......因此该行是所有排序和采样实际发生的地方。但这只是一个猜测。如果任何对 JIT 优化如何工作有更好了解的人都可以对此有所了解，那将很有趣。（您可以通过从预先存在的矩阵中绘制 batch_inds 值来测试它）
但是，您正在为索引排序一小组整数。这是非常快的。分解列的内存读取，以便在基于无序索引进行切片时来回切换它们的顺序，这会给您带来更大的打击。