Discussion: New tactic: `compute_asymptotics`

Stream: general

[Jireh Loreaux]

A discussion thread for #announce > New tactic: `compute_asymptotics`

[Jireh Loreaux]

@Vasilii Nesterov very exciting! A few questions:

1. Is it extensible? Can we teach it to support other functions? If so, how? Does the user write meta code, or just tag lemmas?
2. Have you looked for things throughout Mathlib you can simplify with this tactic? Such a list could be a nice selling point (e.g., provided as a PR targeting your current one).

[Vasilii Nesterov]

1. At this point the tactic isn’t extensible, but I plan to make it so. I think it’s hard to avoid user meta code completely here, though the difficulty of this meta code depends on the properties of the function being added. One of the next steps on my TODO list is to add support for trigonometric functions and the gamma function. I’ve already done sine and cosine, and those took about 400 lines, roughly 50 of which are meta.
2. Yes, I’m thinking about this, but I haven’t searched yet. I guess this tactic is most useful in combinatorics, where it often comes down to computing some asymptotic.

[Jireh Loreaux]

One important reason for the tactic to be extensible is so that you can use it on-the-fly, so I'm glad you have plans in this direction. Otherwise you have to develop all the necessary theory of your special functions so that you can build it into the tactic, and only after you import all of this can you use it.

[Bolton Bailey]

Can or could the tactic be used to obtain concrete bounds for points after which one function exceeds another?

For example, one step of the proof of Bertrand's postulate is to show that

$$\frac{4^n}{2n} > (2n)^{\sqrt{2n} -1} 4^{2n/3}$$

For sufficiently large $n > N$. It's not too important what the $N$ is, but it has to be concrete so that the theorem can be checked manually for smaller values.

[David Ledvinka]

I noticed you wrote TODO: Support different domains including ℕ. If you do get around to this I have a good case where it should be usable to simplify a nasty computation. This comes from the Brownian Motion project so it is not currently in Mathlib but it will be eventually:

import Mathlib

open MeasureTheory ProbabilityTheory Filter
open scoped ENNReal NNReal Topology Asymptotics

variable {T Ω E : Type*} {mΩ : MeasurableSpace Ω}
  {X : T → Ω → E} {c : ℝ≥0∞} {d p q : ℝ} {M β : ℝ≥0} {P : Measure Ω}

variable [PseudoEMetricSpace T] [PseudoEMetricSpace E] [MeasurableSpace E] [BorelSpace E]

noncomputable
def Cp (d p q : ℝ) : ℝ≥0∞ :=
  max (1 / ((2 ^ ((q - d) / p)) - 1) ^ p) (1 / (2 ^ (q - d) - 1))

noncomputable
-- the `max 0 ...` in the blueprint is performed by `ENNReal.ofReal` here
def constL (T : Type*) [PseudoEMetricSpace T] (c : ℝ≥0∞) (d p q β : ℝ) : ℝ≥0∞ :=
  2 ^ (2 * p + 5 * q + 1) * c * (EMetric.diam (.univ : Set T) + 1) ^ (q - d)
  * ∑' (k : ℕ), 2 ^ (k * (β * p - (q - d)))
      * (4 ^ d * (ENNReal.ofReal (Real.logb 2 c.toReal + (k + 2) * d)) ^ q + Cp d p q)

