确定哪些行(或列)在稀疏矩阵中具有值

Question

我需要识别在大型稀疏布尔矩阵中定义了值的行(/列).我想用它来1. view用那些行/列切片(实际)矩阵; 2.切片(/ view)矢量和矩阵,其尺寸与矩阵的边距相同.即结果应该是索引/ Bools的Vector或(最好)迭代器.

我试过了明显的事:

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

但每一种把我的机器上像20ms的,可能是因为他们所分配的1MB左右(至少他们分配cols和rows).这是一个性能关键的功能,所以我希望优化代码.Base代码似乎有一个nzrange稀疏矩阵的迭代器,但我不容易看到如何将它应用于我的情况.

是否有建议的方法这样做？

第二个问题:我还需要对我的稀疏矩阵的视图执行此操作 - 这会是类似的x = view(a,:,:); cols = unique(x.parent.colptr[x.indices[:,2]])还是有专门的功能？稀疏矩阵的视图似乎很棘手(参见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.2ms和2.9ms.它也应该比有问题的方法更快(除了正确性问题).

关于稀疏矩阵的视图,快速解决方案是将视图转换为稀疏矩阵(例如,使用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))

更好的解决方案是根据视图自定义功能.例如,当视图对行和列使用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

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

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