朱莉娅：稀疏矩阵的视图

Question

朱莉娅：稀疏矩阵的视图

performancejuliasparse-matrixbenchmarking

5

我对 Julia 中 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 的内存占用似乎要低得多，但速度比 A*v 慢了 8000 倍。是什么原因造成这些性能差异？

- mattu

2个回答

-1

虽然与问题不直接相关，但是对于所有想要以稀疏格式查看/检查矩阵/数组的人，您可以使用display(YourSparseArray)。

这个答案是受到这个问题相对较多的浏览量的启发，我认为并不是所有人都关心view的性能，而是更感兴趣如何美观地查看稀疏矩阵。

- Dan Doe

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- mbauman · Accepted Answer

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