Bounders for TaylorModels

The following is a blog post about my work at Google Summer of Code 2020 on the project Taylor Models and a guaranteed ODE solver. This project aims to improve some of the methods in the existing package TaylorModels.jl for the Julia Language.

What is a Taylor Model?

]add TaylorModels IntervalArithmetic

# Load the required packages for this notebook

using TaylorModels

using IntervalArithmetic

using Plots

using WebSockets

using WebIO, Interact

A Taylor Model is a tool for obtaining a validated polynomial approximation to a function inline_formula not implemented, where validated refers to the fact that the true function is guaranteed to be represented (or enclosed) by this approximation.

A Taylor Model inline_formula not implemented over the domain inline_formula not implemented for the function inline_formula not implemented is the tuple inline_formula not implemented specified by a inline_formula not implemented-degree polynomial inline_formula not implemented and remainder inline_formula not implemented, such that inline_formula not implemented. In other words, inline_formula not implemented is a rigorous approximation to the function inline_formula not implemented over the domain inline_formula not implemented and inline_formula not implemented provides an enclosure that guarantees that the true function inline_formula not implemented is contained in the set inline_formula not implemented for every inline_formula not implemented in the domain. An important fact about the remainder inline_formula not implemented is that it must satisfy inline_formula not implemented, this means that the exact value of the function inline_formula not implemented is enclosed by inline_formula not implemented.

As the name suggests the polynomial inline_formula not implemented is just the Taylor Series for the function inline_formula not implemented, typically the expansion point of the series is taken to be the midpoint of inline_formula not implemented, this can be generalized to consider other types of series expansions. Notice that the remainder inline_formula not implemented acts as a "tube" around the expansion inline_formula not implemented that "absorbs" all kind of computation errors (truncation, roundoff) to ensure that inline_formula not implemented is contained in inline_formula not implemented.

Remark: It is advisable to run this blog post interactively, in order to be able to exploit the sliders that are included in some cases.

Let's see this in action, as a matter of example we're going to compute the Taylor Model for the function inline_formula not implemented using the package TaylorModels.jl.

D = 0 .. 5 # Taylor Model domain

x0 = mid(D) # expansion point for the Taylor Series

f(x) = exp(x)

xrange = range(D.lo, D.hi, length=100)

@manipulate for order in 2:5

# Let's define the Taylor Model for the univariate function f(x)

  tm = TaylorModel1(order, x0, D) # create a Taylor Model "variable"

  fT = f(tm) # from evaluating this variable on f we get its Taylor Model

  plot(fT, label="Taylor Model, order=$order")

  plot!(xrange, x -> f(x), label="f(x)", lw=3)

  ylims!(0, 150)

end

From the previous animation, we see that as the order of the expansion is increased its remainder inline_formula not implemented gets tighter. This is the expected behavior since there is less error due to the expansion to be included in inline_formula not implemented.

The Taylor Models defines its own arithmetic, that means, every arithmetical operation has a corresponding operation in terms of Taylor Models. So, for example, we can define the product inline_formula not implemented of two Taylor Models for the functions inline_formula not implemented, such that

where inline_formula not implemented corresponds to the resulting Taylor Model of the product inline_formula not implemented. In other words, if we define a inline_formula not implemented Taylor Models for the functions inline_formula not implemented and then multiply them, the resulting Taylor Model is guaranteed to enclose the function inline_formula not implemented, and this holds for all arithmetic operations. In this same way is possible to compute composition of Taylor Models. The details of this Taylor Model arithmetic and composition can be found in [1].

To illustrate this Taylor Model arithmetic, let's compute the following

f(x) = sin(x)

g(x) = cos(x)

h(x) = x^3

s(x) = h(x) * (f(x) + g(x))

D = -1 .. 1

x0 = mid(D)

xrange = range(D.lo, D.hi, length=100)

@manipulate for order in 2:5

  tm = TaylorModel1(order, x0, D)

  fT = f(tm)

  gT = g(tm)

  hT = h(tm)

  sT = hT * (fT + gT)

  plot(sT, label="Taylor Model, order=$order")

  plot!(xrange, x -> s(x), label="s(x)", lw=3)

  ylims!(-3, 3)

end