lemma constL_lt_top (hT : EMetric.diam (Set.univ : Set T) < ∞)
    (hc : c ≠ ∞) (hd_pos : 0 < d) (hp_pos : 0 < p) (hdq_lt : d < q) (hβ_lt : β < (q - d) / p) :
    constL T c d p q β < ∞ := by
  have hq_pos : 0 < q := lt_trans hd_pos hdq_lt
  have hC : Cp d p q ≠ ⊤ := by
    apply max_ne_top <;> apply ENNReal.div_ne_top (by norm_num)
    · apply ne_of_gt
      refine ENNReal.rpow_pos ?_ (by finiteness)
      exact tsub_pos_of_lt (ENNReal.one_lt_rpow (by norm_num) (by bound))
    · exact ne_of_gt <| tsub_pos_of_lt (ENNReal.one_lt_rpow (by norm_num) (by bound))
  have hC_pos : 0 < Cp d p q := by
    apply lt_max_of_lt_right (ENNReal.div_pos (by norm_num) (by finiteness))
  apply ENNReal.mul_lt_top (by finiteness)
  conv =>
    enter [1, 1, _]
    rw [← (ENNReal.ofReal_toReal_eq_iff (a := _ * _)).mpr (by finiteness),
      ENNReal.ofReal_eq_coe_nnreal (by positivity)]
  rw [lt_top_iff_ne_top, ENNReal.tsum_coe_ne_top_iff_summable_coe]
  apply summable_of_ratio_test_tendsto_lt_one (l := 2 ^ (β * p - (q - d)))
  · apply Real.rpow_lt_one_of_one_lt_of_neg (by norm_num)
    simp [← lt_div_iff₀, hp_pos, hβ_lt]
  · filter_upwards with k
    apply ne_of_gt
    simp only [ENNReal.toReal_mul, NNReal.coe_mk]
    apply mul_pos <;> refine ENNReal.toReal_pos_iff.mpr ⟨?_, by finiteness⟩
    · exact ENNReal.rpow_pos (by norm_num) (by norm_num)
    · positivity
  -- Goal becomes `Tendsto` of a function from `ℕ → ℝ`
  simp only [Nat.cast_add, Nat.cast_one, ENNReal.toReal_mul, NNReal.coe_mk, norm_mul,
    Real.norm_eq_abs, ENNReal.abs_toReal, ← div_mul_div_comm, add_mul (b := (1 : ℝ)), one_mul]
  conv => enter [1, _, 2, 1]; rw [ENNReal.toReal_add (by finiteness) (by finiteness)]
  conv => enter [1, _, 2, 2]; rw [ENNReal.toReal_add (by finiteness) (by finiteness)]
  simp only [← ENNReal.toReal_rpow, ENNReal.toReal_ofNat, Nat.ofNat_pos, Real.rpow_add,
    ENNReal.toReal_mul]
  conv => enter [1, _, 1]; rw [mul_div_cancel_left₀ _ (by positivity)]
  conv => enter [3, 1]; rw [← mul_one (_ ^ _)]
  apply Tendsto.const_mul
  conv => enter [1]; change (fun n ↦ _) / (fun n ↦ _)
  rw [← Asymptotics.isEquivalent_iff_tendsto_one]; swap
  · filter_upwards with _
    apply ne_of_gt
    refine lt_of_le_of_lt ?_ <| (add_lt_add_left (ENNReal.toReal_pos (by positivity) hC)) _
    positivity
  refine Asymptotics.IsEquivalent.add_add_of_nonneg
    (by intro _; positivity) (by intro _; positivity) ?_ .refl
  apply Asymptotics.IsEquivalent.mul .refl
  apply Asymptotics.IsEquivalent.rpow (by intro _; positivity)
  have h (k : ℕ) : ∀ᶠ (n : ℕ) in atTop, 0 ≤ Real.logb 2 c.toReal + (n + k + 2) * d := by
    obtain ⟨n₀, hn₀⟩ := exists_nat_gt (- Real.logb 2 c.toReal / d)
    rw [eventually_atTop]
    use n₀
    intro n hn
    grw [hn, add_mul, add_mul, ← le_of_lt ((div_lt_iff₀ hd_pos).mp hn₀), add_assoc, ← add_assoc]
    simp only [add_neg_cancel, zero_add]
    positivity
  apply Asymptotics.IsEquivalent.congr_right; swap
  · filter_upwards [h 0] with n h_nonneg
    rw [← add_zero (n : ℝ), ← Nat.cast_zero, ENNReal.toReal_ofReal h_nonneg]
  apply Asymptotics.IsEquivalent.congr_left; swap
  · filter_upwards [h 1] with n h_nonneg
    rw [← Nat.cast_one, ENNReal.toReal_ofReal h_nonneg]
  apply Asymptotics.IsEquivalent.const_add_of_norm_tendsto_atTop; swap
  · apply Tendsto.comp tendsto_norm_atTop_atTop
    apply tendsto_atTop_add_const_left
    rw [tendsto_mul_const_atTop_of_pos hd_pos]
    repeat apply tendsto_atTop_add_const_right
    exact tendsto_natCast_atTop_iff.mpr tendsto_id
  refine (Asymptotics.IsEquivalent.const_add_of_norm_tendsto_atTop ?_ ?_).symm; swap
  · apply Tendsto.comp tendsto_norm_atTop_atTop
    rw [tendsto_mul_const_atTop_of_pos hd_pos]
    repeat apply tendsto_atTop_add_const_right
    exact tendsto_natCast_atTop_iff.mpr tendsto_id
  refine Asymptotics.IsEquivalent.mul ?_ .refl
  simp only [add_assoc]
  apply Asymptotics.IsEquivalent.add_const_of_norm_tendsto_atTop; swap
  · apply Tendsto.comp tendsto_norm_atTop_atTop
    apply tendsto_atTop_add_const_right
    exact tendsto_natCast_atTop_iff.mpr tendsto_id
  refine (Asymptotics.IsEquivalent.add_const_of_norm_tendsto_atTop .refl ?_).symm
  exact Tendsto.comp tendsto_norm_atTop_atTop (tendsto_natCast_atTop_iff.mpr tendsto_id)

