Theory Weighted_Arithmetic_Geometric_Mean

section ‹The Weighted Arithmetic--Geometric Mean Inequality›
theory Weighted_Arithmetic_Geometric_Mean
  imports Complex_Main
begin

subsection ‹Auxiliary Facts›

lemma root_powr_inverse': "0 < n ⟹ 0 ≤ x ⟹ root n x = x powr (1/n)"
  by (cases "x = 0") (auto simp: root_powr_inverse)

lemma powr_sum_distrib_real_right:
  assumes "a ≠ 0"
  shows   "(∏x∈X. a powr e x :: real) = a powr (∑x∈X. e x)"
  using assms
  by (induction X rule: infinite_finite_induct) (auto simp: powr_add)

lemma powr_sum_distrib_real_left:
  assumes "⋀x. x ∈ X ⟹ a x ≥ 0"
  shows   "(∏x∈X. a x powr e :: real) = (∏x∈X. a x) powr e"
  using assms
  by (induction X rule: infinite_finite_induct)
     (auto simp: powr_mult prod_nonneg)

lemma (in linordered_semidom) prod_mono_strict':
  assumes "i ∈ A"
  assumes "finite A"
  assumes "⋀i. i ∈ A ⟹ 0 ≤ f i ∧ f i ≤ g i"
  assumes "⋀i. i ∈ A ⟹ 0 < g i"
  assumes "f i < g i"
  shows   "prod f A < prod g A"
proof -
  have "prod f A = f i * prod f (A - {i})"
    using assms by (intro prod.remove)
  also have "… ≤ f i * prod g (A - {i})"
    using assms by (intro mult_left_mono prod_mono) auto
  also have "… < g i * prod g (A - {i})"
    using assms by (intro mult_strict_right_mono prod_pos) auto
  also have "… = prod g A"
    using assms by (intro prod.remove [symmetric])
  finally show ?thesis .
qed

lemma prod_ge_pointwise_le_imp_pointwise_eq:
  fixes f :: "'a ⇒ real"
  assumes "finite X"
  assumes ge: "prod f X ≥ prod g X"
  assumes nonneg: "⋀x. x ∈ X ⟹ f x ≥ 0"
  assumes pos: "⋀x. x ∈ X ⟹ g x > 0"
  assumes le: "⋀x. x ∈ X ⟹ f x ≤ g x" and x: "x ∈ X"
  shows   "f x = g x"
