Accuracy of Log-Sum-Exp
This example proves a result of accuracy on the computation of the logarithm of a sum of exponentials. It is presented in details in the paper Formally verified bounds on rounding errors in concrete implementations of logarithm-sum-exponential functions.
Authors: Claude Marché / Paul Bonnot
Topics: Floating-Point Computations
Tools: J3
see also the index (by topic, by tool, by reference, by year)
download ZIP archive
The following is an excerpt of the WhyML code where the main theorem if accuracy is stated.
use int.Int use real.RealInfix use real.Abs use real.FromInt use real.ExpLog use udouble.UDouble (* unbounded doubles *) use sum_real.Sum as SumReal use sum.Sum as USum constant exp_max_value :real axiom exp_max_value_spec: 0.0 <. exp_max_value constant exp_error:real axiom exp_error_bound : 0. <. exp_error <=. 0.5
parameter for the error on exponential
function u_exp (x:udouble) : udouble axiom u_exp_spec : forall x:udouble [to_real (u_exp x)]. abs (to_real x) <=. exp_max_value -> abs (to_real (u_exp x) -. exp (to_real x)) <=. exp (to_real x) *. exp_error
the assumed exponential function on unbounded doubles
lemma u_exp_pos: forall x:udouble [to_real (u_exp x)]. abs (to_real x) <=. exp_max_value -> 0.0 <. to_real (u_exp x) constant log_max_value :real axiom log_max_value_spec: 0.0 <. log_max_value constant log_error:real axiom log_error_bound : 0. <. log_error <=. 1.
parameter for the error on logarithm
function u_log (x:udouble) : udouble axiom u_log_spec : forall x:udouble [to_real (u_log x)]. 0.0 <. to_real x <=. log_max_value -> abs (to_real (u_log x) -. log (to_real x)) <=. abs (log (to_real x)) *. log_error
the assumed logarithm function on unbounded doubles
let ghost function exp_fun (a:int -> udouble) : int -> real = fun i -> exp ((USum.real_fun a) i)
[exp_fun a i = exp (a i)]
let ghost function u_exp_fun (a: int -> udouble) : int -> udouble = fun i -> u_exp (a i)
[u_exp_fun a i = exp (a i)] in unbound floating-point
let ghost function u_sum_of_u_exp (a: int -> udouble) (m n:int) : udouble = USum.u_sum (u_exp_fun a) m n
Sum of approximate exponentials [u_sum_of_u_exp(a,n) = exp a(m) + ... + exp a(n-1)]
let ghost function lse_exact (f : int -> udouble) (size : int) = log (SumReal.sum (exp_fun f) 0 size)
The LSE function itself, exact computation [lse(a,n) = log (exp a(0) + ... + exp a(size-1))]
function u_lse (a: int -> udouble) (size:int) : udouble = u_log (u_sum_of_u_exp a 0 size)
The LSE function, floating-point approximation
let lemma lse_accuracy (a:int -> udouble) (size:int) : unit requires { 1 <= size } requires { from_int (size - 1) <=. 0x1p51 } requires { forall i. 0 <= i < size -> abs (to_real (a i)) <=. exp_max_value } requires { exp exp_max_value *. (1.0 +. exp_error) *. from_int size *. (1.0 +. eps *. from_int (size - 1)) <=. log_max_value } ensures { let err = exp_error +. eps *. from_int (size - 1) *. (1. +. exp_error) in abs (to_real (u_lse a size) -. lse_exact a size) <=. log_error *. abs (lse_exact a size) -. log (1. -. err) *. (1. +. log_error) } = ()
The following is an excerpt of the C code with ACSL specification stating pre-conditions for absence of overflow and post-conditions for accuracy.
/*@ // An axiomatic for floating-point exponential and logarithm
@ // the specifications are given in exp_log.mlw, module ExpLogApprox
@
@ axiomatic ExpLogApprox __attribute__ ((j3_theory ("exp_log.ExpLogApprox"))) {
@
@ // Validity range on exponential, assuming 0.0 < exp_max_value
@ logic real exp_max_value __attribute__((j3_symbol("ExpLogApprox.exp_max_value")));
@
@ // Error parameter on exponential, assuming 0. <=. exp_error <=. 0.5
@ logic real exp_error __attribute__((j3_symbol("ExpLogApprox.exp_error")));
@
@ // Error parameter on logarithm, assuming 0. <=. log_error <=. 1.
@ logic real log_error __attribute__((j3_symbol("ExpLogApprox.log_error")));
@
@ // the approximate exponential, on unbounded doubles
@ // spec: for all x, |u_exp(x) - exp(x)| <= exp(x) * exp_error
@ logic udouble u_exp(udouble x) __attribute__((j3_symbol("ExpLogApprox.u_exp")));
@
@ // Validity range on logarithm, assuming 0.0 < log_max_value
@ logic real log_max_value __attribute__((j3_symbol("ExpLogApprox.log_max_value")));
@
@ // the approximate logarithm, on unbounded doubles
@ // spec: for all x, 0 < x ==> |u_log(x) - log(x)| <= exp(x) * exp_error
@ logic udouble u_log(udouble x) __attribute__((j3_symbol("ExpLogApprox.u_log")));
@ }
@*/
/*@ // An axiomatic for LSE
@
@ axiomatic LSE __attribute__ ((j3_theory ("j3bridge.LSE"))) {
@
@ logic udouble u_sum_of_u_exp(double a[MAX_SIZE], integer m, integer n)
@ __attribute__((j3_symbol("LSE.u_sum_of_u_exp")));
@
@ logic udouble u_lse(double a[MAX_SIZE], integer size)
@ __attribute__((j3_symbol("LSE.u_lse")));
@
@ logic real lse_exact(double a[MAX_SIZE], integer size)
@ __attribute__((j3_symbol("LSE.lse_exact")));
@
@ }
@*/
/*@ requires \abs(x) <= exp_max_value;
@ ensures to_udouble(\result) == u_exp(to_udouble(x));
@ assigns \nothing;
@*/
extern double exp_approx(double x);
/*@ requires 0 < x <= log_max_value;
@ ensures to_udouble(\result) == u_log(to_udouble(x));
@ assigns \nothing;
@*/
extern double log_approx(double x);
/*@ requires 0 < size <= MAX_SIZE;
@ requires \initialized (&a[0..size-1]);
@ // to fit exp_approx pre-condition
@ requires \forall integer i; 0 <= i < size ==> \abs(a[i]) <= exp_max_value;
@ // additional requirement to prevent overflow on addition
@ requires exp_max_value <= 701.0;
@ // the hypothesis (6) of Theorem 3.4
@ requires \exp(exp_max_value) * (1.0 + exp_error) * size * (1.0 + (eps * (size - 1)))
<= log_max_value;
@ // additional requirement to prevent overflow on addition
@ requires log_max_value <= 0x1p1023; // close to max double
@ // the result is equal to the WhyML def of LSE on udouble
@ ensures to_udouble(\result) == u_lse(a, size);
@ // the accuracy property
@ ensures \abs(\result - lse_exact(a, size)) <=
@ log_error * \abs(lse_exact(a,size))
@ - \log(1 - (exp_error + eps * (size - 1) * (1 + exp_error))) * (1 +
@ log_error);
@*/
double log_sum_exp(size_t size) {
int i;
double s = 0.0;
/*@ loop invariant 0 <= i <= size;
@ // to prove the first post-condition
@ loop invariant to_udouble(s) == u_sum_of_u_exp(a, 0, i);
@ // for absence of overflow on the sum
@ loop invariant \forall integer j; 0 <= j < i ==>
@ \abs(to_real(u_exp(to_udouble(a[j])))) <= \exp(exp_max_value) * (1.0 + exp_error) ;
@ // for calling log:
@ loop invariant (i == 0 ? s == 0.0 : 0.0 < s);
@ loop assigns i, s;
@ loop variant (size - i);
@*/
for (i = 0; i < size; i++) {
//@ assert 0.0 <= to_real(u_exp(to_udouble(a[i]))) ;
//@ assert to_real(to_udouble(a[i])) <= exp_max_value ;
//@ assert \exp(to_real(to_udouble(a[i]))) <= \exp(exp_max_value) ;
/*@ assert \abs(to_real(u_exp(to_udouble(a[i]))) - \exp(to_real(to_udouble(a[i]))))
@ <= \exp(exp_max_value) * exp_error ;
@*/
//@ assert to_real(u_exp(to_udouble(a[i]))) <= \exp(exp_max_value) * (1.0+exp_error) ;
//@ assert usum_double_bound(u_sum_of_u_exp(a, 0, i), \exp(exp_max_value) * (1.0 + exp_error), size);
s += exp_approx(a[i]);
}
//@ assert usum_double_bound(to_udouble(s), \exp(exp_max_value) * (1.0 + exp_error), size);
return log_approx(s);
}