overestimation free quadratic expansion

[lucaferranti]

fixes #109

With a quick look at the exponential.jl file, that function isn't used anywhere else anyway?

[mforets]

With a quick look at the exponential.jl file, that function isn't used anywhere else anyway?

indeed, i didn't spot other places where it is used.

otoh, see for example here, method that wants to compute the term 1/2 * δ^2 * A + 1/6 * δ^3 * A². so indeed instead of

IδW = Iδ + 1/2 * δ^2 * A + 1/6 * δ^3 * A²

that method could be improved by using

IδW = Iδ + quadratic_expansion(A, 1/2 * δ^2, 1/6 * δ^3)

[mforets]

@lucaferranti how about that we add the tests from #109 (comment) ?

[lucaferranti]

what about the

quadratic_expansion(A::IntervalMatrix{T}, t) where {T}

from what I can tell, it's the same of the more general one with $\alpha=t$ $\beta=1/2t^2$, is the current implementation of the 2-input methods overestimate free? maybe that could be replaced with

quadratic_expansion(A::IntervalMatrix{T}, t) where {T} = quadratic_expansion(A, t, 1/2*t^2)

[lucaferranti]

julia> t = 1.0
1.0

julia> α = 1.0
1.0

julia> β = 0.5
0.5

julia> A = IntervalMatrix(randn(100, 100) .± abs.(randn(100, 100)));

julia> @btime quadratic_expansion($A, $t);
  135.591 ms (7544198 allocations: 115.27 MiB)

julia> @btime quadratic_expansion($A, $α, $β);
  71.674 ms (3606 allocations: 296.91 KiB)

so maybe we could define

quadratic_expansion(A, t) = quadratic_expansion(A, t, t^2/2)

[mforets]

yes, the code refactoring that you propose sounds good to me

[lucaferranti]

quick question about the docstring, it currently says

### Input
- `A` -- interval matrix
- `t` -- non-negative time value

would it be better to keep it "application agnostic" and say that t is just a real number? Mathematically, I think the algorithm works also if t is negative (at least after the proposed refactoring)

[schillic]

True. Maybe it is better to just write quadratic_expansion(A, t, t^2/2) in the two places where it is used and remove that method altogether.

[lucaferranti]

the function is exported though, is it used in some other packages? According to JuliaHub, it has zero dependents

Just trying to avoid breaking things in other packages, although semantic versioning would take care things don't break by accident

[mforets]

According to JuliaHub, it has zero dependents

hmm i saw that too (here). strange, since IM is a dependency of ReachabilityAnalysis (see here)

[mforets]

about removing quadratic_expansion(A, t), sure, why not. i think this method is only used internally.

the reason to have a special method is that it's the canonical case (for the matrix exp expansion).

[lucaferranti]

just double checking one last time before pulling the trigger, remove or replace quadratic_expansion(A, t)?

[schillic]

As you prefer. I think removing is good because it simplifies the code base.

[lucaferranti]

added the test for quadratic expansion from the issue comment, now this should be ready.

I noticed the tests for exponentials aren't really testing the result is correct, just that the output is an interval matrix, but I guess that can be addressed in a different PR

[schillic]

Yes, we were eager to see results back then 😅

[lucaferranti]

for the release, since this is breaking and there are at least a couple of other breaking changes coming in the upcoming days (moving exponential to ILA, bringing rump multiplication here, etc.), maybe we can wait till Tuesday meeting for the next release?

[schillic]

Just for the record, I went back to the method that was removed here and compared it to the new implementation. The old implementation gives slightly more precise bounds on some random example. Not sure if this is general or a coincidence, and why (could be just accumulated rounding errors due to different arithmetic operations).

julia> A = rand(IntervalMatrix)
2×2 IntervalMatrix{Float64, IntervalArithmetic.Interval{Float64}, Matrix{IntervalArithmetic.Interval{Float64}}}:
 [-1.37152, 0.237901]   [-0.502148, 0.325257]
 [-0.353947, 0.231257]   [1.13864, 1.22323]

julia> B = quadratic_expansion(A, 0.1)  # old implementation
2×2 IntervalMatrix{Float64, IntervalArithmetic.Interval{Float64}, Matrix{IntervalArithmetic.Interval{Float64}}}:
 [-0.128327, 0.0249617]   [-0.0538833, 0.0349019]
 [-0.0379805, 0.0248152]   [0.119766, 0.130693]

julia> C = quadratic_expansion(A, 0.1, 0.5*0.1^2)  # new implementation
2×2 IntervalMatrix{Float64, IntervalArithmetic.Interval{Float64}, Matrix{IntervalArithmetic.Interval{Float64}}}:
 [-0.128327, 0.0249617]   [-0.0538833, 0.0349019]
 [-0.0379805, 0.0248152]   [0.119766, 0.130693]

# they are very similar, which needs to be checked in an awkward way because ≈ does not work for Interval

julia> all(B[i, j].lo ≈ C[i, j].lo for (i, j) in [(1,1), (1,2), (2,1), (2,2)])
true

julia> all(B[i, j].hi ≈ C[i, j].hi for (i, j) in [(1,1), (1,2), (2,1), (2,2)])
true

# how are they formally related?

julia> B ⊆ C
true

julia> C ⊆ B
false

[lucaferranti]

what is the order of magnitude of the overestimation? I'd guess it's due to directed rounding. Are you sure the previous implementation was rigorous?

[lucaferranti]

for example, taking a patch from the previous implementation

            a = inf(aii) * t + inf(aii)^2 * t2d2
            b = sup(aii) * t + sup(aii)^2 * t2d2
            return min(a, b)

those are floating point operations subject to rounding error, the returned minimum could be 1e-15 bigger that the true minimum. If you compute maximum(diam(C)./diam(B)) what do you get. If it is in the order of magnitude of 1e-16 -- 1e14, then it's just due to direct rounding (rounding outwards at each operation). In that case, I wouln't say that the previous one was "more precise", it was simply not rigorous.

[schillic]

Yes, that could very well be the explanation. The difference is in that order. You are right, I should not have said "more precise."