如何加速在Julia中对稀疏数组应用Log（x + 1）的过程

Question

在Julia中，稀疏矩阵仅存储非零元素。

一些函数（例如所有基数的log(x+1)）将零映射为零，因此不需要应用于这些零元素。 （我认为我们会称之为幺半同态。）

我如何利用这个事实来加速操作？

示例代码：

X = sprand(10^4,10^4, 10.0^-5, rand)
function naiveLog2p1(N::SparseMatrixCSC{Float64,Int64})
    log2(1+N) |> sparse
end

运行中：

@time naiveLog2p1(X)

输出是：

elapsed time: 2.580125482 seconds (2289 MB allocated, 6.86% gc time in 3 pauses with 0 full sweep)

第二次运行时（因此预计函数已经被编译）：

elapsed time: 2.499118888 seconds (2288 MB allocated, 8.17% gc time in 3 pauses with 0 full sweep)

稍作更改，可能是因为编译非常简单。

- Frames Catherine White

---

只是想确保您知道内置的 log1p 利用了这一点。julia> @time log1p(X); 经过时间：4.81e-5 秒（分配了 93 kB） - Andreas Noack

@andreasnoackjensen 但是log1p（X）是自然对数 - Frames Catherine White

我可以看出你现在正在寻求所有基础知识，但是你可能也可以按照log1p的定义方法来处理其他基础知识。 - Andreas Noack

Answer 1

根据Julia手册中有关"稀疏矩阵操作"的建议，我将使用findnz()将稀疏矩阵转换为密集矩阵，对值执行对数运算，然后使用sparse()重新构造稀疏矩阵。

function improvedLog2p1(N::SparseMatrixCSC{Float64,Int64})
    I,J,V = findnz(N)
    return sparse(I,J,log2(1+V))
end

@time improvedLog2p1(X)
elapsed time: 0.000553508 seconds (473288 bytes allocated)

- epx

Answer 2

我的解决方案是实际上在数据结构内部进行操作：

mysparselog(N::SparseMatrixCSC) =
    SparseMatrixCSC(N.m, N.n, copy(N.colptr), copy(N.rowval), log2(1+N.nzval))

请注意，如果您想在稀疏矩阵上原地操作，这在实践中可能经常发生，那么这将是一项零内存操作。基准测试显示，这与@Oxinabox的答案类似，因为它在内存操作方面大致相同（尽管该答案实际上没有返回新矩阵，如mean输出所示）：

with warmup times removed:

naiveLog2p1
elapsed time: 1.902405905 seconds (2424151464 bytes allocated, 10.35% gc time)
mean(M) => 0.005568094618997372

mysparselog
elapsed time: 0.022551705 seconds (24071168 bytes allocated)
elapsed time: 0.025841895 seconds (24071168 bytes allocated)
mean(M) => 0.005568094618997372

improvedLog2p1
elapsed time: 0.018682775 seconds (32068160 bytes allocated)
elapsed time: 0.027129497 seconds (32068160 bytes allocated)
mean(M) => 0.004995127985160583

- IainDunning

---

你能澄清一下：“尽管那个答案实际上并没有返回新矩阵”吗？ - Frames Catherine White

Answer 3

你需要的是稀疏非零函数。

nonzeros(A)

返回稀疏矩阵A中结构非零值的向量。这包括在稀疏矩阵中明确存储的零。返回的向量直接指向A的内部非零存储，并且对返回的向量进行任何修改也会改变A。

你可以按以下方式使用它：

function improvedLog2p1(N::SparseMatrixCSC{Float64,Int64})
    M = copy(N)
    ms = nonzeros(M) #Creates a view, 
    ms = log2(1+ms) #changes to ms, change M
    M
end


@time improvedLog2p1(X)

第一次运行的输出为：

 elapsed time: 0.002447847 seconds (157 kB allocated)

第二次运行的输出为：

 0.000102335 seconds (109 kB allocated)

这是速度和内存使用方面提高了4个数量级。

- Frames Catherine White

---

nonzeros 不创建视图，这个函数所做的只是复制 N。 - IainDunning

值得知道。我会进行一些研究并完善我的答案（已核实，在3.0文档中相同）。 - Frames Catherine White

因此，这解释了我的解决方案相对于你的解决方案使用更高的内存（也许略微更快）的原因？ - Frames Catherine White

肯定会使用更高的内存，而且速度也不太清楚。 - IainDunning

只需将该行替换为 for i = 1:length(ms); ms[i] = log2(1+ms[i]); end。 - tholy