LAProof.accuracy_proofs.libvalidsdp

Importing a proof from libValidSDP into LAProof

libValidSDP and LAProof each have proofs about the accuracy of linear-algebra operations in floating-point, but they represent floating-point very differently.

LAProof says that ftype t is an IEEE-754 floating-point number with a number of exponent bits and mantissa bits specified by t, and with all the Infinity and NaN behaviors specified by IEEE-754.
libValidSDP says that a floating pointer number is a real number that satisfies the format predicate, where format is a predicate R→bool. The abstraction in libValidSDP makes some things easier to prove, perhaps; in any case, some very useful things are proved there. But we might not want to use the format abstraction globally, because then we couldn't distinguish infinities from NaNs. Because it is proved (in module libvalidsdp.flocq_float) that the IEEE floats are an instance of a legal format, one can import theorems from libValidSDP into LAProof (though not vice versa).

This module is a demonstration of how to do that. The theorem in libValidSDP is cholesky.lemma_2_1, and it is imported here as LVSDP_lemma_2_1.

In addition, there are proofs relating libValidSDP's definitions of Cholesky decomposition to LAProof's.

From LAProof.accuracy_proofs Require Import preamble common solve_model float_acc_lems.
From libValidSDP Require cholesky flocq_float float_spec float_infnan_spec flocq_float binary_infnan.
From LAProof Require accuracy_proofs.mv_mathcomp.

Definition max_mantissa t : positive := Pos.pow 2 (fprecp t) - 1.

Lemma digits_max_mantissa t: Z.pos (SpecFloat.digits2_pos (max_mantissa t)) = fprec t.

Lemma max_mantissa_bounded t: SpecFloat.bounded (fprec t) (femax t) (max_mantissa t) (femax t - fprec t).

Definition Float_max t := Binary.B754_finite (fprec t) (femax t) false _ _ (max_mantissa_bounded t).

Section WithNaN.

Context {NAN: FPCore.Nans} {t : type}.

Definition default_rel : R :=
  / 2 × Raux.bpow Zaux.radix2 (- fprec t + 1).

Definition default_abs : R :=
  / 2 × Raux.bpow Zaux.radix2 (3 - femax t - fprec t).

Lemma prec_lt_emax: @flocq_float.prec (fprecp t) <? femax t.

Notation F := (ftype t).
Notation eps := (default_rel).
Notation eta := (default_abs).

Lemma default_abs_nonzero: default_abs ≠ 0.

Definition fspec := @flocq_float.flocq_float (fprecp t) (femax t) (fprec_gt_one _) prec_lt_emax.

Lemma fspec_eta_nonzero: float_spec.eta fspec ≠ 0.

Definition iszero {t} (x: ftype t) : bool :=
  match x with Binary.B754_zero _ _ _ ⇒ true | _ ⇒ false end.

Fixpoint fsum_l2r_rec [n: nat] (c : F) : F^n → F :=
  match n with
    | 0%N ⇒ fun _ ⇒ c
    | n'.+1 ⇒
      fun a ⇒ fsum_l2r_rec (BPLUS c (a ord0)) [ffun i ⇒ a (lift ord0 i)]
  end.

Definition fcmsum_l2r [n] (c : F) (x : F^n) : F :=
  fsum_l2r_rec c [ffun i ⇒ BOPP (x i)].

Definition stilde [k] (c : F) (a b : F^k) : F :=
  fcmsum_l2r c [ffun i ⇒ BMULT (a i) (b i) : F].

Definition ytilded [k : nat] (c : F) (a b : F^k) (bk : F) :=
  BDIV (stilde c a b) bk.

Definition ytildes [k : nat] (c : F) (a : F^k):=
  BSQRT (stilde c a a).

Lemma BPLUS_B2R_zero (a : ftype t):
  Binary.is_finite a →
  FT2R (BPLUS a (Zconst t 0)) = FT2R a.

Lemma format_FT2R: ∀ (x: ftype t), is_true (@flocq_float.format (fprecp t) (femax t) (FT2R x)).

Definition LVSDP_NAN : binary_infnan.Nans.
Defined.

Definition mkFS (x: F) : float_spec.FS fspec := float_spec.Build_FS_of (format_FT2R x).

Section A.

Import float_infnan_spec.
Definition origFIS :=
  @binary_infnan.binary_infnan LVSDP_NAN (fprecp t) (femax t)
       (fprec_gt_one t) prec_lt_emax.

