Theory HOL-Analysis.Higher_Differentiability_Multi

section ‹Higher-Order Differentiability in Several Variables›

text ‹
  ‹k›-times Fréchet differentiability and ‹Ck› smoothness for maps between real normed
  vector spaces, and gradients, Hessians and Jacobians for maps between Euclidean spaces.
  Higher derivatives are iterated as directional derivatives
  ‹λy. frechet_derivative f (at y) v›, which keeps the codomain type fixed.  The
  one-variable notions ‹k_times_differentiable_at› and ‹C_k_on› use ‹deriv ^^ k› instead,
  as Taylor expansions and power series require; for ‹f :: real ⇒ real› the two agree
  (‹k_times_Fr_real_iff›, ‹Ck_on_real_iff›), and for functions of one real variable
  ‹C1› agrees with ‹C1_differentiable_on› of HOL-Analysis.
›

theory Higher_Differentiability_Multi
  imports Higher_Order_Derivatives Cartesian_Euclidean_Space
begin

subsection ‹Multi-dimensional ‹k›-times Fréchet differentiability at a point›

text ‹
  Differentiability without continuity: ‹f› is ‹(Suc k)›-times differentiable at ‹x› if
  it is ‹k›-times differentiable near ‹x›, differentiable at ‹x›, and every directional
  derivative ‹λy. frechet_derivative f (at y) v› is ‹k›-times differentiable at ‹x›.
›

primrec k_times_Fr_differentiable_at
  :: "nat ⇒ ('a::real_normed_vector ⇒ 'b::real_normed_vector) ⇒ 'a ⇒ bool"
where
  "k_times_Fr_differentiable_at 0 f x ⟷ True"
| "k_times_Fr_differentiable_at (Suc k) f x ⟷
     (∃A. open A ∧ x ∈ A ∧ (∀y∈A. k_times_Fr_differentiable_at k f y))
   ∧ f differentiable (at x)
   ∧ (∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x)"

text ‹‹1›-times differentiable is Fréchet differentiable.›

lemma one_times_Fr_iff:
  "k_times_Fr_differentiable_at 1 f x ⟷ f differentiable (at x)"
  by auto

text ‹Monotonicity: higher differentiability implies lower.›

lemma k_times_Fr_differentiable_at_mono:
  assumes "m ≤ k" and "k_times_Fr_differentiable_at k f x"
  shows   "k_times_Fr_differentiable_at m f x"
  using assms
proof (induction k arbitrary: m f x)
  case 0
  then have "m = 0" by simp
  then show ?case by simp
next
  case (Suc k)
  note IH = Suc.IH
  note asm = Suc.prems

  show ?case
  proof (cases m)
    case 0
    then show ?thesis by simp
  next
    case (Suc m')
    from asm(1) Suc have m'_le: "m' ≤ k"
      by simp

    from asm(2) obtain A where
      A: "open A"
         "x ∈ A"
         "∀y∈A. k_times_Fr_differentiable_at k f y"
      and fdiff: "f differentiable (at x)"
      and D: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
      by auto

    have A': "∀y∈A. k_times_Fr_differentiable_at m' f y"
      using A(3) IH[OF m'_le] by blast

    have D': "∀v. k_times_Fr_differentiable_at m' (λy. frechet_derivative f (at y) v) x"
      using D IH[OF m'_le] by blast

    show ?thesis
      using Suc A fdiff D' A'
      by (metis IH asm(1,2) le_Suc_eq)
  qed
qed

text ‹Peeling off the top layer.›

lemma k_times_Fr_differentiable_at_SucD:
  assumes "k_times_Fr_differentiable_at (Suc k) f x"
  shows   "k_times_Fr_differentiable_at k f x"
    and   "f differentiable (at x)"
    and   "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
  using assms k_times_Fr_differentiable_at_mono
  by auto

text ‹The derivative field inherits differentiability.›

lemma k_times_Fr_differentiable_at_derivative:
  assumes "k_times_Fr_differentiable_at (Suc k) f x"
  shows   "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
  using assms by simp


subsection ‹Set-wise ‹k›-times Fréchet differentiability›

definition k_times_Fr_differentiable_on
  :: "nat ⇒ ('a::real_normed_vector ⇒ 'b::real_normed_vector) ⇒ 'a set ⇒ bool"
where
  "k_times_Fr_differentiable_on k f S ⟷ (∀x∈S. k_times_Fr_differentiable_at k f x)"

lemma k_times_Fr_differentiable_onI:
  "(⋀x. x ∈ S ⟹ k_times_Fr_differentiable_at k f x) ⟹ k_times_Fr_differentiable_on k f S"
  by (simp add: k_times_Fr_differentiable_on_def)

lemma k_times_Fr_differentiable_onD:
  "k_times_Fr_differentiable_on k f S ⟹ x ∈ S ⟹ k_times_Fr_differentiable_at k f x"
  by (simp add: k_times_Fr_differentiable_on_def)

lemma k_times_Fr_differentiable_on_mono:
  "m ≤ k ⟹ k_times_Fr_differentiable_on k f S ⟹ k_times_Fr_differentiable_on m f S"
  by (simp add: k_times_Fr_differentiable_on_def k_times_Fr_differentiable_at_mono)

lemma k_times_Fr_differentiable_on_subset:
  "S ⊆ T ⟹ k_times_Fr_differentiable_on k f T ⟹ k_times_Fr_differentiable_on k f S"
  by (simp add: k_times_Fr_differentiable_on_def subset_iff)


subsection ‹‹Ck› at a point (with continuity)›

text ‹
  ‹Ck_at k f x›: for ‹k = 0›, ‹f› is continuous at ‹x›; for ‹Suc n›, ‹f› is ‹Cn›
  near ‹x›, differentiable at ‹x›, and every directional derivative is ‹Cn› at ‹x›.
›

primrec Ck_at
  :: "nat ⇒ ('a::real_normed_vector ⇒ 'b::real_normed_vector) ⇒ 'a ⇒ bool"
where
  "Ck_at 0 f x ⟷ continuous (at x) f"
| "Ck_at (Suc k) f x ⟷
     (∃A. open A ∧ x ∈ A ∧ (∀y∈A. Ck_at k f y))
   ∧ f differentiable (at x)
   ∧ (∀v. Ck_at k (λy. frechet_derivative f (at y) v) x)"


subsection ‹‹Ck› on an open set›

text ‹‹Ck› on an open set; the multi-dimensional version of @{const C_k_on}.›

definition Ck_on
  :: "nat ⇒ ('a::real_normed_vector ⇒ 'b::real_normed_vector) ⇒ 'a set ⇒ bool"
where
  "Ck_on k f U ⟷ open U ∧ (∀x∈U. Ck_at k f x)"


subsection ‹Relationships between the three notions›

text ‹‹Ck› implies ‹k›-times differentiable (forgetting continuity).›

lemma Ck_at_imp_k_times_Fr:
  "Ck_at k f x ⟹ k_times_Fr_differentiable_at k f x"
  by (induction k arbitrary: f x) auto

corollary Ck_on_imp_k_times_Fr_on:
  "Ck_on k f U ⟹ k_times_Fr_differentiable_on k f U"
  by (simp add: Ck_on_def k_times_Fr_differentiable_on_def Ck_at_imp_k_times_Fr)


subsection ‹Basic properties of ‹Ck› at a point and on an open set›

lemma Ck_atI_Suc:
  assumes "open A" and "x ∈ A" and "⋀y. y ∈ A ⟹ Ck_at k f y"
    and "f differentiable (at x)"
    and "⋀v. Ck_at k (λy. frechet_derivative f (at y) v) x"
  shows "Ck_at (Suc k) f x"
  unfolding Ck_at.simps(2) using assms by blast

lemma Ck_at_SucE:
  assumes "Ck_at (Suc k) f x"
  obtains A where "open A" and "x ∈ A" and "∀y∈A. Ck_at k f y"
    and "f differentiable (at x)"
    and "∀v. Ck_at k (λy. frechet_derivative f (at y) v) x"
  using assms unfolding Ck_at.simps(2) by blast

lemma Ck_at_imp_continuous: "Ck_at k f x ⟹ continuous (at x) f"
  by (cases k) (auto intro: differentiable_imp_continuous_within)

lemma Ck_at_SucD:
  assumes "Ck_at (Suc k) f x"
  shows "Ck_at k f x"
  using assms
proof (induction k arbitrary: f x)
  case 0
  then show ?case
    by (auto intro: differentiable_imp_continuous_within)
next
  case (Suc k)
  from Suc.prems obtain A where A: "open A" "x ∈ A" "∀y∈A. Ck_at (Suc k) f y"
    and d: "f differentiable (at x)"
    and D: "∀v. Ck_at (Suc k) (λy. frechet_derivative f (at y) v) x"
    by (rule Ck_at_SucE)
  show ?case
  proof (rule Ck_atI_Suc[OF A(1,2) _ d])
    show "Ck_at k f y" if "y ∈ A" for y
      using A(3) that by (blast intro: Suc.IH)
    show "Ck_at k (λy. frechet_derivative f (at y) v) x" for v
      using D by (blast intro: Suc.IH)
  qed
qed

lemma Ck_at_mono:
  assumes "Ck_at k f x" and "m ≤ k"
  shows "Ck_at m f x"
  using assms
proof (induction k)
  case 0
  then show ?case by simp
next
  case (Suc k)
  then show ?case
    using Ck_at_SucD le_Suc_eq by blast
qed

text ‹A transfer principle: ‹Ck_at k f x› only depends on the values of ‹f› on an open
  neighbourhood of ‹x›.›

lemma Ck_at_transfer_open:
  assumes "open U" and "x ∈ U" and "⋀y. y ∈ U ⟹ f y = g y"
    and "Ck_at k f x"
  shows "Ck_at k g x"
  using assms
proof (induction k arbitrary: f g x)
  case 0
  have "∀F y in nhds x. f y = g y"
    using "0.prems"(1-3) by (auto simp: eventually_nhds)
  then show ?case
    using "0.prems"(4) isCont_cong by auto
next
  case (Suc k)
  from Suc.prems(4) obtain A where A: "open A" "x ∈ A" "∀y∈A. Ck_at k f y"
    and d: "f differentiable (at x)"
    and D: "∀v. Ck_at k (λy. frechet_derivative f (at y) v) x"
    by (rule Ck_at_SucE)
  have eq_der: "frechet_derivative f (at y) = frechet_derivative g (at y)" if "y ∈ U" for y
    using Suc.prems using frechet_derivative_transform_within_open that by blast
  show ?case
  proof (rule Ck_atI_Suc)
    show "open (A ∩ U)" and "x ∈ A ∩ U"
      using A(1,2) Suc.prems(1,2) by auto
    show "Ck_at k g y" if "y ∈ A ∩ U" for y
      using that A(3) Suc.IH[where f = f and g = g, OF Suc.prems(1) _ Suc.prems(3)] by blast
    show "g differentiable (at x)"
      using d Suc.prems(1-3) by (rule differentiable_transform_within_open)
    show "Ck_at k (λy. frechet_derivative g (at y) v) x" for v
      by (rule Suc.IH[OF Suc.prems(1,2) _ D[rule_format, of v]]) (simp add: eq_der)
  qed
qed

lemma Ck_onI: "open U ⟹ (⋀x. x ∈ U ⟹ Ck_at k f x) ⟹ Ck_on k f U"
  by (simp add: Ck_on_def)

lemma Ck_on_open: "Ck_on k f U ⟹ open U"
  by (simp add: Ck_on_def)

lemma Ck_onD: "Ck_on k f U ⟹ x ∈ U ⟹ Ck_at k f x"
  by (simp add: Ck_on_def)

lemma Ck_on_mono:
  assumes "Ck_on k f U" and "m ≤ k"
  shows "Ck_on m f U"
  using assms unfolding Ck_on_def by (blast intro: Ck_at_mono)

lemma Ck_on_SucD: "Ck_on (Suc k) f U ⟹ Ck_on k f U"
  by (erule Ck_on_mono) simp

lemma Ck_on_subset:
  assumes "Ck_on k f U" and "open V" and "V ⊆ U"
  shows "Ck_on k f V"
  using assms unfolding Ck_on_def by blast

lemma Ck_on_imp_continuous_on: "Ck_on k f U ⟹ continuous_on U f"
  unfolding Ck_on_def by (auto simp: continuous_on_eq_continuous_at intro: Ck_at_imp_continuous)

lemma Ck_on_cong:
  assumes "U = V" and "⋀x. x ∈ V ⟹ f x = g x"
  shows "Ck_on k f U ⟷ Ck_on k g V"
proof -
  have sym: "⋀x. x ∈ V ⟹ g x = f x"
    using assms(2) by simp
  have "Ck_at k f x ⟷ Ck_at k g x" if V: "open V" "x ∈ V" for x
    using Ck_at_transfer_open[of V x f g k, OF V assms(2)]
      Ck_at_transfer_open[of V x g f k, OF V sym] by blast
  then show ?thesis
    using assms(1) unfolding Ck_on_def by blast
qed

lemma Ck_on_congI:
  assumes "Ck_on k g U" and "⋀x. x ∈ U ⟹ f x = g x"
  shows "Ck_on k f U"
  using Ck_on_cong[of U U f g k, OF refl assms(2)] assms(1) by simp

text ‹On open sets, ‹Ck› satisfies the recursion below.  The closure properties are proved by
  induction along it, following the treatment of ‹higher_differentiable_on› in the AFP
  entry ∗‹Smooth Manifolds› by Immler and Zhan.›

lemma Ck_on_0_iff: "Ck_on 0 f U ⟷ open U ∧ continuous_on U f"
  by (auto simp: Ck_on_def continuous_on_eq_continuous_at)

lemma Ck_on_of_derivatives:
  assumes U: "open U"
    and d: "⋀x. x ∈ U ⟹ f differentiable (at x)"
    and D: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
  shows "Ck_on k f U"
  using d D
proof (induction k arbitrary: f)
  case 0
  then show ?case
    using U by (auto simp: Ck_on_def intro: differentiable_imp_continuous_within)
next
  case (Suc k)
  have Ck: "Ck_on k f U"
    by (rule Suc.IH) (use Suc.prems in ‹auto intro: Ck_on_SucD›)
  show ?case
  proof (rule Ck_onI[OF U])
    fix x assume x: "x ∈ U"
    show "Ck_at (Suc k) f x"
    proof (rule Ck_atI_Suc[OF U x])
      show "Ck_at k f y" if "y ∈ U" for y
        using Ck that by (rule Ck_onD)
      show "f differentiable (at x)"
        using x by (rule Suc.prems(1))
      show "Ck_at k (λy. frechet_derivative f (at y) v) x" for v
        by (rule Ck_at_SucD[OF Ck_onD[OF Suc.prems(2) x]])
    qed
  qed
qed

lemma Ck_on_Suc_iff:
  "Ck_on (Suc k) f U ⟷
     open U ∧ (∀x∈U. f differentiable (at x)) ∧
     (∀v. Ck_on k (λx. frechet_derivative f (at x) v) U)"
proof
  assume "Ck_on (Suc k) f U"
  then show "open U ∧ (∀x∈U. f differentiable (at x)) ∧
      (∀v. Ck_on k (λx. frechet_derivative f (at x) v) U)"
    by (auto simp: Ck_on_def)
next
  assume "open U ∧ (∀x∈U. f differentiable (at x)) ∧
      (∀v. Ck_on k (λx. frechet_derivative f (at x) v) U)"
  then have U: "open U" and d: "⋀x. x ∈ U ⟹ f differentiable (at x)"
    and D: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
    by auto
  have Ck: "Ck_on k f U"
    using U d D by (rule Ck_on_of_derivatives)
  show "Ck_on (Suc k) f U"
  proof (rule Ck_onI[OF U])
    fix x assume x: "x ∈ U"
    show "Ck_at (Suc k) f x"
      by (rule Ck_atI_Suc[OF U x Ck_onD[OF Ck] d[OF x] Ck_onD[OF D x]])
  qed
qed


subsection ‹Bridge to the one-dimensional theory›

text ‹For ‹f :: real ⇒ real›, the two notions of ‹k›-times differentiability agree.›


text ‹Being ‹k›-times differentiable at ‹x› only depends on the function near ‹x›.›

lemma k_times_differentiable_at_transfer_open:
  fixes f g :: "real ⇒ real"
  assumes U: "open U" "x ∈ U"
    and eq: "⋀y. y ∈ U ⟹ f y = g y"
    and Hf: "k_times_differentiable_at k f x"
  shows "k_times_differentiable_at k g x"
  using U eq Hf
proof (induction k arbitrary: f g x U)
  case 0
  then show ?case by simp
next
  case (Suc k)

  text ‹Unfold the definition at ‹x›.›
  from Suc.prems(4) obtain ε where εpos: "ε > 0"
    and ball_f: "⋀y. ¦y - x¦ < ε ⟹ k_times_differentiable_at k f y"
    and der_f: "((deriv ^^ k) f
                   has_derivative (λh. (deriv ^^ Suc k) f x * h)) (at x)"
    by auto

  text ‹Shrink the ball so that it sits inside ‹U›, where ‹f = g›.›
  obtain δ where δpos: "δ > 0" and ballU: "ball x δ ⊆ U"
    using Suc.prems(1,2) open_contains_ball by blast
  define r where "r = min ε δ"
  have rpos: "r > 0" using εpos δpos by (simp add: r_def)

  define B where "B = ball x r"
  have openB: "open B" and xB: "x ∈ B" using rpos by (auto simp: B_def)
  have BsubU: "B ⊆ U"
    using ballU by (auto simp: B_def r_def dist_real_def abs_minus_commute)
  have eqB: "∀y∈B. f y = g y" using BsubU Suc.prems(3) by auto

  text ‹On ‹B› the function ‹f› is ‹k›-times differentiable.›
  have f_on_B: "f k-times_differentiable_on B"
    by (rule k_times_differentiable_onI)
       (auto simp: B_def r_def dist_real_def abs_minus_commute intro!: ball_f)

  text ‹Transfer this to ‹g› and obtain agreement of the lower derivatives.›
  have g_on_B: "g k-times_differentiable_on B"
   and der_agree: "∀y∈B. ∀m<k. ((deriv ^^ m) g
                       has_derivative (*) ((deriv ^^ Suc m) f y)) (at y)"
    using times_differentiable_on_transfer[OF openB f_on_B eqB] by blast+

  text ‹Part 1: the ‹ε›-ball condition for ‹g› (radius ‹r›).›
  have ball_g: "∀y. ¦y - x¦ < r ⟶ k_times_differentiable_at k g y"
    using g_on_B
    by (auto simp: k_times_differentiable_on_def B_def dist_real_def abs_minus_commute)

  text ‹Part 2: the ‹k›-th derivatives of ‹f› and ‹g› coincide on ‹B›.›
  have kth_eq: "∀y∈B. (deriv ^^ k) f y = (deriv ^^ k) g y"
  proof (cases k)
    case 0
    then show ?thesis using eqB by simp
  next
    case (Suc n)
    show ?thesis
    proof
      fix y assume yB: "y ∈ B"
      have "((deriv ^^ n) g has_derivative (*) ((deriv ^^ Suc n) f y)) (at y)"
        using der_agree yB Suc by simp
      hence "deriv ((deriv ^^ n) g) y = (deriv ^^ Suc n) f y"
        by (rule deriv_eq)
      thus "(deriv ^^ k) f y = (deriv ^^ k) g y"
        using Suc by simp
    qed
  qed

  text ‹Transfer the derivative condition at ‹x› from ‹f› to ‹g›.›
  have der_f': "((deriv ^^ k) f has_derivative (*) ((deriv ^^ Suc k) f x)) (at x)"
    using der_f by simp
  have der_g': "((deriv ^^ k) g has_derivative (*) ((deriv ^^ Suc k) f x)) (at x)"
    using has_derivative_transfer_on_open[OF openB xB _ der_f'] kth_eq by blast
  have kSuc_eq: "(deriv ^^ Suc k) g x = (deriv ^^ Suc k) f x"
    using der_g' by (simp add: deriv_eq)

  have der_g: "((deriv ^^ k) g
                  has_derivative (λh. (deriv ^^ Suc k) g x * h)) (at x)"
    using der_g' kSuc_eq by simp

  show "k_times_differentiable_at (Suc k) g x"
    using εpos rpos ball_g der_g by auto
qed

lemma eq_on_open_k_times_differentiable_at:
  fixes f g :: "real ⇒ real"
  assumes U: "open U" "x ∈ U"
    and eq: "⋀y. y ∈ U ⟹ f y = g y"
  shows "k_times_differentiable_at k f x ⟷ k_times_differentiable_at k g x"
  using k_times_differentiable_at_transfer_open[OF U eq]
        k_times_differentiable_at_transfer_open[OF U(1,2), of g f k] eq
  by auto

lemma k_times_differentiable_at_cmult:
  fixes f :: "real ⇒ real"
  shows "k_times_differentiable_at k f x ⟹
         k_times_differentiable_at k (λy. c * f y) x"
  by (rule kth_deriv_cmultE)



lemma k_times_Fr_real_iff:
  fixes f :: "real ⇒ real"
  shows "k_times_Fr_differentiable_at k f x ⟷ k_times_differentiable_at k f x"
proof (induction k arbitrary: f x)
  case 0
  then show ?case by simp
next
  case (Suc k)
  show ?case
  proof
    assume H: "k_times_Fr_differentiable_at (Suc k) f x"

    from H obtain A where
      A: "open A" "x ∈ A" "∀y∈A. k_times_Fr_differentiable_at k f y"
      and df: "f differentiable (at x)"
      and D: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
      unfolding k_times_Fr_differentiable_at.simps(2)
      by blast

    have neigh: "∀y∈A. k_times_differentiable_at k f y"
      using A(3) Suc.IH by blast

    have dk: "k_times_differentiable_at k (deriv f) x"
    proof (cases k)
      case 0
      then show ?thesis by simp
    next
      case (Suc j)

      have H1: "k_times_Fr_differentiable_at (Suc j)
                  (λy. frechet_derivative f (at y) 1) x"
        using D Suc by simp
      then have dk1:
        "k_times_differentiable_at (Suc j)
           (λy. frechet_derivative f (at y) 1) x"
        using Suc.IH Suc by blast

      have eq_deriv: "⋀y. y ∈ A ⟹ frechet_derivative f (at y) 1 = deriv f y"
        using A(3) Suc frechet_derivative_one_eq_deriv k_times_Fr_differentiable_at.simps(2) by blast


      have "k_times_differentiable_at (Suc j)
              (λy. frechet_derivative f (at y) 1) x ⟷
            k_times_differentiable_at (Suc j) (deriv f) x"
        using A(1,2) eq_deriv eq_on_open_k_times_differentiable_at by presburger
      then show ?thesis
        using Suc dk1 by blast
    qed

    show "k_times_differentiable_at (Suc k) f x"
    proof -
      obtain ε where εpos: "ε > 0" and ballA: "ball x ε ⊆ A"
        using A(1,2) open_contains_ball by blast
      have part1: "∀y. ¦y - x¦ < ε ⟶ k_times_differentiable_at k f y"
        using ballA neigh by (auto simp: dist_real_def abs_minus_commute,
                               simp add: dist_norm subsetD)
      have diffk: "(deriv ^^ k) f differentiable (at x)"
      proof (cases k)
        case 0
        then show ?thesis using df by simp
      next
        case (Suc j)
        from dk Suc have d: "(deriv f) (Suc j)-times_differentiable_at x" by simp
        have "((deriv ^^ j) (deriv f)
                 has_real_derivative (deriv ^^ Suc j) (deriv f) x) (at x)"
          using k_times_differentiable_at_le_deriv(2)[OF d lessI] .
        then have "(deriv ^^ j) (deriv f) differentiable (at x)"
          using real_differentiable_def by blast
        then show ?thesis
          using Suc kth_deriv_shift by metis
      qed
      have part2: "((deriv ^^ k) f
                      has_derivative (λh. (deriv ^^ Suc k) f x * h)) (at x)"
      proof -
        from diffk have "((deriv ^^ k) f
                 has_real_derivative deriv ((deriv ^^ k) f) x) (at x)"
          using DERIV_deriv_iff_real_differentiable by blast
        then show ?thesis
          by (simp add: has_field_derivative_def)
      qed
      show ?thesis
        using εpos part1 part2 by auto
    qed
  next
    assume H: "k_times_differentiable_at (Suc k) f x"

    from H obtain ε where εpos: "ε > 0"
      and ball_f: "⋀y. ¦y - x¦ < ε ⟹ k_times_differentiable_at k f y"
      by auto
    define A where "A = ball x ε"
    have A: "open A" "x ∈ A" "∀y∈A. k_times_differentiable_at k f y"
      using εpos ball_f by (auto simp: A_def dist_real_def abs_minus_commute)
    have df: "f differentiable (at x)"
    proof -
      have "f 1-times_differentiable_at x"
        using H k_times_differentiable_at_mono[of 1 "Suc k" f x] by simp
      then show ?thesis
        by (metis one_time_differentiable_at_iff real_differentiable_def)
    qed
    have dk: "k_times_differentiable_at k (deriv f) x"
      using k_times_differentiable_at_derivative[OF H] by simp

    have neigh: "∀y∈A. k_times_Fr_differentiable_at k f y"
      using A(3) Suc.IH by blast

    have D: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
    proof
      fix v
      show "k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
      proof (cases k)
        case 0
        then show ?thesis by simp
      next
        case (Suc j)

        have dk': "k_times_differentiable_at (Suc j) (λy. v * deriv f y) x"
          using dk Suc k_times_differentiable_at_cmult[of "Suc j" "deriv f" x v]
          by blast
        hence dk'': "k_times_differentiable_at (Suc j) (λy. deriv f y * v) x"
          by (simp add: mult.commute)

        have eq_deriv: "⋀y. y ∈ A ⟹ frechet_derivative f (at y) v = deriv f y * v"
          using Suc frechet_derivative_to_deriv neigh by auto

        have "k_times_differentiable_at (Suc j)
                (λy. frechet_derivative f (at y) v) x ⟷
              k_times_differentiable_at (Suc j) (λy. deriv f y * v) x"
          using A(1,2) eq_deriv eq_on_open_k_times_differentiable_at by presburger
        then have "k_times_differentiable_at (Suc j)
                     (λy. frechet_derivative f (at y) v) x"
          using dk'' by blast
        then show ?thesis
          using Suc Suc.IH by blast
      qed
    qed

    show "k_times_Fr_differentiable_at (Suc k) f x"
      unfolding k_times_Fr_differentiable_at.simps(2)
      using A(1,2) neigh df D by blast
  qed