Even if it still needs some preprocessing (simplifications after the goal becomes Tendsto I think you should be able to skip the repetitive computation at the end)

[David Ledvinka]

Another nice thing would be if you could do something similar to @Tomas Skrivan s data_synth tactic #29432 where it could also compute the limit (without you stating what it is) in the case of Tendsto. Not sure if that is feasible with the current algorithm. Perhaps this isn't a priority but Id like to see the day where Lean + Mathlib could compete with Wolfram Alpha/Mathematica as a research tool.

[Vasilii Nesterov]

Yes, different domains are low-hanging fruit. For example, if f : ℕ → ℝ, then we need to "cast back" f to ℝ → ℝ (formally, construct a function g : ℝ → ℝ such that f = g ∘ (Nat.cast : ℕ → ℝ)), compute its limit, and then use docs#Filter.Tendsto.comp.

[Vasilii Nesterov]

And yes, the tactic is able to compute the limit, not just prove that the given one is correct. I'll make an elaborator compute_limit% f source that computes the limit target of f at filter source and returns the proof of Tendsto f source target, so it can be used like have := compute_limit% f source inside the proof.

[Vasilii Nesterov]

Bolton Bailey said:

Can or could the tactic be used to obtain concrete bounds for points after which one function exceeds another?

For example, one step of the proof of Bertrand's postulate is to show that

$$\frac{4^n}{2n} > (2n)^{\sqrt{2n} -1} 4^{2n/3}$$

For sufficiently large $n > N$. It's not too important what the $N$ is, but it has to be concrete so that the theorem can be checked manually for smaller values.

No, I think it's fundamentally out of scope for the current approach of the tactic.

[Vasilii Nesterov]

It can, of course, prove that the RHS is little-o of the LHS, but it gives only the existence of such an N, not a concrete number.

[Manuel Eberl]

I'd been wondering why nobody ported this to Lean yet, considering the seemingly infinite number of people working with Lean. After all, the initial prototype for real_asymp took me only maybe two or three weeks or so back in 2016. In any case, nice work! I'm sure it will be very useful to the Lean community.

Please let me know if you ever figure out how to properly deal with things like $f( x + \exp(-x))$ where $f$ is a function that doesn't have nice properties w.r.t. addition like $\sin$, $\exp$ do (e.g. $\arctan$, $\Gamma$). They show how to handle this in the book by Richardson, Salvy, et al., but they don't really explain why it works and I was never able to make it rigorous. So real_asymp doesn't support this.

