有没有一种有效的方法可以在 Julia 中重用稀疏矩阵的稀疏模式？

Question

我想知道是否有一种有效的方法可以在 Julia 中重用稀疏矩阵的稀疏模式。

假设我有向量

I

、

J

和

V

，我使用

创建一个稀疏矩阵

using SparseArrays
A = sparse(I, J, V)

即使

I

和

J

未排序并且可能出现多次，这似乎也非常有效。

我现在想利用矩阵

A

，或更准确地说是它的稀疏模式，并用新值

W

填充它。下面的代码就是一个例子：

N = 10
I = vec([ 1, 1, 2, 2]' .+ (0:N))
J = vec([ 1, 2, 1, 2]' .+ (0:N))
V = vec([ 1,-1,-1, 1]' .* ones(N+1))
A = sparse(I, J, V)
W = vec([ 1,-1, 1,-1]' .* ones(N+1))
A = sparse(I, J, W)

这里，最后一行应该替换为更有效的内容，因为 A 的稀疏模式保持不变。

我希望这会像这样简单

A.nzval .= 0
@inbounds A[I + size(A, 2) .* (J .- 1)] += W

我还尝试了其他几种（并且更复杂的）方法，例如对

I

和

J

的值进行预排序，然后显式设置

A.nzval

的值。然而，我总是得出这样的结论：最有效的方法仍然是使用

A

从头开始重新创建

A = sparse(I, J, W)

，这显然没有意义，因为此时我们已经有了更多关于

A

的信息。

我如何利用之前创建的矩阵

A

？

Answer 1

这个问题与我们有许多稀疏矩阵需要使用相同的

I

、

J

坐标但不同的值进行初始化的情况相关。由于重复元素坐标的可能性，初始化稀疏矩阵的确切顺序变得复杂。

下面尝试使用

I

、

J

和专门的稀疏矩阵提取“时间表”来创建稀疏矩阵，然后可以更有效地用于不同的值。

如果解释不清楚，代码可能是：

using SparseArrays

function sparse_sched(I,J,n = maximum(I), m = maximum(J))
    D1 = collect(1:length(I))
    Z = sparse(I,J,1:length(I),n, m,(x,y)->(D1[y] = x; x))
    D2 = similar(D1)
    for i in 1:nnz(Z)
        D2[Z.nzval[i]] = i
    end
    for i in eachindex(D2)
        D1[i] == 0 || (D2[i] = D2[D1[i]])
    end
    return D2
end

function apply_sched!(A,S,V,op = +)
    A.nzval .= zero(eltype(A))
    for i in eachindex(V)
        v = A.nzval[S[i]]
        A.nzval[S[i]] = iszero(v) ? V[i] : op(v, V[i])
    end
    return A
end

这些函数可以像下面的代码一样使用：

julia> N = 10;

julia> I = vec([ 1, 1, 2, 2]' .+ (0:N));

julia> J = vec([ 1, 2, 1, 2]' .+ (0:N));

julia> V = vec([ 1,-1,-1, 1]' .* ones(N+1));

julia> W = vec([ 1,-1, 1,-1]' .* ones(N+1));

julia> A = sparse(I, J, V);
julia> sched = sparse_sched(I,J);

julia> apply_sched!(similar(A), sched, W) == sparse(I,J,W)
true

julia> apply_sched!(similar(A), sched, V) == A
true

最后两行表明新的构造方法与标准稀疏矩阵构造相同。

在本地进行基准测试，它似乎提供了 5 倍的速度提升。如果它被大量使用，那么经历这个麻烦可能是值得的。