qed



subsection ‹Basic closure properties›

text ‹These generalise the one-dimensional closure results for @{const C_k_on}.›

lemma k_times_Fr_const:
  "k_times_Fr_differentiable_at k (λ_. c) x"
proof (induction k arbitrary: x c)
  case 0
  then show ?case
    by simp
next
  case (Suc k)
  have "∃A. open A ∧ x ∈ A ∧ (∀y∈A. k_times_Fr_differentiable_at k (λ_. c) y)"
    using Suc.IH by (intro exI[of _ UNIV]) auto
  moreover have "(λ_. c) differentiable (at x)"
    by simp
  moreover have "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative (λ_. c) (at y) v) x"
    by (simp add: Suc)
  ultimately show ?case
    by simp
qed

lemma k_times_Fr_id:
  "k_times_Fr_differentiable_at k (λx. x) x"
proof (induction k arbitrary: x)
  case 0
  then show ?case
    by simp
next
  case (Suc k)
  have nbhd:
    "∃A. open A ∧ x ∈ A ∧ (∀y∈A. k_times_Fr_differentiable_at k (λx. x) y)"
    using Suc.IH by (intro exI[of _ UNIV]) auto
  have diff: "(λx. x) differentiable (at x)"
    by simp
  have derivs: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative (λx. x) (at y) v) x"
    by (simp add: k_times_Fr_const)
  show ?case
    using nbhd diff derivs
    by simp
qed

lemma Ck_at_const:
  "Ck_at k (λ_. c) x"
proof (induction k arbitrary: x c)
  case 0
  then show ?case
    by simp
next
  case (Suc k)
  have nbhd: "∃A. open A ∧ x ∈ A ∧ (∀y∈A. Ck_at k (λ_. c) y)"
    using Suc.IH by (intro exI[of _ UNIV]) auto

  have diff: "(λ_. c) differentiable (at x)"
    by simp

  have derivs: "∀v. Ck_at k (λy. frechet_derivative (λ_. c) (at y) v) x"
    by (simp add: Suc)

  show ?case
    using nbhd diff derivs
    by simp
qed

lemma Ck_on_const:
  "open U ⟹ Ck_on k (λ_. c) U"
  by (simp add: Ck_on_def Ck_at_const)

text ‹
  As in the one-dimensional case, the forward implication is proved by
  induction and the equivalence follows by symmetry of the hypotheses.
›

lemma k_times_Fr_differentiable_at_transfer_open:
  fixes f g :: "'a::real_normed_vector ⇒ 'b::real_normed_vector"
  assumes U: "open U" "x ∈ U"
    and eq: "⋀y. y ∈ U ⟹ f y = g y"
    and Hf: "k_times_Fr_differentiable_at k f x"
  shows "k_times_Fr_differentiable_at k g x"
  using U eq Hf
proof (induction k arbitrary: f g x U)
  case 0
  then show ?case by simp
next
  case (Suc k)

  from Suc.prems(4) obtain A where
    A: "open A" "x ∈ A" "∀y∈A. k_times_Fr_differentiable_at k f y"
    and df: "f differentiable (at x)"
    and Df: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
    unfolding k_times_Fr_differentiable_at.simps(2)
    by blast

  let ?C = "A ∩ U"
  have C: "open ?C" "x ∈ ?C"
    using A Suc.prems by auto

  have neigh: "∀y∈?C. k_times_Fr_differentiable_at k g y"
    by (metis A(3) Int_iff Suc.IH Suc.prems(1,3))
  have evx: "eventually (λy. y ∈ U) (nhds x)"
    using Suc.prems(1,2) by (simp add: eventually_nhds, auto)
  have evx_fg: "eventually (λy. f y = g y) (nhds x)"
    by (rule eventually_mono[OF evx]) (use Suc.prems(3) in auto)


  have dg: "g differentiable (at x)"
    by (metis Suc.prems(1,2,3) df differentiable_transform_within_open)

  have Dg: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative g (at y) v) x"
  proof
    fix v
    show "k_times_Fr_differentiable_at k (λy. frechet_derivative g (at y) v) x"
    proof (cases k)
      case 0
      then show ?thesis by simp
    next
      case (Suc j)

      have eqD:
        "⋀y. y ∈ ?C ⟹ frechet_derivative f (at y) v = frechet_derivative g (at y) v"
      proof -
        fix y
        assume yC: "y ∈ ?C"
        hence yA: "y ∈ A" and yU: "y ∈ U"
          by auto

        have fy: "k_times_Fr_differentiable_at (Suc j) f y"
          using A(3) Suc yA by blast

        hence dfy: "f differentiable (at y)"
          using k_times_Fr_differentiable_at_mono[of 1 "Suc j" f y]
          by (simp add: one_times_Fr_iff)

        have gy: "k_times_Fr_differentiable_at (Suc j) g y"
          using Suc neigh yC by blast

        hence dgy: "g differentiable (at y)"
          using k_times_Fr_differentiable_at_mono[of 1 "Suc j" g y]
          by (simp add: one_times_Fr_iff)

        have evy: "eventually (λz. f z = g z) (nhds y)"
          using Suc.prems(1) yU Suc.prems(3)
          by (simp add: eventually_nhds, auto)

        have "(f has_derivative frechet_derivative f (at y)) (at y)"
          by (simp add: dfy frechet_derivative_works[THEN iffD1])
        then have "(g has_derivative frechet_derivative f (at y)) (at y)"
          using Suc.prems(1,3) has_derivative_transfer_on_open yU by blast
        moreover have "(g has_derivative frechet_derivative g (at y)) (at y)"
          using dgy frechet_derivative_works[THEN iffD1] by blast
        ultimately have "frechet_derivative f (at y) = frechet_derivative g (at y)"
          by (rule has_derivative_unique)
        then show "frechet_derivative f (at y) v = frechet_derivative g (at y) v"
          by simp
      qed

      have "k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
        using Df by blast
      then show ?thesis
        using Suc.IH[OF C(1,2), of
            "λy. frechet_derivative f (at y) v"
            "λy. frechet_derivative g (at y) v"]
          eqD
        by blast
    qed
  qed

  show "k_times_Fr_differentiable_at (Suc k) g x"
    unfolding k_times_Fr_differentiable_at.simps(2)
    using C neigh dg Dg by blast
qed

lemma eq_on_open_k_times_Fr_differentiable_at:
  fixes f g :: "'a::real_normed_vector ⇒ 'b::real_normed_vector"
  assumes U: "open U" "x ∈ U"
    and eq: "⋀y. y ∈ U ⟹ f y = g y"
  shows "k_times_Fr_differentiable_at k f x ⟷ k_times_Fr_differentiable_at k g x"
  using k_times_Fr_differentiable_at_transfer_open[OF U eq]
        k_times_Fr_differentiable_at_transfer_open[OF U(1,2), of g f k] eq
  by auto

lemma k_times_Fr_add:
  assumes "k_times_Fr_differentiable_at k f x"
      and "k_times_Fr_differentiable_at k g x"
  shows "k_times_Fr_differentiable_at k (λy. f y + g y) x"
  using assms
proof (induction k arbitrary: f g x)
  case 0
  then show ?case
    by simp
next
  case (Suc k)
  from Suc.prems(1) obtain A where
    A: "open A" "x ∈ A" "∀y∈A. k_times_Fr_differentiable_at k f y"
    and df: "f differentiable (at x)"
    and Df: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
    unfolding k_times_Fr_differentiable_at.simps(2)
    by blast

  from Suc.prems(2) obtain B where
    B: "open B" "x ∈ B" "∀y∈B. k_times_Fr_differentiable_at k g y"
    and dg: "g differentiable (at x)"
    and Dg: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative g (at y) v) x"
    unfolding k_times_Fr_differentiable_at.simps(2)
    by blast

  let ?C = "A ∩ B"
  have C: "open ?C" "x ∈ ?C"
    using A B by auto

  have neigh: "∀y∈?C. k_times_Fr_differentiable_at k (λz. f z + g z) y"
  proof
    fix y
    assume yC: "y ∈ ?C"
    then have yA: "y ∈ A" and yB: "y ∈ B"
      by auto
    show "k_times_Fr_differentiable_at k (λz. f z + g z) y"
      using Suc.IH A(3)[rule_format, OF yA] B(3)[rule_format, OF yB] by blast
  qed

  have diff: "(λy. f y + g y) differentiable (at x)"
    by (simp add: df dg)
  have Dsum:"∀v. k_times_Fr_differentiable_at k
           (λy. frechet_derivative (λz. f z + g z) (at y) v) x"
  proof
    fix v
    show "k_times_Fr_differentiable_at k
            (λy. frechet_derivative (λz. f z + g z) (at y) v) x"
    proof (cases k)
      case 0
      then show ?thesis
        by simp
    next
      case (Suc j)

      have ksum:
        "k_times_Fr_differentiable_at k
           (λy. frechet_derivative f (at y) v + frechet_derivative g (at y) v) x"
        using Suc.IH Df Dg by blast

      have eqD:
        "⋀y. y ∈ ?C ⟹
          frechet_derivative (λz. f z + g z) (at y) v =
          frechet_derivative f (at y) v + frechet_derivative g (at y) v"
      proof -
        fix y
        assume yC: "y ∈ ?C"
        then have yA: "y ∈ A" and yB: "y ∈ B"
          by auto

        have fy: "k_times_Fr_differentiable_at (Suc j) f y"
          using A(3) Suc yA by blast

        have gy: "k_times_Fr_differentiable_at (Suc j) g y"
          using B(3) Suc yB by blast


        have dfy: "f differentiable (at y)"
          using fy k_times_Fr_differentiable_at_mono[of 1 "Suc j" f y]
          by (simp add: one_times_Fr_iff)

        have dgy: "g differentiable (at y)"
          using gy k_times_Fr_differentiable_at_mono[of 1 "Suc j" g y]
          by (simp add: one_times_Fr_iff)

        have hder: "((λz. f z + g z) has_derivative
             (λh. frechet_derivative f (at y) h + frechet_derivative g (at y) h)) (at y)"
          by (simp add: dfy dgy frechet_derivative_works[THEN iffD1])

        then have "frechet_derivative (λz. f z + g z) (at y) =
              (λh. frechet_derivative f (at y) h + frechet_derivative g (at y) h)"
          using frechet_derivative_at[symmetric] by blast
        then show
          "frechet_derivative (λz. f z + g z) (at y) v =
           frechet_derivative f (at y) v + frechet_derivative g (at y) v"
          by simp
      qed

      have
        "k_times_Fr_differentiable_at k
           (λy. frechet_derivative (λz. f z + g z) (at y) v) x ⟷
         k_times_Fr_differentiable_at k
           (λy. frechet_derivative f (at y) v + frechet_derivative g (at y) v) x"
        by (smt (verit) C(1,2) eqD eq_on_open_k_times_Fr_differentiable_at)
      then show ?thesis
        using ksum by blast
    qed
  qed
  show ?case
    unfolding k_times_Fr_differentiable_at.simps(2)
    using C neigh diff Dsum by blast
qed

lemma k_times_Fr_scaleR:
  assumes "k_times_Fr_differentiable_at k f x"
  shows "k_times_Fr_differentiable_at k (λy. c *R f y) x"
  using assms
proof (induction k arbitrary: f x)
  case 0
  then show ?case
    by simp
