A sparse matrix cannot be indexed by a SVector

General Usage
,
1

As indicated by the topic, I intend to replace some small matrices with the StaticArrays but the global sparse matrix created by SparseArrays.jl cannot be indexed by a SVector. A simple example is given as follows.

using SparseArrays
using StaticArrays

A = spzeros(10, 10)
ind1 = SVector(1, 2, 3, 4)
a = @SArray rand(4, 4)

A[ind1, ind1] += a

@show A

The error is given in REPL as

julia> A[ind1,ind1] += a
ERROR: MethodError: no method matching length(::StaticArrays.StaticIndexing{SVector{4, Int64}})

Closest candidates are:
  length(::Union{Base.KeySet, Base.ValueIterator})
   @ Base abstractdict.jl:58
  length(::Union{SparseArrays.FixedSparseVector{Tv, Ti}, SparseVector{Tv, Ti}} where {Tv, Ti})
   @ SparseArrays C:\lijinze\julia\share\julia\stdlib\v1.9\SparseArrays\src\sparsevector.jl:95
  length(::Union{ZMQ._Message, Base.RefValue{ZMQ._Message}})
   @ ZMQ C:\lijinze\julia\packages\packages\ZMQ\lrABE\src\_message.jl:31
  ...

Stacktrace:
 [1] setindex!(A::SparseMatrixCSC{Float64, Int64}, V::SMatrix{4, 4, Float64, 16}, Ix::SVector{4, Int64}, Jx::SVector{4, Int64})
   @ SparseArrays C:\users\julia\share\julia\stdlib\v1.9\SparseArrays\src\sparsematrix.jl:3170
 [2] top-level scope
   @ REPL[24]:1

Thanks in advance for tackling this problem.

2

The sparse matrix A doesn’t hold any values, except for the default value zero (of default type Float64) which isn’t explicitly stored. Accessing any range of indices gives some dots:

julia> A[ind1, ind1]
4×4 SparseMatrixCSC{Float64, Int64} with 0 stored entries:
  ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅

Accessing via a single index returns the default value 0.0:

julia> A[1]
0.0

You operation can be written as (see the dot before the +):

julia> A[ind1,ind1] .+= a
4×4 view(::SparseMatrixCSC{Float64, Int64}, [1, 2, 3, 4], [1, 2, 3, 4]) with eltype Float64 with indices SOneTo(4)×SOneTo(4):
 0.396905  0.518507  0.0333795  0.428079
 0.105497  0.358759  0.0329777  0.748697
 0.495962  0.948462  0.214173   0.55445
 0.517173  0.866702  0.91867    0.352342

julia> A
10×10 SparseMatrixCSC{Float64, Int64} with 16 stored entries:
 0.396905  0.518507  0.0333795  0.428079   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
 0.105497  0.358759  0.0329777  0.748697   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
 0.495962  0.948462  0.214173   0.55445    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
 0.517173  0.866702  0.91867    0.352342   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅         ⋅         ⋅          ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅         ⋅         ⋅          ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅         ⋅         ⋅          ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅         ⋅         ⋅          ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅         ⋅         ⋅          ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅         ⋅         ⋅          ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅

This is a bit of a lack in the documentation of Sparse Arrays · The Julia Language I would say, where it says:

Arithmetic operations on sparse matrices also work as they do on dense matrices.

The error message points to a forgotten implementation of method length(::StaticArrays.StaticIndexing{SVector{4, Int64}}) somewhere, but I don’t have the nerves to look this up for now, and it wouldn’t help you much for your issue.

3

Thank you very much! It succeeds in addressing my problem.

1 Like
4

This is basically an incorrect return value from to_indices in StaticArrays.jl:

It’s overriding internals Base._getindex but that doesn’t actually work.

StaticIndexing needs to be <: AbstractArray and it needs to define iterate and size (shockingly not getindex but that’s another story…).

Please make an issue for this in StaticArrays.jl

4 Likes
5

From the to_indices docs:

“The returned tuple must only contain either Ints or AbstractArrays of scalar indices that are supported by array A”

StaticArrays.jl is just breaking the interface here.

Issue made here:

3 Likes