# Is working with SubArrays of a SparseMatrix a bad idea?
**URL:** https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000
**Category:** General Usage
**Tags:** question, bug, sparse, views
**Created:** [November 9, 2022, 6:10pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000 "2022-11-09T18:10:45Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [November 9, 2022, 6:10pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/1 "2022-11-09T18:10:45Z")
I’ve been trying to work with `view`s of a `SparseMatrix` for a while now. The goal is to avoid extra copies when you want to operate inplace on a subsection of the sparse matrix. I’ve run into enough issues that I’m starting to think that this is simply a bad idea. The issues fall into two categories: 1) a valid method doesn’t exist for a SubArray of a SparseArray or 2) _some valid method_ exists, but the result for SparseArray is terribly slow which more than negates the advantages of using a sparse array in the first place.
Examples of 1) are [#286](https://github.com/JuliaSparse/SparseArrays.jl/issues/286#issue-1438771133) and [#16](https://github.com/JuliaSparse/SparseArrays.jl/issues/16). Obviously these can be or have been fixed.
There are lots of examples of 2) but a very simple one is that there is no optimized method for `sum` for `SubArray{T, N, AbstractSparseArray}`, which has a major performance impact.
```julia
using SparseArrays, BenchmarkTools
a = sprand(1000, 1000, .01)
aview = view(a, :, :)
amat = Matrix(a)
sum2(x::SubArray{T, N, <:AbstractSparseArray}) where {T, N} = sum(nonzeros(x))
@btime sum($a) # 1.604 μs (0 allocations: 0 bytes)
@btime sum($amat) # 182.875 μs (0 allocations: 0 bytes)
@btime sum($aview) # 4.646 ms (0 allocations: 0 bytes)
@btime sum2($aview) # 1.642 μs (0 allocations: 0 bytes)
```
I don’t intend this to be criticism of the SparseArrays developers - they’ve been awesome. The package itself is quite good and mature. With enough time and developer effort I know all this can be fixed.
However, supporting operations on `SubArray` in Julia in general seems extremely cumbersome for developers since you don’t get things like `SubArrray{T, 1, SomeArray} <: SomeArray{T, 1}` by default. This seems to result in the need for tons of custom methods, type unions, etc that I know have been challenging for package developers to keep up with. E.g. “Why don’t you support operation X efficiently on a view of a reshape of a view of an UpperTriangular of your special array?”
I’m curious to know how others have approached these sorts of issues. I’m also interested in general advice for working with Julia in this context - even if it’s just to keep posting issues until it all works.
---
### Author:  [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)
#### Post date: [November 9, 2022, 7:08pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/2 "2022-11-09T19:08:18Z")
Sparse matrix internal structure is optimized for operations, especially multiplication, rather than arbitrary indexing [Sparse matrix - Wikipedia](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_row_(CSR,_CRS_or_Yale_format)) and [Sparse matrix - Wikipedia](https://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_(CSC_or_CCS)). Creating and/or using a subarray requires such arbitrary indexing. Operations on the subarray can no longer take advantage of the efficient structure of the original sparse matrix. They are likely to pass through an indexing stage. If you have several operations to perform with the subarray, then you should create a new sparse matrix for the subarray.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 9, 2022, 7:25pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/3 "2022-11-09T19:25:32Z")
Forget about the timing! I think there is a correctness bug here:
```julia
julia> sum(@view a[:,1:5])
25.236381449783398
julia> sum(nonzeros(@view a[:,1:5]))
5140.889898161122
```
The bug itself is quite obvious, in `sparsematrix.jl`:
```julia
nonzeros(S::SparseMatrixCSCView) = nonzeros(S.parent)
```
`nonzeros` goes directly to parent without filtering output to view.
(and the next two `nonzeros` specializations for Triangular matrices look fishy as well).
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [November 9, 2022, 7:36pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/4 "2022-11-09T19:36:03Z")
Is this true if you limit yourself to views that include whole columns or groups of adjacent columns?
---
### Author:  [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)
#### Post date: [November 9, 2022, 7:38pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/5 "2022-11-09T19:38:08Z")
That depends on the internals of the view, which I don’t know about, but those could be fast(er).
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [November 9, 2022, 7:38pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/6 "2022-11-09T19:38:39Z")
I see that now too.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 9, 2022, 7:45pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/7 "2022-11-09T19:45:04Z")
The right code should be more along the lines of:
```julia
import SparseArrays: nonzeros
function nonzeros(S::SparseArrays.SparseMatrixCSCView)
s = S.parent.colptr[first(S.indices[2])]
e = S.parent.colptr[last(S.indices[2])+1]
if S.indices[1]==axes(S.parent,1)
return @view S.parent.nzval[s:e-1]
else
res = eltype(S.parent)[]
for i in s:e-1
if S.parent.rowval[i] in S.indices[1]
push!(res, S.parent.nzval[i])
end
end
return res
end
end
```
And as some other `nonzeros` definitions are suspicious, we have for example:
```julia
julia> ut = UpperTriangular(sprand(5,5,0.5))
5×5 UpperTriangular{Float64, SparseMatrixCSC{Float64, Int64}}:
0.0 0.0 0.919615 0.364471 0.0
⋅ 0.268043 0.0 0.254744 0.441155
⋅ ⋅ 0.0 0.0 0.0
⋅ ⋅ ⋅ 0.935346 0.792978
⋅ ⋅ ⋅ ⋅ 0.0
julia> sum(ut)
3.97635045952992
julia> sum(nonzeros(ut))
7.379883839287259
```
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [November 9, 2022, 7:54pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/8 "2022-11-09T19:54:30Z")
The current package does appears to work for 1D views:
> julia\> sum(a[:,1])
> 5.826271341749841
> julia\> sum(view(a, :, 1))
> 5.826271341749841
> julia\> sum(nonzeros(view(a, :, 1)))
> 5.826271341749841
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 9, 2022, 8:02pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/9 "2022-11-09T20:02:38Z")
The problem is the specializations in SparseArray package.
So it would pop-up only on sparse matrix views.
Here is a version for sparse upper triangular matrices of `nonzeros`:
```julia
function nonzeros(S::SparseArrays.UpperTriangular{<:Any,<:SparseArrays.SparseMatrixCSCUnion})
res = eltype(S.data)[]
c = 1
nextc = S.data.colptr[c+1]
for i in eachindex(S.data.nzval)
while i >= nextc
c += 1
nextc = S.data.colptr[c+1]
end
if S.data.rowval[i]<=c
push!(res, S.data.nzval[i])
end
end
return res
end
```
Lower triangular version is identical with `<=c` replaced with `>=c`.
While the Triangular versions are fast. In case of sparse matrices, there should be a debate about forcing triangle structure in the underlying representation (as performance considerations different than dense matrices).
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [November 9, 2022, 8:16pm UTC](https://discourse.julialang.org/t/is-working-with-subarrays-of-a-sparsematrix-a-bad-idea/90000/10 "2022-11-09T20:16:16Z")
Looks related to [https://github.com/JuliaSparse/SparseArrays.jl/issues/64](https://github.com/JuliaSparse/SparseArrays.jl/issues/64)