#90 - Add quadratic expansion with scalars

[mforets]

Closes #90.

[mforets]

This PR addresses #90 by generalizing quadratic_expansion given two scalars alpha and beta. Moreover,

• the new function square is a special case, and experiments show that we can make it faster using some of the ideas in this PR:

n = 2
M = rand(IntervalMatrix, n, n)

@btime square($M)
@btime quadratic_expansion($M, 0.0, 1.0)

  3.552 μs (82 allocations: 3.41 KiB)
  197.505 ns (2 allocations: 160 bytes)
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-0.610303, 0.937016]  [-0.507415, 0.486545]
 [-1.99161, 1.45359]     [0.264239, 2.02786] 

• the function square is more precise for using power ^ when possible (since x^2 cuts negative terms while x*x doesn't, if x is an interval)

• also the speed difference comes from using the default pow (ref. Speed up pow by 4x or more JuliaIntervals/IntervalArithmetic.jl#210 )

• i wonder if B[j, j] = (α + β*A[j, j])*A[j, j] can be handled in a better by inspecting the expression as in the reference..

[mforets]

• i tried Wstar from SUE for discretization calculation #90 on some random examples and the experiments shows that it is more accurate than the naive computation, so that is already an improvement that we can use in the discretization algorithm (ref. this notebook):

Wstar(A, δ) = 1/2 * A * δ^2 + 1/2 * A^2 * δ^3

Wstar (generic function with 1 method)

A = rand(IntervalMatrix) * 100
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-143.198, 51.5174]   [-52.0361, -47.4452]
 [-121.98, -21.5745]  [-105.322, 47.0649]  

δ = 0.01

Ws = Wstar(A, δ)
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-0.0103367, 0.0160024]  [-0.00516672, 0.00409371]
 [-0.0121116, 0.0140785]  [-0.00723271, 0.0110732] 

Wqe = quadratic_expansion(A, 1/2*δ^2, 1/2*δ^3)
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-0.0103367, 0.00707657]  [-0.00516672, 0.00386417]
 [-0.0121116, 0.00905817]  [-0.00723271, 0.00663447]

diam(Ws) - diam(Wqe)
2×2 Array{Float64,2}:
 0.00892577  0.000229537
 0.00502027  0.00443872 

[schillic]

i wonder if B[j, j] = (α + β*A[j, j])*A[j, j] can be handled in a better by inspecting the expression as in the reference..

True, (α + βx) x = αx + βx², so you can replace a multiplication by an addition (which sounds more precise).

[mforets]

Doing (α + βx) x looses the dependency between the arguments, so it may be less precise using interval arithmetic. If there is a general rule, one has to find it by taking the derivative of the expression as in the paper, which is a particular case with alpha, beta being something like t and 1/2t^2.

julia> fprod(α, β, x) = (α + β*x)*x
fprod (generic function with 1 method)

julia> fsum(α, β, x) = α*x + β*x^2
fsum (generic function with 1 method)

julia> y = Interval(-1, 1)
[-1, 1]

julia> fprod(0.0, 1.0, y)
[-1, 1]

julia> fsum(0.0, 1.0, y) # fsum is tighter
[0, 1]

julia> y = Interval(-3, -2)
[-3, -2]

julia> fprod(1.0, 2.0, y) # fprod is tighter
[6, 15]

julia> fsum(1.0, 2.0, y)
[5, 16]