julia - 识别稀疏矩阵中哪些行(或列)具有值

Question

我需要识别在大型稀疏 bool 矩阵中定义了值的行 (/columns)。我想用它来 1. 按那些行/列对矩阵进行切片(实际上是 view )；和 2. 切片 (/view ) 向量和矩阵，其维度与矩阵的边距相同。 IE。结果应该是一个索引向量/ bool 值或(最好)一个迭代器。

我已经尝试过显而易见的:

a = sprand(10000, 10000, 0.01)
cols = unique(a.colptr)
rows = unique(a.rowvals)

但是每个在我的机器上都需要 20 毫秒，可能是因为它们分配了大约 1MB(至少它们分配了 cols 和 rows)。这是在一个性能关键的函数中，所以我希望优化代码。基本代码似乎有一个 nzrange稀疏矩阵的迭代器，但我很难看到如何将其应用于我的案例。

有没有建议的方法来做到这一点？

第二个问题:我还需要对我的稀疏矩阵的 View 执行此操作 - 会像 x = view(a,:,:); cols = unique(x.parent.colptr[x.indices[:,2]]) 那样吗？或者是否有专门的功能？稀疏矩阵的 View 似乎很棘手(参见 https://discourse.julialang.org/t/slow-arithmetic-on-views-of-sparse-matrices/3644 – 不是交叉帖子)

非常感谢!

Answer 1

关于获取稀疏矩阵的非零行和列，以下函数应该非常有效:

nzcols(a::SparseMatrixCSC) = collect(i 
  for i in 1:a.n if a.colptr[i]<a.colptr[i+1])

function nzrows(a::SparseMatrixCSC)
    active = falses(a.m)
    for r in a.rowval
        active[r] = true
    end
    return find(active)
end

对于密度为 0.1 的 10_000x10_000 矩阵，列和行分别需要 0.2 毫秒和 2.9 毫秒。它也应该比所讨论的方法更快(除了正确性问题)。

关于稀疏矩阵的 View ，一个快速的解决方案是将 View 转换为稀疏矩阵(例如使用 b = sparse(view(a,100:199,100:199)) )并使用上面的函数。在代码中:

nzcols(b::SubArray{T,2,P}) where {T,P<:AbstractSparseArray} = nzcols(sparse(b))
nzrows(b::SubArray{T,2,P}) where {T,P<:AbstractSparseArray} = nzrows(sparse(b))

更好的解决方案是根据 View 自定义功能。例如，当 View 对行和列都使用 UnitRanges 时:

# utility predicate returning true if element of sorted v in range r
inrange(v,r) = searchsortedlast(v,last(r))>=searchsortedfirst(v,first(r))

function nzcols(b::SubArray{T,2,P,Tuple{UnitRange{Int64},UnitRange{Int64}}}
  ) where {T,P<:SparseMatrixCSC}
    return collect(i+1-start(b.indexes[2]) 
      for i in b.indexes[2]
      if b.parent.colptr[i]<b.parent.colptr[i+1] && 
        inrange(b.parent.rowval[nzrange(b.parent,i)],b.indexes[1]))
end

function nzrows(b::SubArray{T,2,P,Tuple{UnitRange{Int64},UnitRange{Int64}}}
  ) where {T,P<:SparseMatrixCSC}
    active = falses(length(b.indexes[1]))
    for c in b.indexes[2]
        for r in nzrange(b.parent,c)
            if b.parent.rowval[r] in b.indexes[1]
                active[b.parent.rowval[r]+1-start(b.indexes[1])] = true
            end
        end
    end
    return find(active)
end

它比完整矩阵的版本工作得更快(对于 10,000x10,000 矩阵以上的 100x100 子矩阵，列和行在我的机器上分别需要 16μs 和 12μs，但这些结果不稳定)。

适当的基准测试将使用固定矩阵(或至少固定随机种子)。如果我这样做，我会用这样的基准来编辑这一行。