By the way, while real_asymp is technically extensible, nobody uses that feature. I think I added support for $\Gamma$ and maybe $\text{erf}$ that way (not sure if I ever put that in the AFP), but it's a bit messy and definitely requires user-written code.

[Ruben Van de Velde]

For low values of infinity :)

[Alex Kontorovich]

This is very cool, congrats!!! I think it can shorten a lot of the PNT+ code, much of it is fighting against asymptotic estimates in ways bound can't handle.

[Vasilii Nesterov]

David Ledvinka said:

I noticed you wrote TODO: Support different domains including ℕ. If you do get around to this I have a good case where it should be usable to simplify a nasty computation. This comes from the Brownian Motion project so it is not currently in Mathlib but it will be eventually:

import Mathlib

open MeasureTheory ProbabilityTheory Filter
open scoped ENNReal NNReal Topology Asymptotics

variable {T Ω E : Type*} {mΩ : MeasurableSpace Ω}
  {X : T → Ω → E} {c : ℝ≥0∞} {d p q : ℝ} {M β : ℝ≥0} {P : Measure Ω}

variable [PseudoEMetricSpace T] [PseudoEMetricSpace E] [MeasurableSpace E] [BorelSpace E]

noncomputable
def Cp (d p q : ℝ) : ℝ≥0∞ :=
  max (1 / ((2 ^ ((q - d) / p)) - 1) ^ p) (1 / (2 ^ (q - d) - 1))

noncomputable
-- the `max 0 ...` in the blueprint is performed by `ENNReal.ofReal` here
def constL (T : Type*) [PseudoEMetricSpace T] (c : ℝ≥0∞) (d p q β : ℝ) : ℝ≥0∞ :=
  2 ^ (2 * p + 5 * q + 1) * c * (EMetric.diam (.univ : Set T) + 1) ^ (q - d)
  * ∑' (k : ℕ), 2 ^ (k * (β * p - (q - d)))
      * (4 ^ d * (ENNReal.ofReal (Real.logb 2 c.toReal + (k + 2) * d)) ^ q + Cp d p q)