From this comes the true power of Taylor Models as this enables us to get validated results for a wide variety of calculations by just doing it over Taylor Models.

Bounds

The bound for a model (not just for a Taylor Model) is the enclosure that contains the range of the function inline_formula not implemented over some domain D. So, for example, for an interval-based model, this bound is just the function evaluated in inline_formula not implemented by interval arithmetic. For a Taylor Model, this corresponds to evaluate its polynomial part on interval arithmetic over inline_formula not implemented and then sum its remainder inline_formula not implemented.

# compute the bound for a Taylor Model

f(x) = 1 + x^5 - x^4

D = 0 .. 1

x0 = mid(D)

order = 3

bound_interval = f(D)

tm = TaylorModel1(order, x0, D)

fT = f(tm)

bound_tm = fT(fT.dom - fT.x0)

box_interval = IntervalBox(D, bound_interval)

box_tm = IntervalBox(D, bound_tm)

xrange = range(0, 1, length=100)

plot(xrange, x -> f(x), color=:red, label="f(x)", lw=3)

plot!(box_interval, label="Bound Interval Arithmetic")

plot!(box_tm, label="Bound Taylor Model")

Both types of enclosures contain the range of the function inline_formula not implemented over inline_formula not implemented, but with a great over-estimation of the function range, a way to improve this is to split or mince the original domain inline_formula not implemented in smaller equally large subdomains.

# compute the bound for a Taylor Model

f(x) = 1 + x^5 - x^4

D = mince(0 .. 1, 16) # 16 equally large subdomains

order = 3

p = plot()

interval_boxes = []

tm_boxes = []

for (idx, Di) in enumerate(D)

  x0 = mid(Di)

  bound_interval = f(Di)

  tm = TaylorModel1(order, x0, Di)

  fT = f(tm)

  bound_tm = fT(fT.dom - fT.x0)

  box_interval = IntervalBox(Di, bound_interval)

  box_tm = IntervalBox(Di, bound_tm)

  if idx == 16

    label_interval = "Interval Arithmetic"

    label_tm = "Taylor Model order=$order"

  else

    label_interval = ""

    label_tm = ""

  end

   plot!(box_interval, c=:green, label=label_interval)

   plot!(box_tm, c=:red, label=label_tm)

end

xrange = range(0, 1, length=100)

plot!(xrange, x -> f(x), color=:red, label="f(x)", lw=3)

The last plot shows an important fact: Taylor Models allows for better enclosures of a function.

To improve these bounds we can go further and split the domain even further, but a better (in terms of computational efficiency) choice is to use the so-called range bounders.

Range Bounders for Taylor Models

The main idea of range bounders is to bound the minimum (or maximum) of a Taylor Model over some region in order to tighten its bound. By bounding the minimum (maximum) of the Taylor Model in an interval inline_formula not implemented, then we can just take the lower (upper) bound inline_formula not implemented (inline_formula not implemented) as one of the bound for the Taylor Model as by definition nothing can be below (above) of the minimum (maximum). The kind of range bounders that follows, use the information from the linear and quadratic terms of the Taylor Model in order to estimate this bound inline_formula not implemented and then use it to improve its bound.

Linear Dominated Bounder

The linear dominated bounder (ldb) is an algorithm capable of tightening the bounds of a Taylor Model. To accomplish this, the ldb algorithm uses the linear part of a Taylor Model as a guideline to iteratively reduce the domain. This algorithm is motivated by the fact that far from a local optimum, the linear part of the function typically dominates its behavior, so, we can use exploit it to tightening the bound of the Taylor Model. The algorithm proceeds as follow

1. Re-expand the polynomial inline_formula not implemented at the midpoint inline_formula not implemented of the interval domain inline_formula not implemented.

2. Check the sign of the linear coefficients inline_formula not implemented's and save them for later.

3. Compute the bound of the linear part inline_formula not implemented and of the non-linear part inline_formula not implemented.

4. The minimum of the function is bounded by inline_formula not implemented, define inline_formula not implemented, if inline_formula not implemented is less than some pre-specified threshold, then finish, else

   1. Reduce the i-th domain coordinate inline_formula not implemented by the following criteria

      • If inline_formula not implemented; then

        inline_formula not implemented

      • If inline_formula not implemented; then

        inline_formula not implemented

      • Else; nothing to do, as no linear part.

   2. Go to 1, replacing the original domain inline_formula not implemented with the reduced domain of the previous step.

