Looking for code of sparse matrix * dense matrix

Hi ! I may have found an issue in Julia 1.6 when doing matrix product between a sparse matrix and a full matrix when function zero() and one() and operators + and * have been changed for a particular algebra. I’m searching inside the code source of Julia in julia/stdlib/SparseArrays/src/sparsematrix.jl for this specific matrix product but I did not find it. Can somebody show me the function name. I’ll analyze the code before opening an issue.

Thanks in advance !

julia> A
2×2 SparseMatrixCSC{MP{Float64}, Int64} with 4 stored entries:

# Wrong answer (same for full(A)*A)
julia> A*full(A)

# Good answer
julia> A*A

Can you give a MWE (definition of MP and associated methods?). You can find the method signature with @which, which results in the method *(A::SparseMatrixCSCUnion{TA}, B::AdjOrTransDenseMatrix) where {TA} in stdlib/v1.7/SparseArrays/src/linalg.jl

A debugger is also very useful for these purposes. For example, in VS Code type @enter A*B, then step inside functions until you find the code that you are interested in.

Hi ! Thanks for the link I’ll watch.

Here is the miminal code extracted from my project. I’ll try with the debugger :slight_smile: but it seems to me this problem has been introduced since Julia 1.4.2 because I tested on versions 0.7 1.2.x 1.3.x and got the correct answers.

MP means Max-Plus: a tropical algebra (see here for more information http://www.stanczyk.pro/mpa/download.php?id=max-plus-1.7&ex=pdf )

This code redefines operators + as operator max() and operator * as operator *. We define neutral and absorbing elements zeros() as -Inf and one() as 0. I did not include conversion and promotion (not used in this basic example). The magic in Julia is that all complex algebra functions are automatically working … when there is no regression introduced by newer version :slight_smile:

using LinearAlgebra, SparseArrays

struct MP
    v::Float64
end

Base.zero(::Type{MP}) = MP(typemin(Float64))
Base.zero(x::MP) = zero(typeof(x))

Base.one(::Type{MP}) = MP(zero(Float64))
Base.one(x::MP) = one(typeof(x))

Base.:(+)(x::MP,   y::MP)   = MP(max(x.v, y.v))
Base.:(+)(x::MP,   y::Real) = MP(max(x.v, y))
Base.:(+)(x::Real, y::MP)   = MP(max(x, y.v))

Base.:(*)(x::MP,   y::MP)   = MP(x.v + y.v)
Base.:(*)(x::MP,   y::Real) = MP(x.v + y)
Base.:(*)(x::Real, y::MP)   = MP(x + y.v)

#Not needed for this example:
#MP(x::MP) = MP(x.v)
#Base.promote_rule(::Type{MP}, ::Type{U}) where {U} = MP
#Base.convert(::MP, x::Number) = MP(x)

This code gives the good result:

julia> A=[MP(1) MP(2); MP(3) MP(4)]
julia> A*A
2×2 Array{MP,2}:
 MP(5.0)  MP(6.0)
 MP(7.0)  MP(8.0)

This one does not give the correct answer since Julia 1.4.2 (1.5.x, 1.6.0):

julia> sparse(A)*A
2×2 Matrix{MP}:
 MP(6.0)  MP(7.0)
 MP(8.0)  MP(9.0)

When commenting

#Base.:(+)(x::MP,   y::Real) = MP(max(x.v, y))
#Base.:(+)(x::Real, y::MP)   = MP(max(x,   y.v))
#Base.:(*)(x::MP,   y::Real) = MP(x.v + y)
#Base.:(*)(x::Real, y::MP)   = MP(x   + y.v)

I’ll have this interesting error:

ERROR: MethodError: no method matching *(::MP, ::Bool)
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
  *(::StridedArray{P, N} where N, ::Real) where P<:Dates.Period at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Dates/src/deprecated.jl:44
  *(::Union{SparseVector{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false} where var"#s832"<:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, SubArray{Tv, 1, var"#s832", Tuple{Base.Slice{Base.OneTo{Int64}}}, false} where var"#s832"<:AbstractSparseVector{Tv, Ti}} where {Tv, Ti}, ::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsevector.jl:1448
  ...
Stacktrace:
 [1] mul!(C::Matrix{MP}, A::SparseMatrixCSC{MP, Int64}, B::Matrix{MP}, α::Bool, β::Bool)
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:41
 [2] *(A::SparseMatrixCSC{MP, Int64}, B::Matrix{MP})
   @ SparseArrays /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:53
 [3] top-level scope
   @ REPL[12]:1

This is not good since we are not supposed to call *(::MP, ::Bool) bool is not for Max-Plus algebra. Note that we may want to write MP(1.0) + 3.0 just for the user to avoid typing too many characters 3.0 meaning MP(3.0)) but MP(1.0) + false has nonsense.

By the way this is a second regression I found. Identity matrix that used to be created with the eye() function was made like this:

[ one() zero() zero()
  zero() one() zero()
  zero() zero() one()]

Since eye() no longer exist, replaced in Julia 1.x (I guess) by I which is not working because AGAIN with introduction boolean instead of zero() and one()

julia> I
UniformScaling{Bool}
true*I

and therefore Identity are no longer working in Max-Plus. I succeeded to make a workaround in my package but I have to find a real fix. See Matrix identity not calling one() ? · Issue #33036 · JuliaLang/julia · GitHub

Ah, the problem is that sparse matrices assumes that false will work as a strong zero (false+x=x, and false*x=false) for your type. It would be a bit weird, but defining
MP(x::Bool)= x ? one(MP) : zero(MP)
should fix it.

It does not seem to fix. Tomorrow, I’ll try to diff the julia code between 1.3.x and 1.4.x and see if something odd has changed.

Regarding your specific question (locate the source code for the dense times sparse matrix product):
https://github.com/JuliaLang/julia/blob/e3adc399e6bf69d0801c243f254ff8624f50d3c4/stdlib/SparseArrays/src/linalg.jl#L87-L115

Regarding the broader issue (implement the max-plus algebra): I would recommend against implementing Number <-> MP conversion and mixed arithmetic. On a more philosophical level, the problem with mixed arithmetic is that it is not clear whether MP(1) + 2 should be evaluated according to standard or max-plus arithmetic. In more practical terms, the issue is that the Julia standard library implicitly assumes standard arithmetic in at least a couple of places. This leads to silent errors if you do implement mixed arithmetic, whereas if you don’t then you get an error message which is much easier to detect and address.

1 Like

@attersi Thank you for your search, I did not get time to dive into code since hard for me to understand and missing comments.

For the moment since sparse * sparse is working I made a fallback for dense matrix * sparse where I convert the dense to sparse. This is allowing me to get the same results as in previous Julia versions. Once my work done i’ll get back investigating on this regression.

Concerning Number <-> MP: contamination is fine since working with MP number make you dive in (max,+) algebra and no more in classic algebra. For classic addition: MP.v + scalar is okay