lemma constL_lt_top (hT : EMetric.diam (Set.univ : Set T) < ∞)
    (hc : c ≠ ∞) (hd_pos : 0 < d) (hp_pos : 0 < p) (hdq_lt : d < q) (hβ_lt : β < (q - d) / p) :
    constL T c d p q β < ∞ := by
  have hq_pos : 0 < q := lt_trans hd_pos hdq_lt
  have hC : Cp d p q ≠ ⊤ := by
    apply max_ne_top <;> apply ENNReal.div_ne_top (by norm_num)
    · apply ne_of_gt
      refine ENNReal.rpow_pos ?_ (by finiteness)
      exact tsub_pos_of_lt (ENNReal.one_lt_rpow (by norm_num) (by bound))
    · exact ne_of_gt <| tsub_pos_of_lt (ENNReal.one_lt_rpow (by norm_num) (by bound))
  have hC_pos : 0 < Cp d p q := by
    apply lt_max_of_lt_right (ENNReal.div_pos (by norm_num) (by finiteness))
  apply ENNReal.mul_lt_top (by finiteness)
  conv =>
    enter [1, 1, _]
    rw [← (ENNReal.ofReal_toReal_eq_iff (a := _ * _)).mpr (by finiteness),
      ENNReal.ofReal_eq_coe_nnreal (by positivity)]
  rw [lt_top_iff_ne_top, ENNReal.tsum_coe_ne_top_iff_summable_coe]
  apply summable_of_ratio_test_tendsto_lt_one (l := 2 ^ (β * p - (q - d)))
  · apply Real.rpow_lt_one_of_one_lt_of_neg (by norm_num)
    simp [← lt_div_iff₀, hp_pos, hβ_lt]
  · filter_upwards with k
    apply ne_of_gt
    simp only [ENNReal.toReal_mul, NNReal.coe_mk]
    apply mul_pos <;> refine ENNReal.toReal_pos_iff.mpr ⟨?_, by finiteness⟩
    · exact ENNReal.rpow_pos (by norm_num) (by norm_num)
    · positivity
  -- Goal becomes `Tendsto` of a function from `ℕ → ℝ`
  simp only [Nat.cast_add, Nat.cast_one, ENNReal.toReal_mul, NNReal.coe_mk, norm_mul,
    Real.norm_eq_abs, ENNReal.abs_toReal, ← div_mul_div_comm, add_mul (b := (1 : ℝ)), one_mul]
  conv => enter [1, _, 2, 1]; rw [ENNReal.toReal_add (by finiteness) (by finiteness)]
  conv => enter [1, _, 2, 2]; rw [ENNReal.toReal_add (by finiteness) (by finiteness)]
  simp only [← ENNReal.toReal_rpow, ENNReal.toReal_ofNat, Nat.ofNat_pos, Real.rpow_add,
    ENNReal.toReal_mul]
  conv => enter [1, _, 1]; rw [mul_div_cancel_left₀ _ (by positivity)]
  conv => enter [3, 1]; rw [← mul_one (_ ^ _)]
  apply Tendsto.const_mul
  conv => enter [1]; change (fun n ↦ _) / (fun n ↦ _)
  rw [← Asymptotics.isEquivalent_iff_tendsto_one]; swap
  · filter_upwards with _
    apply ne_of_gt
    refine lt_of_le_of_lt ?_ <| (add_lt_add_left (ENNReal.toReal_pos (by positivity) hC)) _
    positivity
  refine Asymptotics.IsEquivalent.add_add_of_nonneg
    (by intro _; positivity) (by intro _; positivity) ?_ .refl
  apply Asymptotics.IsEquivalent.mul .refl
  apply Asymptotics.IsEquivalent.rpow (by intro _; positivity)
  have h (k : ℕ) : ∀ᶠ (n : ℕ) in atTop, 0 ≤ Real.logb 2 c.toReal + (n + k + 2) * d := by
    obtain ⟨n₀, hn₀⟩ := exists_nat_gt (- Real.logb 2 c.toReal / d)
    rw [eventually_atTop]
    use n₀
    intro n hn
    grw [hn, add_mul, add_mul, ← le_of_lt ((div_lt_iff₀ hd_pos).mp hn₀), add_assoc, ← add_assoc]
    simp only [add_neg_cancel, zero_add]
    positivity
  apply Asymptotics.IsEquivalent.congr_right; swap
  · filter_upwards [h 0] with n h_nonneg
    rw [← add_zero (n : ℝ), ← Nat.cast_zero, ENNReal.toReal_ofReal h_nonneg]
  apply Asymptotics.IsEquivalent.congr_left; swap
  · filter_upwards [h 1] with n h_nonneg
    rw [← Nat.cast_one, ENNReal.toReal_ofReal h_nonneg]
  apply Asymptotics.IsEquivalent.const_add_of_norm_tendsto_atTop; swap
  · apply Tendsto.comp tendsto_norm_atTop_atTop
    apply tendsto_atTop_add_const_left
    rw [tendsto_mul_const_atTop_of_pos hd_pos]
    repeat apply tendsto_atTop_add_const_right
    exact tendsto_natCast_atTop_iff.mpr tendsto_id
  refine (Asymptotics.IsEquivalent.const_add_of_norm_tendsto_atTop ?_ ?_).symm; swap
  · apply Tendsto.comp tendsto_norm_atTop_atTop
    rw [tendsto_mul_const_atTop_of_pos hd_pos]
    repeat apply tendsto_atTop_add_const_right
    exact tendsto_natCast_atTop_iff.mpr tendsto_id
  refine Asymptotics.IsEquivalent.mul ?_ .refl
  simp only [add_assoc]
  apply Asymptotics.IsEquivalent.add_const_of_norm_tendsto_atTop; swap
  · apply Tendsto.comp tendsto_norm_atTop_atTop
    apply tendsto_atTop_add_const_right
    exact tendsto_natCast_atTop_iff.mpr tendsto_id
  refine (Asymptotics.IsEquivalent.add_const_of_norm_tendsto_atTop .refl ?_).symm
  exact Tendsto.comp tendsto_norm_atTop_atTop (tendsto_natCast_atTop_iff.mpr tendsto_id)

