Julia：稀疏矩阵的视图

Question

我对 view 如何作用于 Julia 中的稀疏矩阵感到很困惑：

using LinearAlgebra, SparseArrays, BenchmarkTools
v = SparseVector(1000, [212,554,873], [.3, .4, .3]);
A = sparse(diagm(rand(1000)));  # same effect observed for non-diag matrix
B = view(A, :, :);

julia> @btime A*v;
  1.536 μs (4 allocations: 23.84 KiB)

julia> @btime B*v;
  11.557 ms (5 allocations: 288 bytes)

B*v 的内存占用似乎要低得多，但比 A*v 慢 8000 倍。发生了什么以及导致这些性能差异的原因是什么？

Answer 1

2021 年 6 月更新：现在实现了下面提到的稀疏视图缺少的专用算法，因此性能很多更合理天（Julia 1.6+）：

julia> @btime A*v;
  2.063 μs (4 allocations: 23.84 KiB)

julia> @btime B*v;
  2.836 μs (9 allocations: 25.30 KiB)

你可以看到，是的，确实，这是对 *:

使用稀疏特化

julia> @which B*v
*(A::Union{SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 2, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}, I}, L} where {I<:AbstractUnitRange, var"#s832"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, L}, LowerTriangular{Tv, var"#s831"} where var"#s831"<:Union{SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 2, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}, I}, L} where {I<:AbstractUnitRange, var"#s832"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, L}}, UpperTriangular{Tv, var"#s832"} where var"#s832"<:Union{SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 2, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}, I}, L} where {I<:AbstractUnitRange, var"#s832"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, L}}} where {Tv, Ti}, B::AbstractSparseVector{Tv, Ti} where {Tv, Ti}) in SparseArrays at /Users/mbauman/Julia/release-1.6/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:163

以前，该方法没有实现，这意味着发生了以下情况：

旧答案：不是慢 8 倍，而是 8000x。原因是因为 Julia 使用多重分派来使用专门的算法，这些算法可以利用矩阵和向量的稀疏存储来 完全避免 在它知道只会为零的数组部分上工作。您可以看到使用 @which:

调用了哪个算法

julia> @which A*v
*(A::SparseArrays.AbstractSparseMatrixCSC, x::AbstractSparseArray{Tv,Ti,1} where Ti where Tv) in SparseArrays at /Users/mbauman/Julia/master/usr/share/julia/stdlib/v1.4/SparseArrays/src/sparsevector.jl:1722

julia> @which B*v
*(A::AbstractArray{T,2}, x::AbstractArray{S,1}) where {T, S} in LinearAlgebra at /Users/mbauman/Julia/master/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/matmul.jl:50

前者使用高度专业化的稀疏实现，而后者使用稍微更通用的接口，也可以支持视图。现在，理想情况下，我们会检测像 view(A, :, :) 这样的微不足道的情况，并将它们专门化以有效地相同，但请注意，在一般情况下，可能无法保留矩阵的稀疏性和结构：

julia> view(A, ones(Int, 1000), ones(Int, 1000))
1000×1000 view(::SparseMatrixCSC{Float64,Int64}, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) with eltype Float64:
 0.306159  0.306159  0.306159  0.306159  …  0.306159  0.306159  0.306159
 0.306159  0.306159  0.306159  0.306159     0.306159  0.306159  0.306159
 0.306159  0.306159  0.306159  0.306159     0.306159  0.306159  0.306159
 ⋮                                       ⋱
 0.306159  0.306159  0.306159  0.306159     0.306159  0.306159  0.306159
 0.306159  0.306159  0.306159  0.306159     0.306159  0.306159  0.306159
 0.306159  0.306159  0.306159  0.306159  …  0.306159  0.306159  0.306159

Answer 2

与问题没有直接关系，但对于所有想要以稀疏格式获取 look/view/examine 和 matrix/array 的人，您可以使用 display(YourSparseArray).

这个答案的动机是这个问题的观点相对较多，我认为这些观点并不都关心 view 的表现，而是对如何看到稀疏的矩阵很好。

Julia：稀疏矩阵的视图

Julia: view of sparse matrix

performance

benchmarking

sparse-matrix

julia