Юлия: вид разреженной матрицы - PullRequest
Новая
       36

Юлия: вид разреженной матрицы

3 голосов
/ 04 ноября 2019

Я довольно смущен тем, как 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. Что происходит и чем обусловлены эти различия в производительности?

1 Ответ

5 голосов
/ 04 ноября 2019

Это не в 8 раз медленнее, а 8000x медленнее. Причина в том, что Джулия использует множественную диспетчеризацию для использования специализированных алгоритмов, которые могут использовать редкое хранилище матрицы и вектора, чтобы полностью избегать работы с секциями массива, которые, как он знает, будут просто равны нулю. Вы можете увидеть, какой алгоритм вызывается с помощью @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
...