next
  case (Suc k)
  from Suc.prems obtain A where
    A: "open A" "x ∈ A" "∀y∈A. k_times_Fr_differentiable_at k f y"
    and df: "f differentiable (at x)"
    and Df: "∀v. k_times_Fr_differentiable_at k (λy. frechet_derivative f (at y) v) x"
    unfolding k_times_Fr_differentiable_at.simps(2)
    by blast

  have neigh: "∀y∈A. k_times_Fr_differentiable_at k (λz. c *R f z) y"
    using A(3) Suc.IH by blast

  have diff: "(λy. c *R f y) differentiable (at x)"
    by (simp add: df)

  have Dscale: "∀v. k_times_Fr_differentiable_at k
           (λy. frechet_derivative (λz. c *R f z) (at y) v) x"
  proof
    fix v
    show "k_times_Fr_differentiable_at k
            (λy. frechet_derivative (λz. c *R f z) (at y) v) x"
    proof (cases k)
      case 0
      then show ?thesis
        by simp
    next
      case (Suc j)

      have kscaled:
        "k_times_Fr_differentiable_at k
           (λy. c *R frechet_derivative f (at y) v) x"
        using Suc.IH Df by blast

      have eqD:
        "⋀y. y ∈ A ⟹
          frechet_derivative (λz. c *R f z) (at y) v =
          c *R frechet_derivative f (at y) v"
      proof -
        fix y
        assume yA: "y ∈ A"

        have fy: "k_times_Fr_differentiable_at (Suc j) f y"
          using A(3) Suc yA by blast

        hence dfy: "f differentiable (at y)"
          using k_times_Fr_differentiable_at_mono[of 1 "Suc j" f y]
          by (simp add: one_times_Fr_iff)

        have hder: "((λz. c *R f z) has_derivative (λh. c *R frechet_derivative f (at y) h)) (at y)"
          by (simp add: dfy frechet_derivative_works[THEN iffD1] has_derivative_scaleR_right)

        have "frechet_derivative (λz. c *R f z) (at y) = (λh. c *R frechet_derivative f (at y) h)"
          by (metis frechet_derivative_at hder)
        then show "frechet_derivative (λz. c *R f z) (at y) v =  c *R frechet_derivative f (at y) v"
          by simp
      qed

      have  "k_times_Fr_differentiable_at k
           (λy. frechet_derivative (λz. c *R f z) (at y) v) x ⟷
         k_times_Fr_differentiable_at k
           (λy. c *R frechet_derivative f (at y) v) x"
        by (smt (verit) A(1,2) eqD eq_on_open_k_times_Fr_differentiable_at)
      then show ?thesis
        using kscaled by blast
    qed
  qed

  show ?case
    unfolding k_times_Fr_differentiable_at.simps(2)
    using A(1,2) neigh diff Dscale by blast
qed

lemma Ck_on_add:
  assumes "Ck_on k f U" and "Ck_on k g U"
  shows "Ck_on k (λy. f y + g y) U"
  using assms
proof (induction k arbitrary: f g)
  case 0
  then show ?case
    by (auto simp: Ck_on_0_iff intro: continuous_intros)
next
  case (Suc k)
  from Suc.prems have U: "open U"
    and df: "⋀x. x ∈ U ⟹ f differentiable (at x)"
    and dg: "⋀x. x ∈ U ⟹ g differentiable (at x)"
    and Df: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
    and Dg: "⋀v. Ck_on k (λx. frechet_derivative g (at x) v) U"
    by (auto simp: Ck_on_Suc_iff)
  show ?case
    unfolding Ck_on_Suc_iff
    using U df dg Suc.IH[OF Df Dg]
    by (auto simp: frechet_derivative_add cong: Ck_on_cong)
qed

text ‹A bounded bilinear operation preserves ‹Ck›; this covers products, scalar
  multiplication and inner products.›

lemma Ck_on_bilinear:
  assumes P: "bounded_bilinear P" and "Ck_on k f U" and "Ck_on k g U"
  shows "Ck_on k (λy. P (f y) (g y)) U"
  using assms(2,3)
proof (induction k arbitrary: f g)
  case 0
  then show ?case
    by (auto simp: Ck_on_0_iff intro: bounded_bilinear.continuous_on[OF P])
next
  case (Suc k)
  from Suc.prems have U: "open U"
    and df: "⋀x. x ∈ U ⟹ f differentiable (at x)"
    and dg: "⋀x. x ∈ U ⟹ g differentiable (at x)"
    and Df: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
    and Dg: "⋀v. Ck_on k (λx. frechet_derivative g (at x) v) U"
    by (auto simp: Ck_on_Suc_iff)
  have fk: "Ck_on k f U" and gk: "Ck_on k g U"
    using Suc.prems by (auto intro: Ck_on_SucD)
  have dP: "(λy. P (f y) (g y)) differentiable (at x)" if "x ∈ U" for x
    using bounded_bilinear.FDERIV[OF P frechet_derivative_works[THEN iffD1, OF df[OF that]]
        frechet_derivative_works[THEN iffD1, OF dg[OF that]]]
    unfolding differentiable_def by blast
  show ?case
    unfolding Ck_on_Suc_iff
    using U dP Ck_on_add[OF Suc.IH[OF fk Dg] Suc.IH[OF Df gk]]
    by (auto simp: frechet_derivative_bilinear[OF P] df dg cong: Ck_on_cong)
qed

lemma Ck_on_sum:
  fixes F :: "'i ⇒ 'a::real_normed_vector ⇒ 'b::real_normed_vector"
  assumes fin: "finite I"
      and ne: "I ≠ {}"
      and Ck: "⋀i. i ∈ I ⟹ Ck_on k (F i) U"
  shows "Ck_on k (λy. ∑i∈I. F i y) U"
  using fin ne Ck
proof (induction rule: finite_induct)
  case empty
  then show ?case by simp
next
  case (insert i I)
  have Ci: "Ck_on k (F i) U"
    using insert.prems by simp
  show ?case
  proof (cases "I = {}")
    case True
    then show ?thesis
      using Ci insert.hyps by simp
  next
    case False
    have CI: "Ck_on k (λy. ∑j∈I. F j y) U"
      using insert.IH[OF False] insert.prems by blast
    show ?thesis
      using Ck_on_add[OF Ci CI] insert.hyps by simp
  qed
qed

lemma Ck_on_scaleR:
  assumes "Ck_on k f U"
  shows "Ck_on k (λy. c *R f y) U"
  using Ck_on_bilinear[OF bounded_bilinear_scaleR Ck_on_const[OF Ck_on_open[OF assms]] assms] .

lemma Ck_on_id:
  "open U ⟹ Ck_on k (λx. x) U"
proof (induction k)
  case 0
  then show ?case by (simp add: Ck_on_0_iff)
next
  case (Suc k)
  then show ?case by (simp add: Ck_on_Suc_iff Ck_on_const)
qed

lemma Ck_on_neg:
  assumes "Ck_on k f U"
  shows "Ck_on k (λy. - f y) U"
proof -
  have "Ck_on k (λy. (-1) *R f y) U"
    by (rule Ck_on_scaleR[OF assms])
  thus ?thesis by simp
qed

lemma Ck_on_sub:
  assumes "Ck_on k f U" and "Ck_on k g U"
  shows "Ck_on k (λy. f y - g y) U"
proof -
  have "Ck_on k (λy. f y + (- g y)) U"
    by (rule Ck_on_add[OF assms(1) Ck_on_neg[OF assms(2)]])
  thus ?thesis by simp
qed

lemma Ck_on_mult:
  fixes f g :: "'a::real_normed_vector ⇒ 'b::real_normed_algebra"
  assumes "Ck_on k f U" and "Ck_on k g U"
  shows "Ck_on k (λy. f y * g y) U"
  using Ck_on_bilinear[OF bounded_bilinear_mult assms] .

lemma Ck_on_pow:
  fixes f :: "'a::real_normed_vector ⇒ real"
  assumes "Ck_on k f U"
  shows "Ck_on k (λy. (f y) ^ n) U"
proof (induction n)
  case 0
  have "open U" using assms by (simp add: Ck_on_def)
  then show ?case
    using Ck_on_const by simp
next
  case (Suc n)
  have "Ck_on k (λy. f y * (f y) ^ n) U"
    by (rule Ck_on_mult[OF assms Suc])
  thus ?case by (simp add: power_Suc2)
qed

lemma Ck_on_inverse:
  fixes f :: "'a::real_normed_vector ⇒ 'b::real_normed_field"
  assumes "Ck_on k f U" and "⋀y. y ∈ U ⟹ f y ≠ 0"
  shows "Ck_on k (λy. inverse (f y)) U"
  using assms
proof (induction k arbitrary: f)
  case 0
  then show ?case
    by (auto simp: Ck_on_0_iff intro!: continuous_on_inverse)
next
  case (Suc k)
  from Suc.prems(1) have U: "open U"
    and df: "⋀x. x ∈ U ⟹ f differentiable (at x)"
    and Df: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
    by (auto simp: Ck_on_Suc_iff)
  have inv: "Ck_on k (λy. inverse (f y)) U"
    using Ck_on_SucD[OF Suc.prems(1)] Suc.prems(2) by (rule Suc.IH)
  have D: "Ck_on k (λx. - (inverse (f x) * frechet_derivative f (at x) v * inverse (f x))) U"
    for v
    by (intro Ck_on_neg Ck_on_mult inv Df)
  show ?case
    unfolding Ck_on_Suc_iff
    using U df Suc.prems(2) D
    by (auto simp: frechet_derivative_inverse cong: Ck_on_cong)
qed

lemma Ck_on_divide:
  fixes f g :: "'a::real_normed_vector ⇒ real"
  assumes "Ck_on k f U" and "Ck_on k g U" and "⋀y. y ∈ U ⟹ g y ≠ 0"
  shows "Ck_on k (λy. f y / g y) U"
proof -
  have inv_g: "Ck_on k (λy. inverse (g y)) U"
    by (rule Ck_on_inverse[OF assms(2,3)])
  have "Ck_on k (λy. f y * inverse (g y)) U"
    by (rule Ck_on_mult[OF assms(1) inv_g])
  thus ?thesis by (simp add: divide_inverse)
qed

lemma Ck_on_inner:
  fixes f g :: "'a::real_normed_vector ⇒ 'b::real_inner"
  assumes "Ck_on k f U" and "Ck_on k g U"
  shows "Ck_on k (λy. f y ∙ g y) U"
  using Ck_on_bilinear[OF bounded_bilinear_inner assms] .

lemma Ck_on_norm_sq:
  fixes f :: "'a::real_normed_vector ⇒ 'b::real_inner"
  assumes "Ck_on k f U"
  shows "Ck_on k (λy. (norm (f y))2) U"
proof -
  have "Ck_on k (λy. f y ∙ f y) U"
    by (rule Ck_on_inner[OF assms assms])
  thus ?thesis by (simp add: dot_square_norm)
qed

lemma Ck_on_compose:
  fixes f :: "'a::real_normed_vector ⇒ 'b::euclidean_space"
    and g :: "'b ⇒ 'c::real_normed_vector"
  assumes "Ck_on k g V" and "Ck_on k f U" and fUV: "⋀y. y ∈ U ⟹ f y ∈ V"
  shows "Ck_on k (λy. g (f y)) U"
  using assms(1,2)
proof (induction k arbitrary: g)
  case 0
  then show ?case
    using fUV by (auto simp: Ck_on_0_iff intro: continuous_on_compose2)
next
  case (Suc k)
  from Suc.prems(1) have dg: "⋀y. y ∈ V ⟹ g differentiable (at y)"
    and Dg: "⋀w. Ck_on k (λy. frechet_derivative g (at y) w) V"
    by (auto simp: Ck_on_Suc_iff)
  from Suc.prems(2) have U: "open U"
    and df: "⋀x. x ∈ U ⟹ f differentiable (at x)"
    and Df: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
    by (auto simp: Ck_on_Suc_iff)
  have fk: "Ck_on k f U"
    using Suc.prems(2) by (rule Ck_on_SucD)
  have D: "Ck_on k (λx. ∑i∈Basis. (frechet_derivative f (at x) v ∙ i) *R
      frechet_derivative g (at (f x)) i) U" for v
  proof (rule Ck_on_sum[OF finite_Basis nonempty_Basis])
    fix i :: 'b
    have "Ck_on k (λx. frechet_derivative f (at x) v ∙ i) U"
      by (rule Ck_on_inner[OF Df Ck_on_const[OF U]])
    moreover have "Ck_on k (λx. frechet_derivative g (at (f x)) i) U"
      by (rule Suc.IH[OF Dg fk])
    ultimately show "Ck_on k (λx. (frechet_derivative f (at x) v ∙ i) *R
        frechet_derivative g (at (f x)) i) U"
      by (rule Ck_on_bilinear[OF bounded_bilinear_scaleR])
  qed
  show ?case
    unfolding Ck_on_Suc_iff
    using U df dg fUV D
    by (auto simp: frechet_derivative_compose_euclidean cong: Ck_on_cong
        intro: differentiable_chain_at[unfolded o_def])
qed

lemma Ck_on_Pair:
  assumes "Ck_on k f U" and "Ck_on k g U"
  shows "Ck_on k (λy. (f y, g y)) U"
  using assms
proof (induction k arbitrary: f g)
  case 0
  then show ?case
    by (auto simp: Ck_on_0_iff intro: continuous_intros)
next
  case (Suc k)
  from Suc.prems have U: "open U"
    and df: "⋀x. x ∈ U ⟹ f differentiable (at x)"
    and dg: "⋀x. x ∈ U ⟹ g differentiable (at x)"
    and Df: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
    and Dg: "⋀v. Ck_on k (λx. frechet_derivative g (at x) v) U"
    by (auto simp: Ck_on_Suc_iff)
  show ?case
    unfolding Ck_on_Suc_iff
    using U df dg Suc.IH[OF Df Dg]
    by (auto simp: frechet_derivative_Pair cong: Ck_on_cong)
qed

text ‹A bounded linear map is its own derivative, hence ‹Ck› for every ‹k›.›

lemma Ck_at_bounded_linear:
  fixes T :: "'b::real_normed_vector ⇒ 'c::real_normed_vector"
  assumes T: "bounded_linear T"
  shows "Ck_at k T x"
proof (induction k arbitrary: x)
  case 0
  have "continuous (at x) T"
    by (rule bounded_linear.continuous[OF T continuous_ident])
  thus ?case by simp
next
  case (Suc k)
  have der: "(T has_derivative T) (at y)" for y
    using bounded_linear.has_derivative[OF T has_derivative_ident] by simp
  show ?case
    unfolding Ck_at.simps(2)
  proof (intro conjI allI)
    show "∃A. open A ∧ x ∈ A ∧ (∀y∈A. Ck_at k T y)"
      by (intro exI[where x = UNIV] conjI) (auto simp: Suc.IH)
    show "T differentiable (at x)" using der unfolding differentiable_def by blast
    fix v
    have "frechet_derivative T (at y) v = T v" for y
      using frechet_derivative_at[OF der] by simp
    hence "(λy. frechet_derivative T (at y) v) = (λy. T v)" by (rule ext)
    thus "Ck_at k (λy. frechet_derivative T (at y) v) x"
      by (simp add: Ck_at_const)
  qed
qed

lemma Ck_on_bounded_linear:
  fixes T :: "'b::real_normed_vector ⇒ 'c::real_normed_vector"
  assumes "bounded_linear T" and "open V"
  shows "Ck_on k T V"
  using assms by (simp add: Ck_on_def Ck_at_bounded_linear)

lemma Ck_on_bounded_linear_compose:
  fixes T :: "'b::euclidean_space ⇒ 'c::real_normed_vector"
    and f :: "'a::real_normed_vector ⇒ 'b"
  assumes T: "bounded_linear T" and f: "Ck_on k f U"
  shows "Ck_on k (λy. T (f y)) U"
  by (rule Ck_on_compose[OF Ck_on_bounded_linear[OF T open_UNIV] f]) simp

text ‹For ‹f :: real ⇒ real›, const‹Ck_on› agrees with the one-dimensional notion
  const‹C_k_on›.›

lemma Ck_on_real_iff:
  fixes f :: "real ⇒ real"
  shows "Ck_on k f U ⟷ C_k_on k f U"
proof (induction k arbitrary: f)
  case 0
  show ?case
    by (simp add: Ck_on_0_iff C0_on_def)
next
  case (Suc k)
  show ?case
  proof
    assume "Ck_on (Suc k) f U"
    then have U: "open U" and d: "⋀x. x ∈ U ⟹ f differentiable (at x)"
      and D: "⋀v. Ck_on k (λx. frechet_derivative f (at x) v) U"
      by (auto simp: Ck_on_Suc_iff)
    have "Ck_on k (deriv f) U"
      using D[of 1] by (rule Ck_on_congI) (simp add: d frechet_derivative_one_eq_deriv)
    then have "C_k_on k (deriv f) U"
      by (simp add: Suc.IH)
    moreover have "f differentiable_on U"
      using U d by (simp add: differentiable_on_eq_differentiable_at)
    ultimately show "C_k_on (Suc k) f U"
      by (simp add: C_k_on_Suc_iff)
  next
    assume "C_k_on (Suc k) f U"
    then have diff: "f differentiable_on U" and "C_k_on k (deriv f) U"
      by (simp_all add: C_k_on_Suc_iff)
    then have Ck: "Ck_on k (deriv f) U"
      by (simp add: Suc.IH)
    have U: "open U"
      using Ck by (rule Ck_on_open)
    have d: "⋀x. x ∈ U ⟹ f differentiable (at x)"
      using diff U by (simp add: differentiable_on_eq_differentiable_at)
    have "Ck_on k (λx. frechet_derivative f (at x) v) U" for v
    proof (rule Ck_on_congI)
      show "Ck_on k (λx. v * deriv f x) U"
        using Ck_on_mult[OF Ck_on_const[OF U] Ck] .
      show "frechet_derivative f (at x) v = v * deriv f x" if "x ∈ U" for x
        using d[OF that] by (rule frechet_derivative_to_deriv)
    qed
    then show "Ck_on (Suc k) f U"
      using U d by (simp add: Ck_on_Suc_iff)
  qed
qed

text ‹For functions of one real variable, ‹C1› agrees with @{const C1_differentiable_on}.›

lemma Ck_on_1_iff_C1_differentiable_on:
  fixes f :: "real ⇒ 'b::real_normed_vector"
  assumes "open U"
  shows "Ck_on 1 f U ⟷ f C1_differentiable_on U"
proof -
  have fd: "frechet_derivative f (at x) v = v *R vector_derivative f (at x)"
    if "f differentiable (at x)" for x v
  proof -
    have "(f has_derivative (λh. h *R vector_derivative f (at x))) (at x)"
      using vector_derivative_works[THEN iffD1, OF that] by (simp only: has_vector_derivative_def)
    then have "frechet_derivative f (at x) = (λh. h *R vector_derivative f (at x))"
      by (rule frechet_derivative_at[symmetric])
    then show ?thesis
      by simp
  qed
  have "Ck_on 1 f U ⟷ (∀x∈U. f differentiable (at x)) ∧
      (∀v. continuous_on U (λx. frechet_derivative f (at x) v))"
    using assms by (simp add: One_nat_def Ck_on_Suc_iff Ck_on_0_iff)
  also have "… ⟷ (∀x∈U. f differentiable (at x)) ∧
      continuous_on U (λx. vector_derivative f (at x))"
  proof (rule conj_cong[OF refl])
    assume d: "∀x∈U. f differentiable (at x)"
    show "(∀v. continuous_on U (λx. frechet_derivative f (at x) v)) ⟷
        continuous_on U (λx. vector_derivative f (at x))"
    proof
      assume "∀v. continuous_on U (λx. frechet_derivative f (at x) v)"
      then have "continuous_on U (λx. frechet_derivative f (at x) 1)"
        by blast
      then show "continuous_on U (λx. vector_derivative f (at x))"
        by (rule continuous_on_eq) (simp add: d fd)
    next
      assume c: "continuous_on U (λx. vector_derivative f (at x))"
      show "∀v. continuous_on U (λx. frechet_derivative f (at x) v)"
      proof
        fix v
        have "continuous_on U (λx. v *R vector_derivative f (at x))"
          by (intro continuous_intros c)
        then show "continuous_on U (λx. frechet_derivative f (at x) v)"
          by (rule continuous_on_eq) (simp add: d fd)
      qed
    qed
  qed
  also have "… ⟷ f C1_differentiable_on U"
    by (simp add: C1_differentiable_on_eq)
  finally show ?thesis .
