如何在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)

变化很小，大概是因为编译起来很简单。

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)

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

Answer 3

您要查找的是稀疏非零函数。

nonzeros(A)

返回稀疏矩阵A中结构非零值的一个向量。这包括显式存储在稀疏矩阵返回的矢量直接指向内部非零A的存储，以及对返回的载体的任何修改都将发生突变CCD_ 8。

你可以如下使用：

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个数量级。