Sparse array filled with constant in Julia

Question

I have a matrix, M, with very simple form: M[i,j]=a if i==j and M[i,j]=b everywhere else. M is very large (>10000) so I cannot initialize it or store it without using some sort of sparse matrix. It seems that the best way to store this matrix would be to use some sort of sparse matrix where non-entries are set to b instead of zero.

I have tried

M = spdiagm((a-b), N, N) .+ b

But doing it this way stores M as a 10 by 10 matrix with 100 entries which seems to mean that there is no compression.

Is there a better way to initialize this matrix?

Answer 1

It depends on what you want to do with your matrix, but if you are only interested in performing matrix-vector products, you can use the FillArrays.jl package in conjunction with the LinearMaps.jl package:

julia> using LinearMaps, FillArrays

julia> n=100 # can be bigger
100

julia> a=2.0
2.0

julia> b=3.0
3.0

julia> M=(a-b)*LinearMap(Eye(n,n))+b*LinearMap(Ones(n,n))
100×100 LinearMaps.LinearCombination{Float64} with 2 maps:
  100×100 LinearMaps.ScaledMap{Float64} with scale: -1.0 of
    100×100 LinearMaps.WrappedMap{Float64} of
      100×100 Eye{Float64}
  100×100 LinearMaps.ScaledMap{Float64} with scale: 3.0 of
    100×100 LinearMaps.WrappedMap{Float64} of
      100×100 Ones{Float64}

Then you can perform operations like

julia> x=rand(n);

julia> M*x;

julia> x'*M;
    
julia> M*M
100×100 LinearMaps.CompositeMap{Float64} with 2 maps:
  100×100 LinearMaps.LinearCombination{Float64} with 2 maps:
    100×100 LinearMaps.ScaledMap{Float64} with scale: -1.0 of
      100×100 LinearMaps.WrappedMap{Float64} of
        100×100 Eye{Float64}
    100×100 LinearMaps.ScaledMap{Float64} with scale: 3.0 of
      100×100 LinearMaps.WrappedMap{Float64} of
        100×100 Ones{Float64}
  100×100 LinearMaps.LinearCombination{Float64} with 2 maps:
    100×100 LinearMaps.ScaledMap{Float64} with scale: -1.0 of
      100×100 LinearMaps.WrappedMap{Float64} of
        100×100 Eye{Float64}
    100×100 LinearMaps.ScaledMap{Float64} with scale: 3.0 of
      100×100 LinearMaps.WrappedMap{Float64} of
        100×100 Ones{Float64}

The matrices are stored using a "lazy" approach and not expanded in the memory, hence n can be big. However one must take care that M is not an AbstractMatrix :

julia> supertype(typeof(M))
LinearMap{Float64}