qed


subsection ‹Gradient for ‹realn ⇒ real››

definition grad_fun :: "(real^'n::finite ⇒ real) ⇒ real^'n ⇒ real^'n"
  ("∇")
  where "∇ f x = (THE g :: real^'n. GDERIV f x :> g)"

lemma grad_fun_eq:
  assumes "GDERIV f x :> g"
  shows "∇ f x = g"
  unfolding grad_fun_def using assms gradient_unique
  by (metis the_equality)

lemma grad_fun_satisfies_GDERIV:
  assumes "GDERIV f x :> g"
  shows "GDERIV f x :> ∇ f x"
  using assms grad_fun_eq by blast

lemma frechet_eq_inner_gradient:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "(f has_derivative L) (at x)" and "GDERIV f x :> ∇ f x"
  shows "L v = v ∙ ∇ f x"
  using assms has_derivative_unique gderiv_def by blast

subsection ‹Hessian for ‹realn ⇒ real››

text ‹
  The multi-dimensional Hessian: the Fréchet derivative of the gradient,
  represented as a matrix.
›

definition has_hessian ::
    "(real^'n::finite ⇒ real) ⇒ real^'n ⇒ real^'n^'n ⇒ bool"
    ("(HESS (_)/ (_)/ :> (_))" [1000, 1000, 60] 60)
  where "HESS f x :> H ⟷ (∇ f has_derivative (λv. H *v v)) (at x)"

lemma hessian_unique:
  "HESS f x :> H ⟹ HESS f x :> H' ⟹ H = H'"
  unfolding has_hessian_def
  by (metis has_derivative_unique matrix_eq)

definition hess_fun :: "(real^'n::finite ⇒ real) ⇒ real^'n ⇒ real^'n^'n"
  ("∇2")
  where "∇2 f x = (THE H :: real^'n^'n. HESS f x :> H)"

lemma hess_fun_eq:
  assumes "HESS f x :> H"
  shows "∇2 f x = H"
  unfolding hess_fun_def using assms hessian_unique
  by (metis the_equality)

lemma hessian_eq_jacobian_of_gradient:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "HESS f x :> H"
  shows "H = matrix (frechet_derivative (∇ f) (at x))"
  by (metis assms frechet_derivative_at[symmetric] has_hessian_def matrix_of_matrix_vector_mul)

text ‹
  The Hessian entries are iterated partial derivatives:
  ‹(∇2 f x) $ i $ j = ∂j (∂i f) (x)›.
›

lemma hessian_eq_double_nabla:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "HESS f x :> ∇2 f x"
  shows "∀i j. ∇2 f x $ i $ j = (∇ (λy. ∇ f y $ i)) x $ j"
