Every degree-at-most-four scalar real polynomial minimal graph is affine, in every dimension

/- Polynomial minimal graphs of degree at most four are affine.
   Submission body for Math30Catalog.source28.
   The challenge supplies its imports and the surrounding Bounty namespace. -/


/- Component: Compatibility -/
section ComponentCompatibility

/-
The extensionality proofs below are adapted from Mathlib/Algebra/MvPolynomial/Funext.lean.
Copyright (c) 2020 Johan Commelin. All rights reserved.
Released under Apache 2.0 license as described in Mathlib's LICENSE.
Authors: Johan Commelin.

The evaluation compatibility proofs are adapted from
Mathlib/Algebra/MvPolynomial/Polynomial.lean.
Copyright (c) 2023 Kim Morrison. All rights reserved.
Released under Apache 2.0 license as described in Mathlib's LICENSE.
Authors: Kim Morrison.

The homogeneous derivative and Euler proofs are adapted from
Mathlib/RingTheory/MvPolynomial/EulerIdentity.lean.
Copyright (c) 2024 Junyan Xu. All rights reserved.
Released under Apache 2.0 license as described in Mathlib's LICENSE.
Authors: Junyan Xu.
-/

namespace QuarticCompatibility

open MvPolynomial

theorem polynomial_eval_eval₂ {R S σ : Type*} [CommSemiring R] [CommSemiring S]
    {x : S} (f : R →+* Polynomial S) (g : σ → Polynomial S) (p : MvPolynomial σ R) :
    Polynomial.eval x (eval₂ f g p) =
      eval₂ ((Polynomial.evalRingHom x).comp f) (fun s => Polynomial.eval x (g s)) p := by
  apply induction_on p
  · simp
  · intro p q hp hq
    simp [hp, hq]
  · intro p n hp
    simp [hp]

theorem eval_polynomial_eval_finSuccEquiv {R : Type*} {n : ℕ} {x : Fin n → R}
    [CommSemiring R] (f : MvPolynomial (Fin (n + 1)) R) (q : MvPolynomial (Fin n) R) :
    (eval x) (Polynomial.eval q (finSuccEquiv R n f)) = eval (Fin.cases (eval x q) x) f := by
  simp only [finSuccEquiv_apply, coe_eval₂Hom, polynomial_eval_eval₂, eval_eval₂]
  conv in RingHom.comp _ _ =>
    refine @RingHom.ext _ _ _ _ _ (RingHom.id _) fun r => ?_
    simp
  simp only [eval₂_id]
  congr
  funext i
  refine Fin.cases (by simp) (by simp) i

section Extensionality

variable {R : Type*} [CommRing R] [IsDomain R]

theorem polynomial_funext_fin {n : ℕ} {p : MvPolynomial (Fin n) R}
    (s : Fin n → Set R) (hs : ∀ i, (s i).Infinite)
    (h : ∀ x ∈ Set.pi .univ s, eval x p = 0) : p = 0 := by
  induction n with
  | zero =>
    have hp := h 0 finZeroElim
    rw [eq_C_of_isEmpty p, eval_C] at hp
    rw [eq_C_of_isEmpty p, hp, map_zero]
  | succ n ih =>
    apply (finSuccEquiv R n).injective
    rw [map_zero]
    apply Polynomial.eq_zero_of_infinite_isRoot
    apply ((hs 0).image (C_injective ..).injOn).mono
    rintro _ ⟨r, hr, rfl⟩
    refine ih (s ·.succ) (fun _ ↦ hs _) fun x hx ↦ ?_
    rw [eval_polynomial_eval_finSuccEquiv]
    exact h _ fun i _ ↦ i.cases (by simpa [eval_C] using hr) (by simpa using hx)

