Matrix identity not calling one() ?

[Lecrapouille]

Hi, I'm trying to implement the max-plus semiring algebra in Julia 1.0.3 like the follow code:

struct MP{T} <: Number val::T end
Base.:+(a::MP, b::MP) = MP(max(a.val, b.val))
Base.:*(a::MP, b::MP) = MP(a.val + b.val)
Base.show(io::IO, k::MP) = show(io, k.val)
Base.convert(::MP{T}, x) where T = MP(T(x))
Base.zero(::MP{T}) where T = MP(typemin(T))
Base.one(::MP{T}) where T = MP(zero(T))
Base.zero(::Type{MP{T}}) where T = MP(typemin(T))
Base.one(::Type{MP{T}}) where T = MP(zero(T))
mparray(A::Array) = map(MP, A)

The following code gives the expected behavior :

julia> Base.zero(MP{Float64})
-Inf

julia> Base.one(MP{Float64})
0.0

But the following code does not gives the expected behavior :

julia> Matrix{MP{Float64}}(I, 2, 2)
2×2 Array{MP{Float64},2}:
  1.0  -Inf
 -Inf   1.0

What I was expected:

julia> Matrix{MP{Float64}}(I, 2, 2)
2×2 Array{MP{Float64},2}:
  0.0  -Inf
 -Inf   0.0

I think Matrix{MP{Float64}}(I, 2, 2) creates a matrix initialized by Base.zero(MP{Float64}) calling MP(typemin(Float64)) giving -Inf but the I does not create the diagonal initialized by one(MP{Float64}) which would called MP(zero(Float64)) giving 0.0.

I followed this link JuliaLang/julia#30298 so fortunatly I can fix by writting code like this:

julia> Matrix{MP{Float64}}(0I, 2, 2)
2×2 Array{MP{Float64},2}:
  0.0  -Inf
 -Inf   0.0

Is it because I use Bool that I have this behavior ? And how to replace code 0I by one(MP{Float64})I ? I really think that Julia shall implement I calling one() and without giving Type::T calls one(Bool).

I'm starting learning julia! Thx!

[fredrikekre]

I is a UniformScaling, which may have any value on the diagonal (not just 1) so we just convert that value to the correct type. In this case we call MP{Float64}(true), see https://github.com/JuliaLang/julia/blob/cd16f6e74cdb23f9986a387aef29b2f0631f5653/stdlib/LinearAlgebra/src/uniformscaling.jl#L386.
Consider using https://discourse.julialang.org/ for further questions.

[Lecrapouille]

Thank but it was not a discourse question, what you said prooved me it's a bug (in the way not mathematicaly consistent). As I said A is created with zero() which is correct. But s.λ was not init with one(T) therefore it does not contain the correct value: s.λ shall be one(MP{Float64}) = 0.0 = false. Look:

julia> UniformScaling{MP{Float64}}(one(MP{Float64}))
UniformScaling{MP{Float64}}
0.0*I

and the doc says: An object of type UniformScaling, representing an identity matrix of any size.
So I thinking this way I = [one(T) zero(T) ....; zero(T) one(T) zero(T) ... ]

[Lecrapouille]

PS: May I rephrase to be clear. When you said I is a UniformScaling, which may have any value on the diagonal (not just 1) I agree with you for:

• 0.0I implies UniformScaling{MP{Float64}}(0.0)
• 1.0I implies UniformScaling{MP{Float64}}(1.0)

But I disagree when:

• I implies UniformScaling{MP{Float64}}(1.0) it shall call UniformScaling{MP{Float64}}(one(MP{Float64})) because the doc says representing an identity matrix

[Lecrapouille]

@fredrikekre, @andreasnoack I'm sorry to disturb you again (this will be my last comment, up to you to take into account) but reading a little more the code and comments I comfort myself that: https://github.com/JuliaLang/julia/blob/cd16f6e74cdb23f9986a387aef29b2f0631f5653/stdlib/LinearAlgebra/src/uniformscaling.jl#L49 is not good idea.

By the way: codes I(3) and (0.7*I)(3) in the doc/comment are not working. I have ERROR: MethodError: objects of type UniformScaling{Bool} are not callable. Missing a constructor ?

I know max-plus algebra is disturbing: we can confuse between MP(1) = 1 (coming from MP(true)) and one() = MP(zero(T)) = 0 named %1 which is neutral for operators + and * in this algebra. But also zero() = MP(typemin(T)) = -Inf named %0 neutral for opertor+.

So ideally it should be const I = UniformScaling(one()) which is not possible because one() is template of T. In my case I can fix in other way by doing something awful like this:

# Copy constructor
MP{T}(x::MP) where T = MP{T}(x.val)

# and similar to uniformscaling.jl:
const J=UniformScaling(one(MP{Float64}))

# Now working
Matrix{MP{Float64}}(J, 2,2)
2×2 Array{MP{Float64},2}:
  0.0  -Inf
 -Inf   0.0

My idea would be: is it not possible to change I for a lambda function (or composition function such as methods(∘) ?) something like I = UniformScaling∘one ? But I don't know to formulate it in Julia.