Even if it still needs some preprocessing (simplifications after the goal becomes Tendsto I think you should be able to skip the repetitive computation at the end)

I've added support for other domains and codomains of functions and tried your example. Some manual work with casts is still needed, but it got a little better.

import Mathlib

open MeasureTheory ProbabilityTheory Filter
open scoped ENNReal NNReal Topology Asymptotics

universe u

noncomputable
def Cp (d p q : ℝ) : ℝ≥0∞ :=
  max (1 / ((2 ^ ((q - d) / p)) - 1) ^ p) (1 / (2 ^ (q - d) - 1))

set_option maxHeartbeats 0 in
theorem constL_lt_top.extracted_1_10 {T : Type u} {c : ℝ≥0∞} {d p q : ℝ} {β : ℝ≥0} [inst : PseudoEMetricSpace T]
  (hT : EMetric.diam (@Set.univ T) < ⊤) (hc : c ≠ ⊤) (hd_pos : 0 < d) (hp_pos : 0 < p) (hdq_lt : d < q)
  (hβ_lt : ↑β < (q - d) / p) (hq_pos : 0 < q) (hC : Cp d p q ≠ ⊤) (hC_pos : 0 < Cp d p q) :
  Tendsto
    (fun n =>
      ‖NNReal.toReal ⟨(2 ^ (↑(n + 1) * (↑β * p - (q - d))) *
                  (4 ^ d * ENNReal.ofReal (Real.logb 2 c.toReal + (↑(n + 1) + 2) * d) ^ q + Cp d p q)).toReal,
              sorry⟩‖ /
        ‖NNReal.toReal ⟨(2 ^ (↑n * (↑β * p - (q - d))) *
                  (4 ^ d * ENNReal.ofReal (Real.logb 2 c.toReal + (↑n + 2) * d) ^ q + Cp d p q)).toReal,
              sorry⟩‖)
    atTop (𝓝 (2 ^ (↑β * p - (q - d)))) := by
  -- fighting with casts begins
  generalize Cp d p q = C at *
  cases C with
  | top => simp at hC
  | coe C =>
  simp
  simp [← ENNReal.toReal_rpow]
  have h1 (n) : (4 ^ d * ENNReal.ofReal (Real.logb 2 c.toReal + (n + 1 + 2) * d) ^ q + ↑C).toReal =
      4 ^ d * (Real.logb 2 c.toReal + (n + 1 + 2) * d) ^ q + C := by
    sorry
  have h2 : ∀ n, (4 ^ d * ENNReal.ofReal (Real.logb 2 c.toReal + (n + 2) * d) ^ q + ↑C).toReal =
      4 ^ d * (Real.logb 2 c.toReal + (n + 2) * d) ^ q + C := by
    sorry
  simp_rw [h1, h2]
  -- fighting with casts ends
  -- The tactic asks us to find the sign of `Real.log 2 * (β * p - (q - d))`
  have : 0 < Real.log 2 * (β * p - (q - d)) := by
    sorry
  compute_asymptotics
  -- the goal now is to prove
  -- 4 ^ d * d ^ q *
  --   ((d ^ q)⁻¹ * (4 ^ d)⁻¹ *
  --     (Real.exp
  --         (-(Real.log 2 * (↑β * p - (q - d)) *
  --             (Real.log 2 * (↑β * p - (q - d)) / (Real.log 2 * (↑β * p - (q - d)))))))⁻¹) =
  -- 2 ^ (↑β * p - (q - d))
  -- (LHS = the limit the tactic found, RHS = the desired limit)
  sorry