Add centered forms for scalar functions

Author: koehlerson

Project.toml

@@ -3,7 +3,10 @@ uuid = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2"
 version = "0.4.2"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
 IntervalArithmetic = "0.15, 0.16, 0.17"

src/IntervalOptimisation.jl

@@ -3,6 +3,7 @@ module IntervalOptimisation
 export minimise, maximise,
        minimize, maximize
 export HeapedVector, SortedVector
+export mean_value_form_scalar, third_order_taylor_form_scalar
 
 include("StrategyBase.jl")
 using .StrategyBase
@@ -13,10 +14,11 @@ using .SortedVectors
 include("HeapedVectors.jl")
 using .HeapedVectors
 
-using IntervalArithmetic
+using IntervalArithmetic, IntervalRootFinding
+using LinearAlgebra
 
 include("optimise.jl")
-
+include("centered_forms.jl") 
 
 const minimize = minimise
 const maximize = maximise

src/centered_forms.jl

@@ -0,0 +1,33 @@
+import ForwardDiff: gradient, jacobian, hessian
+
+gradient(f, X::IntervalBox) = gradient(f, X.v)
+# jacobian(f, X::IntervalBox) = jacobian(f, X.v)
+hessian(f, X::IntervalBox) = hessian(f, X.v)
+
+"""
+Calculate the mean-value form of a function ``f:\\mathbb{R}^n \\to \\mathbb{R}``
+using the gradient ``\nabla f``;
+this gives a second-order approximation.
+"""
+function mean_value_form_scalar(f, X)
+    m = IntervalBox(mid(X))
+
+    return f(m) + gradient(f, X.v) ⋅ (X - m)
+end
+
+mean_value_form_scalar(f) = X -> mean_value_form_scalar(f, X)
+
+
+"""
+Calculate a third-order Taylor form of ``f:\\mathbb{R}^n \\to \\mathbb{R}`` using the Hessian.
+"""
+function third_order_taylor_form_scalar(f, X)
+    m = IntervalBox(mid(X))
+
+    H = hessian(f, X)
+    δ = X - m
+
+    return f(m) + gradient(f, m) ⋅ δ + 0.5 * sum(δ[i]*H[i,j]*δ[j] for i in 1:length(X) for j in 1:length(X))
+end
+
+third_order_taylor_form_scalar(f) = X -> third_order_taylor_form_scalar(f, X)

src/optimise.jl

@@ -1,6 +1,9 @@
 numeric_type(::Interval{T}) where {T} = T
 numeric_type(::IntervalBox{N, T}) where {N, T} = T
 
+interval_mid(X::Interval) = Interval(mid(X))
+interval_mid(X::IntervalBox) = IntervalBox(mid(X))
+
 """
     minimise(f, X, structure = SortedVector, tol=1e-3)
     or
@@ -37,7 +40,7 @@ function minimise(f, X::T; structure = HeapedVector, tol=1e-3) where {T}
         end
 
         # find candidate for upper bound of global minimum by just evaluating a point in the interval:
-        m = sup(f(Interval.(mid.(X))))   # evaluate at midpoint of current interval
+        m = sup(f(interval_mid(X)))   # evaluate at midpoint of current interval
 
         if m < global_min
             global_min = m