Definition the_FIS : float_infnan_spec.Float_infnan_spec :=
           {| FIS := origFIS.(FIS);
              FIS0 := origFIS.(FIS0);
              FIS1 := origFIS.(FIS1);
              finite0 := origFIS.(finite0);
              finite1 := origFIS.(finite1);
              fis := fspec;
              m := origFIS.(m);
              m_ge_2 := origFIS.(m_ge_2);
              FIS2FS := mkFS;
              FIS2FS_spec := origFIS.(FIS2FS_spec);
              FIS2FS0 := origFIS.(FIS2FS0);
              FIS2FS1 := origFIS.(FIS2FS1);
              firnd := origFIS.(firnd);
              firnd_spec := origFIS.(firnd_spec);
              firnd_spec_f := origFIS.(firnd_spec_f);
              fiopp := _;
              fiopp_spec := origFIS.(fiopp_spec);
              fiopp_spec_f1 := origFIS.(fiopp_spec_f1);
              fiopp_spec_f := origFIS.(fiopp_spec_f);
              fiplus := _;
              fiplus_spec := origFIS.(fiplus_spec);
              fiplus_spec_fl := origFIS.(fiplus_spec_fl);
              fiplus_spec_fr := origFIS.(fiplus_spec_fr);
              fiplus_spec_f := origFIS.(fiplus_spec_f);
              fiminus := _;
              fiminus_spec := origFIS.(fiminus_spec);
              fiminus_spec_fl := origFIS.(fiminus_spec_fl);
              fiminus_spec_fr := origFIS.(fiminus_spec_fr);
              fiminus_spec_f := origFIS.(fiminus_spec_f);
              fimult:= _;
              fimult_spec := origFIS.(fimult_spec);
              fimult_spec_fl := origFIS.(fimult_spec_fl);
              fimult_spec_fr := origFIS.(fimult_spec_fr);
              fimult_spec_f := origFIS.(fimult_spec_f);
              fidiv := _;
              fidiv_spec := origFIS.(fidiv_spec);
              fidiv_spec_fl := origFIS.(fidiv_spec_fl);
              fidiv_spec_f := origFIS.(fidiv_spec_f);
              fisqrt := _;
              fisqrt_spec := origFIS.(fisqrt_spec);
              fisqrt_spec_f1 := origFIS.(fisqrt_spec_f1);
              fisqrt_spec_f := origFIS.(fisqrt_spec_f);
              ficompare := origFIS.(ficompare);
              ficompare_spec := origFIS.(ficompare_spec);
              ficompare_spec_eq := origFIS.(ficompare_spec_eq);
              ficompare_spec_eq_f := origFIS.(ficompare_spec_eq_f);
          |}.

Definition F' := the_FIS.(FIS).

Lemma FS_val_mkFS: ∀ x: F', float_spec.FS_val (mkFS x) = (FT2R x ).

Lemma FS_val_fplus: ∀ x y: ftype t,
  is_true (the_FIS.(finite) (the_FIS.(fiplus) x y)) →
  float_spec.FS_val (float_spec.fplus (mkFS x) (mkFS y)) = FT2R (BPLUS x y).

Lemma FS_val_fplus': ∀ x y: ftype t,
  Binary.is_finite (BPLUS x y) →
  float_spec.FS_val (float_spec.fplus (mkFS x) (mkFS y)) = FT2R (BPLUS x y).

Lemma FS_val_fmult: ∀ x y: ftype t,
  is_true (the_FIS.(finite) (the_FIS.(fimult) x y)) →
  float_spec.FS_val (float_spec.fmult (mkFS x) (mkFS y)) = FT2R (BMULT x y).

Lemma FS_val_fmult': ∀ x y: ftype t,
  Binary.is_finite (BMULT x y) →
  float_spec.FS_val (float_spec.fmult (mkFS x) (mkFS y)) = FT2R (BMULT x y).

Lemma FS_val_fopp: ∀ x: ftype t,
  is_true (the_FIS.(finite) (the_FIS.(fiopp) x)) →
  float_spec.FS_val (float_spec.fopp (mkFS x)) = FT2R (BOPP x).

Lemma FS_val_fopp': ∀ x: ftype t,
  Binary.is_finite (BOPP x) →
  float_spec.FS_val (float_spec.fopp (mkFS x)) = FT2R (BOPP x).

