import Mathlib.Algebra.Order.BigOperators.Group.Finset
import Mathlib.Algebra.BigOperators.Ring.Finset
import Mathlib.Data.Real.Basic
import Mathlib.Tactic.FieldSimp
import Mathlib.Tactic.Ring
import Mathlib.Tactic.Linarith

/-!
Broyden (1965), DOI 10.1090/S0025-5718-1965-0198670-6.
The direct rank-one secant update and its least-change characterization.
Rectangular real matrices are allowed. No convergence or invertibility claim.
-/
namespace PaperToLean.Broyden
noncomputable section
open scoped BigOperators
variable {m n : ℕ}
abbrev Vec (n : ℕ) := Fin n → ℝ
abbrev Mat (m n : ℕ) := Fin m → Vec n
def dot (x y : Vec n) : ℝ := ∑ i, x i * y i
def distance (C B : Mat m n) : ℝ := ∑ j, ∑ i, (C j i - B j i) ^ 2
def update (B : Mat m n) (s : Vec n) (y : Vec m) : Mat m n :=
  fun j i => B j i + (y j - dot (B j) s) / dot s s * s i

theorem action (B : Mat m n) (s v : Vec n) (y : Vec m) (j : Fin m) :
    dot (update B s y j) v = dot (B j) v +
      (y j - dot (B j) s) / dot s s * dot s v := by
  simp only [dot, update, add_mul, Finset.sum_add_distrib, Finset.mul_sum]
  congr 1
  apply Finset.sum_congr rfl
  intro i _
  ring

theorem secant (B : Mat m n) (s : Vec n) (y : Vec m) (hs : dot s s ≠ 0) :
    ∀ j, dot (update B s y j) s = y j := by
  intro j
  rw [action]
  field_simp
  ring

theorem unchanged_on_orthogonal (B : Mat m n) (s v : Vec n) (y : Vec m)
    (hv : dot s v = 0) : ∀ j, dot (update B s y j) v = dot (B j) v := by
  intro j
  simp [action, hv]

theorem distance_nonneg (C B : Mat m n) : 0 ≤ distance C B := by
  exact Finset.sum_nonneg (fun _ _ => Finset.sum_nonneg (fun _ _ => sq_nonneg _))

theorem distance_zero_iff (C B : Mat m n) : distance C B = 0 ↔ C = B := by
  constructor
  · intro h
    have hj := (Finset.sum_eq_zero_iff_of_nonneg
      (fun j (_ : j ∈ Finset.univ) =>
        Finset.sum_nonneg (fun i (_ : i ∈ Finset.univ) => sq_nonneg (C j i - B j i)))).mp h
    funext j i
    have hi := (Finset.sum_eq_zero_iff_of_nonneg
      (fun i (_ : i ∈ Finset.univ) => sq_nonneg (C j i - B j i))).mp
        (hj j (Finset.mem_univ j))
    exact sub_eq_zero.mp (sq_eq_zero_iff.mp (hi i (Finset.mem_univ i)))
  · rintro rfl
    simp [distance]

/-- The exact gap, stronger than the usual least-change inequality. -/
theorem pythagorean (B C : Mat m n) (s : Vec n) (y : Vec m)
    (hs : dot s s ≠ 0) (hC : ∀ j, dot (C j) s = y j) :
    distance C B = distance (update B s y) B + distance C (update B s y) := by
  have row (b c s : Vec n) (t : ℝ) :
      (∑ i, (c i - b i) ^ 2) = (∑ i, (t * s i) ^ 2) +
      (∑ i, (c i - (b i + t * s i)) ^ 2) +
      2 * t * (dot c s - dot b s - t * dot s s) := by
    simp only [dot, mul_sub, Finset.mul_sum, ← Finset.sum_sub_distrib,
      ← Finset.sum_add_distrib]
    apply Finset.sum_congr rfl
    intro i _
    ring
  unfold distance
  rw [← Finset.sum_add_distrib]
  apply Finset.sum_congr rfl
  intro j _
  have hzero : dot (C j) s - dot (B j) s -
      (y j - dot (B j) s) / dot s s * dot s s = 0 := by
    rw [hC j]
    field_simp
    ring
  have h := row (B j) (C j) s ((y j - dot (B j) s) / dot s s)
  simpa only [hzero, mul_zero, add_zero, update, add_sub_cancel_left] using h

theorem least_change (B C : Mat m n) (s : Vec n) (y : Vec m)
    (hs : dot s s ≠ 0) (hC : ∀ j, dot (C j) s = y j) :
    distance (update B s y) B ≤ distance C B := by
  rw [pythagorean B C s y hs hC]
  exact le_add_of_nonneg_right (distance_nonneg C (update B s y))

theorem unique_minimizer (B C : Mat m n) (s : Vec n) (y : Vec m)
    (hs : dot s s ≠ 0) (hC : ∀ j, dot (C j) s = y j)
    (hmin : distance C B ≤ distance (update B s y) B) : C = update B s y := by
  have hp := pythagorean B C s y hs hC
  have hn := distance_nonneg C (update B s y)
  apply (distance_zero_iff C (update B s y)).mp
  linarith

#print axioms action
#print axioms secant
#print axioms unchanged_on_orthogonal
#print axioms distance_nonneg
#print axioms distance_zero_iff
#print axioms pythagorean
#print axioms least_change
#print axioms unique_minimizer
end
end PaperToLean.Broyden