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

Question

我需要在大型稀疏布尔矩阵中识别具有定义值的行（/列）。我想用它来 1. 将矩阵切片（实际上是 view）rows/columns；和 2. slice (/view) 向量和矩阵，其维度与矩阵的边距相同。 IE。结果可能应该是一个 Vector of indices / Bools 或者（最好）一个迭代器。

我试过显而易见的方法：

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

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

是否有建议的方法？

第二个问题：我还需要对我的稀疏矩阵的视图执行此操作 - 是否类似于 x = view(a,:,:); cols = unique(x.parent.colptr[x.indices[:,2]]) 或者是否有专门的功能？稀疏矩阵的视图似乎很棘手（cf https://discourse.julialang.org/t/slow-arithmetic-on-views-of-sparse-matrices/3644 – 不是交叉post）

非常感谢！

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 毫秒。它也应该比有问题的方法更快（除了正确性问题）。

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

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

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

Answer 2

如果索引不是范围，转换为稀疏矩阵的后备方法有效，但这里是索引的版本，它是向量。如果索引是混合的，则需要制作另一组版本。相当重复，但这是 Julia 的强项，当版本完成后，代码将使用调用者中的类型正确选择优化方法，而不需要太多努力。

function sortedintersecting(v1, v2)
    i,j = start(v1), start(v2)
    while i <= length(v1) && j <= length(v2)
        if v1[i] == v2[j] return true
        elseif v1[i] > v2[j] j += 1
        else i += 1
        end
    end
    return false
end

function nzcols(b::SubArray{T,2,P,Tuple{Vector{Int64},Vector{Int64}}}
  ) where {T,P<:SparseMatrixCSC}
    brows = sort(unique(b.indexes[1]))
    return [k 
      for (k,i) in enumerate(b.indexes[2])
      if b.parent.colptr[i]<b.parent.colptr[i+1] && 
        sortedintersecting(brows,b.parent.rowval[nzrange(b.parent,i)])]
end

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

-- 附录 --

由于注意到 nzrows 用于 Vector{Int} 索引有点慢，这是通过将 findin 替换为利用排序的版本来提高其速度的尝试：

function findin2(inds,v,w)
    i,j = start(v),start(w)
    res = Vector{Int}()
    while i<=length(v) && j<=length(w)
        if v[i]==w[j]
            push!(res,inds[i])
            i += 1
        elseif (v[i]<w[j]) i += 1
        else j += 1
        end
    end
    return res
end

function nzrows(b::SubArray{T,2,P,Tuple{Vector{Int64},Vector{Int64}}}
  ) where {T,P<:SparseMatrixCSC}
    active = falses(length(b.indexes[1]))
    inds = sortperm(b.indexes[1])
    brows = (b.indexes[1])[inds] 
    for c in b.indexes[2]
      active[findin2(inds,brows,b.parent.rowval[nzrange(b.parent,c)])] = true
    end
    return find(active)
end

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

Identify which rows (or columns) have values in sparse Matrix

sparse-matrix

julia