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 倍，而是慢了 8000 倍。原因是 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