Thx!

[StefanKarpinski]

It’s unclear what your objection is here. Can you clarify what the actual problem is?

[StefanKarpinski]

We’re fairly familiar with unusual algebras like max-plus: https://dspace.mit.edu/bitstream/handle/1721.1/115964/JuliaSemiring_HPEC2013_Paper%281%29.pdf

[andreasnoack]

The inconsistency is that in Matrix{T}(I, m, n) we use zero(T) when constructing the zeros but convert(T, I.λ) when constructing the diagonal. It would probably be more consistent to use convert(T, zero(I.λ)) for the off diagonals.

It might not be @Lecrapouille's preferred solution but it would be more consistent. Notice, that you are actually telling Julia to convert true to "one" when you define Base.convert(::MP{T}, x) where T = MP(T(x)) so maybe you'll need to reconsider that definition.

By the way: codes I(3) and (0.7*I)(3) in the doc/comment are not working. I have ERROR: MethodError: objects of type UniformScaling{Bool} are not callable. Missing a constructor ?

The docs you are reading must be for a more recent version of Julia. The operations were only introduced in the just released version 1.2.

[fredrikekre]

Alternatively, one(T) * I.λ maybe?

[Lecrapouille]

@StefanKarpinski ok I try to clarify

Identity matrices for max-plus (for example the case of a 2x2 matrix) are:

  0.0  -Inf
 -Inf   0.0

Julia 0.4.7, eye(MP{Float64},2,2) returns the correct matrix because it uses correctly one() and zero(). See my enclosed file, code source MaxPlus.jl.txt I made at this time and inspired by the PDF you linked.

But Julia >= 1.03 (>= 0.7 ?) returns an incorrect answer:

Matrix{MP{Float64}}(I, 2, 2) =
  1.0  -Inf
 -Inf   1.0

So this is clearly a regression.

From an algebra point of view:

Matrix{T}(I, 2, 2) = [one(T); zero(T); zero(T); one(T)]

And with your API:

Matrix{T}(nI, 2, 2) = n*Matrix{T}(I, 2, 2) = [n*one(T); n*zero(T); n*zero(T); n*one(T)]

In the current Julia implementation zero(T) are called correctly but one(T) is not called correctly.
The purpose of one() and zero() is not purely for programmers for initializing their data. They are here to define neutral/absorbing element for operators +, * when defining your own algebra.

Thanks to @fredrikekre who gave me the line of code which looks incorrect to me (still point of view of the algebra):

const I = UniformScaling(true)

Now, let suppose to replace it by:

const I = UniformScaling(one())

Unfortunatly, this line does not compile because of missing T. But just let suppose it compiles.

For boolean: UniformScaling(one()) -> UniformScaling(true) and we obtain the same I, right ?
Idem for Int and Float: UniformScaling(one()) -> UniformScaling(1) and UniformScaling(1.0)

So finally the code:

Matrix{T}(I, 2, 2) = [one(T); zero(T); zero(T); one(T)]

will produce for Bool, Int, Float correct values: [true false; false true], [1.0 0.0; 0.0 1.0]
but also correct values for maxplux [0.0 -Inf; -Inf 0.0]

For Matrix{T}(nI, 2, 2) = n*Matrix{T}(I, 2, 2) = [n*one(T); n*zero(T); n*zero(T); n*one(T)]. Example with n=2: Float: [2.0 0.0; 0.0 2.0] MP: [2.0 -Inf; -Inf 2.0] Because in max-plus operator* does the addition (2+0). This operator is correctly overloaded in Julia (as well as operator+ computing the maximun).

@andreasnoack, @fredrikekre concerning your proposal of fixes concerning:

function Matrix{T}(s::UniformScaling, dims::Dims{2}) where {T}
    A = zeros(T, dims)
    v = T(s.λ)
    for i in diagind(dims...)
        @inbounds A[i] = v
    end
    return A
end

@fredrikekre one(T) * I.λ is not correct because it will return 1.0 (instead of expected I.λ = one(T) = 0.0) because I.λ was init with the wrong value. I.λ = true => MP(true) = MP(1) = 1.0 and operator* does the addition operation => 0.0 + 1.0 = 1.0.

@andreasnoack convert(T, zero(I.λ)) is only correct for MP because zero(I.λ) will return false whatever the value of λ (wich reach my idea) => 0.0. Unfortunately not working for Int, Float. The matrix will be fully zeroed.

My proposal

I.λ shall contain the correct value and therefore init with one(T) of the correct type. Using Bool by default is not good. So the problem is how to write a compilable version of const I = UniformScaling(one()) with the type of Matrix{T}? I'm too newbie to have the answer. I had the idea to change I into a composition function returning one() and passing the template T to one (this is easy to do in C++).

Anyway the fix for my case is the follow code:

const mpI = UniformScaling(one(MP{Float64}).val)
Matrix{MP{Float64}}(mpI, 2,2)

[Lecrapouille]

Close since this code is fixing it:

MP(x::Bool) = x ? one(MP{Float64}) : zero(MP{Float64})