proof (intro allI)
  fix i j
  have row_grad: "GDERIV (λy. ∇ f y $ i) x :> (∇2 f x) $ i"
  proof -
    have H: "(∇ f has_derivative (*v) (∇2 f x)) (at x)"
      using assms unfolding has_hessian_def by simp
    have Hcomp: "((λy. ∇ f y ∙ axis i 1) has_derivative
         (λv. ((*v) (∇2 f x)) v ∙ axis i 1)) (at x within UNIV)"
      using H by (subst (asm) has_derivative_componentwise_within[where S = UNIV],
                  auto simp: Basis_vec_def)
    have comp_fun:  "(λy. ∇ f y ∙ axis i 1) = (λy. ∇ f y $ i)"
      by (rule ext, simp add: cart_eq_inner_axis)
    have comp_deriv: "(λv. ((*v) (∇2 f x)) v ∙ axis i 1) = (λv. v ∙ ((∇2 f x) $ i))"
      by (rule ext, simp add: inner_axis' inner_commute matrix_vector_mul_component)
    from Hcomp show ?thesis
      unfolding gderiv_def by (simp add: comp_fun comp_deriv)
  qed
  hence "∇ (λy. ∇ f y $ i) x = (∇2 f x) $ i"
    by (rule grad_fun_eq)
  then show "∇2 f x $ i $ j = (∇ (λy. ∇ f y $ i)) x $ j"
    by simp
qed


subsection ‹Connecting ‹Ck› to the Hessian›

text ‹Consequences of ‹Ck› for gradients and Hessians.›

lemma Ck_2_imp_gradient_exists:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "Ck_on 2 f U" and "x ∈ U"
  shows "∃g. GDERIV f x :> g"
proof -
  from assms have "Ck_at 2 f x"
    by (simp add: Ck_on_def)
  then have "f differentiable (at x)"
    by (metis Ck_at.simps(2) Suc_1)
  then show ?thesis
    by (rule Fr_diff_imp_gradient_exists)
qed

lemma Ck_2_imp_hessian_exists:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "Ck_on 2 f U" and "x ∈ U"
  shows "HESS f x :> ∇2 f x"
proof -
  from assms have C2: "Ck_at 2 f x"
    by (simp add: Ck_on_def)

  then obtain A where
    A: "open A" "x ∈ A" "∀y∈A. Ck_at 1 f y"
    and diffx: "f differentiable (at x)"
    and D: "∀v. Ck_at 1 (λy. frechet_derivative f (at y) v) x"
    by (metis Ck_at.simps(2) Suc_1)

  let ?H = "(χ i. ∇ (λy. ∇ f y $ i) x)"

  have H_wit: "HESS f x :> ?H"
  proof (unfold has_hessian_def)
    have comp: "∀i∈Basis. ((λy. ∇ f y ∙ i) has_derivative (λv. ((*v) ?H) v ∙ i)) (at x)"
    proof clarify
      fix b :: "real^'n"
      assume b: "b ∈ Basis"
      then obtain i where i: "b = axis i 1"
        by (auto simp: Basis_vec_def)

      let ?Fi = "(λy. frechet_derivative f (at y) (axis i 1))"
      let ?Gi = "(λy. ∇ f y $ i)"

      have Fi_C1: "Ck_at 1 ?Fi x"
        using D by blast
      hence Fi_diff: "?Fi differentiable (at x)"
        by simp

      have eqA: "⋀y. y ∈ A ⟹ ?Fi y = ?Gi y"
        by (metis (lifting) A(3) Ck_at.simps(2) Fr_diff_imp_gradient_exists Suc_eq_plus1 add_0
            frechet_derivative_at grad_fun_eq gderiv_def inner_axis' inner_real_def lambda_one)


      have ev_eq: "eventually (λy. ?Fi y = ?Gi y) (nhds x)"
      proof -
        have "∃S. open S ∧ x ∈ S ∧ (∀y∈S. ?Fi y = ?Gi y)"
          using A eqA by blast
        then show ?thesis
          by (simp add: eventually_nhds)
      qed

      have Gi_diff: "?Gi differentiable (at x)"
        by (metis (no_types, lifting) A(1,2) Fi_diff differentiable_transform_within_open eqA)


      from Fr_diff_imp_gradient_exists[OF Gi_diff]
      obtain gi where gi: "GDERIV ?Gi x :> gi"
        by blast

      have gradGi: "GDERIV ?Gi x :> ∇ ?Gi x"
        using gi by (rule grad_fun_satisfies_GDERIV)

      have dGi: "(?Gi has_derivative (λv. v ∙ (?H $ i))) (at x)"
        using gradGi unfolding gderiv_def by simp

      have "((λy. ∇ f y ∙ b) has_derivative (λv. ((*v) ?H) v ∙ b)) (at x)"
        by (metis (no_types, lifting) ext cart_eq_inner_axis dGi i
            inner_commute matrix_vector_mul_component)
      then show "((λy. ∇ f y ∙ b) has_derivative (λv. ((*v) ?H) v ∙ b)) (at x)".
    qed

    then show "(∇ f has_derivative (*v) ?H) (at x)"
      using has_derivative_componentwise_within by blast
  qed
  show ?thesis
    using H_wit hess_fun_eq by fastforce
qed

lemma Ck_2_imp_hessian_continuous:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "Ck_on 2 f U"
  shows "continuous_on U (∇2 f)"
proof -
  have openU: "open U"
    using assms by (simp add: Ck_on_def)

  have comp_cont: "⋀x i j. x ∈ U ⟹ continuous (at x) (λy. ∇2 f y $ i $ j)"
  proof -
    fix x i j
    assume xU: "x ∈ U"

    from assms xU have C2x: "Ck_at 2 f x"
      by (simp add: Ck_on_def)

    from C2x obtain A where
      A: "open A" "x ∈ A" "∀y∈A. Ck_at 1 f y"
      and diffx: "f differentiable (at x)"
      and Dx: "∀v. Ck_at 1 (λy. frechet_derivative f (at y) v) x"
      by (metis Ck_at.simps(2) Suc_1)

    let ?Fi = "(λy. frechet_derivative f (at y) (axis i 1))"
    let ?K  = "(λy. frechet_derivative ?Fi (at y) (axis j 1))"
    let ?Hc = "(λy. ∇2 f y $ i $ j)"

    have Fi_C1: "Ck_at 1 ?Fi x"
      using Dx by simp

    have K_cont: "continuous (at x) ?K"
      using Fi_C1 by simp

    have eq_on_U: "⋀y. y ∈ U ⟹ frechet_derivative ?Fi (at y) (axis j 1) = ∇2 f y $ i $ j"
    proof -
      fix y
      assume yU: "y ∈ U"

      from assms yU have C2y: "Ck_at 2 f y"
        by (simp add: Ck_on_def)

      have dy: "f differentiable (at y)"
        by (metis C2y Ck_at.simps(2) Suc_1)

      have Fi_C1_y: "Ck_at 1 ?Fi y"
        using C2y by (metis Ck_at.simps(2) Suc_1)

      have Fi_diff_y: "?Fi differentiable (at y)"
        using Fi_C1_y by simp

      let ?Gi = "(λz. ∇ f z $ i)"

      have FG_eq_on_U: "⋀z. z ∈ U ⟹ ?Fi z = ?Gi z"
      proof -
        fix z
        assume zU: "z ∈ U"

        from assms zU have C2z: "Ck_at 2 f z"
          by (simp add: Ck_on_def)

        have dz: "f differentiable (at z)"
          by (metis C2z Ck_at.simps(2) Suc_1)

        from Fr_diff_imp_gradient_exists[OF dz]
        obtain g where g: "GDERIV f z :> g"
          by blast
        have g_eq: "∇ f z = g"
          using g by (rule grad_fun_eq)
        have "(f has_derivative (λv. v ∙ g)) (at z)"
          using g unfolding gderiv_def by simp
        hence fd_eq: "frechet_derivative f (at z) = (λv. v ∙ g)"
          by (metis frechet_derivative_at)
        show "?Fi z = ?Gi z"
          by (simp add: fd_eq g_eq inner_axis')
      qed

      have ev_FG: "eventually (λz. ?Fi z = ?Gi z) (nhds y)"
        using FG_eq_on_U eventually_nhds openU yU by blast


      have Gi_diff_y: "?Gi differentiable (at y)"
        by (metis (no_types, lifting) FG_eq_on_U Fi_diff_y differentiable_transform_within_open openU yU)


      then have fd_Fi_Gi: "frechet_derivative ?Fi (at y) = frechet_derivative ?Gi (at y)"
        by (smt (verit, best) FG_eq_on_U frechet_derivative_transform_within_open openU yU)

      from Fr_diff_imp_gradient_exists[OF Gi_diff_y]
      obtain gi where gi: "GDERIV ?Gi y :> gi"
        by blast

      have gi_eq: "∇ ?Gi y = gi"
        using gi by (rule grad_fun_eq)

      have "(?Gi has_derivative (λv. v ∙ gi)) (at y)"
        using gi unfolding gderiv_def by simp
      hence fd_Gi: "frechet_derivative ?Gi (at y) = (λv. v ∙ gi)"
        by (metis frechet_derivative_at)

      have fd_Gi_axis: "frechet_derivative ?Gi (at y) (axis j 1) = ∇ ?Gi y $ j"
        by (metis cart_eq_inner_axis fd_Gi gi_eq inner_commute)
      have Hess_y: "HESS f y :> ∇2 f y"
        using assms yU by (rule Ck_2_imp_hessian_exists)

      have hess_eq: "∇2 f y $ i $ j = ∇ ?Gi y $ j"
        using hessian_eq_double_nabla[OF Hess_y] by simp

      show "frechet_derivative ?Fi (at y) (axis j 1) = ∇2 f y $ i $ j"
        using fd_Fi_Gi fd_Gi_axis hess_eq by simp
    qed
    have ev_eq: "eventually (λy. ?K y = ?Hc y) (nhds x)"
      using eq_on_U eventually_nhds openU xU by blast
    show "continuous (at x) ?Hc"
      using K_cont ev_eq isCont_cong by fastforce
  qed

  show ?thesis
    unfolding continuous_on
  proof
    fix x
    assume xU: "x ∈ U"

    have isCont_H: "isCont (λy. χ i. χ j. ∇2 f y $ i $ j) x"
      unfolding isCont_def
    proof (rule tendsto_vec_lambda)
      fix i
      show "((λy. χ j. ∇2 f y $ i $ j) ⤏ (χ j. ∇2 f x $ i $ j)) (at x)"
      proof (rule tendsto_vec_lambda)
        fix j
        from comp_cont[OF xU, of i j]
        show "((λy. ∇2 f y $ i $ j) ⤏ ∇2 f x $ i $ j) (at x)"
          unfolding isCont_def by simp
      qed
    qed
    then have "continuous (at x) (∇2 f)"
      by simp
    then show "(∇2 f ⤏ ∇2 f x) (at x within U)"
      by (metis at_within_open continuous_within openU xU)
  qed
qed

text ‹The proof evaluates the second difference of ‹f› along ‹ei› and ‹ej› by the mean
  value theorem in both orders and lets the increments tend to ‹0›.›

lemma mixed_coordinate_second_derivative_eq:
  fixes f :: "real^'n::finite ⇒ real"
  assumes openU: "open U"
      and xU: "x ∈ U"
      and C2: "Ck_on 2 f U"
  shows "(∇ (λy. ∇ f y $ i)) x $ j = (∇ (λy. ∇ f y $ j)) x $ i"
proof -
  (* Notation *)
  let ?ei = "axis i 1 :: real^'n"
  let ?ej = "axis j 1 :: real^'n"

  (* Work inside a ball *)
  obtain r where r_pos: "r > 0" and rU: "ball x r ⊆ U"
    using openU xU by (meson open_contains_ball)

  define δ where "δ = r / 4"
  have δ_pos: "δ > 0" using r_pos by (simp add: δ_def)

  (* Any point x + s·ei + t·ej with |s|,|t| < δ lies in U. *)
  have inU: "⟦ ¦s¦ < δ; ¦t¦ < δ ⟧ ⟹ x + s *R ?ei + t *R ?ej ∈ U" for s t
  proof -
    assume s_bd: "¦s¦ < δ" and t_bd: "¦t¦ < δ"
    have "norm (s *R ?ei + t *R ?ej) ≤ ¦s¦ + ¦t¦"
      by (simp add: norm_triangle_le)
    also have "… < δ + δ" using s_bd t_bd by linarith
    also have "… = r / 2" by (simp add: δ_def)
    also have "… < r" using r_pos by linarith
    finally show "x + s *R ?ei + t *R ?ej ∈ U"
      by (metis (no_types, lifting) add.assoc basic_trans_rules(31)
          dist_0_norm dist_add_cancel group_cancel.rule0 mem_ball rU)
  qed

  (* Names for partial derivatives *)
  (* Ps(s,t) = ∂if at x + s·ei + t·ej *)
  define Ps where "Ps s t = ∇ f (x + s *R ?ei + t *R ?ej) $ i" for s t
  (* Qt(s,t) = ∂jf at x + s·ei + t·ej *)
  define Qt where "Qt s t = ∇ f (x + s *R ?ei + t *R ?ej) $ j" for s t

  (* Basic differentiability facts *)
  have diff_at: "f differentiable (at z)" if "z ∈ U" for z
    by (metis Ck_at.simps(2) Ck_on_def Suc_1 C2 that)

  have grad_exists: "GDERIV f z :> ∇ f z" if "z ∈ U" for z
    using Fr_diff_imp_gradient_exists[OF diff_at[OF that]]
      grad_fun_satisfies_GDERIV by blast

  have Hess_exists: "HESS f z :> ∇2 f z" if "z ∈ U" for z
    using C2 that by (rule Ck_2_imp_hessian_exists)

  have hcont: "continuous_on U (∇2 f)"
    using C2 by (rule Ck_2_imp_hessian_continuous)

  (* Row-gradient lemma *)
  (* GDERIV (λy. ∇ f y $ k) z :> (∇2f z) $ k  for z ∈ U *)
  have row_grad: "GDERIV (λy. ∇ f y $ k) z :> (∇2 f z) $ k"
    if zU: "z ∈ U" for z k
  proof -
    have H: "(∇ f has_derivative (*v) (∇2 f z)) (at z)"
      using Hess_exists[OF zU] unfolding has_hessian_def .
    have "((λy. ∇ f y ∙ axis k 1) has_derivative
         (λv. ((*v) (∇2 f z)) v ∙ axis k 1)) (at z within UNIV)"
      using H by (subst (asm) has_derivative_componentwise_within[where S = UNIV],
                  auto simp: Basis_vec_def)
    thus ?thesis
      unfolding gderiv_def
      by (simp add: inner_axis' inner_commute matrix_vector_mul_component)
  qed

  (* Derivatives of the slice maps *)
  (* ∂/∂s [Ps(s,t)] = (∇2f)$i$i  and  ∂/∂t [Ps(s,t)] = (∇2f)$i$j *)

  have Ps_has_deriv_t:
    "((λt'. Ps s t') has_real_derivative (∇2 f (x + s *R ?ei + t *R ?ej)) $ i $ j)
       (at t)"
    if s_bd: "¦s¦ < δ" and t_bd: "¦t¦ < δ" for s t
  proof -
    let ?z = "x + s *R ?ei + t *R ?ej"
    have zU: "?z ∈ U" using inU[OF s_bd t_bd] .
    have rg: "GDERIV (λy. ∇ f y $ i) ?z :> (∇2 f ?z) $ i"
      by (rule row_grad[OF zU])
    have fd: "((λy. ∇ f y $ i) has_derivative (λv. v ∙ ((∇2 f ?z) $ i))) (at ?z)"
      using rg unfolding gderiv_def .
    have lin: "((λt'. x + s *R ?ei + t' *R ?ej) has_derivative (λdt. dt *R ?ej)) (at t)"
      by (intro derivative_eq_intros) auto
    have chain:
      "((λt'. ∇ f (x + s *R ?ei + t' *R ?ej) $ i) has_derivative
         (λdt. (dt *R ?ej) ∙ ((∇2 f ?z) $ i))) (at t)"
      using has_derivative_compose[OF lin fd] by (simp add: o_def)
    have "(λdt. (dt *R ?ej) ∙ ((∇2 f ?z) $ i))
        = (λdt. dt * ((∇2 f ?z) $ i $ j))"
      by (rule ext, simp add: inner_axis' mult.commute)
    thus ?thesis
      using chain unfolding Ps_def has_field_derivative_def
      by (simp add: mult_commute_abs)
  qed

  have Qt_has_deriv_s:
    "((λs'. Qt s' t) has_real_derivative (∇2 f (x + s *R ?ei + t *R ?ej)) $ j $ i)
       (at s)"
    if s_bd: "¦s¦ < δ" and t_bd: "¦t¦ < δ" for s t
  proof -
    let ?z = "x + s *R ?ei + t *R ?ej"
    have zU: "?z ∈ U" using inU[OF s_bd t_bd] .
    have rg: "GDERIV (λy. ∇ f y $ j) ?z :> (∇2 f ?z) $ j"
      by (rule row_grad[OF zU])
    have fd: "((λy. ∇ f y $ j) has_derivative (λv. v ∙ ((∇2 f ?z) $ j))) (at ?z)"
      using rg unfolding gderiv_def .
    have lin: "((λs'. x + s' *R ?ei + t *R ?ej) has_derivative (λds. ds *R ?ei)) (at s)"
      by (intro derivative_eq_intros) auto
    have chain:
      "((λs'. ∇ f (x + s' *R ?ei + t *R ?ej) $ j) has_derivative
         (λds. (ds *R ?ei) ∙ ((∇2 f ?z) $ j))) (at s)"
      using has_derivative_compose[OF lin fd] by (simp add: o_def)
    have "(λds. (ds *R ?ei) ∙ ((∇2 f ?z) $ j))
        = (λds. ds * ((∇2 f ?z) $ j $ i))"
      by (rule ext, simp add: inner_axis' mult.commute)
    thus ?thesis
      using chain unfolding Qt_def has_field_derivative_def
      by (metis (no_types, lifting) ext mult.commute)
  qed

  (* Similarly for ∂/∂s [Φ(s,t)] and ∂/∂t [Φ(s,t)] *)
  have Phi_has_deriv_s: "((λs'. f (x + s' *R ?ei + t *R ?ej)) has_real_derivative Ps s t) (at s)"
    if s_bd: "¦s¦ < δ" and t_bd: "¦t¦ < δ" for s t
  proof -
    let ?z = "x + s *R ?ei + t *R ?ej"
    have zU: "?z ∈ U" using inU[OF s_bd t_bd] .
    have fd: "(f has_derivative (λv. v ∙ ∇ f ?z)) (at ?z)"
      using grad_exists[OF zU] unfolding gderiv_def .
    have lin: "((λs'. x + s' *R ?ei + t *R ?ej) has_derivative (λds. ds *R ?ei)) (at s)"
      by (intro derivative_eq_intros) auto
    have chain: "((λs'. f (x + s' *R ?ei + t *R ?ej)) has_derivative
         (λds. (ds *R ?ei) ∙ ∇ f ?z)) (at s)"
      using has_derivative_compose[OF lin fd] by (simp add: o_def)
    have "(λds. (ds *R ?ei) ∙ ∇ f ?z) = (λds. ds * (∇ f ?z $ i))"
      by (rule ext, simp add: inner_axis' mult.commute)
    thus ?thesis
      using chain unfolding Ps_def has_field_derivative_def
      by (metis (full_types, lifting) ext mult.commute)
  qed

  have Phi_has_deriv_t: "((λt'. f (x + s *R ?ei + t' *R ?ej)) has_real_derivative Qt s t) (at t)"
    if s_bd: "¦s¦ < δ" and t_bd: "¦t¦ < δ" for s t
  proof -
    let ?z = "x + s *R ?ei + t *R ?ej"
    have zU: "?z ∈ U" using inU[OF s_bd t_bd] .
    have fd: "(f has_derivative (λv. v ∙ ∇ f ?z)) (at ?z)"
      using grad_exists[OF zU] unfolding gderiv_def .
    have lin: "((λt'. x + s *R ?ei + t' *R ?ej) has_derivative (λdt. dt *R ?ej)) (at t)"
      by (intro derivative_eq_intros) auto
    have chain: "((λt'. f (x + s *R ?ei + t' *R ?ej)) has_derivative
         (λdt. (dt *R ?ej) ∙ ∇ f ?z)) (at t)"
      using has_derivative_compose[OF lin fd] by (simp add: o_def)
    have "(λdt. (dt *R ?ej) ∙ ∇ f ?z) = (λdt. dt * (∇ f ?z $ j))"
      by (rule ext, simp add: inner_axis' mult.commute)
    thus ?thesis
      using chain unfolding Qt_def has_field_derivative_def
      by (metis (full_types, lifting) ext mult.commute)
  qed


  (* Continuity of the relevant Hessian entries *)


  (* For the ε-δ argument we only need: *)
  have Hij_cont_at_0:
    "∀ε>0. ∃δ'>0. ∀s t. ¦s¦ < δ' ∧ ¦t¦ < δ' ⟶
       ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ i $ j - (∇2 f x) $ i $ j¦ < ε"
  proof -
    have cont_comp: "isCont (λp. (∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ i $ j) (0,0)"
    proof -
      have cont_hij: "continuous_on U (λz. (∇2 f z) $ i $ j)"
        using hcont by (simp add: continuous_on_component)
      have isCont_hij: "isCont (λz. (∇2 f z) $ i $ j) x"
        using cont_hij openU xU continuous_on_eq_continuous_at by blast
      have isCont_slice: "isCont (λp. x + fst p *R ?ei + snd p *R ?ej) (0::real, 0::real)"
        by (intro continuous_intros)
      have at_zero: "(λp. x + fst p *R ?ei + snd p *R ?ej) (0::real, 0::real) = x"
        by simp
      have "isCont (λz. (∇2 f z) $ i $ j)
              ((λp. x + fst p *R ?ei + snd p *R ?ej) (0, 0))"
        using isCont_hij by (simp add: at_zero)
      thus ?thesis
        by (rule isCont_o2[OF isCont_slice])
    qed
    show ?thesis
    proof (intro allI impI)
      fix ε :: real
      assume eps: "ε > 0"

      (* Step 1: unfold isCont to tendsto, then to eventually_at *)
      from cont_comp
      have "((λp. (∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ i $ j) ⤏
              (∇2 f x) $ i $ j) (at (0,0))"
        unfolding isCont_def by simp

      (* Step 2: instantiate tendsto_iff at ε *)
      from this[unfolded tendsto_iff] eps
      have "eventually (λp. dist ((∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ i $ j)
                                  ((∇2 f x) $ i $ j) < ε) (at (0,0))"
        by simp

      (* Step 3: unfold eventually_at to get r'' with the p ≠ (0,0) guard *)
      then obtain r'' where r''_pos: "r'' > 0"
        and r''_bd: "∀p. p ≠ (0::real, 0::real) ∧ dist p (0,0) < r'' ⟶
             dist ((∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ i $ j)
                  ((∇2 f x) $ i $ j) < ε"
        unfolding eventually_at by auto

      (* Step 4: extend to ALL p by case-splitting on p = (0,0) *)
      define δ' where "δ' = min δ (r'' / 2)"
      have "δ' > 0" using δ_pos r''_pos by (simp add: δ'_def)
      moreover have "∀s t. ¦s¦ < δ' ∧ ¦t¦ < δ' ⟶
        ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ i $ j -
         (∇2 f x) $ i $ j¦ < ε"
      proof (intro allI impI)
        fix s t assume st: "¦s¦ < δ' ∧ ¦t¦ < δ'"
        show "¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ i $ j -
               (∇2 f x) $ i $ j¦ < ε"
        proof (cases "s = 0 ∧ t = 0")
          case True
          then show ?thesis using eps by simp
        next
          case False
          then have "(s, t) ≠ (0::real, 0::real)" by auto
          moreover have "dist (s,t) (0::real, 0::real) < r''"
          proof -
            have "dist (s,t) (0::real, 0::real) ≤ ¦s¦ + ¦t¦"
              using sqrt_sum_squares_le_sum_abs by (simp add: dist_Pair_Pair)
            also have "… < r''" using st by (simp add: δ'_def)
            finally show ?thesis.
          qed
          ultimately have "dist ((∇2 f (x + s *R ?ei + t *R ?ej)) $ i $ j)
                                ((∇2 f x) $ i $ j) < ε"
            using r''_bd by auto
          thus ?thesis by (simp add: dist_real_def)
        qed
      qed
      ultimately show "∃δ'>0. ∀s t. ¦s¦ < δ' ∧ ¦t¦ < δ' ⟶
        ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ i $ j -
         (∇2 f x) $ i $ j¦ < ε"
        by blast
    qed
  qed

  have Hji_cont_at_0:
  "∀ε>0. ∃δ'>0. ∀s t. ¦s¦ < δ' ∧ ¦t¦ < δ' ⟶
     ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ j $ i - (∇2 f x) $ j $ i¦ < ε"
  proof -
    have cont_comp: "isCont (λp. (∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ j $ i) (0,0)"
    proof -
      have cont_hji: "continuous_on U (λz. (∇2 f z) $ j $ i)"
        using hcont by (simp add: continuous_on_component)
      have isCont_hji: "isCont (λz. (∇2 f z) $ j $ i) x"
        using cont_hji openU xU continuous_on_eq_continuous_at by blast
      have isCont_slice: "isCont (λp. x + fst p *R ?ei + snd p *R ?ej) (0::real, 0::real)"
        by (intro continuous_intros)
      have at_zero: "(λp. x + fst p *R ?ei + snd p *R ?ej) (0::real, 0::real) = x"
        by simp
      have "isCont (λz. (∇2 f z) $ j $ i)
              ((λp. x + fst p *R ?ei + snd p *R ?ej) (0, 0))"
        using isCont_hji by (simp add: at_zero)
      thus ?thesis
        by (rule isCont_o2[OF isCont_slice])
    qed
    show ?thesis
    proof (intro allI impI)
      fix ε :: real
      assume eps: "ε > 0"

      from cont_comp
      have "((λp. (∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ j $ i) ⤏
              (∇2 f x) $ j $ i) (at (0,0))"
        unfolding isCont_def by simp

      from this[unfolded tendsto_iff] eps
      have "eventually (λp. dist ((∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ j $ i)
                                  ((∇2 f x) $ j $ i) < ε) (at (0,0))"
        by simp

      then obtain r'' where r''_pos: "r'' > 0"
        and r''_bd: "∀p. p ≠ (0::real, 0::real) ∧ dist p (0,0) < r'' ⟶
             dist ((∇2 f (x + fst p *R ?ei + snd p *R ?ej)) $ j $ i)
                  ((∇2 f x) $ j $ i) < ε"
        unfolding eventually_at by auto

      define δ' where "δ' = min δ (r'' / 2)"
      have "δ' > 0"
        using δ_pos r''_pos by (simp add: δ'_def)
      moreover have "∀s t. ¦s¦ < δ' ∧ ¦t¦ < δ' ⟶
        ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ j $ i -
         (∇2 f x) $ j $ i¦ < ε"
      proof (intro allI impI)
        fix s t
        assume st: "¦s¦ < δ' ∧ ¦t¦ < δ'"
        show "¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ j $ i -
               (∇2 f x) $ j $ i¦ < ε"
        proof (cases "s = 0 ∧ t = 0")
          case True
          then show ?thesis
            using eps by simp
        next
          case False
          then have "(s, t) ≠ (0::real, 0::real)"
            by auto
          moreover have "dist (s,t) (0::real, 0::real) < r''"
          proof -
            have "dist (s,t) (0::real, 0::real) ≤ ¦s¦ + ¦t¦"
              using sqrt_sum_squares_le_sum_abs by (simp add: dist_Pair_Pair)
            also have "… < r''"
              using st by (simp add: δ'_def)
            finally show ?thesis .
          qed
          ultimately have "dist ((∇2 f (x + s *R ?ei + t *R ?ej)) $ j $ i)
                                ((∇2 f x) $ j $ i) < ε"
            using r''_bd by auto
          thus ?thesis
            by (simp add: dist_real_def)
        qed
      qed
      ultimately show "∃δ'>0. ∀s t. ¦s¦ < δ' ∧ ¦t¦ < δ' ⟶
        ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ j $ i -
         (∇2 f x) $ j $ i¦ < ε"
        by blast
    qed
  qed


  (* The rectangle increment *)
  define Δ where
    "Δ h k = f (x + h *R ?ei + k *R ?ej)
           - f (x + h *R ?ei)
           - f (x + k *R ?ej)
           + f x" for h k

  (* MVT, direction 1: differentiate in s first, then t *)
  have dir1:
    "∃ξ η. ¦ξ¦ ≤ ¦h¦ ∧ ¦η¦ ≤ ¦k¦ ∧
            Δ h k = h * k * (∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j"
    if h_pos: "h > 0" and k_pos: "k > 0"
       and h_bd: "h < δ" and k_bd: "k < δ"
    for h k
  proof -
    (* g(s) = f(x + s·ei + k·ej) - f(x + s·ei) *)
    define g where "g s = f (x + s *R ?ei + k *R ?ej) - f (x + s *R ?ei)" for s

    have g_deriv: "(g has_real_derivative (Ps s k - Ps s 0)) (at s)"
      if "¦s¦ < δ" for s
    proof -
      have "((λs'. f (x + s' *R ?ei + k *R ?ej)) has_real_derivative Ps s k) (at s)"
        using Phi_has_deriv_s[of s k] that k_bd k_pos by linarith
      moreover have "((λs'. f (x + s' *R ?ei + 0 *R ?ej)) has_real_derivative Ps s 0) (at s)"
        using Phi_has_deriv_s[of s 0] that δ_pos by auto
      ultimately show ?thesis
        unfolding g_def
        by (subst derivative_eq_intros, simp_all)
    qed

    (* Apply MVT to g on [0,h] *)
    have g_deriv_on_seg: "⋀x. 0 ≤ x ⟹ x ≤ h ⟹ (g has_real_derivative (Ps x k - Ps x 0)) (at x)"
    proof -
      fix x :: real
      assume x0: "0 ≤ x"
      assume xh: "x ≤ h"
      have "¦x¦ = x"
        using x0 by simp
      also have "... ≤ h"
        using xh by simp
      also have "... < δ"
        using h_bd by simp
      finally have "¦x¦ < δ" .
      thus "(g has_real_derivative (Ps x k - Ps x 0)) (at x)"
        by (rule g_deriv)
    qed

    have g_diff: "∃ξ. 0 < ξ ∧ ξ < h ∧ Δ h k = h * (Ps ξ k - Ps ξ 0)"
    proof -
      obtain ξ where ξ:
        "0 < ξ" "ξ < h"
        "g h - g 0 = (h - 0) * (Ps ξ k - Ps ξ 0)"
        using MVT2[of 0 h g "λx. Ps x k - Ps x 0"]
          h_pos g_deriv_on_seg
        by blast
      have "g h - g 0 = Δ h k"
        by (simp add: g_def Δ_def)
      with ξ show ?thesis
        by auto
    qed
    then obtain ξ where ξ_pos: "0 < ξ" and ξ_lt: "ξ < h"
      and eq1: "Δ h k = h * (Ps ξ k - Ps ξ 0)" by blast

    (* Now apply MVT to p(t) = Ps(ξ,t) on [0,k] *)
    define p where "p t = Ps ξ t" for t

    have p_deriv: "(p has_real_derivative (∇2 f (x + ξ *R ?ei + t *R ?ej)) $ i $ j) (at t)"
      if "¦t¦ < δ" for t
      using Ps_has_deriv_t[of ξ t] ξ_lt h_bd that
      unfolding p_def
      using ξ_pos by argo


    (* MVT application to p on [0,k] *)
    have p_deriv_on_seg: "⋀t. 0 ≤ t ⟹ t ≤ k ⟹
       (p has_real_derivative (∇2 f (x + ξ *R ?ei + t *R ?ej)) $ i $ j) (at t)"
    proof -
      fix t :: real
      assume t0: "0 ≤ t"
      assume tk: "t ≤ k"
      have "¦t¦ = t"
        using t0 by simp
      also have "... ≤ k"
        using tk by simp
      also have "... < δ"
        using k_bd by simp
      finally have "¦t¦ < δ".
      thus "(p has_real_derivative (∇2 f (x + ξ *R ?ei + t *R ?ej)) $ i $ j) (at t)"
        by (rule p_deriv)
    qed

    have p_diff: "∃η. 0 < η ∧ η < k ∧
        Ps ξ k - Ps ξ 0 = k * (∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j"
    proof -
      obtain η where η:
        "0 < η"
        "η < k"
        "p k - p 0 = (k - 0) * ((∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j)"
        using MVT2[of 0 k p "λt. (∇2 f (x + ξ *R ?ei + t *R ?ej)) $ i $ j"] k_pos p_deriv_on_seg
        by blast
      have "p k - p 0 = Ps ξ k - Ps ξ 0"
        unfolding p_def by simp
      with η show ?thesis
        by auto
    qed
    then obtain η where η_pos: "0 < η" and η_lt: "η < k"
      and eq2: "Ps ξ k - Ps ξ 0 = k * (∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j"
      by blast

    have "Δ h k = h * (k * (∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j)"
      using eq1 eq2 by simp
    hence "Δ h k = h * k * (∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j"
      by (simp add: mult.assoc)
    moreover have "¦ξ¦ ≤ ¦h¦" using ξ_pos ξ_lt h_pos by linarith
    moreover have "¦η¦ ≤ ¦k¦" using η_pos η_lt k_pos by linarith
    ultimately show ?thesis by blast
  qed

  (* MVT, direction 2: differentiate in t first, then s *)
  have dir2:
    "∃ξ' η'. ¦ξ'¦ ≤ ¦h¦ ∧ ¦η'¦ ≤ ¦k¦ ∧
              Δ h k = h * k * (∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i"
    if h_pos: "h > 0" and k_pos: "k > 0"
       and h_bd: "h < δ" and k_bd: "k < δ"
    for h k
  proof -
    (* g̃(t) = f(x + h·ei + t·ej) - f(x + t·ej) *)
    define g' where "g' t = f (x + h *R ?ei + t *R ?ej) - f (x + t *R ?ej)" for t

    have g'_deriv: "(g' has_real_derivative (Qt h t - Qt 0 t)) (at t)"
      if "¦t¦ < δ" for t
    proof -
      have "((λt'. f (x + h *R ?ei + t' *R ?ej)) has_real_derivative Qt h t) (at t)"
        using Phi_has_deriv_t[of h t] h_bd that
        using h_pos by linarith
      moreover have "((λt'. f (x + 0 *R ?ei + t' *R ?ej)) has_real_derivative Qt 0 t) (at t)"
        using Phi_has_deriv_t[of 0 t] δ_pos that by auto
      ultimately show ?thesis
        unfolding g'_def by (subst derivative_eq_intros, simp_all)
    qed

    (* MVT on g̃ over [0,k] *)
    have g'_diff: "∃η'. 0 < η' ∧ η' < k ∧ Δ h k = k * (Qt h η' - Qt 0 η')"
    proof -
      have "g' k - g' 0 = Δ h k"
        by (simp add: g'_def Δ_def)
      moreover have g'_deriv_on_seg:
        "⋀t. 0 ≤ t ⟹ t ≤ k ⟹ (g' has_real_derivative (Qt h t - Qt 0 t)) (at t)"
      proof -
        fix t :: real
        assume t0: "0 ≤ t"
        assume tk: "t ≤ k"
        have "¦t¦ = t"
          using t0 by simp
        also have "... ≤ k"
          using tk by simp
        also have "... < δ"
          using k_bd by simp
        finally have "¦t¦ < δ".
        thus "(g' has_real_derivative (Qt h t - Qt 0 t)) (at t)"
          by (rule g'_deriv)
      qed
      moreover obtain η' where "0 < η'" "η' < k"
        and "g' k - g' 0 = k * (Qt h η' - Qt 0 η')"
        using MVT2[of 0 k g' "λt. Qt h t - Qt 0 t"] k_pos g'_deriv_on_seg
        by auto
      ultimately show ?thesis
        by auto
    qed

    then obtain η' where η'_pos: "0 < η'" and η'_lt: "η' < k"
      and eq1': "Δ h k = k * (Qt h η' - Qt 0 η')" by blast

    (* MVT on q(s) = Qt(s, η') over [0,h] *)
    define q where "q s = Qt s η'" for s

    have q_deriv_on_seg:
      "⋀s. 0 ≤ s ⟹ s ≤ h ⟹
        (q has_real_derivative (∇2 f (x + s *R ?ei + η' *R ?ej)) $ j $ i) (at s)"
    proof -
      fix s :: real
      assume s0: "0 ≤ s"
      assume sh: "s ≤ h"
      have "¦s¦ = s"
        using s0 by simp
      also have "... ≤ h"
        using sh by simp
      also have "... < δ"
        using h_bd by simp
      finally have "¦s¦ < δ" .
      thus "(q has_real_derivative (∇2 f (x + s *R ?ei + η' *R ?ej)) $ j $ i) (at s)"
        using Qt_has_deriv_s η'_lt η'_pos ‹q ≡ λs. Qt s η'› k_bd by fastforce
    qed

    have q_diff: "∃ξ'. 0 < ξ' ∧ ξ' < h ∧
        Qt h η' - Qt 0 η' = h * (∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i"
    proof -
      obtain ξ' where ξ':
        "0 < ξ'"
        "ξ' < h"
        "q h - q 0 = (h - 0) * ((∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i)"
        using MVT2[of 0 h q "λs. (∇2 f (x + s *R ?ei + η' *R ?ej)) $ j $ i"] h_pos q_deriv_on_seg
        by blast
      have "q h - q 0 = Qt h η' - Qt 0 η'"
        unfolding q_def by simp
      with ξ' show ?thesis
        by auto
    qed
    then obtain ξ' where ξ'_pos: "0 < ξ'" and ξ'_lt: "ξ' < h"
      and eq2': "Qt h η' - Qt 0 η' = h * (∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i"
      by blast

    have "Δ h k = k * (h * (∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i)"
      using eq1' eq2' by simp
    hence "Δ h k = h * k * (∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i"
      by (simp add: mult.commute mult.assoc)
    moreover have "¦ξ'¦ ≤ ¦h¦" using ξ'_pos ξ'_lt h_pos by linarith
    moreover have "¦η'¦ ≤ ¦k¦" using η'_pos η'_lt k_pos by linarith
    ultimately show ?thesis by blast
  qed

  (* Combine: equality of Hessian entries *)
  have "(∇2 f x) $ i $ j = (∇2 f x) $ j $ i"
  proof (rule ccontr)
    assume neq: "(∇2 f x) $ i $ j ≠ (∇2 f x) $ j $ i"

    define ε where "ε = ¦(∇2 f x) $ i $ j - (∇2 f x) $ j $ i¦ / 3"
    then have ε_pos: "ε > 0" using neq by simp

    (* By continuity, get δ1 for the (i,j) entry and δ2 for the (j,i) entry *)
    obtain δ1 where δ1_pos: "δ1 > 0"
      and δ1_bd: "∀s t. ¦s¦ < δ1 ∧ ¦t¦ < δ1 ⟶
        ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ i $ j - (∇2 f x) $ i $ j¦ < ε"
      using Hij_cont_at_0 ε_pos by blast

    obtain δ2 where δ2_pos: "δ2 > 0"
      and δ2_bd: "∀s t. ¦s¦ < δ2 ∧ ¦t¦ < δ2 ⟶
        ¦(∇2 f (x + s *R ?ei + t *R ?ej)) $ j $ i - (∇2 f x) $ j $ i¦ < ε"
      using Hji_cont_at_0 ε_pos by blast

    define δ3 where "δ3 = min δ (min δ1 δ2)"
    have δ3_pos: "δ3 > 0" using δ_pos δ1_pos δ2_pos by (simp add: δ3_def)

    (* Pick concrete h, k *)
    define h where "h = δ3 / 2"
    define k where "k = δ3 / 2"
    have h_pos: "h > 0" and k_pos: "k > 0"
      using δ3_pos by (simp_all add: h_def k_def)
    have h_bd: "h < δ" and k_bd: "k < δ"
      using δ3_pos by (simp_all add: h_def k_def δ3_def, auto)

    (* Apply dir1 and dir2 *)
    obtain ξ η where ξ_bd: "¦ξ¦ ≤ h" and η_bd: "¦η¦ ≤ k"
      and eq_ij: "Δ h k = h * k * (∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j"
      using dir1[OF h_pos k_pos h_bd k_bd] h_pos k_pos by auto

    obtain ξ' η' where ξ'_bd: "¦ξ'¦ ≤ h" and η'_bd: "¦η'¦ ≤ k"
      and eq_ji: "Δ h k = h * k * (∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i"
      using dir2[OF h_pos k_pos h_bd k_bd] h_pos k_pos by auto

    (* Both ξ,η and ξ',η' are within δ1 and δ2 bounds *)
    have "¦ξ¦ < δ1" and "¦η¦ < δ1"
      using ξ_bd η_bd δ3_pos by (simp_all add: h_def k_def δ3_def)
    hence close_ij:
      "¦(∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j - (∇2 f x) $ i $ j¦ < ε"
      using δ1_bd by blast

    have "¦ξ'¦ < δ2" and "¦η'¦ < δ2"
      using ξ'_bd η'_bd δ3_pos by (simp_all add: h_def k_def δ3_def)
    hence close_ji:
      "¦(∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i - (∇2 f x) $ j $ i¦ < ε"
      using δ2_bd by blast

    (* From eq_ij and eq_ji, since h*k > 0 we can cancel: *)
    have "(∇2 f (x + ξ *R ?ei + η *R ?ej)) $ i $ j =
          (∇2 f (x + ξ' *R ?ei + η' *R ?ej)) $ j $ i"
      using eq_ij eq_ji h_pos k_pos by simp

    (* Triangle inequality gives contradiction *)
    hence "¦(∇2 f x) $ i $ j - (∇2 f x) $ j $ i¦ < 2 * ε"
      using close_ij close_ji by linarith
    hence "¦(∇2 f x) $ i $ j - (∇2 f x) $ j $ i¦
            < 2 * ¦(∇2 f x) $ i $ j - (∇2 f x) $ j $ i¦ / 3"
      by (simp add: ε_def)
    moreover have "¦(∇2 f x) $ i $ j - (∇2 f x) $ j $ i¦ > 0"
      using neq by simp
    ultimately show False
      by (simp add: field_simps)
  qed

  (* Translate to gradient notation *)
  have Hx: "HESS f x :> ∇2 f x"
    using Hess_exists xU by blast

  have rowi: "(∇2 f x) $ i $ j = (∇ (λy. ∇ f y $ i)) x $ j"
    using hessian_eq_double_nabla[OF Hx] by simp
  have rowj: "(∇2 f x) $ j $ i = (∇ (λy. ∇ f y $ j)) x $ i"
    using hessian_eq_double_nabla[OF Hx] by simp

  show ?thesis
    using ‹(∇2 f x) $ i $ j = (∇2 f x) $ j $ i› rowi rowj by simp
qed



theorem clairaut_hessian_symmetric:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "open U"
      and "x ∈ U"
      and "Ck_on 2 f U"
  shows "transpose (∇2 f x) = ∇2 f x"
proof -
  have H: "HESS f x :> ∇2 f x"
    using assms(2,3) by (subst Ck_2_imp_hessian_exists, simp_all)

  have sym_entries: "∀i j. ∇2 f x $ i $ j = ∇2 f x $ j $ i"
  proof (intro allI)
    fix i j
    have ij: "∇2 f x $ i $ j = (∇ (λy. ∇ f y $ i)) x $ j"
      using hessian_eq_double_nabla[OF H] by simp
    have ji: "∇2 f x $ j $ i = (∇ (λy. ∇ f y $ j)) x $ i"
      using hessian_eq_double_nabla[OF H] by simp
    have mix: "(∇ (λy. ∇ f y $ i)) x $ j = (∇ (λy. ∇ f y $ j)) x $ i"
      by (rule mixed_coordinate_second_derivative_eq[OF assms])
    show "∇2 f x $ i $ j = ∇2 f x $ j $ i"
      using ij ji mix by simp
  qed
  then show ?thesis
    by (simp add: Finite_Cartesian_Product.transpose_def)
qed

text ‹Equivalently, all mixed partials commute.›

corollary mixed_partials_commute:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "open U" and "x ∈ U" and "Ck_on 2 f U"
  shows "∇2 f x $ i $ j = ∇2 f x $ j $ i"
  using clairaut_hessian_symmetric[OF assms]
  by (metis (no_types, lifting) Finite_Cartesian_Product.transpose_def vec_lambda_beta)

subsection ‹Basic algebra of gradients›

lemma grad_fun_add:
  fixes f g :: "real^'n::finite ⇒ real"
  assumes "∃gf. GDERIV f x :> gf"
      and "∃gg. GDERIV g x :> gg"
  shows "∇ (λy. f y + g y) x = ∇ f x + ∇ g x"
proof -
  have Gf: "GDERIV f x :> ∇ f x"
    using assms(1) by (blast intro: grad_fun_satisfies_GDERIV)
  have Gg: "GDERIV g x :> ∇ g x"
    using assms(2) by (blast intro: grad_fun_satisfies_GDERIV)
  have "GDERIV (λy. f y + g y) x :> ∇ f x + ∇ g x"
    by (rule GDERIV_add[OF Gf Gg])
  thus ?thesis
    by (rule grad_fun_eq)
qed

lemma grad_fun_scaleR:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "∃gf. GDERIV f x :> gf"
  shows "∇ (λy. c * f y) x = c *R ∇ f x"
proof -
  have Gf: "GDERIV f x :> ∇ f x"
    using assms by (blast intro: grad_fun_satisfies_GDERIV)
  have "GDERIV (λy. c * f y) x :> c *R ∇ f x"
    by (rule GDERIV_cmult[OF Gf])
  thus ?thesis
    by (rule grad_fun_eq)
qed

lemma grad_fun_neg:
  fixes f :: "real^'n::finite ⇒ real"
  assumes "∃gf. GDERIV f x :> gf"
  shows "∇ (λy. - f y) x = - ∇ f x"
proof -
  have "∇ (λy. (-1) * f y) x = (-1) *R ∇ f x"
    by (rule grad_fun_scaleR[OF assms])
  thus ?thesis by simp
qed

lemma grad_fun_sub:
  fixes f g :: "real^'n::finite ⇒ real"
  assumes "∃gf. GDERIV f x :> gf"
      and "∃gg. GDERIV g x :> gg"
  shows "∇ (λy. f y - g y) x = ∇ f x - ∇ g x"
proof -
  have "∇ (λy. f y + (- g y)) x = ∇ f x + ∇ (λy. - g y) x"
    using assms grad_fun_add GDERIV_minus by blast
  also have "∇ (λy. - g y) x = - ∇ g x"
    by (rule grad_fun_neg[OF assms(2)])
  finally show ?thesis by simp
qed


subsection ‹Constants and affine maps›

lemma grad_fun_const:
  fixes c :: real
  shows "∇ (λ_. c) x = 0"
  by (rule grad_fun_eq[OF GDERIV_const])

lemma grad_fun_affine:
  fixes a :: real and b :: "real^'n::finite"
  shows "∇ (λx. a + x ∙ b) x = b"
  by (rule grad_fun_eq[OF GDERIV_affine])

lemma grad_fun_sum:
  fixes F :: "'i ⇒ real^'n::finite ⇒ real"
  assumes exG: "⋀i. i ∈ I ⟹ ∃g. GDERIV (F i) x :> g"
  shows "∇ (λy. ∑i∈I. F i y) x = (∑i∈I. ∇ (F i) x)"
proof -
  have G: "⋀i. i ∈ I ⟹ GDERIV (F i) x :> ∇ (F i) x"
    using exG by (blast intro: grad_fun_satisfies_GDERIV)
  have "GDERIV (λy. ∑i∈I. F i y) x :> (∑i∈I. ∇ (F i) x)"
    by (rule GDERIV_sum[OF G])
  thus ?thesis
    by (rule grad_fun_eq)
qed


subsection ‹Hessian: constants and affine maps›

lemma HESS_const_zero:
  fixes c :: real
  shows "HESS (λ_. c) x :> 0"
  unfolding has_hessian_def
  by (metis (no_types, lifting) ext grad_fun_const has_derivative_const matrix_vector_mult_0)

lemma HESS_affine_zero:
  fixes a :: real and b :: "real^'n::finite"
  shows "HESS (λx. a + x ∙ b) x :> 0"
  unfolding has_hessian_def
  by (metis (no_types, lifting) ext grad_fun_affine has_derivative_const matrix_vector_mult_0)

lemma hessian_const_zero:
  fixes c :: real
  shows "∇2 (λ_. c) x = 0"
  using HESS_const_zero by (metis hess_fun_eq)

lemma hessian_affine_zero:
  fixes a :: real and b :: "real^'n::finite"
  shows "∇2 (λx. a + x ∙ b) x = 0"
  using HESS_affine_zero by (metis hess_fun_eq)


subsection ‹Coordinate formulas›

lemma HESS_row_gradient:
  fixes f :: "real^'n::finite ⇒ real"
  assumes H: "HESS f x :> Hx"
  shows "GDERIV (λy. ∇ f y $ i) x :> Hx $ i"
proof -
  have Hd: "(∇ f has_derivative (*v) Hx) (at x)"
    using H unfolding has_hessian_def by simp
  have Hcomp:
    "((λy. ∇ f y ∙ axis i 1) has_derivative
       (λv. ((*v) Hx) v ∙ axis i 1)) (at x within UNIV)"
    using Hd
    by (subst (asm) has_derivative_componentwise_within[where S = UNIV])
       (auto simp: Basis_vec_def)
  have comp_fun:
    "(λy. ∇ f y ∙ axis i 1) = (λy. ∇ f y $ i)"
    by (rule ext) (simp add: cart_eq_inner_axis)
  have comp_deriv:
    "(λv. ((*v) Hx) v ∙ axis i 1) = (λv. v ∙ (Hx $ i))"
    by (metis (no_types) cart_eq_inner_axis inner_commute matrix_vector_mul_component)
  show ?thesis
    using Hcomp
    unfolding gderiv_def
    by (simp add: comp_fun comp_deriv)
qed

lemma HESS_row_eq:
  fixes f :: "real^'n::finite ⇒ real"
  assumes H: "HESS f x :> Hx"
  shows "∇ (λy. ∇ f y $ i) x = Hx $ i"
  by (rule grad_fun_eq[OF HESS_row_gradient[OF H]])

lemma HESS_component_eq:
  fixes f :: "real^'n::finite ⇒ real"
  assumes H: "HESS f x :> Hx"
  shows "Hx $ i $ j = (∇ (λy. ∇ f y $ i)) x $ j"
  using HESS_row_eq[OF H, of i] by simp


subsection ‹Hessian algebra at the predicate level›

lemma HESS_add:
  fixes f g :: "real^'n::finite ⇒ real"
  assumes Hf: "HESS f x :> Hf'" and Hg: "HESS g x :> Hg'"
      and eq: "⋀y. y ∈ A ⟹ ∇ (λz. f z + g z) y = ∇ f y + ∇ g y"
      and Aop: "open A" and xA: "x ∈ A"
  shows "HESS (λy. f y + g y) x :> Hf' + Hg'"
proof -
  have dsum: "((λy. ∇ f y + ∇ g y) has_derivative
               (λv. Hf' *v v + Hg' *v v)) (at x)"
    using has_derivative_add
      Hf[unfolded has_hessian_def] Hg[unfolded has_hessian_def] by blast
  have dtrans: "((λy. ∇ (λz. f z + g z) y) has_derivative
                 (λv. Hf' *v v + Hg' *v v)) (at x)"
    by (smt (verit, best) Aop dsum eq has_derivative_transfer_on_open xA)

  have "⋀v. (Hf' + Hg') *v v = Hf' *v v + Hg' *v v"
    by (simp add: matrix_vector_mult_def vec_eq_iff sum.distrib distrib_right)
  thus ?thesis
    unfolding has_hessian_def using dtrans by presburger
qed

lemma HESS_scaleR:
  fixes f :: "real^'n::finite ⇒ real"
  assumes Hf: "HESS f x :> Hf'"
      and eq: "⋀y. y ∈ A ⟹ ∇ (λz. c * f z) y = c *R ∇ f y"
      and Aop: "open A" and xA: "x ∈ A"
  shows "HESS (λy. c * f y) x :> c *R Hf'"
proof -
  have dscale: "((λy. c *R ∇ f y) has_derivative
                 (λv. c *R (Hf' *v v))) (at x)"
    using Hf[unfolded has_hessian_def]
    by (intro has_derivative_scaleR_right)
  have dtrans: "((λy. ∇ (λz. c * f z) y) has_derivative
                 (λv. c *R (Hf' *v v))) (at x)"
    using Aop dscale xA by (force simp: eq has_derivative_transform_within_open)
  have "⋀v. (c *R Hf') *v v = c *R (Hf' *v v)"
    by (simp add: matrix_vector_mult_def vec_eq_iff scaleR_sum_right,
        simp add: sum_distrib_left vector_space_over_itself.scale_scale)
  thus ?thesis
    unfolding has_hessian_def using dtrans by presburger
qed


subsection ‹Linearity of the Hessian on ‹C2› maps›

lemma hessian_add_on_C2:
  fixes f g :: "real^'n::finite ⇒ real"
  assumes Cf: "Ck_on 2 f U"
      and Cg: "Ck_on 2 g U"
      and xU: "x ∈ U"
  shows "∇2 (λy. f y + g y) x = ∇2 f x + ∇2 g x"
proof (rule vec_eq_iff[THEN iffD2], intro allI)
  fix i
  have openU: "open U"
    using Cf by (simp add: Ck_on_def)
  have Hf: "HESS f x :> ∇2 f x"
    using Cf xU by (rule Ck_2_imp_hessian_exists)
  have Hg: "HESS g x :> ∇2 g x"
    using Cg xU by (rule Ck_2_imp_hessian_exists)
  have Hfg: "HESS (λy. f y + g y) x :> ∇2 (λy. f y + g y) x"
    using Ck_on_add[OF Cf Cg] xU by (rule Ck_2_imp_hessian_exists)
  let ?φ = "λy. ∇ (λz. f z + g z) y $ i"
  let ?ψ = "λy. ∇ f y $ i + ∇ g y $ i"
  have eqU: "⋀y. y ∈ U ⟹ ?φ y = ?ψ y"
  proof -
    fix y assume yU: "y ∈ U"
    have Gf: "GDERIV f y :> ∇ f y"
      using Ck_2_imp_gradient_exists[OF Cf yU]
      by (blast intro: grad_fun_satisfies_GDERIV)
    have Gg: "GDERIV g y :> ∇ g y"
      using Ck_2_imp_gradient_exists[OF Cg yU]
      by (blast intro: grad_fun_satisfies_GDERIV)
    have "GDERIV (λz. f z + g z) y :> ∇ f y + ∇ g y"
      by (rule GDERIV_add[OF Gf Gg])
    hence "∇ (λz. f z + g z) y = ∇ f y + ∇ g y"
      by (rule grad_fun_eq)
    thus "?φ y = ?ψ y" by simp
  qed
  have Grow_f: "GDERIV (λy. ∇ f y $ i) x :> (∇2 f x) $ i"
    by (rule HESS_row_gradient[OF Hf])
  have Grow_g: "GDERIV (λy. ∇ g y $ i) x :> (∇2 g x) $ i"
    by (rule HESS_row_gradient[OF Hg])
  have Gψ: "GDERIV ?ψ x :> ((∇2 f x + ∇2 g x) $ i)"
    using GDERIV_add[OF Grow_f Grow_g] by simp
  have Dψ: "(?ψ has_derivative (λv. v ∙ ((∇2 f x + ∇2 g x) $ i))) (at x)"
    using Gψ unfolding gderiv_def by simp
  have Dφ: "(?φ has_derivative (λv. v ∙ ((∇2 f x + ∇2 g x) $ i))) (at x)"
    by (smt (verit, best) Dψ eqU has_derivative_transform_within_open openU xU)
  have Gφ: "GDERIV ?φ x :> ((∇2 f x + ∇2 g x) $ i)"
    using Dφ unfolding gderiv_def by simp
  show "∇2 (λy. f y + g y) x $ i = (∇2 f x + ∇2 g x) $ i"
    using Gφ HESS_row_eq Hfg grad_fun_eq by fastforce
qed

lemma hessian_scaleR_on_C2:
  fixes f :: "real^'n::finite ⇒ real"
  assumes Cf: "Ck_on 2 f U"
      and xU: "x ∈ U"
  shows "∇2 (λy. c * f y) x = c *R ∇2 f x"
proof (rule vec_eq_iff[THEN iffD2], intro allI)
  fix i
  have openU: "open U"
    using Cf by (simp add: Ck_on_def)
  have Hf: "HESS f x :> ∇2 f x"
    using Cf xU by (rule Ck_2_imp_hessian_exists)
  have Hcf: "HESS (λy. c * f y) x :> ∇2 (λy. c * f y) x"
    using Ck_on_scaleR[OF Cf] xU by (subst Ck_2_imp_hessian_exists, auto)
  let ?φ = "λy. ∇ (λz. c * f z) y $ i"
  let ?ψ = "λy. c * (∇ f y $ i)"
  have eqU: "⋀y. y ∈ U ⟹ ?φ y = ?ψ y"
  proof -
    fix y assume yU: "y ∈ U"
    have Gf: "GDERIV f y :> ∇ f y"
      using Ck_2_imp_gradient_exists[OF Cf yU]
      by (blast intro: grad_fun_satisfies_GDERIV)
    have "GDERIV (λz. c * f z) y :> c *R ∇ f y"
      by (rule GDERIV_cmult[OF Gf])
    hence "∇ (λz. c * f z) y = c *R ∇ f y"
      by (rule grad_fun_eq)
    thus "?φ y = ?ψ y" by simp
  qed
  have Grow_f: "GDERIV (λy. ∇ f y $ i) x :> (∇2 f x) $ i"
    by (rule HESS_row_gradient[OF Hf])
  have Gψ: "GDERIV ?ψ x :> c *R ((∇2 f x) $ i)"
    using GDERIV_cmult[OF Grow_f] by simp
  have Dψ: "(?ψ has_derivative (λv. v ∙ (c *R ((∇2 f x) $ i)))) (at x)"
    using Gψ unfolding gderiv_def by simp
  have Dφ: "(?φ has_derivative (λv. v ∙ (c *R ((∇2 f x) $ i)))) (at x)"
    using Dψ openU xU by (fastforce simp: eqU has_derivative_transform_within_open)
  have Gφ: "GDERIV ?φ x :> c *R ((∇2 f x) $ i)"
    using Dφ unfolding gderiv_def by simp
  show "∇2 (λy. c * f y) x $ i = (c *R ∇2 f x) $ i"
    using Gφ HESS_row_eq Hcf grad_fun_eq by fastforce
qed

lemma hessian_sub_on_C2:
  fixes f g :: "real^'n::finite ⇒ real"
  assumes Cf: "Ck_on 2 f U"
      and Cg: "Ck_on 2 g U"
      and xU: "x ∈ U"
  shows "∇2 (λy. f y - g y) x = ∇2 f x - ∇2 g x"
proof -
  have "∇2 (λy. f y + (-1) * g y) x = ∇2 f x + (-1) *R ∇2 g x"
  proof (subst hessian_add_on_C2)
    show "Ck_on 2 f U"
      by (rule Cf)
    show "Ck_on 2 (λy. (-1) * g y) U"
      using Ck_on_scaleR[OF Cg] by (metis ext real_scaleR_def)
    show "x ∈ U"
      by (rule xU)
    show "∇2 f x + ∇2 (λy. - 1 * g y) x = ∇2 f x + - 1 *R ∇2 g x"
      by (metis Cg hessian_scaleR_on_C2 xU)
  qed
  thus ?thesis by simp
qed

lemma hessian_sum_on_C2:
  fixes F :: "'i ⇒ real^'n::finite ⇒ real"
  assumes fin: "finite I"
      and C2: "⋀i. i ∈ I ⟹ Ck_on 2 (F i) U"
      and xU: "x ∈ U"
  shows "∇2 (λy. ∑i∈I. F i y) x = (∑i∈I. ∇2 (F i) x)"
  using fin C2
proof (induction rule: finite_induct)
  case empty
  show ?case by (simp add: hessian_const_zero)
next
  case (insert i I)
  have Ci: "Ck_on 2 (F i) U"
    using insert.prems by simp
  have openU: "open U"
    using Ci by (simp add: Ck_on_def)
  have C2_I: "⋀j. j ∈ I ⟹ Ck_on 2 (F j) U"
    using insert.prems by simp
  have CI: "Ck_on 2 (λy. ∑j∈I. F j y) U"
  proof (cases "I = {}")
    case True
    then show ?thesis
      using Ck_on_const[OF openU] by simp
  next
    case False
    then show ?thesis
      using Ck_on_sum[OF insert.hyps(1) False C2_I]
      by presburger
  qed
  have IH: "∇2 (λy. ∑j∈I. F j y) x = (∑j∈I. ∇2 (F j) x)"
    using insert.IH C2_I by blast
  have "∇2 (λy. ∑j∈insert i I. F j y) x
        = ∇2 (λy. F i y + (∑j∈I. F j y)) x"
    by (simp add: insert.hyps(1,2))
  also have "… = ∇2 (F i) x + ∇2 (λy. ∑j∈I. F j y) x"
    by (rule hessian_add_on_C2[OF Ci CI xU])
  also have "… = ∇2 (F i) x + (∑j∈I. ∇2 (F j) x)"
    by (simp add: IH)
  also have "… = (∑j∈insert i I. ∇2 (F j) x)"
    using insert.hyps by simp
  finally show ?case.
qed

lemma second_directional_derivative_eq_hessian_quadratic_form:
  fixes f :: "real^'n::finite ⇒ real"
  assumes C2: "Ck_on 2 f U"
      and xU: "x ∈ U"
  shows "frechet_derivative (λy. frechet_derivative f (at y) v) (at x) v
       = v ∙ ((∇2 f x) *v v)"
proof -
  have openU: "open U"
    using C2 by (simp add: Ck_on_def)

  have H: "HESS f x :> ∇2 f x"
    using C2 xU by (rule Ck_2_imp_hessian_exists)

  have eqU: "⋀y. y ∈ U ⟹ frechet_derivative f (at y) v = v ∙ ∇ f y"
  proof -
    fix y
    assume yU: "y ∈ U"

    from Ck_2_imp_gradient_exists[OF C2 yU]
    obtain g where g: "GDERIV f y :> g"
      by blast

    have Gy: "GDERIV f y :> ∇ f y"
      using g by (rule grad_fun_satisfies_GDERIV)

    have "(f has_derivative (λw. w ∙ ∇ f y)) (at y)"
      using Gy unfolding gderiv_def by simp
    hence "frechet_derivative f (at y) = (λw. w ∙ ∇ f y)"
      by (subst frechet_derivative_at, auto)

    thus "frechet_derivative f (at y) v = v ∙ ∇ f y"
      by simp
  qed
  have "∃A. open A ∧ x ∈ A ∧ (∀y∈A. frechet_derivative f (at y) v = v ∙ ∇ f y)"
    using openU xU eqU by blast
  then have ev_eq: "eventually (λy. frechet_derivative f (at y) v = v ∙ ∇ f y) (nhds x)"
    by (simp add: eventually_nhds)
  have Dgrad: "(∇ f has_derivative (*v) (∇2 f x)) (at x)"
    using H unfolding has_hessian_def by simp
  have Dcomp: "((λy. v ∙ ∇ f y) has_derivative (λh. v ∙ (((*v) (∇2 f x)) h))) (at x)"
    using Dgrad by (auto intro!: derivative_eq_intros)
  have Dfd: "((λy. frechet_derivative f (at y) v) has_derivative
       (λh. v ∙ (((*v) (∇2 f x)) h))) (at x)"
    by (metis (no_types, lifting) Dcomp eqU has_derivative_transform_within_open openU xU)
  have FD: "frechet_derivative (λy. frechet_derivative f (at y) v) (at x)
    = (λh. v ∙ (((*v) (∇2 f x)) h))"
    by (metis Dfd frechet_derivative_at)
  show ?thesis
    by (simp add: FD)
qed

subsection ‹Outer product of vectors›

definition outer_prod :: "real^'n ⇒ real^'n ⇒ real^'n^'n"
  where "outer_prod a b = (χ i j. a $ i * b $ j)"

lemma outer_prod_component [simp]:
  "outer_prod a b $ i $ j = a $ i * b $ j"
  by (simp add: outer_prod_def)

lemma outer_prod_row:
  "outer_prod a b $ i = (a $ i) *R b"
  by (simp add: vec_eq_iff outer_prod_def)

lemma outer_prod_commute:
  "transpose (outer_prod a b) = outer_prod b a"
  by (simp add: vec_eq_iff transpose_def outer_prod_def mult.commute)

lemma outer_prod_add_left:
  "outer_prod (a + b) c = outer_prod a c + outer_prod b c"
  by (simp add: vec_eq_iff outer_prod_def distrib_right)

lemma outer_prod_add_right:
  "outer_prod a (b + c) = outer_prod a b + outer_prod a c"
  by (simp add: vec_eq_iff outer_prod_def distrib_left)

lemma outer_prod_scaleR_left:
  "outer_prod (c *R a) b = c *R outer_prod a b"
  by (simp add: vec_eq_iff outer_prod_def)

lemma outer_prod_scaleR_right:
  "outer_prod a (c *R b) = c *R outer_prod a b"
  by (simp add: vec_eq_iff outer_prod_def)

lemma outer_prod_zero_left [simp]:
  "outer_prod 0 b = 0"
  by (simp add: vec_eq_iff outer_prod_def)

lemma outer_prod_zero_right [simp]:
  "outer_prod a 0 = 0"
  by (simp add: vec_eq_iff outer_prod_def)

lemma outer_prod_mult_vec:
  "outer_prod a b *v v = (b ∙ v) *R a"
  by (simp add: matrix_vector_mul_component outer_prod_row vec_eq_iff)


subsection ‹Gradient product rule›

lemma grad_fun_mult:
  fixes f g :: "real^'n::finite ⇒ real"
  assumes "∃gf. GDERIV f x :> gf"
      and "∃gg. GDERIV g x :> gg"
  shows "∇ (λy. f y * g y) x = f x *R ∇ g x + g x *R ∇ f x"
proof -
  have Gf: "GDERIV f x :> ∇ f x"
    using assms(1) by (blast intro: grad_fun_satisfies_GDERIV)
  have Gg: "GDERIV g x :> ∇ g x"
    using assms(2) by (blast intro: grad_fun_satisfies_GDERIV)
  have "GDERIV (λy. f y * g y) x :> f x *R ∇ g x + g x *R ∇ f x"
    by (rule GDERIV_mult[OF Gf Gg])
  thus ?thesis
    by (rule grad_fun_eq)
qed


subsection ‹Hessian product rule›

text ‹
  ‹∇2(fg) = f ∇2g + g ∇2f + ∇f ⊗ ∇g + ∇g ⊗ ∇f› for ‹C2› functions
  ‹f, g : ℝn → ℝ›, where ‹⊗› is the outer product @{const outer_prod}.
›

lemma hessian_mult_on_C2:
  fixes f g :: "real^'n::finite ⇒ real"
  assumes Cf: "Ck_on 2 f U"
      and Cg: "Ck_on 2 g U"
      and xU: "x ∈ U"
  shows "∇2 (λy. f y * g y) x =
           f x *R ∇2 g x + g x *R ∇2 f x
         + outer_prod (∇ f x) (∇ g x)
         + outer_prod (∇ g x) (∇ f x)"
proof (rule vec_eq_iff[THEN iffD2], intro allI)
  fix i

  have openU: "open U"
    using Cf by (simp add: Ck_on_def)

  (* Hessians exist *)
  have Hf: "HESS f x :> ∇2 f x"
    using Cf xU by (rule Ck_2_imp_hessian_exists)
  have Hg: "HESS g x :> ∇2 g x"
    using Cg xU by (rule Ck_2_imp_hessian_exists)

  have Cfh: "Ck_on 2 (λy. f y * g y) U"
    by (simp add: Cf Cg Ck_on_mult)

  have Hfg: "HESS (λy. f y * g y) x :> ∇2 (λy. f y * g y) x"
    using Cfh xU by (rule Ck_2_imp_hessian_exists)

  (* Gradient existence on U *)
  have Gf_at: "⋀y. y ∈ U ⟹ GDERIV f y :> ∇ f y"
    using Ck_2_imp_gradient_exists[OF Cf]
    by (blast intro: grad_fun_satisfies_GDERIV)
  have Gg_at: "⋀y. y ∈ U ⟹ GDERIV g y :> ∇ g y"
    using Ck_2_imp_gradient_exists[OF Cg]
    by (blast intro: grad_fun_satisfies_GDERIV)

  (* Row gradients of Hessians *)
  have Hf_row: "GDERIV (λy. ∇ f y $ i) x :> (∇2 f x) $ i"
    by (rule HESS_row_gradient[OF Hf])
  have Hg_row: "GDERIV (λy. ∇ g y $ i) x :> (∇2 g x) $ i"
    by (rule HESS_row_gradient[OF Hg])

  (* The i-th component of ∇(fg) *)
  (* On U: ∇(fg)(y) $ i = f(y) * ∇g(y) $ i + g(y) * ∇f(y) $ i *)
  let ?φ = "λy. ∇ (λz. f z * g z) y $ i"
  let ?ψ = "λy. f y * (∇ g y $ i) + g y * (∇ f y $ i)"

  have eqU: "⋀y. y ∈ U ⟹ ?φ y = ?ψ y"
  proof -
    fix y assume yU: "y ∈ U"
    have "GDERIV (λz. f z * g z) y :> f y *R ∇ g y + g y *R ∇ f y"
      by (rule GDERIV_mult[OF Gf_at[OF yU] Gg_at[OF yU]])
    hence "∇ (λz. f z * g z) y = f y *R ∇ g y + g y *R ∇ f y"
      by (rule grad_fun_eq)
    thus "?φ y = ?ψ y" by simp
  qed

  (* Gradient of y ↦ f(y) * ∇g(y)$i *)
  have Gf_x: "GDERIV f x :> ∇ f x"
    using Gf_at[OF xU] .
  have Gg_x: "GDERIV g x :> ∇ g x"
    using Gg_at[OF xU] .

  have G_term1: "GDERIV (λy. f y * (∇ g y $ i)) x :>
                   f x *R (∇2 g x) $ i + (∇ g x $ i) *R ∇ f x"
    by (rule GDERIV_mult[OF Gf_x Hg_row])

  (* Gradient of y ↦ g(y) * ∇f(y)$i *)
  have G_term2: "GDERIV (λy. g y * (∇ f y $ i)) x :>
                   g x *R (∇2 f x) $ i + (∇ f x $ i) *R ∇ g x"
    by (rule GDERIV_mult[OF Gg_x Hf_row])

  (* Gradient of ψ by addition *)
  have Gψ: "GDERIV ?ψ x :>
               (f x *R (∇2 g x) $ i + (∇ g x $ i) *R ∇ f x)
             + (g x *R (∇2 f x) $ i + (∇ f x $ i) *R ∇ g x)"
    by (rule GDERIV_add[OF G_term1 G_term2])

  (* Transfer from ψ to φ using agreement on U *)
  have Dψ: "(?ψ has_derivative
      (λv. v ∙ ((f x *R (∇2 g x) $ i + (∇ g x $ i) *R ∇ f x)
              + (g x *R (∇2 f x) $ i + (∇ f x $ i) *R ∇ g x)))) (at x)"
    using Gψ unfolding gderiv_def by simp

  have Dφ: "(?φ has_derivative
      (λv. v ∙ ((f x *R (∇2 g x) $ i + (∇ g x $ i) *R ∇ f x)
              + (g x *R (∇2 f x) $ i + (∇ f x $ i) *R ∇ g x)))) (at x)"
    using Dψ openU xU by (fastforce simp: eqU has_derivative_transform_within_open)

  have Gφ: "GDERIV ?φ x :>
               (f x *R (∇2 g x) $ i + (∇ g x $ i) *R ∇ f x)
             + (g x *R (∇2 f x) $ i + (∇ f x $ i) *R ∇ g x)"
    using Dφ unfolding gderiv_def by simp

  (* Assemble the row *)
  have row_eq: "∇ ?φ x =
      (f x *R (∇2 g x) $ i + (∇ g x $ i) *R ∇ f x)
    + (g x *R (∇2 f x) $ i + (∇ f x $ i) *R ∇ g x)"
    by (rule grad_fun_eq[OF Gφ])

  have lhs: "∇2 (λy. f y * g y) x $ i = ∇ ?φ x"
    using HESS_row_eq[OF Hfg] by simp

  (* Express the RHS in terms of the target matrix *)
  let ?M = "f x *R ∇2 g x + g x *R ∇2 f x
          + outer_prod (∇ f x) (∇ g x)
          + outer_prod (∇ g x) (∇ f x)"

  have rhs: "?M $ i =
      (f x *R (∇2 g x) $ i + (∇ g x $ i) *R ∇ f x)
    + (g x *R (∇2 f x) $ i + (∇ f x $ i) *R ∇ g x)"
    by (simp add: vec_eq_iff outer_prod_row algebra_simps)

  show "∇2 (λy. f y * g y) x $ i = ?M $ i"
    using lhs row_eq rhs by simp
qed

lemma grad_fun_compose:
  fixes g :: "real^'m::finite ⇒ real"
    and F :: "real^'n::finite ⇒ real^'m"
  assumes "∃gg. GDERIV g (F x) :> gg"
      and "F differentiable (at x)"
  shows "∇ (λy. g (F y)) x = transpose (jacobian F (at x)) *v ∇ g (F x)"
proof -
  have Gg: "GDERIV g (F x) :> ∇ g (F x)"
    using assms(1) by (blast intro: grad_fun_satisfies_GDERIV)
  have "GDERIV (λy. g (F y)) x :> transpose (jacobian F (at x)) *v ∇ g (F x)"
    by (rule GDERIV_compose'[OF Gg assms(2)])
  thus ?thesis
    by (rule grad_fun_eq)
qed


subsection ‹Component closure›

lemma Ck_on_component:
  fixes F :: "'a::real_normed_vector ⇒ real^'m::finite"
  assumes "Ck_on k F U"
  shows "Ck_on k (λx. F x $ r) U"
  by (rule Ck_on_bounded_linear_compose[OF bounded_linear_vec_nth assms])


text ‹
  For ‹C2› maps ‹g› and ‹F›:
  ‹∇2(g ∘ F)(x) = JT ** ∇2g(F x) ** J + Σr (∇g(F x) $ r) *R ∇2Fr(x)›,
  where ‹J = jacobian F (at x)› and ‹Fr y = F y $ r›.
›

lemma hessian_compose_on_C2:
  fixes g :: "real^'m::finite ⇒ real"
    and F :: "real^'n::finite ⇒ real^'m"
  assumes Cg: "Ck_on 2 g V"
      and CF: "Ck_on 2 F U"
      and FUV: "⋀y. y ∈ U ⟹ F y ∈ V"
      and xU: "x ∈ U"
  shows "∇2 (λy. g (F y)) x =
           transpose (jacobian F (at x)) ** ∇2 g (F x) ** jacobian F (at x)
         + (∑r∈UNIV. (∇ g (F x) $ r) *R ∇2 (λy. F y $ r) x)"
         (is "?LHS = ?RHS")
proof (rule vec_eq_iff[THEN iffD2], intro allI)
  fix i :: 'n

  have openU: "open U" using CF by (simp add: Ck_on_def)
  have openV: "open V" using Cg by (simp add: Ck_on_def)

  (* C2 closure: g ∘ F is C2 on U *)
  have CgF: "Ck_on 2 (λy. g (F y)) U"
    using Ck_on_compose[OF Cg CF FUV] .

  (* Component C2 *)
  have CF_r: "⋀r. Ck_on 2 (λy. F y $ r) U"
    using CF by (rule Ck_on_component)

  (* Hessians exist *)
  have HgF: "HESS (λy. g (F y)) x :> ∇2 (λy. g (F y)) x"
    using CgF xU by (rule Ck_2_imp_hessian_exists)
  have Hg: "HESS g (F x) :> ∇2 g (F x)"
    using Cg FUV[OF xU] by (rule Ck_2_imp_hessian_exists)
  have HF_r: "⋀r. HESS (λy. F y $ r) x :> ∇2 (λy. F y $ r) x"
    using CF_r xU by (rule Ck_2_imp_hessian_exists)

  (* Differentiability of F on U *)
  have F_diff: "⋀y. y ∈ U ⟹ F differentiable (at y)"
  proof -
    fix y assume "y ∈ U"
    then have "Ck_at 2 F y"
      using CF by (simp add: Ck_on_def)
    thus "F differentiable (at y)"
      by (metis Ck_at.simps(2) Suc_1)
  qed

  (* Gradient existence *)
  have Gg_at: "⋀z. z ∈ V ⟹ GDERIV g z :> ∇ g z"
    using Ck_2_imp_gradient_exists[OF Cg]
    by (blast intro: grad_fun_satisfies_GDERIV)

  have GF_r_at: "⋀r y. y ∈ U ⟹ GDERIV (λy. F y $ r) y :> ∇ (λy. F y $ r) y"
    using Ck_2_imp_gradient_exists[OF CF_r]
    by (blast intro: grad_fun_satisfies_GDERIV)

  (* Row gradients of component Hessians *)
  have HF_r_row: "⋀r. GDERIV (λy. ∇ (λz. F z $ r) y $ i) x :> (∇2 (λy. F y $ r) x) $ i"
    by (rule HESS_row_gradient[OF HF_r])

  (* Row gradient of the Hessian of g *)
  have Hg_row: "⋀r. GDERIV (λz. ∇ g z $ r) (F x) :> (∇2 g (F x)) $ r"
    by (rule HESS_row_gradient[OF Hg])

  (* On U, ∇(g ∘ F)(y) $ i = Σ_r (∇(F_r)(y) $ i) * (∇g(F(y)) $ r) *)
  let ?φ = "λy. ∇ (λz. g (F z)) y $ i"
  let ?ψ = "λy. ∑r∈UNIV. ∇ (λz. F z $ r) y $ i * ∇ g (F y) $ r"

  have eqU: "⋀y. y ∈ U ⟹ ?φ y = ?ψ y"
  proof -
    fix y :: "real^'n"
    assume yU: "y ∈ U"

    have Fy_V: "F y ∈ V" using FUV[OF yU] .
    have Gy: "GDERIV g (F y) :> ∇ g (F y)"
      using Gg_at[OF Fy_V] .
    have Fy_diff: "F differentiable (at y)"
      using F_diff[OF yU] .

    have grad_comp: "∇ (λz. g (F z)) y = transpose (jacobian F (at y)) *v ∇ g (F y)"
      by (rule grad_fun_compose[where g=g and F=F], blast intro: Gy, rule Fy_diff)

    have "?φ y = (transpose (jacobian F (at y)) *v ∇ g (F y)) $ i"
      using grad_comp by simp
    also have "… = (∑r∈UNIV. transpose (jacobian F (at y)) $ i $ r * ∇ g (F y) $ r)"
      by (simp add: matrix_vector_mult_def)
    also have "… = (∑r∈UNIV. jacobian F (at y) $ r $ i * ∇ g (F y) $ r)"
      by (simp add: transpose_def)
    also have "… = ?ψ y"
    proof (rule sum.cong[OF refl])
      fix r :: 'm
      assume "r ∈ UNIV"

      have Fr_diff: "(λz. F z $ r) differentiable (at y)"
      proof -
        have FD: "(F has_derivative frechet_derivative F (at y)) (at y)"
          using Fy_diff frechet_derivative_works[THEN iffD1] by blast

        have Hcomp:"((λz. F z ∙ axis r 1) has_derivative
                     (λh. frechet_derivative F (at y) h ∙ axis r 1)) (at y within UNIV)"
          using FD
          by (subst (asm) has_derivative_componentwise_within[where S = UNIV],
              auto simp: Basis_vec_def)

        have comp_fun: "(λz. F z ∙ axis r 1) = (λz. F z $ r)"
          by (rule ext) (simp add: cart_eq_inner_axis)
        have coord_D: "((λz. F z $ r) has_derivative
                        (λh. frechet_derivative F (at y) h $ r)) (at y)"
          using Hcomp by (simp add: inner_axis)
        then show ?thesis
          unfolding differentiable_def by blast
      qed

      show "jacobian F (at y) $ r $ i * ∇ g (F y) $ r =  ∇ (λz. F z $ r) y $ i * ∇ g (F y) $ r"
      proof -
        have Jcomp: "jacobian F (at y) $ r $ i = frechet_derivative (λz. F z $ r) (at y) (axis i 1)"
          using jacobian_component[OF Fy_diff, of r i] by simp

        have GFr: "GDERIV (λz. F z $ r) y :> ∇ (λz. F z $ r) y"
          using Fr_diff_imp_gradient_exists[OF Fr_diff]
          by (blast intro: grad_fun_satisfies_GDERIV)

        have DFr: "((λz. F z $ r) has_derivative (λh. h ∙ ∇ (λz. F z $ r) y)) (at y)"
          using GFr unfolding gderiv_def by simp

        have FD_eq:  "frechet_derivative (λz. F z $ r) (at y) = (λh. h ∙ ∇ (λz. F z $ r) y)"
          by (subst frechet_derivative_at[OF DFr], simp)

        have "frechet_derivative (λz. F z $ r) (at y) (axis i 1)
              = axis i 1 ∙ ∇ (λz. F z $ r) y"
          by (simp add: FD_eq)
        also have "... = ∇ (λz. F z $ r) y $ i"
          using inner_commute by (simp add: cart_eq_inner_axis, auto)
        finally show ?thesis
          using Jcomp by simp
      qed
      then have "jacobian F (at y) $ r $ i = ∇ (λz. F z $ r) y $ i"
        using jacobian_component[OF Fy_diff]
        by (metis (mono_tags, lifting) Fr_diff Fr_diff_imp_gradient_exists cart_eq_inner_axis
            frechet_derivative_at grad_fun_eq gderiv_def inner_commute)
      qed
      thus "∇ (λz. g (F z)) y $ i = (∑r∈UNIV. ∇ (λz. F z $ r) y $ i * ∇ g (F y) $ r)"
        using calculation by presburger
  qed

  (* Differentiate ψ at x using GDERIV_sum, GDERIV_mult *)
  (* ψ y = Σ_r a_r y * b_r y, where a_r y = ∇(F_r)(y) $ i and b_r y = ∇g(F y) $ r *)

  have Ga_r: "⋀r. GDERIV (λy. ∇ (λz. F z $ r) y $ i) x :> (∇2 (λy. F y $ r) x) $ i"
    using HF_r_row .

  have Gb_r: "⋀r. GDERIV (λy. ∇ g (F y) $ r) x :> transpose (jacobian F (at x)) *v ((∇2 g (F x)) $ r)"
  proof -
    fix r :: 'm
    have "GDERIV (λz. ∇ g z $ r) (F x) :> (∇2 g (F x)) $ r"
      using Hg_row .
    thus "GDERIV (λy. ∇ g (F y) $ r) x :> transpose (jacobian F (at x)) *v ((∇2 g (F x)) $ r)"
      by (rule GDERIV_compose'[OF _ F_diff[OF xU]])
  qed

  have G_term_r: "⋀r. GDERIV (λy. ∇ (λz. F z $ r) y $ i * ∇ g (F y) $ r) x :>
      ∇ (λz. F z $ r) x $ i *R (transpose (jacobian F (at x)) *v ((∇2 g (F x)) $ r))
    + ∇ g (F x) $ r *R ((∇2 (λy. F y $ r) x) $ i)"
    by (rule GDERIV_mult[OF Ga_r Gb_r])

  define G_r where "G_r r =
      ∇ (λz. F z $ r) x $ i *R (transpose (jacobian F (at x)) *v ((∇2 g (F x)) $ r))
    + ∇ g (F x) $ r *R ((∇2 (λy. F y $ r) x) $ i)" for r

  have Gψ: "GDERIV ?ψ x :> (∑r∈UNIV. G_r r)"
      using G_term_r by (simp add: G_r_def GDERIV_sum)

  (* Transfer from ψ to φ *)
  have Dψ: "(?ψ has_derivative (λv. v ∙ (∑r∈UNIV. G_r r))) (at x)"
    using Gψ unfolding gderiv_def .

  have Dφ: "(?φ has_derivative (λv. v ∙ (∑r∈UNIV. G_r r))) (at x)"
    using Dψ openU xU by (force simp: eqU has_derivative_transform_within_open)

  have Gφ: "GDERIV ?φ x :> (∑r∈UNIV. G_r r)"
    using Dφ unfolding gderiv_def by simp

  have row_eq: "∇ ?φ x = (∑r∈UNIV. G_r r)"
    by (rule grad_fun_eq[OF Gφ])

  have lhs: "?LHS $ i = ∇ ?φ x"
    using HESS_row_eq[OF HgF] by simp

  (* Match the RHS row *)

  have sum_second: "(∑r∈UNIV. ∇ g (F x) $ r *R ((∇2 (λy. F y $ r) x) $ i))
                  = (∑r∈UNIV. (∇ g (F x) $ r) *R ∇2 (λy. F y $ r) x) $ i"
    by simp

  (* matrix algebra: the key identity relating the sum to J^T H J *)
  have sum_first: "(∑r∈UNIV. ∇ (λz. F z $ r) x $ i *R
                      (transpose (jacobian F (at x)) *v ((∇2 g (F x)) $ r)))
                 = (transpose (jacobian F (at x)) ** ∇2 g (F x) ** jacobian F (at x)) $ i"
  proof -
    have grad_jac: "∇ (λz. F z $ r) x $ i = jacobian F (at x) $ r $ i" for r
    proof -
      have Fr_diff: "(λz. F z $ r) differentiable (at x)"
        using F_diff[OF xU] by (metis CF_r Ck_at.simps(2) Ck_on_def Suc_1 xU)
      have GFr: "GDERIV (λz. F z $ r) x :> ∇ (λz. F z $ r) x"
        using Fr_diff_imp_gradient_exists[OF Fr_diff]
        by (blast intro: grad_fun_satisfies_GDERIV)
      have FD_eq: "frechet_derivative (λz. F z $ r) (at x) = (λh. h ∙ ∇ (λz. F z $ r) x)"
        using GFr unfolding gderiv_def  by (metis frechet_derivative_at)
      have "jacobian F (at x) $ r $ i = frechet_derivative (λz. F z $ r) (at x) (axis i 1)"
        using jacobian_component[OF F_diff[OF xU]].
      also have "… = ∇ (λz. F z $ r) x $ i"
        by (simp add: FD_eq cart_eq_inner_axis inner_commute,
            metis (no_types, lifting) ext FD_eq cart_eq_inner_axis)
      finally show ?thesis by simp
    qed
    then have "(∑r∈UNIV. ∇ (λz. F z $ r) x $ i *R  (transpose (jacobian F (at x)) *v ((∇2 g (F x)) $ r)))
        = (∑r∈UNIV. jacobian F (at x) $ r $ i *R   (transpose (jacobian F (at x)) *v ((∇2 g (F x)) $ r)))"
      by simp
    also have "… = (transpose (jacobian F (at x)) ** ∇2 g (F x) ** jacobian F (at x)) $ i"
      by (rule row_transpose_mult_both[symmetric])
    finally show ?thesis.
  qed
  then have "(∑r∈UNIV. G_r r) = ?RHS $ i"
    unfolding G_r_def sum.distrib vector_add_component
    using sum_second by (rule arg_cong2[where f = "(+)"])
  then show "?LHS $ i = ?RHS $ i"
    using lhs row_eq by simp
qed


subsection ‹Affine composition (special case)›

text ‹For affine ‹F y = A *v y + b› the chain rule becomes
  ‹∇2(g ∘ F)(x) = AT ** ∇2g(A *v x + b) ** A›.›

lemma hessian_affine_compose_on_C2:
  fixes g :: "real^'m::finite ⇒ real"
    and A :: "real^'n^'m"
    and b :: "real^'m"
  assumes Cg: "Ck_on 2 g V"
      and sub: "⋀y. y ∈ U ⟹ A *v y + b ∈ V"
      and oU: "open U"
      and xU: "x ∈ U"
  shows "∇2 (λy. g (A *v y + b)) x = transpose A ** ∇2 g (A *v x + b) ** A"
proof -
  define F where "F y = A *v y + b" for y
  have bl: "bounded_linear ((*v) A)"
    by simp
  have CF: "Ck_on 2 F U"
    unfolding F_def using Ck_on_add[OF Ck_on_bounded_linear[OF bl oU] Ck_on_const[OF oU]] .
  have F_diff: "⋀y. F differentiable (at y)"
    unfolding F_def by (simp add: bounded_linear_imp_differentiable)
  have jac_eq: "jacobian F (at y) = A" for y
    unfolding jacobian_def F_def by (metis bl bounded_linear_imp_has_derivative
              frechet_derivative_at has_derivative_add_const matrix_of_matrix_vector_mul)
  have comp_hess_zero: "∇2 (λy. F y $ r) x = 0" for r
  proof -
    have fn_eq: "(λy. F y $ r) = (λy. b $ r + y ∙ (A $ r))"
    proof (rule ext)
      fix y :: "real^'n"
      have "F y $ r = (A *v y + b) $ r"
        by (simp add: F_def)
      also have "… = (∑j∈UNIV. A $ r $ j * y $ j) + b $ r"
        by (simp add: matrix_vector_mult_def)
      also have "… = b $ r + y ∙ (A $ r)"
        by (simp add: inner_vec_def, meson mult.commute)
      finally show "F y $ r = b $ r + y ∙ (A $ r)".
    qed
    have "HESS (λy. b $ r + y ∙ (A $ r)) x :> 0"
      by (rule HESS_affine_zero)
    hence "HESS (λy. F y $ r) x :> 0"
      by (simp add: fn_eq)
    thus ?thesis
      by (metis hess_fun_eq)
  qed
  then have grad_zero_sum: "(∑r∈UNIV. ∇ g (F x) $ r *R ∇2 (λy. F y $ r) x) = 0"
    by simp
  have "∇2 (λy. g (F y)) x =
         transpose (jacobian F (at x)) ** ∇2 g (F x) ** jacobian F (at x)
       + (∑r∈UNIV. ∇ g (F x) $ r *R ∇2 (λy. F y $ r) x)"
    by (rule hessian_compose_on_C2[OF Cg CF _ xU], simp add: F_def sub)
  also have "… = transpose A ** ∇2 g (F x) ** A + 0"
    by (simp add: jac_eq grad_zero_sum)
  also have "… = transpose A ** ∇2 g (A *v x + b) ** A"
    by (simp add: F_def)
  ultimately show ?thesis by (simp add: F_def)
qed

subsection ‹Summary of the hierarchy›

text ‹
  const‹Ck_at› implies const‹k_times_Fr_differentiable_at›, and const‹Ck_on› implies
  const‹k_times_Fr_differentiable_on›.  For ‹f :: real ⇒ real›,
  const‹k_times_Fr_differentiable_at› agrees with const‹k_times_differentiable_at›,
  and const‹Ck_on› with const‹C_k_on›.
›

end