Lemma FS_val_fdiv: ∀ x y: ftype t,
  is_true (the_FIS.(finite) (the_FIS.(fidiv) x y)) →
  is_true (the_FIS.(finite) y) →
  float_spec.FS_val (float_spec.fdiv (mkFS x) (mkFS y)) = FT2R (BDIV x y).

Lemma FS_val_fdiv': ∀ x y: ftype t,
  Binary.is_finite (BDIV x y) →
  Binary.is_finite y →
  float_spec.FS_val (float_spec.fdiv (mkFS x) (mkFS y)) = FT2R (BDIV x y).

Lemma FS_val_fsqrt': ∀ x: ftype t,
  Binary.is_finite (BSQRT x) →
  float_spec.FS_val (float_spec.fsqrt (mkFS x)) = FT2R (BSQRT x).

Lemma FT2R_congr: ∀ x y : ftype t, feq x y → FT2R x = FT2R y.

Lemma FT2R_feq: ∀ x y: ftype t, Binary.is_finite x → Binary.is_finite y → FT2R x = FT2R y → feq x y.

Lemma FT2R_BDIV_congr: ∀ x x' y y': ftype t,
  Binary.is_finite x → Binary.is_finite x' → Binary.is_finite y → Binary.is_finite y' →
  FT2R x = FT2R x' →
  FT2R y = FT2R y' →
  FT2R (BDIV x y) = FT2R (BDIV x' y').

End A.

Lemma default_abs_eq: eta = float_spec.eta fspec.

Lemma default_rel_eq: eps = float_spec.eps fspec.

Lemma FS_val_ext: ∀ {format} x y,
  @float_spec.FS_val format x = float_spec.FS_val y → x = y.

Lemma BDIV_finite_e: ∀ (x y: ftype t) (H: Binary.is_finite (BDIV x y)), Binary.is_finite x.

Lemma BSQRT_finite_e: ∀ (x: ftype t) (H: Binary.is_finite (BSQRT x)), Binary.is_finite x.

Demonstrations and applications

Everything abovbe is preamble. Now we demonstrate how to relate libValidSDP definitions to corresponding LAProof definitions, and how to import libValidSDP theorems into LAProof.

Lemma fsum_l2r_rec_finite_e: ∀ k (c: ftype t) (a: ftype t ^ k.+1),
  Binary.is_finite (fsum_l2r_rec c a) →
  Binary.is_finite c
  ∧ (∀ i, Binary.is_finite (a i))
  ∧ Binary.is_finite (fsum_l2r_rec (BPLUS c (a ord0)) [ffun i ⇒ a (rshift 1 i)]).

Lemma fsum_l2r_rec_finite_e1: ∀ k (c: ftype t) (a: ftype t ^ k),
  Binary.is_finite (fsum_l2r_rec c a) →
  Binary.is_finite c
  ∧ (∀ i, Binary.is_finite (a i)).

Lemma LVSDP_fcmsum_eq:
∀ [k] (c: F) (a: F^k)
   (FIN: Binary.is_finite (fcmsum_l2r c a)),
   mkFS (fcmsum_l2r c a) = fcmsum.fcmsum_l2r (mkFS c) [ffun i ⇒ float_spec.FS_val (mkFS (a i))].

Lemma LVSDP_stilde_eq: ∀ [k] (a b : F ^ k) (c: F),
    Binary.is_finite (stilde c a b) →
    cholesky.stilde (mkFS c) [ffun i ⇒ mkFS (fun_of_fin a i)] [ffun i ⇒ mkFS (fun_of_fin b i)] = mkFS (stilde c a b).

Lemma LVSDP_ytilded_eq: ∀ [k] (a b : F ^ k) (c bk: F),
    Binary.is_finite bk →
    Binary.is_finite (ytilded c a b bk) →
  float_spec.FS_val (cholesky.ytilded (mkFS c) [ffun i ⇒ mkFS (fun_of_fin a i)] [ffun i ⇒ mkFS (fun_of_fin b i)] (mkFS bk)) =
  FT2R (ytilded c a b bk).

Lemma LVSDP_ytildes_eq: ∀ [k] (a : F ^ k) (c: F),
    Binary.is_finite (ytildes c a) →
  float_spec.FS_val (cholesky.ytildes (mkFS c) [ffun i ⇒ mkFS (fun_of_fin a i)] ) =
  FT2R (ytildes c a).

