Theory Taylor_Models
theory Taylor_Models
imports
"Horner_Eval"
"Polynomial_Expression_Additional"
"Taylor_Models_Misc"
"HOL-Decision_Procs.Approximation"
"HOL-Library.Function_Algebras"
"HOL-Library.Set_Algebras"
"Affine_Arithmetic.Straight_Line_Program"
"Affine_Arithmetic.Affine_Approximation"
begin
text ‹TODO: get rid of float poly/float inteval and use real poly/real interval
and data refinement?›
section ‹Multivariate Taylor Models›
subsection ‹Computing interval bounds on arithmetic expressions›
text ‹This is a wrapper around the "approx" function.
It computes range bounds on floatarith expressions.›
fun compute_bound_fa :: "nat ⇒ floatarith ⇒ float interval list ⇒ float interval option"
where "compute_bound_fa prec f I = approx prec f (map Some I)"
lemma compute_bound_fa_correct:
"interpret_floatarith f i ∈⇩r ivl"
if "compute_bound_fa prec f I = Some ivl"
"i all_in I"
for i::"real list"
proof-
have bounded: "bounded_by i (map Some I)"
using that(2)
unfolding bounded_by_def
by (auto simp: bounds_of_interval_eq_lower_upper set_of_eq)
from that have Some: "approx prec f (map Some I) = Some ivl"
by (auto simp: lower_Interval upper_Interval min_def split: option.splits if_splits)
from approx[OF bounded Some]
show ?thesis by (auto simp: set_of_eq)
qed
subsection ‹Definition of Taylor models and notion of rangeity›
text ‹Taylor models are a pair of a polynomial and an absolute error bound.›
datatype taylor_model = TaylorModel (tm_poly: "float poly") (tm_bound: "float interval")
text ‹Taylor model for a real valuation of variables›
primrec insertion :: "(nat ⇒ 'a) ⇒ 'a poly ⇒ 'a::{plus,zero,minus,uminus,times,one,power}"
where
"insertion bs (C c) = c"
| "insertion bs (poly.Bound n) = bs n"
| "insertion bs (Neg a) = - insertion bs a"
| "insertion bs (poly.Add a b) = insertion bs a + insertion bs b"
| "insertion bs (Sub a b) = insertion bs a - insertion bs b"
| "insertion bs (Mul a b) = insertion bs a * insertion bs b"
| "insertion bs (Pw t n) = insertion bs t ^ n"
| "insertion bs (CN c n p) = insertion bs c + (bs n) * insertion bs p"
definition range_tm :: "(nat ⇒ real) ⇒ taylor_model ⇒ real interval" where
"range_tm e tm = interval_of (insertion e (tm_poly tm)) + real_interval (tm_bound tm)"
lemma Ipoly_num_params_cong: "Ipoly xs p = Ipoly ys p"
if "⋀i. i < num_params p ⟹ xs ! i = ys ! i"
using that
by (induction p; auto)
lemma insertion_num_params_cong: "insertion e p = insertion f p"
if "⋀i. i < num_params p ⟹ e i = f i"
using that
by (induction p; auto)
lemma insertion_eq_IpolyI: "insertion xs p = Ipoly ys p"
if "⋀i. i < num_params p ⟹ xs i = ys ! i"
using that
by (induction p; auto)
lemma Ipoly_eq_insertionI: "Ipoly ys p = insertion xs p"
if "⋀i. i < num_params p ⟹ xs i = ys ! i"
using that
by (induction p; auto)
lemma range_tmI:
"x ∈⇩i range_tm e tm"
if x: "x ∈⇩i interval_of (insertion e ((tm_poly tm))) + real_interval (tm_bound tm)"
for e::"nat⇒real"
by (auto simp: range_tm_def x)
lemma range_tmD:
"x ∈⇩i interval_of (insertion e (tm_poly tm)) + real_interval (tm_bound tm)"
if "x ∈⇩i range_tm e tm"
for e::"nat⇒real"
using that
by (auto simp: range_tm_def)
subsection ‹Interval bounds for Taylor models›
text ‹Bound a polynomial by simply approximating it with interval arguments.›
fun compute_bound_poly :: "nat ⇒ float interval poly ⇒ (float interval list) ⇒ (float interval list) ⇒ float interval" where
"compute_bound_poly prec (poly.C f) I a = f"
| "compute_bound_poly prec (poly.Bound n) I a = round_interval prec (I ! n - (a ! n))"
| "compute_bound_poly prec (poly.Add p q) I a =
round_interval prec (compute_bound_poly prec p I a + compute_bound_poly prec q I a)"
| "compute_bound_poly prec (poly.Sub p q) I a =
round_interval prec (compute_bound_poly prec p I a - compute_bound_poly prec q I a)"
| "compute_bound_poly prec (poly.Mul p q) I a =
mult_float_interval prec (compute_bound_poly prec p I a) (compute_bound_poly prec q I a)"
| "compute_bound_poly prec (poly.Neg p) I a = -compute_bound_poly prec p I a"
| "compute_bound_poly prec (poly.Pw p n) I a = power_float_interval prec n (compute_bound_poly prec p I a)"
| "compute_bound_poly prec (poly.CN p n q) I a =
round_interval prec (compute_bound_poly prec p I a +
mult_float_interval prec (round_interval prec (I ! n - (a ! n))) (compute_bound_poly prec q I a))"
text ‹Bounds on Taylor models are simply a bound on its polynomial, widened by the approximation error.›
fun compute_bound_tm :: "nat ⇒ float interval list ⇒ float interval list ⇒ taylor_model ⇒ float interval"
where "compute_bound_tm prec I a (TaylorModel p e) = compute_bound_poly prec p I a + e"
lemma compute_bound_tm_def:
"compute_bound_tm prec I a tm = compute_bound_poly prec (tm_poly tm) I a + (tm_bound tm)"
by (cases tm) auto
lemma real_of_float_in_real_interval_of[intro, simp]: "real_of_float x ∈⇩r X" if "x ∈⇩i X"
using that
by (auto simp: set_of_eq)
lemma in_set_of_round_interval[intro, simp]:
"x ∈⇩r round_interval prec X" if "x ∈⇩r X"
using round_ivl_correct[of X prec] that
by (auto simp: set_of_eq)
lemma in_set_real_minus_interval[intro, simp]:
"x - y ∈⇩r X - Y" if "x ∈⇩r X" "y ∈⇩r Y"
using that
by (auto simp: set_of_eq)
lemma real_interval_plus: "real_interval (a + b) = real_interval a + real_interval b"
by (simp add: interval_eqI)
lemma real_interval_uminus: "real_interval (- b) = - real_interval b"
by (simp add: interval_eqI)
lemma real_interval_of: "real_interval (interval_of b) = interval_of b"
by (simp add: interval_eqI)
lemma real_interval_minus: "real_interval (a - b) = real_interval a - real_interval b"
using real_interval_plus[of a "-b"] real_interval_uminus[of b]
by (auto simp: interval_eq_iff)
lemma in_set_real_plus_interval[intro, simp]:
"x + y ∈⇩r X + Y" if "x ∈⇩r X" "y ∈⇩r Y"
using that
by (auto simp: set_of_eq)
lemma in_set_neg_plus_interval[intro, simp]:
"- y ∈⇩r - Y" if "y ∈⇩r Y"
using that
by (auto simp: set_of_eq)
lemma in_set_real_times_interval[intro, simp]:
"x * y ∈⇩r X * Y" if "x ∈⇩r X" "y ∈⇩r Y"
using that
by (auto simp: real_interval_times intro!: times_in_intervalI)
lemma real_interval_one: "real_interval 1 = 1"
by (simp add: interval_eqI)
lemma real_interval_zero: "real_interval 0 = 0"
by (simp add: interval_eqI)
lemma real_interval_power: "real_interval (a ^ b) = real_interval a ^ b"
by (induction b arbitrary: a)
(auto simp: real_interval_times real_interval_one)
lemma in_set_real_power_interval[intro, simp]:
"x ^ n ∈⇩r X ^ n" if "x ∈⇩r X"
using that
by (auto simp: real_interval_power intro!: set_of_power_mono)
lemma power_float_interval_real_interval[intro, simp]:
"x ^ n ∈⇩r power_float_interval prec n X" if "x ∈⇩r X"
by (auto simp: real_interval_power that intro!: power_float_intervalI)
lemma in_set_mult_float_interval[intro, simp]:
"x * y ∈⇩r mult_float_interval prec X Y" if "x ∈⇩r X" "y ∈⇩r Y"
using mult_float_interval[of X Y] in_set_real_times_interval[OF that] that(1) that(2)
by blast
lemma in_set_real_minus_swapI: "e i ∈⇩r I ! i - a ! i"
if "x - e i ∈⇩r a ! i" "x ∈⇩r I ! i"
using that
by (auto simp: set_of_eq)
definition develops_at_within::"(nat ⇒ real) ⇒ float interval list ⇒ float interval list ⇒ bool"
where "develops_at_within e a I ⟷ (a all_subset I) ∧ (∀i < length I. e i ∈⇩r I ! i - a ! i)"
lemma develops_at_withinI:
assumes all_in: "a all_subset I"
assumes e: "⋀i. i < length I ⟹ e i ∈⇩r I ! i - a ! i"
shows "develops_at_within e a I"
using assms by (auto simp: develops_at_within_def)
lemma develops_at_withinD:
assumes "develops_at_within e a I"
shows "a all_subset I"
"⋀i. i < length I ⟹ e i ∈⇩r I ! i - a ! i"
using assms by (auto simp: develops_at_within_def)
lemma compute_bound_poly_correct:
fixes p::"float poly"
assumes "num_params p ≤ length I"
assumes dev: "develops_at_within e a I"
shows "insertion e (p::real poly) ∈⇩r compute_bound_poly prec (map_poly interval_of p) I a"
using assms(1)
proof (induction p)
case (C x)
then show ?case by auto
next
case (Bound i)
then show ?case
using dev
by (auto simp: develops_at_within_def)
next
case (Add p1 p2)
then show ?case by force
next
case (Sub p1 p2)
then show ?case by force
next
case (Mul p1 p2)
then show ?case by force
next
case (Neg p)
then show ?case by force
next
case (Pw p x2a)
then show ?case by force
next
case (CN p1 i p2)
then show ?case
using dev
by (auto simp: develops_at_within_def)
qed
lemma compute_bound_tm_correct:
fixes I :: "float interval list" and f :: "real list ⇒ real"
assumes n: "num_params (tm_poly t) ≤ length I"
assumes dev: "develops_at_within e a I"
assumes x0: "x0 ∈⇩i range_tm e t"
shows "x0 ∈⇩r compute_bound_tm prec I a t"
proof -
let ?I = "insertion e (tm_poly t)"
have "x0 = ?I + (x0 - ?I)" by simp
also have "… ∈⇩r compute_bound_tm prec I a t"
unfolding compute_bound_tm_def
apply (rule in_set_real_plus_interval)
apply (rule compute_bound_poly_correct)
apply (rule assms)
apply (rule dev)
using range_tmD[OF x0]
by (auto simp: set_of_eq)
finally show "x0 ∈⇩r compute_bound_tm prec I a t" .
qed
lemma compute_bound_tm_correct_subset:
fixes I :: "float interval list" and f :: "real list ⇒ real"
assumes n: "num_params (tm_poly t) ≤ length I"
assumes dev: "develops_at_within e a I"
shows "set_of (range_tm e t) ⊆ set_of (real_interval (compute_bound_tm prec I a t))"
using assms
by (auto intro!: compute_bound_tm_correct)
lemma compute_bound_poly_mono:
assumes "num_params p ≤ length I"
assumes mem: "I all_subset J" "a all_subset I"
shows "set_of (compute_bound_poly prec p I a) ⊆ set_of (compute_bound_poly prec p J a)"
using assms(1)
proof (induction p arbitrary: a)
case (C x)
then show ?case by auto
next
case (Bound x)
then show ?case using mem
by (simp add: round_interval_mono set_of_sub_inc)
next
case (Add p1 p2)
then show ?case using mem
by (simp add: round_interval_mono set_of_add_inc)
next
case (Sub p1 p2)
then show ?case using mem
by (simp add: round_interval_mono set_of_sub_inc)
next
case (Mul p1 p2)
then show ?case using mem
by (simp add: round_interval_mono mult_float_interval_mono')
next
case (Neg p)
then show ?case using mem
by (simp add: round_interval_mono set_of_neg_inc)
next
case (Pw p x2a)
then show ?case using mem
by (simp add: power_float_interval_mono)
next
case (CN p1 x2a p2)
then show ?case using mem
by (simp add: round_interval_mono mult_float_interval_mono'
set_of_add_inc set_of_sub_inc)
qed
lemma compute_bound_tm_mono:
fixes I :: "float interval list" and f :: "real list ⇒ real"
assumes "num_params (tm_poly t) ≤ length I"
assumes "I all_subset J"
assumes "a all_subset I"
shows "set_of (compute_bound_tm prec I a t) ⊆ set_of (compute_bound_tm prec J a t)"
apply (simp add: compute_bound_tm_def)
apply (rule set_of_add_inc_left)
apply (rule compute_bound_poly_mono)
using assms
by (auto simp: set_of_eq)
subsection ‹Computing taylor models for basic, univariate functions›
definition tm_const :: "float ⇒ taylor_model"
where "tm_const c = TaylorModel (poly.C c) 0"
context includes floatarith_notation begin
definition tm_pi :: "nat ⇒ taylor_model"
where "tm_pi prec = (
let pi_ivl = the (compute_bound_fa prec Pi [])
in TaylorModel (poly.C (mid pi_ivl)) (centered pi_ivl)
)"
lemma zero_real_interval[intro,simp]: "0 ∈⇩r 0"
by (auto simp: set_of_eq)
lemma range_TM_tm_const[simp]: "range_tm e (tm_const c) = interval_of c"
by (auto simp: range_tm_def real_interval_zero tm_const_def)
lemma num_params_tm_const[simp]: "num_params (tm_poly (tm_const c)) = 0"
by (auto simp: tm_const_def)
lemma num_params_tm_pi[simp]: "num_params (tm_poly (tm_pi prec)) = 0"
by (auto simp: tm_pi_def Let_def)
lemma range_tm_tm_pi: "pi ∈⇩i range_tm e (tm_pi prec)"
proof-
have "⋀prec. real_of_float (lb_pi prec) ≤ real_of_float (ub_pi prec)"
using iffD1[OF atLeastAtMost_iff, OF pi_boundaries]
using order_trans by auto
then obtain ivl_pi where ivl_pi_def: "compute_bound_fa prec Pi [] = Some ivl_pi"
by (simp add: approx.simps)
show ?thesis
unfolding range_tm_def Let_def
using compute_bound_fa_correct[OF ivl_pi_def, of "[]"]
by (auto simp: set_of_eq Let_def centered_def ivl_pi_def tm_pi_def
simp del: compute_bound_fa.simps)
qed
subsubsection ‹Derivations of floatarith expressions›
text ‹Compute the nth derivative of a floatarith expression›
fun deriv :: "nat ⇒ floatarith ⇒ nat ⇒ floatarith"
where "deriv v f 0 = f"
| "deriv v f (Suc n) = DERIV_floatarith v (deriv v f n)"
lemma isDERIV_DERIV_floatarith:
assumes "isDERIV v f vs"
shows "isDERIV v (DERIV_floatarith v f) vs"
using assms
proof(induction f)
case (Power f m)
then show ?case
by (cases m) (auto simp: isDERIV_Power)
qed (simp_all add: numeral_eq_Suc add_nonneg_eq_0_iff )
lemma isDERIV_is_analytic:
"isDERIV i (Taylor_Models.deriv i f n) xs"
if "isDERIV i f xs"
using isDERIV_DERIV_floatarith that
by(induction n) auto
lemma deriv_correct:
assumes "isDERIV i f (xs[i:=t])" "i < length xs"
shows "((λx. interpret_floatarith (deriv i f n) (xs[i:=x])) has_real_derivative interpret_floatarith (deriv i f (Suc n)) (xs[i:=t]))
(at t within S)"
apply(simp)
apply (rule has_field_derivative_at_within)
apply(rule DERIV_floatarith)
apply fact
apply (rule isDERIV_is_analytic)
apply fact
done
text ‹Faster derivation for univariate functions, producing smaller terms and thus less over-approximation.›
text ‹TODO: Extend to Arctan, Log!›