As an example of this algorithm let's improve the bound of the Taylor Model for the function inline_formula not implemented over the domain inline_formula not implemented. To this end, we're going to mince inline_formula not implemented in 16 equally large subdomains and compute a Taylor Model on every one of them,

f(x) = 1 + x^5 - x^4

order = 5

subdomains = mince(0 .. 1, 16)

xrange = range(0, 1, length=100)

boxes_tm = []

boxes_ldb = []

for i in 1:length(subdomains)

  p = plot(legend=false)

  plot!(xrange, x -> f(x), c=:red,)

  D = subdomains[i]

  x0 = mid(D)

  tm = TaylorModel1(order, x0, D)

  fT = f(tm)

  bound_tm = fT(fT.dom - fT.x0)

  bound_ldb = linear_dominated_bounder(fT)

  box_tm = IntervalBox(D, bound_tm)

  box_ldb = IntervalBox(D, bound_ldb)

  push!(boxes_tm, box_tm)

  push!(boxes_ldb, box_ldb)

end

@manipulate for i in 1:length(subdomains)

  box_ldb = boxes_ldb[i]

  p = plot()

  for box in boxes_tm

    plot!(box, label="", fill=0.)

  end

  plot!(box_ldb, label="TM5 + Linear dominated bounder")

  plot!(xrange, x -> f(x), label="f(x)", lw=3)

end

As we see, there are some subdomains with a great over-estimation of the bound, the last one, and the one that has the optimum.

Quadratic Fast Bounder

The idea behind the quadratic fast bounder (qfb) is to bound the second-order terms of the Taylor Series expansion, this becomes noticeably important when there is a local optimum, as the linear part vanishes and therefore the ldb algorithm can't properly work.

The key idea of this algorithm boils down to the following decomposition of a Taylor Model

Then, the lower bound for this model is just

where inline_formula not implemented denotes the lower bound and inline_formula not implemented is a factor introduced for convenience. For this algorithm, inline_formula not implemented is chosen to be

with inline_formula not implemented being the hessian of inline_formula not implemented. Note that, near a local minimizer this hessian becomes positive definite, that means inline_formula not implemented for every vector inline_formula not implemented, therefore inline_formula not implemented trivially. Better yet, if we choose inline_formula not implemented as

for some inline_formula not implemented in the domain, then inline_formula not implemented has linear terms and for a careful selection of inline_formula not implemented, inline_formula not implemented can contain no linear terms and still have inline_formula not implemented. By this selection, we're bounding not just the second-order terms but the first ones too. It turns out that this inline_formula not implemented must be chosen so that is a solution to the following system of linear equations

where inline_formula not implemented is the vector of linear terms of inline_formula not implemented.

Let's use this range bounder for the same function as in the ldb example

f(x) = 1 + x^5 - x^4

order = 5

subdomains = mince(0 .. 1, 16)

xrange = range(0, 1, length=100)

boxes_tm = []

boxes_qfb = []

for i in 1:length(subdomains)

  p = plot(legend=false)

  plot!(xrange, x -> f(x), c=:red,)

  D = subdomains[i]

  x0 = mid(D)

  tm = TaylorModel1(order, x0, D)

  fT = f(tm)

  bound_tm = fT(fT.dom - fT.x0)

  bound_qfb = quadratic_fast_bounder(fT)

  box_tm = IntervalBox(D, bound_tm)

  box_qfb = IntervalBox(D, bound_qfb)

  push!(boxes_tm, box_tm)

  push!(boxes_qfb, box_qfb)

end

@manipulate for i in 1:length(subdomains)

  box_ldb = boxes_qfb[i]

  p = plot()

  for box in boxes_tm

    plot!(box, label="", fill=0.)

  end

  plot!(box_ldb, label="TM5 + Quadratic fast bounder")

  plot!(xrange, x -> f(x), label="f(x)", lw=3)

end

Bounds comparison for distinct functions

f(x) = x^3 * sin(x)

D = -3 .. 3

n = 16

subdomains = mince(D, n)

