Unexpected negative lower bound for product of positive interval row vector and positive interval vector

[ChangheYang1999]

Hi, I found a surprising result when multiplying a positive interval row vector by a positive interval vector.

All entries of both vectors are strictly positive intervals, so the resulting dot product should also have a positive lower bound. However, b * a returns an interval with a negative lower bound.

Minimal working example:

using Pkg
using IntervalArithmetic

Pkg.status("IntervalArithmetic")

a = Interval{Float64}[
    interval(77.0713, 645.097),
    interval(59.9388, 282.904),
    interval(3.8025, 6.30973)
]

b = Interval{Float64}[
    interval(2.35903e14, 3.16785e15)  interval(2.36564e14, 1.38742e15)  interval(3.89381e12, 1.68618e13)
]

println(a)
println(b)
println(size(a), " ", size(b))
println(b * a)

bb = vec(b)
s = zero(eltype(a))
for i in eachindex(a)
    s += bb[i] * a[i]
end
println("manual dot = ", s)

On my side, the package version is:

IntervalArithmetic v1.0.9

The output of b * a is:

Interval{Float64}[[-9.28652e+17, 2.43619e+18]_com]

This is surprising because all entries are positive. The expected lower bound should be positive. A direct endpoint estimate gives approximately:

[3.24e16, 2.436e18]

The upper bound agrees with the matrix multiplication result, but the lower bound is much too wide and becomes negative.

I understand that interval matrix multiplication may return an enclosure rather than the tightest possible interval. However, since this is a dot product of positive intervals, I wonder whether this behavior is expected, or whether the matrix multiplication path is using a method that loses sign information too aggressively.

Is this intended behavior in v1.0.9, or could it be a bug / missed optimization in interval matrix-vector multiplication?

Thanks!

[schillic]

Smaller example:

julia> a = Interval{Float64}[interval(2, 30)];

julia> b = Interval{Float64}[interval(2, 3);;];

julia> b * a
1-element Vector{Interval{Float64}}:
 [-10.0, 90.0]_com

[lbenet]

Following on @schillic minimal example:

julia> b[1,1] * a[1]
[4.0, 90.0]_com

That is, the product is ok, which I interpret as a problem related to b being a matrix and a a vector...

[OlivierHnt]

I don't recall the exact statement of Rump's theorem for his algorithm, but that should be (I hope 😅) within the expected error bounds. Since v1.0, this algo is the default algorithm since the naive approach does not scale well and becomes impractical for large matrix products.

I believe it is mentioned in the docs somewhere that you can use configure to use different algorithms:

julia> IntervalArithmetic.configure(matmul=:slow)
Configuration options:
  - bound type: Float64
  - flavor: set_based
  - interval rounding: correct
  - power mode: fast
  - matrix multiplication mode: slow

julia> a = Interval{Float64}[interval(2, 30)];

julia> b = Interval{Float64}[interval(2, 3);;];

julia> b * a
1-element Vector{Interval{Float64}}:
 [4.0, 90.0]_com

EDIT: I closed the issue, but feel free to re-open it if I overlooked something of course 🙂 !

EDIT 2: I suppose we could

• check if the operands have all positive entries and cut-off negative values. So that would return [0, 90] instead of [-10, 90]...
• and / or we could change the matrix product behaviour when the size of the matrices are small...

[dpsanders]

If you are literally just needing a dot product, you can do

using LinearAlgebra: dot

dot(vec(b), a)

[schillic]

I still think this is wrong. Even with intervals, a negative number should never pop up when you multiply and add.

[OlivierHnt]

Sure the algorithm could be refined, and an easy way is to get rid of the negative parts (as I suggested above).

[schillic]

I meant more that this looks like a bug.
To be fair, I do not know how the algorithm works. I would just expect that it does not use subtractions. But maybe it does.

[OlivierHnt]

The algorithm involves taking the midpoint-radius form of the matrices and performing several matrix multiplications in floating-point arithmetic with directed rounding, before (in a nutshell) returning the result C = [mC - rC, mC + rC].
So there is an inherent overestimation, which I believe is what we are looking at but it would be nice to check this (here is a reference, https://www.tuhh.de/ti3/paper/rump/Ru11a.pdf, IIRC we use Algorithm 4.9).

[ChangheYang1999]

The algorithm involves taking the midpoint-radius form of the matrices and performing several matrix multiplications in floating-point arithmetic with directed rounding, before (in a nutshell) returning the result C = [mC - rC, mC + rC]. So there is an inherent overestimation, which I believe is what we are looking at but it would be nice to check this (here is a reference, https://www.tuhh.de/ti3/paper/rump/Ru11a.pdf, IIRC we use Algorithm 4.9).

Thanks. That makes sense to me.