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 
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.
See the discussion in
https://github.com/JuliaLang/julia/issues/10410
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).
Maybe you could store 1./A, so that your Infs are 0s?
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.
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?
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))
@StefanKarpinski
this for example
@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
One can think of <:AbstractMatrix types as
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
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 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.
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
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?
MWE, i hope that there is a way to solve this , 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.
I seem to recall that our suggestion was the other way around: define a custom array type that wraps around a
SparseMatrixCSC, not aSparseMatrixCSCthat 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).