theorem polynomial_funext_set {σ : Type*} {p q : MvPolynomial σ R}
    (s : σ → Set R) (hs : ∀ i, (s i).Infinite)
    (h : ∀ x ∈ Set.pi .univ s, eval x p = eval x q) : p = q := by
  suffices ∀ p, (∀ x ∈ Set.pi .univ s, eval x p = 0) → p = 0 by
    rw [← sub_eq_zero, this (p - q)]
    intro x hx
    simp_rw [map_sub, h x hx, sub_self]
  intro p h
  obtain ⟨n, f, hf, p, rfl⟩ := exists_fin_rename p
  suffices p = 0 by rw [this, map_zero]
  refine polynomial_funext_fin (s ∘ f) (fun _ ↦ hs _) fun x hx ↦ ?_
  choose g hg using fun i ↦ (hs i).nonempty
  convert! h (Function.extend f x g) fun i _ ↦ ?_
  · simp only [eval, eval₂Hom_rename, Function.extend_comp hf]
  obtain ⟨i, rfl⟩ | nex := em (∃ x, f x = i)
  · rw [hf.extend_apply]; exact hx _ ⟨⟩
  · simp_rw [Function.extend, dite_eq_right nex, hg]

theorem polynomial_funext [Infinite R] {σ : Type*} {p q : MvPolynomial σ R}
    (h : ∀ x : σ → R, eval x p = eval x q) : p = q :=
  polynomial_funext_set _ (fun _ ↦ Set.infinite_univ) fun _ _ ↦ h _

end Extensionality

section Homogeneity

open Finsupp

variable {R σ M : Type*} [CommSemiring R] {φ : MvPolynomial σ R}