Lemma LVSDP_lemma_2_1 k (a b : F^k) (c bk : F)
   (Hbk : ¬iszero bk)
   (FINbk: Binary.is_finite bk)
   (FINyt: Binary.is_finite (ytilded c a b bk)):
  Rabs (FT2R bk × FT2R (ytilded c a b bk) - (FT2R c - \sum_i (FT2R (a i) × FT2R (b i))%Re))
  < INR k.+1 × eps × (Rabs (FT2R bk × FT2R (ytilded c a b bk)) + \sum_i Rabs (FT2R(a i) × FT2R(b i)))
    + (1 + INR k.+1 × eps) × (INR k.+1 + Rabs (FT2R bk)) × eta.

Import mv_mathcomp.

Definition cholesky_bound (n: nat) := FT2R (Float_max t) - (eps × INR(2*(n-1)) + INR(n+1)×FT2R(Float_max t)).

Lemma cholesky_jik_spec_backwards_finite:
  ∀ n (A R: 'M[F]_n),
  A^T = A →
  cholesky_jik_spec A R →
  (∀ i j: 'I_n, i ≤ j → Binary.is_finite (R i j)) →
  (∀ i j, Binary.is_finite (A i j)).

Lemma ord_enum_S: ∀ n, ord_enum (S n) = (@inord n 0) :: (map (@inord n \o S \o (@nat_of_ord n)) (ord_enum n)).

Lemma Forall_take_ord_enum: ∀ [n] (u: 'I_n), Forall (fun x: 'I_n ⇒ is_true (x< u)) (take u (ord_enum n)).

Lemma stilde_subtract_loop: ∀ [n] (c: F) (R: 'M_n.+1) (i j: 'I_n.+1) (Hij: (i≤j)%N),
  feq (stilde c [ffun k : 'I_i ⇒ R (inord (nat_of_ord k)) i] [ffun k ⇒ R (inord (nat_of_ord k)) j])
  (subtract_loop_jik c R i j i).

Lemma ytilded_subtract_loop: ∀ n (A R: 'M[F]_n.+1) (i j: 'I_n.+1),
(∀ i j: 'I_n.+1, (i≤j)%N → Binary.is_finite (R i j)) →
   (i<j)%N →
  feq (ytilded (A i j) [ffun k: 'I_i ⇒ R (inord (nat_of_ord k)) i] [ffun k ⇒ R (inord (nat_of_ord k)) j] (R i i)) (BDIV (subtract_loop_jik (A i j) R i j i) (R i i)).

Add Parametric Morphism: BSQRT
with signature @feq t ==> @feq t
as BSQRT_mor.

Lemma ytildes_subtract_loop: ∀ n (A R: 'M[F]_n.+1) (i: 'I_n.+1),
  feq (ytildes (A i i) [ffun k: 'I_i ⇒ R (inord (nat_of_ord k)) i]) (BSQRT (subtract_loop_jik (A i i) R i i i)).

Lemma LVSDP_cholesky_spec: ∀ n (A R: 'M[F]_n.+1),
  A^T = A →
  (∀ i j: 'I_n.+1, (i ≤ j)%N → Binary.is_finite (R i j)) →
  cholesky_jik_spec A R →
  libValidSDP.cholesky.cholesky_spec (map_mx mkFS A) (map_mx mkFS R).

Lemma BSQRT_nonneg: ∀ x:F,
   match BSQRT x with Binary.B754_finite _ _ false _ _ _ ⇒ true
                      | Binary.B754_zero _ _ _ ⇒ true
                      | Binary.B754_infinity _ _ false ⇒ true
                      | Binary.B754_nan _ _ _ _ _ ⇒ true
                      | _ ⇒ false
   end = true.

Main Cholesky equivalence theorem

The predicate cholesky_success A R states that matrix R is the result of running the inner-product-form Cholesky decomposition of matrix A, and the resulting diagonal is all positive. Furthermore, we have proofs that if the C function densemat_cfactor returns 1, then cholesky_success holds. LibValidSDP.cholesky.cholesky_success is a similar predicate (but in libValidSDP's characterization of floating point). LibValidSDP has accuracy theorems, for example if cholesky_success holds, then R^T R is close to A.

This theorem LVSDP_cholesky_success states that LAProof's version of cholesky_success implies LibValidSDP's version, and therefore we can use libValidSDP's accuracy theorems.

Lemma LVDSP_cholesky_success: ∀ [n] (A R: 'M[F]_n.+1),
  A^T = A →
  cholesky_success A R →
  libValidSDP.cholesky.cholesky_success (map_mx mkFS A) (map_mx mkFS R).

End WithNaN.