Sparse with default value "Inf"

Hi guys i’m working on a problem which requires solving many linear assignment problems (LAP) and i wonder if it’s possible to declare a sparse matrix object in Julia which defaults to “Inf” rather than 0, as this would turn into a massive optimization in storage for the nature of my problem.

So, is there any way to declare a Sparse Array which defaults to some other value?

Thank you :slight_smile:

You could create your own type to do this.

Perhaps for prototyping, you could just use a Dict to store which values are not -Inf, so that the bookkeeping is easier than CSC — if I understand correctly, you will not be needing any linear algebra anyway.

1 Like

See the discussion in

In particular, as suggested in that issue, you can easily create your own type wrapping SparseArray that returns a different default (or use some other data structure as suggested above).

3 Likes

Maybe you could store 1./A, so that your Infs are 0s?

3 Likes

Even the approach I described there wouldn’t work with a default value of Inf since you can’t add anything to Inf and get a finite value.

1 Like

i’m working with Dict’s in the end.
it is quite sad though, as it is pretty common to use a different default for sparse matrixes.

Really? What sparse matrix library supports this?

3 Likes

You can also define your own type that redefines zero and operations (this idea is from @andreasnoack originally I think), e.g. something like:

struct SemiRing{T,A,M,Z,O} <: Number
        val::T
end

Base.:+(x::SemiRing{T, A, M, Z, O}, y::SemiRing{T, A, M, Z, O}) where {T, A, M, Z, O} = SemiRing{T, A, M, Z, O}(A(x.val, y.val))
Base.:*(x::SemiRing{T, A, M, Z, O}, y::SemiRing{T, A, M, Z, O}) where {T, A, M, Z, O} = SemiRing{T, A, M, Z, O}(M(x.val, y.val))
Base.zero(::Type{SemiRing{T, A, M, Z, O}}) where {T, A, M, Z, O} = SemiRing{T, A, M, Z, O}(Z)
Base.zero(::SemiRing{T, A, M, Z, O}) where {T, A, M, Z, O} = SemiRing{T, A, M, Z, O}(Z)
Base.one(::Type{SemiRing{T, A, M, Z, O}}) where {T, A, M, Z, O} = SemiRing{T, A, M, Z, O}(O)
Base.one(::SemiRing{T, A, M, Z, O}) where {T, A, M, Z, O} = SemiRing{T, A, M, Z, O}(O)

then you use it as

S = SemiRing{Float64,max,+,-Inf,0.0}
using SparseArrays
x=sparse(1:10, 1:10, ones(S,10))
1 Like

@StefanKarpinski
this for example
https://reference.wolfram.com/language/ref/SparseArray.html

@abraunst
Thank you so much, i already prototyped most of the code needed, and it is using dictionaries, for now i’m focusing on correctness, as soon as i finish prototyping the whole algorithm i’ll try to make it work on Sparse SemiRings

1 Like

One can think of <:AbstractMatrix types as

  1. a rectangular table of data
  2. a matrix (ie as understood in linear algebra)

It is my understanding that the SparseArrays library targets the latter, so I am curious what the application is for this use case. Do you have some examples and references?

If it is not about linear algebra, just a kind of compressed storage, then it would make sense to define a subtype of AbstractArray, along the lines of various packages in

1 Like

Do you have some examples and references?

The applications in linear programming are endless, for instance, my problem is solving many very large linear assignment problems (with some modifications w.r.t the various hungarian methods already at disposal) which have many blocking costs (typical size: 200000x200000 matrix with ~ 99% of the values == Inf),

If it is not about linear algebra, just a kind of compressed storage, then it would make sense to define a subtype of AbstractArray , along the lines of various packages in

I would not do any matmul or factorization on this matrices, the things i would use are the broadcasting machinery, sliced views, indexing, masked indexing.

that’s all, if i had the time and the knowledge, i would make package, and maybe i will try to do it after i finish working out the actual algorithms :slight_smile: for my work
thank you all

I am not an expert in those algorithms, but I would assume that they make some use of the sparsity information directly, so it might make sense to just provide that in the interface.

Eg given a Ax \le b constraint, A could be provided as (row, column, value) triplets for finite entries with the understanding that it is \infty everywhere else, instead of going to the somewhat roundabout solution of wrapping that information in a matrix type, then extracting it.

2 Likes

Okay guys, i’m reworking my code now to use Sparse Matrices, containing a custom type, which allows me to redefine zero and get past the original issue,
Thank you again @abraunst

Now my other issue, is that if i broadcast some calculation involving the elements which are not stored, even though this calculation produces a “0” this zero is stored in the sparse matrix, i do not want to use dropzeros!, is there some way to avoid storing “0” after a broadcast?

because i know for sure that all the operations i allow on this custom sparse type, will never affect "0"s

by the way i want to thank everybody, every answer has slowly steered me in the right direction, thank you guys

2 Likes

Can you produce a MWE? For instance, following my example

julia> S = SemiRing{Float64,max,+,-Inf,0.0}
SemiRing{Float64,max,+,-Inf,0.0}
julia> x=sparse(1:10, 1:10, ones(S,10));
julia> (t->S(-Inf)).(x)
10×10 SparseMatrixCSC{Any,Int64} with 0 stored entries

seems to work as expected.

EDIT: the eltype being Any seems a bug?

1 Like

MWE, i hope that there is a way to solve this :slight_smile:, thank you again

struct Weight
    w::Float64
end
import Base.-
Base.:-(x::Weight, y::T) where {T<:Number} = Weight(x.w-y)
Base.:-(y::T,x::Weight) where {T<:Number} = Weight(y-x.w)
import Base.convert
Base.convert(::Type{Weight}, y::T) where {T<:Number} = Weight(y)
import Base.zero
Base.zero(::Type{Weight}) = Weight(Inf)
Base.zero(::Weight) = Weight(Inf)


#TESTS
using SparseArrays
A=spzeros(Weight,10,10)
@show typeof(A)
A.=Inf #Works
snumzeros=size(A.nzval)
@show snumzeros
@show typeof(A)
A.-=0   #ALLOCATES UNNECESSARY ZEROS
snumzeros=size(A.nzval)

@show snumzeros
@show typeof(A)
;
typeof(A) = SparseMatrixCSC{Weight,Int64}
snumzeros = (0,)
typeof(A) = SparseMatrixCSC{Weight,Int64}
snumzeros = (100,)
typeof(A) = SparseMatrixCSC{Weight,Int64}

Would MappedArrays.jl work for this?

Unfortunately it wont help

I seem to recall that our suggestion was the other way around: define a custom array type that wraps around a SparseMatrixCSC, not a SparseMatrixCSC that contains a custom number type.

zero is supposed to be the additive identity. Defining a method that fundamentally changes the meaning of this function will break all generic code that uses zero, and will be a continuing headache.

2 Likes

I seem to recall that our suggestion was the other way around: define a custom array type that wraps around a SparseMatrixCSC , not a SparseMatrixCSC that contains a custom number type.

i see the deep uglyness in asserting W(\infty)=0, unfortunately i do not know how to start doing what you suggested.

I still want a Sparse Matrix, which is sparsified w.r.t. a custom value.

I think it is not wrong for a user-defined type (in which the addition can be redefined, e.g. as min).

1 Like