Matrix Multiplication of Intervals Raises Error

[chelseas]

julia> bnds = [Interval(-1., 1.); Interval(-1., 1.)]
2-element Vector{Interval{Float64}}:
 Interval{Float64}(Interval{Float64}(-1.0, 1.0))
 Interval{Float64}(Interval{Float64}(-1.0, 1.0))
julia> [0. 1.]
1×2 Matrix{Float64}:
 0.0  1.0
julia> [0. 1.] * bnds
ERROR: MethodError: no method matching zero(::MinkowskiSum{Float64, LinearMap{…}, LinearMap{…}})
The function `zero` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  zero(::Type{Union{}}, Any...)
   @ Base number.jl:310
  zero(::Type{CDDLib.GMPRationalMut})
   @ CDDLib ~/.julia/packages/CDDLib/OpQiu/src/mytype.jl:50
  zero(::Type{Dates.DateTime})
   @ Dates ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Dates/src/types.jl:458
  ...

Stacktrace:
 [1] _generic_matvecmul!(C::Vector{…}, tA::Char, A::Matrix{…}, B::Vector{…}, _add::LinearAlgebra.MulAddMul{…})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:845
 [2] generic_matvecmul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:787 [inlined]
 [3] _mul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:73 [inlined]
 [4] mul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:70 [inlined]
 [5] mul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:253 [inlined]
 [6] *(A::Matrix{Float64}, x::Vector{Interval{Float64}})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:60
 [7] top-level scope
   @ REPL[97]:1
Some type information was truncated. Use `show(err)` to see complete types.

And if I try to assert the shapes (should be correct before anyway, 1x2 * 2x1) I still get an error:

[0. 1.] * reshape(bnds, 2, 1)
ERROR: MethodError: no method matching zero(::MinkowskiSum{Float64, LinearMap{…}, LinearMap{…}})
The function `zero` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  zero(::Type{Union{}}, Any...)
   @ Base number.jl:310
  zero(::Type{CDDLib.GMPRationalMut})
   @ CDDLib ~/.julia/packages/CDDLib/OpQiu/src/mytype.jl:50
  zero(::Type{Dates.DateTime})
   @ Dates ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/Dates/src/types.jl:458
  ...

Stacktrace:
 [1] _generic_matmatmul!(C::Matrix{…}, A::Matrix{…}, B::Matrix{…}, _add::LinearAlgebra.MulAddMul{…})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:918
 [2] generic_matmatmul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:868 [inlined]
 [3] _mul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:287 [inlined]
 [4] mul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:285 [inlined]
 [5] mul!
   @ ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:253 [inlined]
 [6] *(A::Matrix{Float64}, B::Matrix{Interval{Float64}})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:114
 [7] top-level scope
   @ REPL[99]:1
Some type information was truncated. Use `show(err)` to see complete types.

I feel like this is basic functionality that should work, no?

[OlivierHnt]

Thank you for reporting this issue :). I cannot reproduce it on the latest version of IntervalArithmetic ; are you sure you are using v0.22.21?

[chelseas]

I tried this and it didn't work either:

import Base.zero
zero(::Interval{Float64}) = Interval(0.0, 0.0)

No I am using "0.21.2" through LazySets and they have a restriction on the version :/// :

⌅ [d1acc4aa] IntervalArithmetic v0.21.2 (<v0.22.21): LazySets

[chelseas]

For now I will commit a little type piracy

import Base.*
function *(A::Matrix{Float64}, B::Matrix{Interval{Float64}})
    C = Matrix{Interval}(undef, size(A, 1), size(B, 2))
    for i in 1:size(A, 1)
        for j in 1:size(B, 2)
            s = Interval(0.0, 0.0)
            for k in 1:size(A, 2)
                s += A[i, k] * B[k, j]
            end
            C[i, j] = s
        end
    end
    return C
end

but this should be resolved in a better way in the future.

[schillic]

Your Interval is a LazySets.Interval, not an IntervalArithmetic.Interval. Hence I think this is not an issue of this library.

When using IntervalArithmetic.Interval, the code even works with old versions of IntervalArithmetic (e.g., v0.21.2).

julia> using IntervalArithmetic

julia> bnds = [interval(-1., 1.); interval(-1., 1.)]
2-element Vector{Interval{Float64}}:
 [-1.0, 1.0]
 [-1.0, 1.0]

julia> [0. 1.] * bnds
1-element Vector{Interval{Float64}}:
 [-1.0, 1.0]

[schillic]

Also note that bnds is a Vector, not a Matrix, so I am not sure whether your * method would resolve the issue.

Generally, consider using https://github.com/JuliaReach/IntervalMatrices.jl for calculations with interval matrices.

[OlivierHnt]

Thanks @schillic for clarifying this. I close this issue then; if anything comes up, feel free to re-open it or open a new one.

[chelseas]

@schillic Is the LazySets.Interval not a wrapper of IntervalArithmetic.Interval?

[schillic]

Yes it is, but it does not act like a Real but like a set. Hence, for instance, a different * method is used:

julia> @which [0. 1.] * bnds  # IntervalArithmetic.Interval; note the `S<:Real` for the eltype of `bnds`
*(A::Union{LinearAlgebra.Adjoint{<:Any, <:StridedMatrix{T}}, LinearAlgebra.Transpose{<:Any, <:StridedMatrix{T}}, StridedMatrix{T}}, x::StridedVector{S}) where {T<:Union{Float32, Float64, ComplexF64, ComplexF32}, S<:Real}
     @ LinearAlgebra .../share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:53

julia> @which [0. 1.] * bnds  # LazySets.Interval
*(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T, S}
     @ LinearAlgebra .../share/julia/stdlib/v1.11/LinearAlgebra/src/matmul.jl:58

But that is actually not the issue here. Both methods would yield the same error.

The issue is that the generic matrix-vector multiplication tries to initialize the solution vector, which requires the zero element of addition. Hence the call to zero in the OP, which results in an error because it is not defined for sets.

More generally, LazySets does not define what a matrix multiplied with a vector of sets should be. And I am not convinced it would be a good idea to define it (but feel free to open an issue in LazySets to discuss it).