朱莉娅：稀疏矩阵的视图

Question

我view对朱莉娅对稀疏矩阵的行为感到困惑不解：

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似乎占用的内存少得多，但比慢8000倍A*v。这是怎么回事，是什么导致这些性能差异？

Answer 1

它不会慢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