xrange = range(D.lo, D.hi, length=100)

boxes_interval = [IntervalBox(d, f(d)) for d in subdomains]

@manipulate for i in eachindex(subdomains), order in 3:5, bounder in ("nobounder", "ldb", "qfb")

  d = subdomains[i]

  x0 = mid(d)

  plot(boxes_interval, label="Interval Arithmetic", fill=0., legend=:top)

  tm = TaylorModel1(order, x0, d)

  fT = f(tm)

  if bounder == "nobounder"

    bound = fT(fT.dom - fT.x0)

    box = IntervalBox(d, bound)

    plot!(box, label="Naive Taylor Model")

  elseif bounder == "ldb"

    bound = linear_dominated_bounder(fT)

    box = IntervalBox(d, bound)

    plot!(box, label="Linear dominated bounder")

  elseif bounder == "qfb"

    bound = quadratic_fast_bounder(fT)

    box = IntervalBox(d, bound)

    plot!(box, label="Quadratic fast bounder")

  end

  plot!(xrange, x -> f(x), label="f(x)", c=:red, lw=2)

end

f(x) = x^2 * sin(x^2 - 3) + x

D = -3 .. 3

n = 16

subdomains = mince(D, n)

xrange = range(D.lo, D.hi, length=100)

boxes_interval = [IntervalBox(d, f(d)) for d in subdomains]

@manipulate for i in eachindex(subdomains), order in 3:5, bounder in ("nobounder", "ldb", "qfb")

  d = subdomains[i]

  x0 = mid(d)

  plot(boxes_interval, label="Interval Arithmetic", fill=0., legend=:top)

  tm = TaylorModel1(order, x0, d)

  fT = f(tm)

  if bounder == "nobounder"

    bound = fT(fT.dom - fT.x0)

    box = IntervalBox(d, bound)

    plot!(box, label="Naive Taylor Model")

  elseif bounder == "ldb"

    bound = linear_dominated_bounder(fT)

    box = IntervalBox(d, bound)

    plot!(box, label="Linear dominated bounder")

  elseif bounder == "qfb"

    bound = quadratic_fast_bounder(fT)

    box = IntervalBox(d, bound)

    plot!(box, label="Quadratic fast bounder")

  end

  plot!(xrange, x -> f(x), label="f(x)", c=:red, lw=2)

end

f(x) = x^2 * sin(x^2 - 3) + cos(3 - x^3) * x

D = -3 .. 3

n = 16

subdomains = mince(D, n)

xrange = range(D.lo, D.hi, length=200)

boxes_interval = [IntervalBox(d, f(d)) for d in subdomains]

@manipulate for i in eachindex(subdomains), order in 3:5, bounder in ("nobounder", "ldb", "qfb")

  d = subdomains[i]

  x0 = mid(d)

  plot(boxes_interval, label="Interval Arithmetic", fill=0., legend=:top)

  tm = TaylorModel1(order, x0, d)

  fT = f(tm)

  if bounder == "nobounder"

    bound = fT(fT.dom - fT.x0)

    box = IntervalBox(d, bound)

    plot!(box, label="Naive Taylor Model")

  elseif bounder == "ldb"

    bound = linear_dominated_bounder(fT)

    box = IntervalBox(d, bound)

    plot!(box, label="Linear dominated bounder")

  elseif bounder == "qfb"

    bound = quadratic_fast_bounder(fT)

    box = IntervalBox(d, bound)

    plot!(box, label="Quadratic fast bounder")

  end

  plot!(xrange, x -> f(x), label="f(x)", c=:red, lw=2)

end

References

[1] Mioara Joldes, Rigorous Polynomial Approximations and Applications, PhD Thesis, École Normale Supérieure de Lyon, ENS-Lyon (2011).

[2] Kyoko Makino, Rigorous Analysis of Nonlinear Motion in Particle Accelerators, PhD Thesis, Michigan State University (1991).

[3] Kyoko Makino & Martin Berz, Verified Global Optimization with Taylor Model based Range Bounders, Transactions on Computers (2005).

[4] Martin Berz, Kyoko Makino, Youn-Kyung Kim, Long-term stability of the Tevatron by verified global optimization, Nuclear Instruments and Methods in Physics Research A (2006).