theorem weighted_homogeneous_pderiv [AddCancelCommMonoid M] {w : σ → M}
    {n n' : M} {i : σ} (h : φ.IsWeightedHomogeneous w n) (h' : n' + w i = n) :
    (pderiv i φ).IsWeightedHomogeneous w n' := by
  rw [← mem_weightedHomogeneousSubmodule, weightedHomogeneousSubmodule_eq_finsupp_supported,
    AddMonoidAlgebra.supported_eq_span_single] at h
  refine Submodule.span_induction ?_ ?_ (fun p q _ _ hp hq ↦ ?_) (fun r p _ h ↦ ?_) h
  · rintro _ ⟨m, hm, rfl⟩
    simp_rw [single_eq_monomial, pderiv_monomial, one_mul]
    by_cases hi : m i = 0
    · rw [hi, Nat.cast_zero, monomial_zero]; apply isWeightedHomogeneous_zero
    convert! isWeightedHomogeneous_monomial ..
    rw [← add_right_cancel_iff (a := w i), h', ← hm, weight_sub_single_add hi]
  · rw [map_zero]; apply isWeightedHomogeneous_zero
  · rw [map_add]; exact hp.add hq
  · rw [(pderiv i).map_smul]; exact (weightedHomogeneousSubmodule ..).smul_mem _ h

theorem homogeneous_pderiv {n : ℕ} {i : σ} (h : φ.IsHomogeneous n) :
    (pderiv i φ).IsHomogeneous (n - 1) := by
  obtain _ | n := n
  · rw [← totalDegree_zero_iff_isHomogeneous, totalDegree_eq_zero_iff_eq_C] at h
    rw [h, pderiv_C]; apply isHomogeneous_zero
  · exact weighted_homogeneous_pderiv h rfl

variable [Fintype σ] {n : ℕ}

open Finset in
theorem weighted_sum_X_mul_pderiv {w : σ → ℕ} (h : φ.IsWeightedHomogeneous w n) :
    ∑ i : σ, w i • (X i * pderiv i φ) = n • φ := by
  rw [← mem_weightedHomogeneousSubmodule, weightedHomogeneousSubmodule_eq_finsupp_supported,
    AddMonoidAlgebra.supported_eq_span_single] at h
  refine Submodule.span_induction ?_ ?_ (fun p q _ _ hp hq ↦ ?_) (fun r p _ h ↦ ?_) h
  · rintro _ ⟨m, hm, rfl⟩
    simp_rw [single_eq_monomial, X_mul_pderiv_monomial, smul_smul, ← sum_smul, mul_comm (w _)]
    congr
    rwa [Set.mem_ofPred, weight_apply, sum_fintype] at hm
    intro; apply zero_smul
  · simp
  · simp_rw [map_add, left_distrib, smul_add, sum_add_distrib, hp, hq]
  · simp_rw [(pderiv _).map_smul, nsmul_eq_mul, mul_smul_comm, ← Finset.smul_sum,
      ← nsmul_eq_mul, h]

theorem sum_X_mul_pderiv (h : φ.IsHomogeneous n) :
    ∑ i : σ, X i * pderiv i φ = n • φ := by
  simp_rw [← weighted_sum_X_mul_pderiv h, Pi.one_apply, one_smul]

end Homogeneity

end QuarticCompatibility
end ComponentCompatibility


/- Component: MatrixAlgebra -/
section ComponentMatrixAlgebra

/-!
Algebraic invariance of the minimal-graph numerator under a matrix whose rows
are orthonormal.  The identities hold over any commutative ring, including a
multivariate polynomial ring.
-/

namespace QuarticMinimal.MatrixAlgebra

open Matrix

variable {R : Type*} [CommRing R]
variable {m n : Type*} [Fintype m] [Fintype n] [DecidableEq m]

theorem dot_transpose_mulVec (A : Matrix m n R) (hA : A * Aᵀ = 1) (v w : m → R) :
    (Aᵀ *ᵥ v) ⬝ᵥ (Aᵀ *ᵥ w) = v ⬝ᵥ w := by
  rw [dotProduct_transpose_mulVec, mulVec_mulVec, hA, one_mulVec, dotProduct_comm]

theorem trace_conjugate (A : Matrix m n R) (hA : A * Aᵀ = 1) (H : Matrix m m R) :
    trace (Aᵀ * H * A) = trace H := by
  rw [trace_mul_cycle, hA, one_mul]

theorem conjugate_mulVec (A : Matrix m n R) (hA : A * Aᵀ = 1)
    (H : Matrix m m R) (v : m → R) :
    (Aᵀ * H * A) *ᵥ (Aᵀ *ᵥ v) = Aᵀ *ᵥ (H *ᵥ v) := by
  calc
    (Aᵀ * H * A) *ᵥ (Aᵀ *ᵥ v) = (Aᵀ * H * (A * Aᵀ)) *ᵥ v := by
      rw [mulVec_mulVec, Matrix.mul_assoc]
    _ = (Aᵀ * H) *ᵥ v := by rw [hA, Matrix.mul_one]
    _ = Aᵀ *ᵥ (H *ᵥ v) := (mulVec_mulVec v Aᵀ H).symm

theorem contraction_conjugate (A : Matrix m n R) (hA : A * Aᵀ = 1)
    (H : Matrix m m R) (v : m → R) :
    (Aᵀ *ᵥ v) ⬝ᵥ ((Aᵀ * H * A) *ᵥ (Aᵀ *ᵥ v)) = v ⬝ᵥ (H *ᵥ v) := by
  rw [conjugate_mulVec A hA, dot_transpose_mulVec A hA]

theorem minimal_expression_invariant (A : Matrix m n R) (hA : A * Aᵀ = 1)
    (H : Matrix m m R) (v : m → R) :
    (1 + (Aᵀ *ᵥ v) ⬝ᵥ (Aᵀ *ᵥ v)) * trace (Aᵀ * H * A) -
      (Aᵀ *ᵥ v) ⬝ᵥ ((Aᵀ * H * A) *ᵥ (Aᵀ *ᵥ v)) =
    (1 + v ⬝ᵥ v) * trace H - v ⬝ᵥ (H *ᵥ v) := by
  rw [dot_transpose_mulVec A hA, trace_conjugate A hA, contraction_conjugate A hA]

omit [Fintype m] in
theorem map_mul_transpose {S : Type*} [CommRing S] (f : R →+* S)
    (A : Matrix m n R) (hA : A * Aᵀ = 1) :
    A.map f * (A.map f)ᵀ = 1 := by
  have hm := congrArg (fun M : Matrix m m R => M.map f) hA
  simpa [Matrix.transpose_map] using hm

omit [Fintype m] in
theorem algebraMap_mul_transpose [Algebra ℝ R]
    (A : Matrix m n ℝ) (hA : A * Aᵀ = 1) :
    A.map (algebraMap ℝ R) * (A.map (algebraMap ℝ R))ᵀ = 1 :=
  map_mul_transpose (algebraMap ℝ R) A hA

end QuarticMinimal.MatrixAlgebra
end ComponentMatrixAlgebra


/- Component: PolynomialBasics -/
section ComponentPolynomialBasics

noncomputable section
open scoped BigOperators
open MvPolynomial

namespace QuarticPolynomialBasics

variable {σ : Type*}

lemma exponent_decompose_one (m : σ →₀ ℕ) (h : 1 ≤ m.degree) :
    ∃ (i : σ) (r : σ →₀ ℕ), r + Finsupp.single i 1 = m := by
  classical
  have hm : m ≠ 0 := by intro hm; simp [hm] at h
  obtain ⟨i, hi⟩ : ∃ i, m i ≠ 0 := by
    by_contra hh
    push Not at hh
    exact hm (Finsupp.ext hh)
  exact ⟨i, m - Finsupp.single i 1, Finsupp.sub_add_single_one_cancel hi⟩

lemma exponent_decompose_two (m : σ →₀ ℕ) (h : 2 ≤ m.degree) :
    ∃ (i j : σ) (r : σ →₀ ℕ), r + Finsupp.single i 1 + Finsupp.single j 1 = m := by
  obtain ⟨j, s, hs⟩ := exponent_decompose_one m (by omega)
  have hsdeg : 1 ≤ s.degree := by
    have hh := congrArg Finsupp.degree hs
    simp only [map_add, Finsupp.degree_single] at hh
    omega
  obtain ⟨i, r, hr⟩ := exponent_decompose_one s hsdeg
  exact ⟨i, j, r, by rw [hr, hs]⟩

lemma coeff_second_partial (P : MvPolynomial σ ℝ) (i j : σ) (r : σ →₀ ℕ) :
    (pderiv i (pderiv j P)).coeff r =
      P.coeff (r + Finsupp.single i 1 + Finsupp.single j 1) *
        (((r + Finsupp.single i 1 : σ →₀ ℕ) j : ℝ) + 1) * (r i + 1) := by
  rw [coeff_pderiv, coeff_pderiv]

lemma second_partial_comm (P : MvPolynomial σ ℝ) (i j : σ) :
    pderiv i (pderiv j P) = pderiv j (pderiv i P) := by
  classical
  ext r
  rw [coeff_second_partial, coeff_second_partial]
  by_cases h : i = j
  · subst j; rfl
  · simp [h, Ne.symm h, add_right_comm, mul_right_comm]

lemma totalDegree_le_one_of_second_partials_zero (P : MvPolynomial σ ℝ)
    (h : ∀ i j, pderiv i (pderiv j P) = 0) : P.totalDegree ≤ 1 := by
  classical
  apply Finset.sup_le
  intro m hm
  change m.degree ≤ 1
  by_contra hh
  obtain ⟨i, j, r, hr⟩ := exponent_decompose_two m (by omega)
  have hc := congrArg (fun Q : MvPolynomial σ ℝ ↦ Q.coeff r) (h i j)
  rw [coeff_second_partial, hr, AddMonoidAlgebra.coeff_zero] at hc
  have hn : P.coeff m ≠ 0 := mem_support_iff.mp hm
  have ha : (((r + Finsupp.single i 1 : σ →₀ ℕ) j : ℕ) : ℝ) + 1 ≠ 0 := by positivity
  have hb : ((r i : ℕ) : ℝ) + 1 ≠ 0 := by positivity
  exact (mul_ne_zero (mul_ne_zero hn ha) hb) hc

lemma second_partials_zero_of_totalDegree_le_one (P : MvPolynomial σ ℝ)
    (h : P.totalDegree ≤ 1) (i j : σ) : pderiv i (pderiv j P) = 0 := by
  classical
  ext r
  rw [coeff_second_partial, AddMonoidAlgebra.coeff_zero]
  have hc : P.coeff (r + Finsupp.single i 1 + Finsupp.single j 1) = 0 := by
    apply coeff_eq_zero_of_totalDegree_lt
    change P.totalDegree < (r + Finsupp.single i 1 + Finsupp.single j 1).degree
    simp only [map_add, Finsupp.degree_single]
    omega
  simp [hc]

lemma totalDegree_le_one_iff_second_partials_zero (P : MvPolynomial σ ℝ) :
    P.totalDegree ≤ 1 ↔ ∀ i j, pderiv i (pderiv j P) = 0 :=
  ⟨second_partials_zero_of_totalDegree_le_one P,
    totalDegree_le_one_of_second_partials_zero P⟩

lemma totalDegree_pderiv_le_of_totalDegree_le (P : MvPolynomial σ ℝ)
    (i : σ) (d : ℕ) (h : P.totalDegree ≤ d + 1) :
    (pderiv i P).totalDegree ≤ d := by
  classical
  apply Finset.sup_le
  intro m hm
  change m.degree ≤ d
  have hc : P.coeff (m + Finsupp.single i 1) ≠ 0 := by
    have hh := mem_support_iff.mp hm
    rw [coeff_pderiv] at hh
    exact (mul_ne_zero_iff.mp hh).1
  have hd := le_totalDegree (mem_support_iff.mpr hc)
  change (m + Finsupp.single i 1).degree ≤ P.totalDegree at hd
  simp only [map_add, Finsupp.degree_single] at hd
  omega

lemma totalDegree_pderiv_le (P : MvPolynomial σ ℝ) (i : σ) :
    (pderiv i P).totalDegree ≤ P.totalDegree - 1 :=
  totalDegree_pderiv_le_of_totalDegree_le P i _ (by omega)

lemma pderiv_homogeneousComponent (P : MvPolynomial σ ℝ) (i : σ) (d : ℕ) :
    pderiv i (homogeneousComponent (d + 1) P) =
      homogeneousComponent d (pderiv i P) := by
  classical
  ext m
  simp only [coeff_pderiv, coeff_homogeneousComponent, map_add, Finsupp.degree_single]
  by_cases h : m.degree = d <;> simp [h]

lemma isHomogeneous_pderiv_succ (P : MvPolynomial σ ℝ) (i : σ) (d : ℕ)
    (h : P.IsHomogeneous (d + 1)) : (pderiv i P).IsHomogeneous d := by
  have hh := pderiv_homogeneousComponent P i d
  rw [homogeneousComponent_eq_self h] at hh
  rw [hh]
  exact homogeneousComponent_isHomogeneous _ _

lemma totalDegree_sub_homogeneousComponent_le (P : MvPolynomial σ ℝ)
    (d : ℕ) (h : P.totalDegree ≤ d) :
    (P - homogeneousComponent d P).totalDegree ≤ d - 1 := by
  classical
  apply Finset.sup_le
  intro m hm
  change m.degree ≤ d - 1
  have hc := mem_support_iff.mp hm
  simp only [coeff_sub, coeff_homogeneousComponent] at hc
  split_ifs at hc with hd
  · exact False.elim (hc (sub_self _))
  · simp only [sub_zero] at hc
    have hh := le_totalDegree (mem_support_iff.mpr hc)
    change m.degree ≤ P.totalDegree at hh
    omega

lemma top_homogeneousComponent_ne_zero (P : MvPolynomial σ ℝ)
    (h : 0 < P.totalDegree) : homogeneousComponent P.totalDegree P ≠ 0 := by
  intro hh
  have hdeg := totalDegree_sub_homogeneousComponent_le P P.totalDegree le_rfl
  rw [hh, sub_zero] at hdeg
  omega

lemma totalDegree_sub_top_homogeneousComponent_lt (P : MvPolynomial σ ℝ)
    (h : 0 < P.totalDegree) :
    (P - homogeneousComponent P.totalDegree P).totalDegree < P.totalDegree := by
  have hdeg := totalDegree_sub_homogeneousComponent_le P P.totalDegree le_rfl
  omega

lemma homogeneousComponent_mul_top (P Q : MvPolynomial σ ℝ) (d e : ℕ)
    (hP : P.totalDegree ≤ d) (hQ : Q.totalDegree ≤ e) :
    homogeneousComponent (d + e) (P * Q) =
      homogeneousComponent d P * homogeneousComponent e Q := by
  by_cases hd : d = 0
  · subst d
    have hp0 : P.totalDegree = 0 := Nat.eq_zero_of_le_zero hP
    rw [(totalDegree_eq_zero_iff_eq_C.mp hp0)]
    simp [homogeneousComponent_C_mul, homogeneousComponent_zero]
  by_cases he : e = 0
  · subst e
    have hq0 : Q.totalDegree = 0 := Nat.eq_zero_of_le_zero hQ
    rw [(totalDegree_eq_zero_iff_eq_C.mp hq0)]
    rw [mul_comm P (C (Q.coeff 0)), homogeneousComponent_C_mul]
    simp [homogeneousComponent_zero, mul_comm]
  have hdp := totalDegree_sub_homogeneousComponent_le P d hP
  have hdq := totalDegree_sub_homogeneousComponent_le Q e hQ
  have htp : (homogeneousComponent d P).totalDegree ≤ d :=
    (homogeneousComponent_isHomogeneous d P).totalDegree_le
  have hz1 : homogeneousComponent (d + e) ((P - homogeneousComponent d P) * Q) = 0 := by
    apply homogeneousComponent_eq_zero
    have hh := totalDegree_mul (P - homogeneousComponent d P) Q
    omega
  have hz2 : homogeneousComponent (d + e)
      (homogeneousComponent d P * (Q - homogeneousComponent e Q)) = 0 := by
    apply homogeneousComponent_eq_zero
    have hh := totalDegree_mul (homogeneousComponent d P) (Q - homogeneousComponent e Q)
    omega
  have heq : P * Q = homogeneousComponent d P * homogeneousComponent e Q +
      ((P - homogeneousComponent d P) * Q +
        homogeneousComponent d P * (Q - homogeneousComponent e Q)) := by ring
  rw [heq, map_add, map_add, hz1, hz2, add_zero, add_zero]
  exact homogeneousComponent_eq_self
    ((homogeneousComponent_isHomogeneous d P).mul (homogeneousComponent_isHomogeneous e Q))

/-- The polynomial numerator of the scalar minimal graph equation. This is
definitionally equal to the corresponding operator in the challenge. -/
def minimalGraphOperator {n : ℕ} (P : MvPolynomial (Fin n) ℝ) :
    MvPolynomial (Fin n) ℝ :=
  (1 + ∑ i : Fin n, (pderiv i P) ^ 2) *
      (∑ i : Fin n, pderiv i (pderiv i P)) -
    ∑ i : Fin n, ∑ j : Fin n,
      pderiv i P * pderiv j P * pderiv i (pderiv j P)

lemma minimalGraphOperator_zero_variables (P : MvPolynomial (Fin 0) ℝ) :
    minimalGraphOperator P = 0 := by
  simp [minimalGraphOperator]

lemma totalDegree_zero_variables (P : MvPolynomial (Fin 0) ℝ) :
    P.totalDegree = 0 := by
  rw [MvPolynomial.eq_C_of_isEmpty P, totalDegree_C]

lemma minimalGraphOperator_one_variable (P : MvPolynomial (Fin 1) ℝ) :
    minimalGraphOperator P = pderiv 0 (pderiv 0 P) := by
  simp only [minimalGraphOperator, Fin.sum_univ_one]
  ring

lemma affine_of_minimalGraphOperator_eq_zero_one_variable
    (P : MvPolynomial (Fin 1) ℝ) (h : minimalGraphOperator P = 0) :
    P.totalDegree ≤ 1 := by
  apply totalDegree_le_one_of_second_partials_zero
  intro i j
  have hi : i = 0 := Subsingleton.elim _ _
  have hj : j = 0 := Subsingleton.elim _ _
  subst i; subst j
  rwa [minimalGraphOperator_one_variable] at h

lemma minimalGraphOperator_eq_zero_of_totalDegree_le_one
    {n : ℕ} (P : MvPolynomial (Fin n) ℝ) (h : P.totalDegree ≤ 1) :
    minimalGraphOperator P = 0 := by
  simp [minimalGraphOperator, second_partials_zero_of_totalDegree_le_one P h]

def polynomialLaplacian {n : ℕ} (P : MvPolynomial (Fin n) ℝ) :
    MvPolynomial (Fin n) ℝ := ∑ i : Fin n, pderiv i (pderiv i P)

def cubicOperator {n : ℕ} (P : MvPolynomial (Fin n) ℝ) :
    MvPolynomial (Fin n) ℝ :=
  (∑ i : Fin n, (pderiv i P) ^ 2) * polynomialLaplacian P -
    ∑ i : Fin n, ∑ j : Fin n,
      pderiv i P * pderiv j P * pderiv i (pderiv j P)

lemma minimalGraphOperator_eq_laplacian_add_cubic {n : ℕ}
    (P : MvPolynomial (Fin n) ℝ) :
    minimalGraphOperator P = polynomialLaplacian P + cubicOperator P := by
  unfold minimalGraphOperator cubicOperator polynomialLaplacian
  ring

lemma totalDegree_second_partial_le (P : MvPolynomial σ ℝ) (d : ℕ)
    (h : P.totalDegree ≤ d + 2) (i j : σ) :
    (pderiv i (pderiv j P)).totalDegree ≤ d := by
  apply totalDegree_pderiv_le_of_totalDegree_le
  exact totalDegree_pderiv_le_of_totalDegree_le P j (d + 1) (by omega)

lemma homogeneousComponent_second_partial (P : MvPolynomial σ ℝ) (d : ℕ)
    (i j : σ) : homogeneousComponent d (pderiv i (pderiv j P)) =
      pderiv i (pderiv j (homogeneousComponent (d + 2) P)) := by
  rw [← pderiv_homogeneousComponent, ← pderiv_homogeneousComponent]

lemma homogeneousComponent_derivative_triple (P : MvPolynomial σ ℝ) (d : ℕ)
    (h : P.totalDegree ≤ d + 2) (i j k l : σ) :
    homogeneousComponent (3 * d + 2)
      (pderiv i P * pderiv j P * pderiv k (pderiv l P)) =
      pderiv i (homogeneousComponent (d + 2) P) *
        pderiv j (homogeneousComponent (d + 2) P) *
        pderiv k (pderiv l (homogeneousComponent (d + 2) P)) := by
  have hgi : (pderiv i P).totalDegree ≤ d + 1 :=
    totalDegree_pderiv_le_of_totalDegree_le P i (d + 1) (by omega)
  have hgj : (pderiv j P).totalDegree ≤ d + 1 :=
    totalDegree_pderiv_le_of_totalDegree_le P j (d + 1) (by omega)
  have hpair : (pderiv i P * pderiv j P).totalDegree ≤ (d + 1) + (d + 1) :=
    (totalDegree_mul _ _).trans (Nat.add_le_add hgi hgj)
  have hkl := totalDegree_second_partial_le P d h k l
  rw [show 3 * d + 2 = ((d + 1) + (d + 1)) + d by omega,
    homogeneousComponent_mul_top _ _ _ _ hpair hkl,
    homogeneousComponent_mul_top _ _ _ _ hgi hgj,
    ← pderiv_homogeneousComponent, ← pderiv_homogeneousComponent,
    homogeneousComponent_second_partial]

lemma homogeneousComponent_cubicOperator_top {n : ℕ}
    (P : MvPolynomial (Fin n) ℝ) (d : ℕ) (h : P.totalDegree ≤ d + 2) :
    homogeneousComponent (3 * d + 2) (cubicOperator P) =
      cubicOperator (homogeneousComponent (d + 2) P) := by
  unfold cubicOperator polynomialLaplacian
  simp only [Finset.sum_mul, Finset.mul_sum, map_sub, map_sum, pow_two]
  simp_rw [homogeneousComponent_derivative_triple P d h]

lemma homogeneousComponent_minimalGraphOperator_top {n : ℕ}
    (P : MvPolynomial (Fin n) ℝ) (d : ℕ) (h : P.totalDegree ≤ d + 2) :
    homogeneousComponent (3 * d + 2) (minimalGraphOperator P) =
      cubicOperator (homogeneousComponent (d + 2) P) := by
  rw [minimalGraphOperator_eq_laplacian_add_cubic, map_add,
    homogeneousComponent_cubicOperator_top P d h]
  have hlap : (polynomialLaplacian P).totalDegree ≤ d := by
    apply totalDegree_finsetSum_le
    intro i hi
    exact totalDegree_second_partial_le P d h i i
  have hz : homogeneousComponent (3 * d + 2) (polynomialLaplacian P) = 0 := by
    apply homogeneousComponent_eq_zero
    omega