proof (rule ccontr)
  assume "f x ≠ g x"
  with le[of x] and x have "f x < g x"
    by auto
  hence "prod f X < prod g X"
    using x and le and nonneg and pos and ‹finite X› 
    by (intro prod_mono_strict') auto
  with ge show False
    by simp
qed

lemma powr_right_real_eq_iff:
  assumes "a ≥ (0 :: real)"
  shows   "a powr x = a powr y ⟷ a = 0 ∨ a = 1 ∨ x = y"
  using assms by (auto simp: powr_def)

lemma powr_left_real_eq_iff:
  assumes "a ≥ (0 :: real)" "b ≥ 0" "x ≠ 0"
  shows   "a powr x = b powr x ⟷ a = b"
  using assms by (auto simp: powr_def)

lemma exp_real_eq_one_plus_iff:
  fixes x :: real
  shows "exp x = 1 + x ⟷ x = 0"
proof (cases "x = 0")
  case False
  define f :: "real ⇒ real" where "f = (λx. exp x - 1 - x)"
  have deriv: "(f has_field_derivative (exp x - 1)) (at x)" for x
    by (auto simp: f_def intro!: derivative_eq_intros)

  have "∃z. z>min x 0 ∧ z < max x 0 ∧ f (max x 0) - f (min x 0) = 
            (max x 0 - min x 0) * (exp z - 1)"
    using MVT2[of "min x 0" "max x 0" f "λx. exp x - 1"] deriv False
    by (auto simp: min_def max_def)
  then obtain z where "z ∈ {min x 0<..<max x 0}"
     "f (max x 0) - f (min x 0) = (max x 0 - min x 0) * (exp z - 1)"
    by (auto simp: f_def)
  thus ?thesis using False
    by (cases x "0 :: real" rule: linorder_cases) (auto simp: f_def)
qed auto


subsection ‹The Inequality›

text ‹
  We first prove the equality under the assumption that all the $a_i$ and $w_i$ are positive.
›
lemma weighted_arithmetic_geometric_mean_pos:
  fixes a w :: "'a ⇒ real"
  assumes "finite X"
  assumes pos1: "⋀x. x ∈ X ⟹ a x > 0"
  assumes pos2: "⋀x. x ∈ X ⟹ w x > 0"
  assumes sum_weights: "(∑x∈X. w x) = 1"
  shows   "(∏x∈X. a x powr w x) ≤ (∑x∈X. w x * a x)"
proof -
  note nonneg1 = less_imp_le[OF pos1]
  note nonneg2 = less_imp_le[OF pos2]
  define A where "A = (∑x∈X. w x * a x)"
  define r where "r = (λx. a x / A - 1)"
  from sum_weights have "X ≠ {}" by auto
  hence "A ≠ 0"
    unfolding A_def using nonneg1 nonneg2 pos1 pos2 ‹finite X›
    by (subst sum_nonneg_eq_0_iff) force+
  moreover from nonneg1 nonneg2 have "A ≥ 0"
    by (auto simp: A_def intro!: sum_nonneg)
  ultimately have "A > 0" by simp

  have "(∏x∈X. (1 + r x) powr w x) = (∏x∈X. (a x / A) powr w x)"
    by (simp add: r_def)
  also have "… = (∏x∈X. a x powr w x) / (∏x∈X. A powr w x)"
    unfolding prod_dividef [symmetric]
    using assms pos2 ‹A > 0› by (intro prod.cong powr_divide) (auto intro: less_imp_le)
  also have "(∏x∈X. A powr w x) = exp ((∑x∈X. w x) * ln A)"
    using ‹A > 0› and ‹finite X› by (simp add: powr_def exp_sum sum_distrib_right)
  also have "(∑x∈X. w x) = 1" by fact
  also have "exp (1 * ln A) = A"
    using ‹A > 0› by simp
  finally have lhs: "(∏x∈X. (1 + r x) powr w x) = (∏x∈X. a x powr w x) / A" .

  have "(∏x∈X. exp (w x * r x)) = exp (∑x∈X. w x * r x)"
    using ‹finite X› by (simp add: exp_sum)
  also have "(∑x∈X. w x * r x) = (∑x∈X. a x * w x) / A - 1"
    using ‹A > 0› by (simp add: r_def algebra_simps sum_subtractf sum_divide_distrib sum_weights)
  also have "(∑x∈X. a x * w x) / A = 1"
    using ‹A > 0› by (simp add: A_def mult.commute)
  finally have rhs: "(∏x∈X. exp (w x * r x)) = 1" by simp

  have "(∏x∈X. a x powr w x) / A = (∏x∈X. (1 + r x) powr w x)"
    by (fact lhs [symmetric])
  also have "(∏x∈X. (1 + r x) powr w x) ≤ (∏x∈X. exp (w x * r x))"
  proof (intro prod_mono conjI)
    fix x assume x: "x ∈ X"
    have "1 + r x ≤ exp (r x)"
      by (rule exp_ge_add_one_self)
    hence "(1 + r x) powr w x ≤ exp (r x) powr w x"
      using nonneg1[of x] nonneg2[of x] x ‹A > 0›
      by (intro powr_mono2) (auto simp: r_def field_simps)
    also have "… = exp (w x * r x)"
      by (simp add: powr_def)
    finally show "(1 + r x) powr w x ≤ exp (w x * r x)" .
  qed auto
  also have "(∏x∈X. exp (w x * r x)) = 1" by (fact rhs)
  finally show "(∏x∈X. a x powr w x) ≤ A"
    using ‹A > 0› by (simp add: field_simps)
qed

text ‹
  We can now relax the positivity assumptions to non-negativity: if one of the $a_i$ is
  zero, the theorem becomes trivial (note that $0^0 = 0$ by convention for the real-valued
  power operator \<^term>‹(powr) :: real ⇒ real ⇒ real›).

  Otherwise, we can simply remove all the indices that have weight 0 and apply the
  above auxiliary version of the theorem.
›
theorem weighted_arithmetic_geometric_mean:
  fixes a w :: "'a ⇒ real"
  assumes "finite X"
  assumes nonneg1: "⋀x. x ∈ X ⟹ a x ≥ 0"
  assumes nonneg2: "⋀x. x ∈ X ⟹ w x ≥ 0"
  assumes sum_weights: "(∑x∈X. w x) = 1"
  shows   "(∏x∈X. a x powr w x) ≤ (∑x∈X. w x * a x)"
proof (cases "∃x∈X. a x = 0")
  case True
  hence "(∏x∈X. a x powr w x) = 0"
    using ‹finite X› by simp
  also have "… ≤ (∑x∈X. w x * a x)"
    by (intro sum_nonneg mult_nonneg_nonneg assms)
  finally show ?thesis .
next
  case False
  have "(∑x∈X-{x. w x = 0}. w x) = (∑x∈X. w x)"
    by (intro sum.mono_neutral_left assms) auto
  also have "… = 1" by fact
  finally have sum_weights': "(∑x∈X-{x. w x = 0}. w x) = 1" .

  have "(∏x∈X. a x powr w x) = (∏x∈X-{x. w x = 0}. a x powr w x)"
    using ‹finite X› False by (intro prod.mono_neutral_right) auto
  also have "… ≤ (∑x∈X-{x. w x = 0}. w x * a x)" using assms False
    by (intro weighted_arithmetic_geometric_mean_pos sum_weights')
       (auto simp: order.strict_iff_order)
  also have "… = (∑x∈X. w x * a x)"
    using ‹finite X› by (intro sum.mono_neutral_left) auto
  finally show ?thesis .
qed

text ‹
  We can derive the regular arithmetic/geometric mean inequality from this by simply
  setting all the weights to $\frac{1}{n}$:
›
corollary arithmetic_geometric_mean:
  fixes a :: "'a ⇒ real"
  assumes "finite X"
  defines "n ≡ card X"
  assumes nonneg: "⋀x. x ∈ X ⟹ a x ≥ 0"
  shows   "root n (∏x∈X. a x) ≤ (∑x∈X. a x) / n"
proof (cases "X = {}")
  case False
  with assms have n: "n > 0"
    by auto
  have "(∏x∈X. a x powr (1 / n)) ≤ (∑x∈X. (1 / n) * a x)"
    using n assms by (intro weighted_arithmetic_geometric_mean) auto
  also have "(∏x∈X. a x powr (1 / n)) = (∏x∈X. a x) powr (1 / n)"
    using nonneg by (subst powr_sum_distrib_real_left) auto
  also have "… = root n (∏x∈X. a x)"
    using ‹n > 0› nonneg by (subst root_powr_inverse') (auto simp: prod_nonneg)
  also have "(∑x∈X. (1 / n) * a x) = (∑x∈X. a x) / n"
    by (subst sum_distrib_left [symmetric]) auto
  finally show ?thesis .
qed (auto simp: n_def)


subsection ‹The Equality Case›

text ‹
  Next, we show that weighted arithmetic and geometric mean are equal if and only if all the 
  $a_i$ are equal.

  We first prove the more difficult direction as a lemmas and again first assume positivity
  of all $a_i$ and $w_i$ and will relax this somewhat later.
›
lemma weighted_arithmetic_geometric_mean_eq_iff_pos:
  fixes a w :: "'a ⇒ real"
  assumes "finite X"
  assumes pos1: "⋀x. x ∈ X ⟹ a x > 0"
  assumes pos2: "⋀x. x ∈ X ⟹ w x > 0"
  assumes sum_weights: "(∑x∈X. w x) = 1"
  assumes eq: "(∏x∈X. a x powr w x) = (∑x∈X. w x * a x)"
  shows   "∀x∈X. ∀y∈X. a x = a y"
proof -
  note nonneg1 = less_imp_le[OF pos1]
  note nonneg2 = less_imp_le[OF pos2]
  define A where "A = (∑x∈X. w x * a x)"
  define r where "r = (λx. a x / A - 1)"
  from sum_weights have "X ≠ {}" by auto
  hence "A ≠ 0"
    unfolding A_def using nonneg1 nonneg2 pos1 pos2 ‹finite X›
    by (subst sum_nonneg_eq_0_iff) force+
  moreover from nonneg1 nonneg2 have "A ≥ 0"
    by (auto simp: A_def intro!: sum_nonneg)
  ultimately have "A > 0" by simp

  have r_ge: "r x ≥ -1" if x: "x ∈ X" for x
    using ‹A > 0› pos1[OF x] by (auto simp: r_def field_simps)

  have "(∏x∈X. (1 + r x) powr w x) = (∏x∈X. (a x / A) powr w x)"
    by (simp add: r_def)
  also have "… = (∏x∈X. a x powr w x) / (∏x∈X. A powr w x)"
    unfolding prod_dividef [symmetric]
    using assms pos2 ‹A > 0› by (intro prod.cong powr_divide) (auto intro: less_imp_le)
  also have "(∏x∈X. A powr w x) = exp ((∑x∈X. w x) * ln A)"
    using ‹A > 0› and ‹finite X› by (simp add: powr_def exp_sum sum_distrib_right)
  also have "(∑x∈X. w x) = 1" by fact
  also have "exp (1 * ln A) = A"
    using ‹A > 0› by simp
  finally have lhs: "(∏x∈X. (1 + r x) powr w x) = (∏x∈X. a x powr w x) / A" .

  have "(∏x∈X. exp (w x * r x)) = exp (∑x∈X. w x * r x)"
    using ‹finite X› by (simp add: exp_sum)
  also have "(∑x∈X. w x * r x) = (∑x∈X. a x * w x) / A - 1"
    using ‹A > 0› by (simp add: r_def algebra_simps sum_subtractf sum_divide_distrib sum_weights)
  also have "(∑x∈X. a x * w x) / A = 1"
    using ‹A > 0› by (simp add: A_def mult.commute)
  finally have rhs: "(∏x∈X. exp (w x * r x)) = 1" by simp

  have "a x = A" if x: "x ∈ X" for x
  proof -
    have "(1 + r x) powr w x = exp (w x * r x)"
    proof (rule prod_ge_pointwise_le_imp_pointwise_eq
             [where f = "λx. (1 + r x) powr w x" and g = "λx. exp (w x * r x)"])
      show "(1 + r x) powr w x ≤ exp (w x * r x)" if x: "x ∈ X" for x
      proof -
        have "1 + r x ≤ exp (r x)"
          by (rule exp_ge_add_one_self)
        hence "(1 + r x) powr w x ≤ exp (r x) powr w x"
          using nonneg1[of x] nonneg2[of x] x ‹A > 0›
          by (intro powr_mono2) (auto simp: r_def field_simps)
        also have "… = exp (w x * r x)"
          by (simp add: powr_def)
        finally show "(1 + r x) powr w x ≤ exp (w x * r x)" .
      qed
    next
      show "(∏x∈X. (1 + r x) powr w x) ≥ (∏x∈X. exp (w x * r x))"
      proof -
        have "(∏x∈X. (1 + r x) powr w x) = (∏x∈X. a x powr w x) / A"
          by (fact lhs)
        also have "… = 1"
          using ‹A ≠ 0› by (simp add: eq A_def)
        also have "… = (∏x∈X. exp (w x * r x))"
          by (simp add: rhs)
        finally show ?thesis by simp
      qed
    qed (use x ‹finite X› in auto)

    also have "exp (w x * r x) = exp (r x) powr w x"
      by (simp add: powr_def)
    finally have "1 + r x = exp (r x)"
      using x pos2[of x] r_ge[of x] by (subst (asm) powr_left_real_eq_iff) auto
    hence "r x = 0"
      using exp_real_eq_one_plus_iff[of "r x"] by auto
    hence "a x = A"
      using ‹A > 0› by (simp add: r_def field_simps)
    thus ?thesis
      by (simp add: )
  qed
  thus "∀x∈X. ∀y∈X. a x = a y"
    by auto
qed

text ‹
  We can now show the full theorem and relax the positivity condition on the $a_i$ to
  non-negativity. This is possible because if some $a_i$ is zero and the two means
  coincide, then the product is obviously 0, but the sum can only be 0 if \<^term>‹all›
  the $a_i$ are 0.
›
theorem weighted_arithmetic_geometric_mean_eq_iff:
  fixes a w :: "'a ⇒ real"
  assumes "finite X"
  assumes nonneg1: "⋀x. x ∈ X ⟹ a x ≥ 0"
  assumes pos2:    "⋀x. x ∈ X ⟹ w x > 0"
  assumes sum_weights: "(∑x∈X. w x) = 1"
  shows   "(∏x∈X. a x powr w x) = (∑x∈X. w x * a x) ⟷ X ≠ {} ∧ (∀x∈X. ∀y∈X. a x = a y)"
proof
  assume *: "X ≠ {} ∧ (∀x∈X. ∀y∈X. a x = a y)"
  from * have "X ≠ {}"
    by blast

  from * obtain c where c: "⋀x. x ∈ X ⟹ a x = c" "c ≥ 0"
  proof (cases "X = {}")
    case False
    then obtain x where "x ∈ X" by blast
    thus ?thesis using * that[of "a x"] nonneg1[of x] by metis
  next
    case True
    thus ?thesis
      using that[of 1] by auto
  qed

  have "(∏x∈X. a x powr w x) = (∏x∈X. c powr w x)"
    by (simp add: c)
  also have "… = c"
    using assms c ‹X ≠ {}› by (cases "c = 0") (auto simp: powr_sum_distrib_real_right)
  also have "… = (∑x∈X. w x * a x)"
    using sum_weights by (simp add: c(1) flip: sum_distrib_left sum_distrib_right)
  finally show "(∏x∈X. a x powr w x) = (∑x∈X. w x * a x)" .
next
  assume *: "(∏x∈X. a x powr w x) = (∑x∈X. w x * a x)"
  have "X ≠ {}"
    using * by auto
  moreover have "(∀x∈X. ∀y∈X. a x = a y)"
  proof (cases "∃x∈X. a x = 0")
    case False
    with nonneg1 have pos1: "∀x∈X. a x > 0"
      by force
    thus ?thesis
      using weighted_arithmetic_geometric_mean_eq_iff_pos[of X a w] assms *
      by blast
  next
    case True
    hence "(∏x∈X. a x powr w x) = 0"
      using assms by auto
    with * have "(∑x∈X. w x * a x) = 0"
      by auto
    also have "?this ⟷ (∀x∈X. w x * a x = 0)"
      using assms by (intro sum_nonneg_eq_0_iff mult_nonneg_nonneg) (auto intro: less_imp_le)
    finally have "(∀x∈X. a x = 0)"
      using pos2 by force
    thus ?thesis
      by auto
  qed
  ultimately show "X ≠ {} ∧ (∀x∈X. ∀y∈X. a x = a y)"
    by blast
qed

text ‹
  Again, we derive a version for the unweighted arithmetic/geometric mean.
›
corollary arithmetic_geometric_mean_eq_iff:
  fixes a :: "'a ⇒ real"
  assumes "finite X"
  defines "n ≡ card X"
  assumes nonneg: "⋀x. x ∈ X ⟹ a x ≥ 0"
  shows   "root n (∏x∈X. a x) = (∑x∈X. a x) / n ⟷ (∀x∈X. ∀y∈X. a x = a y)"
proof (cases "X = {}")
  case False
  with assms have "n > 0"
    by auto
  have "(∏x∈X. a x powr (1 / n)) = (∑x∈X. (1 / n) * a x) ⟷
          X ≠ {} ∧ (∀x∈X. ∀y∈X. a x = a y)"
    using assms False by (intro weighted_arithmetic_geometric_mean_eq_iff) auto
  also have "(∏x∈X. a x powr (1 / n)) = (∏x∈X. a x) powr (1 / n)"
    using nonneg by (subst powr_sum_distrib_real_left) auto
  also have "… = root n (∏x∈X. a x)"
    using ‹n > 0› nonneg by (subst root_powr_inverse') (auto simp: prod_nonneg)
  also have "(∑x∈X. (1 / n) * a x) = (∑x∈X. a x) / n"
    by (subst sum_distrib_left [symmetric]) auto
  finally show ?thesis using False by auto
qed (auto simp: n_def)


subsection ‹The Binary Version›

text ‹
  For convenience, we also derive versions for only two numbers:
›
corollary weighted_arithmetic_geometric_mean_binary:
  fixes w1 w2 x1 x2 :: real
  assumes "x1 ≥ 0" "x2 ≥ 0" "w1 ≥ 0" "w2 ≥ 0" "w1 + w2 = 1"
  shows   "x1 powr w1 * x2 powr w2 ≤ w1 * x1 + w2 * x2"
proof -
  let ?a = "λb. if b then x1 else x2"
  let ?w = "λb. if b then w1 else w2"
  from assms have "(∏x∈UNIV. ?a x powr ?w x) ≤ (∑x∈UNIV. ?w x * ?a x)" 
    by (intro weighted_arithmetic_geometric_mean) (auto simp add: UNIV_bool)
  thus ?thesis by (simp add: UNIV_bool add_ac mult_ac)
qed

corollary weighted_arithmetic_geometric_mean_eq_iff_binary:
  fixes w1 w2 x1 x2 :: real
  assumes "x1 ≥ 0" "x2 ≥ 0" "w1 > 0" "w2 > 0" "w1 + w2 = 1"
  shows   "x1 powr w1 * x2 powr w2 = w1 * x1 + w2 * x2 ⟷ x1 = x2"
proof -
  let ?a = "λb. if b then x1 else x2"
  let ?w = "λb. if b then w1 else w2"
  from assms have "(∏x∈UNIV. ?a x powr ?w x) = (∑x∈UNIV. ?w x * ?a x)
                      ⟷ (UNIV :: bool set) ≠ {} ∧ (∀x∈UNIV. ∀y∈UNIV. ?a x = ?a y)" 
    by (intro weighted_arithmetic_geometric_mean_eq_iff) (auto simp add: UNIV_bool)
  thus ?thesis by (auto simp: UNIV_bool add_ac mult_ac)
qed

corollary arithmetic_geometric_mean_binary:
  fixes x1 x2 :: real
  assumes "x1 ≥ 0" "x2 ≥ 0"
  shows   "sqrt (x1 * x2) ≤ (x1 + x2) / 2"
  using weighted_arithmetic_geometric_mean_binary[of x1 x2 "1/2" "1/2"] assms
  by (simp add: powr_half_sqrt field_simps real_sqrt_mult)

corollary arithmetic_geometric_mean_eq_iff_binary:
  fixes x1 x2 :: real
  assumes "x1 ≥ 0" "x2 ≥ 0"
  shows   "sqrt (x1 * x2) = (x1 + x2) / 2 ⟷ x1 = x2"
  using weighted_arithmetic_geometric_mean_eq_iff_binary[of x1 x2 "1/2" "1/2"] assms
  by (simp add: powr_half_sqrt field_simps real_sqrt_mult)

end