feat: fitting sphere to data case study

CvxLean/Command/Solve/Float/Coeffs.lean

@@ -3,6 +3,19 @@ import CvxLean.Lib.Cones.All
 import CvxLean.Command.Solve.Float.ProblemData
 import CvxLean.Command.Solve.Float.RealToFloat
 
+
+/-!
+# Extract coefficients from problem to generate problem data
+
+TODO
+
+## TODO
+
+* This is probably a big source of inefficency for the `solve` command. We should come up with
+  a better way to extract the numerical values from the Lean expressions.
+* A first step is to not `unrollVectors` and turn thos expressions into floats directly.
+-/
+
 namespace CvxLean
 
 open Lean Meta Elab Tactic
@@ -173,13 +186,24 @@ unsafe def unrollVectors (constraints : Expr) : MetaM (Array Expr) := do
           res := res.push (← mkAppM ``Real.expCone #[ai, bi, ci])
     -- Vector second-order cone.
     | .app (.app (.app (.app (.app (.const ``Real.Vec.soCone _)
-        exprN@(.app (.const ``Fin _) n)) (.app (.const ``Fin _) m)) finTypeN) t) X =>
+        exprN@(.app (.const ``Fin _) _n)) (.app (.const ``Fin _) m)) finTypeN) t) X =>
         let m : Nat ← evalExpr Nat (mkConst ``Nat) m
         for i in [:m] do
           let idxExpr ← mkFinIdxExpr i m
           let ti := mkApp t idxExpr
           let Xi := mkApp X idxExpr
           res := res.push (mkAppN (mkConst ``Real.soCone) #[exprN, finTypeN, ti, Xi])
+    -- Vector rotated second-order cone.
+    -- Vector second-order cone.
+    | .app (.app (.app (.app (.app (.app (.const ``Real.Vec.rotatedSoCone _)
+        exprN@(.app (.const ``Fin _) _n)) (.app (.const ``Fin _) m)) finTypeN) v) w) X =>
+        let m : Nat ← evalExpr Nat (mkConst ``Nat) m
+        for i in [:m] do
+          let idxExpr ← mkFinIdxExpr i m
+          let vi := mkApp v idxExpr
+          let wi := mkApp w idxExpr
+          let Xi := mkApp X idxExpr
+          res := res.push (mkAppN (mkConst ``Real.rotatedSoCone) #[exprN, finTypeN, vi, wi, Xi])
     | _ =>
         res := res.push c
 
@@ -254,7 +278,10 @@ unsafe def determineCoeffsFromExpr (minExpr : Meta.MinimizationExpr) :
     let mut idx := 0
     for c in cs do
       trace[Meta.debug] "Coeffs going through constraint {c}."
+      let mut isTrivial := false
       match Expr.consumeMData c with
+      | .const ``True _ => do
+          isTrivial := true
       | .app (.const ``Real.zeroCone _) e => do
           let e ← realToFloat e
           let res ← determineScalarCoeffsAux e p floatDomain
@@ -331,7 +358,8 @@ unsafe def determineCoeffsFromExpr (minExpr : Meta.MinimizationExpr) :
               idx := idx + 1
       | _ => throwError "No match: {c}."
       -- New group, add idx.
-      sections := sections.push idx
+      if !isTrivial then
+        sections := sections.push idx
     return (data, sections)
 
   let (objectiveDataA, objectiveDataB) := objectiveData

CvxLean/Command/Solve/Float/RealToFloat.lean

@@ -24,15 +24,25 @@ partial def realToFloat (e : Expr) : MetaM Expr := do
   for translation in translations do
     let (mvars, _, pattern) ← lambdaMetaTelescope translation.real
     if ← isDefEq pattern e then
-    -- TODO: Search for conditions.
+      -- TODO: Search for conditions.
       let args ← mvars.mapM instantiateMVars
       return mkAppNBeta translation.float args
     else
       trace[Meta.debug] "`real-to-float` error: no match for \n{pattern} \n{e}"
   match e with
   | Expr.app a b => return mkApp (← realToFloat a) (← realToFloat b)
-  | Expr.lam n ty b d => return mkLambda n d (← realToFloat ty) (← realToFloat b)
-  | Expr.forallE n ty b d => return mkForall n d (← realToFloat ty) (← realToFloat b)
+  | Expr.lam n ty b d => do
+      withLocalDecl n d (← realToFloat ty) fun fvar => do
+        let b := b.instantiate1 fvar
+        let bF ← realToFloat b
+        mkLambdaFVars #[fvar] bF
+      -- return mkLambda n d (← realToFloat ty) (← realToFloat b)
+  | Expr.forallE n ty b d => do
+      withLocalDecl n d (← realToFloat ty) fun fvar => do
+        let b := b.instantiate1 fvar
+        let bF ← realToFloat b
+        mkForallFVars #[fvar] bF
+      -- return mkForall n d (← realToFloat ty) (← realToFloat b)
   | Expr.mdata m e => return mkMData m (← realToFloat e)
   | Expr.letE n ty t b _ => return mkLet n (← realToFloat ty) (← realToFloat t) (← realToFloat b)
   | Expr.proj typeName idx struct => return mkProj typeName idx (← realToFloat struct)
@@ -183,12 +193,18 @@ addRealToFloat (i) : @instHDiv Real i :=
 addRealToFloat (i) : @HPow.hPow Real Nat Real i :=
   fun f n => Float.pow f (Float.ofNat n)
 
-addRealToFloat (i) : @HPow.hPow Real Real Real i :=
+addRealToFloat (i) : @HPow.hPow.{0, 0, 0} Real Real Real i :=
   fun f n => Float.pow f n
 
-addRealToFloat (i) : @instHPow Real i :=
+addRealToFloat (β) (i) : @instHPow Real β i :=
   @HPow.mk Float Float Float Float.pow
 
+addRealToFloat (n) (i) : @HPow.hPow (Fin n → Real) Real (Fin n → Real) i :=
+  fun (x : Fin n → Float) (p : Float) (i : Fin n) => Float.pow (x i) p
+
+addRealToFloat (n) (β) (i) : @instHPow (Fin n → Real) β i :=
+  @HPow.mk (Fin n → Float) Float (Fin n → Float) (fun x p i => Float.pow (x i) p)
+
 addRealToFloat (i) : @LE.le Real i :=
   Float.le
 
@@ -210,6 +226,12 @@ addRealToFloat : @Real.sqrt :=
 addRealToFloat : @Real.log :=
   Float.log
 
+def Float.norm {n : ℕ} (x : Fin n → Float) : Float :=
+  Float.sqrt (Vec.Computable.sum (fun i => (Float.pow (x i) 2)))
+
+addRealToFloat (n) (i) : @Norm.norm.{0} (Fin n → ℝ) i :=
+  @Float.norm n
+
 addRealToFloat (i) : @OfScientific.ofScientific Real i :=
   Float.ofScientific
 
@@ -260,9 +282,12 @@ addRealToFloat (n) (i) (hn) : @Vec.sum.{0} ℝ i (Fin n) hn :=
 addRealToFloat (n) (i) (hn) : @Matrix.sum (Fin n) Real hn i :=
   @Matrix.Computable.sum n
 
-addRealToFloat (n) : @Vec.cumsum n :=
+addRealToFloat (n) (i) : @Vec.cumsum.{0} ℝ i n :=
   @Vec.Computable.cumsum n
 
+addRealToFloat : @Vec.norm :=
+  @Vec.Computable.norm
+
 addRealToFloat (n) (i1) (i2) (i3) : @Matrix.dotProduct (Fin n) ℝ i1 i2 i3 :=
   @Matrix.Computable.dotProduct n
 

CvxLean/Examples/All.lean

@@ -1,5 +1,5 @@
 import CvxLean.Examples.CovarianceEstimation
 import CvxLean.Examples.FittingSphere
 import CvxLean.Examples.HypersonicShapeDesign
-import CvxLean.Examples.OptimalVehicleSpeed
 import CvxLean.Examples.TrussDesign
+import CvxLean.Examples.VehicleSpeedScheduling

CvxLean/Examples/FittingSphere.lean

@@ -2,9 +2,68 @@ import CvxLean
 
 noncomputable section
 
-namespace FittingSphere
+open CvxLean Minimization Real BigOperators Matrix Finset
+
+section LeastSquares
+
+def leastSquares {n : ℕ} (a : Fin n → ℝ) :=
+  optimization (x : ℝ)
+    minimize (∑ i, ((a i - x) ^ 2) : ℝ)
+
+@[reducible]
+def mean {n : ℕ} (a : Fin n → ℝ) : ℝ := (1 / n) * ∑ i, (a i)
+
+/-- It is useful to rewrite the sum of squares in the following way to prove
+`leastSquares_optimal_eq_mean`, following Marty Cohen's answer in
+https://math.stackexchange.com/questions/2554243. -/
+lemma leastSquares_alt_objFun {n : ℕ} (hn : 0 < n) (a : Fin n → ℝ) (x : ℝ) :
+  (∑ i, ((a i - x) ^ 2)) = n * ((x - mean a) ^ 2 + (mean (a ^ 2) - (mean a) ^ 2)) := by
+  calc
+  -- 1) Σ (aᵢ - x)² = Σ (aᵢ² - 2aᵢx + x²)
+  _ = ∑ i, ((a i) ^ 2 - 2 * (a i) * x + (x ^ 2)) := by
+    congr; funext i; simp; ring
+  -- 2) ... = Σ aᵢ² - 2xΣ aᵢ + nx²
+  _ = ∑ i, ((a i) ^ 2) - 2 * x * ∑ i, (a i) + n * (x ^ 2) := by
+    rw [sum_add_distrib, sum_sub_distrib, ← sum_mul, ← mul_sum]; simp [sum_const]; ring
+  -- 3) ... = n{a²} - 2xn{a} + nx²
+  _ = n * mean (a ^ 2) - 2 * x * n * mean a + n * (x ^ 2) := by
+    simp [mean]; field_simp; ring
+  -- 4) ... = n((x - {a})² + ({a²} - {a}²))
+  _ = n * ((x - mean a) ^ 2 + (mean (a ^ 2) - (mean a) ^ 2)) := by
+    simp [mean]; field_simp; ring
+
+/-- Key result about least squares: `x* = mean a`. -/
+lemma leastSquares_optimal_eq_mean {n : ℕ} (hn : 0 < n) (a : Fin n → ℝ) (x : ℝ)
+  (h : (leastSquares a).optimal x) : x = mean a := by
+  simp [optimal, feasible, leastSquares] at h
+  replace h : ∀ y, (x - mean a) ^ 2 ≤ (y - mean a) ^ 2  := by
+    intros y
+    have hy := h y
+    have h_rw_x := leastSquares_alt_objFun hn a x
+    have h_rw_y := leastSquares_alt_objFun hn a y
+    simp only [rpow_two] at h_rw_x h_rw_y ⊢
+    rwa [h_rw_x, h_rw_y, mul_le_mul_left (by positivity), add_le_add_iff_right] at hy
+  have hmean := h (mean a)
+  simp at hmean
+  have hz := le_antisymm hmean (sq_nonneg _)
+  rwa [sq_eq_zero_iff, sub_eq_zero] at hz
+
+def Vec.leastSquares {n : ℕ} (a : Fin n → ℝ) :=
+  optimization (x : ℝ)
+    minimize (Vec.sum ((a - Vec.const n x) ^ 2) : ℝ)
+
+/-- Same as `leastSquares_optimal_eq_mean` in vector notation. -/
+lemma vec_leastSquares_optimal_eq_mean {n : ℕ} (hn : 0 < n) (a : Fin n → ℝ) (x : ℝ)
+  (h : (Vec.leastSquares a).optimal x) : x = mean a := by
+  apply leastSquares_optimal_eq_mean hn a
+  simp [Vec.leastSquares, leastSquares, optimal, feasible] at h ⊢
+  intros y
+  simp only [Vec.sum, Pi.pow_apply, Pi.sub_apply, Vec.const, rpow_two] at h
+  exact h y
 
-open CvxLean Minimization Real BigOperators Matrix
+end LeastSquares
+
+namespace FittingSphere
 
 -- Dimension.
 variable (n : ℕ)
@@ -19,7 +78,144 @@ def fittingSphere :=
   optimization (c : Fin n → ℝ) (r : ℝ)
     minimize (∑ i, (‖(x i) - c‖ ^ 2 - r ^ 2) ^ 2 : ℝ)
     subject to
-      _ : True
+      h₁ : 0 < r
+
+instance : ChangeOfVariables fun (ct : (Fin n → ℝ) × ℝ) => (ct.1, sqrt (ct.2 + ‖ct.1‖ ^ 2)) :=
+  { inv := fun (c, r) => (c, r ^ 2 - ‖c‖ ^ 2),
+    condition := fun (_, t) => 0 ≤ t,
+    property := fun ⟨c, t⟩ h => by simp [sqrt_sq h] }
+
+set_option trace.Meta.debug true
+
+equivalence' eqv/fittingSphereT (n m : ℕ) (x : Fin m → Fin n → ℝ) : fittingSphere n m x := by
+  -- Change of variables.
+  equivalence_step =>
+    apply ChangeOfVariables.toEquivalence
+      (fun (ct : (Fin n → ℝ) × ℝ) => (ct.1, sqrt (ct.2 + ‖ct.1‖ ^ 2)))
+    . rintro _ h; exact le_of_lt h
+  rename_vars [c, t]
+  -- Clean up.
+  conv_constr h₁ => dsimp
+  conv_obj => dsimp
+  -- Rewrite objective.
+  equivalence_step =>
+    apply Equivalence.rewrite_objFun
+      (g := fun (ct : (Fin n → ℝ) × ℝ) =>
+        Vec.sum (((Vec.norm x) ^ 2 - 2 * (Matrix.mulVec x ct.1) - Vec.const m ct.2) ^ 2))
+    . rintro ⟨c, t⟩ h
+      dsimp at h ⊢; simp [Vec.sum, Vec.norm, Vec.const]
+      congr; funext i; congr 1;
+      rw [@norm_sub_sq ℝ (Fin n → ℝ) _ (PiLp.normedAddCommGroup _ _) (PiLp.innerProductSpace _)]
+      rw [sq_sqrt (rpow_two _ ▸ le_of_lt (sqrt_pos.mp <| h))]
+      simp [mulVec, inner, dotProduct]
+  rename_vars [c, t]
+
+#print fittingSphereT
+
+relaxation rel/fittingSphereConvex (n m : ℕ) (x : Fin m → Fin n → ℝ) : fittingSphereT n m x := by
+  relaxation_step =>
+    apply Relaxation.weaken_constraint (cs' := fun _ => True)
+    . rintro ⟨c, t⟩ _; trivial
+
+/-- If the squared error is zero, then `aᵢ = x`. -/
+lemma vec_squared_norm_error_eq_zero_iff {n m : ℕ} (a : Fin m → Fin n → ℝ) (x : Fin n → ℝ) :
+    ∑ i, ‖a i - x‖ ^ 2 = 0 ↔ ∀ i, a i = x := by
+  simp [rpow_two]
+  rw [sum_eq_zero_iff_of_nonneg (fun _ _ => sq_nonneg _)]
+  constructor
+  . intros h i
+    have hi := h i (by simp)
+    rw [sq_eq_zero_iff, @norm_eq_zero _ (PiLp.normedAddCommGroup _ _).toNormedAddGroup] at hi
+    rwa [sub_eq_zero] at hi
+  . intros h i _
+    rw [sq_eq_zero_iff, @norm_eq_zero _ (PiLp.normedAddCommGroup _ _).toNormedAddGroup, sub_eq_zero]
+    exact h i
+
+/-- This tells us that solving the relaxed problem is sufficient for optimal points if the solution
+is non-trivial. -/
+lemma optimal_relaxed_implies_optimal (hm : 0 < m) (c : Fin n → ℝ) (t : ℝ)
+  (h_nontrivial : x ≠ Vec.const m c)
+  (h_opt : (fittingSphereConvex n m x).optimal (c, t)) : (fittingSphereT n m x).optimal (c, t) := by
+  simp [fittingSphereT, fittingSphereConvex, optimal, feasible] at h_opt ⊢
+  constructor
+  . let a := Vec.norm x ^ 2 - 2 * mulVec x c
+    have h_ls : optimal (Vec.leastSquares a) t := by
+      refine ⟨trivial, ?_⟩
+      intros y _
+      simp [objFun, Vec.leastSquares]
+      exact h_opt c y
+    -- Apply key result about least squares to `a` and `t`.
+    have ht_eq := vec_leastSquares_optimal_eq_mean hm a t h_ls
+    have hc2_eq : ‖c‖ ^ 2 = (1 / m) * ∑ i : Fin m, ‖c‖ ^ 2 := by
+      simp [sum_const]
+      field_simp; ring
+    have ht : t + ‖c‖ ^ 2 = (1 / m) * ∑ i, ‖(x i) - c‖ ^ 2 := by
+      rw [ht_eq]; dsimp [mean]
+      rw [hc2_eq, mul_sum, mul_sum, mul_sum, ← sum_add_distrib]
+      congr; funext i; rw [← mul_add]
+      congr; simp [Vec.norm]
+      rw [@norm_sub_sq ℝ (Fin n → ℝ) _ (PiLp.normedAddCommGroup _ _) (PiLp.innerProductSpace _)]
+      congr
+    -- We use the result to establish that `t + ‖c‖ ^ 2` is non-negative.
+    have h_tc2_nonneg : 0 ≤ t + ‖c‖ ^ 2 := by
+      rw [ht]
+      apply mul_nonneg (by norm_num)
+      apply sum_nonneg
+      intros i _
+      rw [rpow_two]
+      exact sq_nonneg _
+    cases (lt_or_eq_of_le h_tc2_nonneg) with
+    | inl h_tc2_lt_zero =>
+        -- If it is positive, we are done.
+        convert h_tc2_lt_zero; simp
+    | inr h_tc2_eq_zero =>
+        -- Otherwise, it contradicts the non-triviality assumption.
+        exfalso
+        rw [ht, zero_eq_mul] at h_tc2_eq_zero
+        rcases h_tc2_eq_zero with (hc | h_sum_eq_zero)
+        . simp at hc; linarith
+        rw [vec_squared_norm_error_eq_zero_iff] at h_sum_eq_zero
+        apply h_nontrivial
+        funext i
+        exact h_sum_eq_zero i
+  . intros c' x' _
+    exact h_opt c' x'
+
+#print fittingSphereConvex
+
+-- We proceed to solve the problem on a concrete example.
+-- https://github.com/cvxgrp/cvxbook_additional_exercises/blob/main/python/sphere_fit_data.py
+
+@[optimization_param]
+def nₚ := 2
+
+@[optimization_param]
+def mₚ := 10
+
+@[optimization_param]
+def xₚ : Fin mₚ → Fin nₚ → ℝ := Matrix.transpose <| ![
+  ![1.824183228637652032e+00, 1.349093690455489103e+00, 6.966316403935147727e-01,
+    7.599387854623529392e-01, 2.388321695850912363e+00, 8.651370608981923116e-01,
+    1.863922545015865406e+00, 7.099743941474848663e-01, 6.005484882320809570e-01,
+    4.561429569892232472e-01],
+  ![-9.644136284187876385e-01, 1.069547315003422927e+00, 6.733229334437943470e-01,
+    7.788072961810316164e-01, -9.467465278344706636e-01, -8.591303443863639311e-01,
+    1.279527420871080956e+00, 5.314829019311283487e-01, 6.975676079749143499e-02,
+    -4.641873429414754559e-01]]
+
+-- We use the `solve` command on the data above.
+
+set_option maxHeartbeats 1000000
+
+solve fittingSphereConvex nₚ mₚ xₚ
+
+-- Finally, we recover the solution to the original problem.
+
+def sol := eqv.backward_map nₚ mₚ xₚ.float fittingSphereConvex.solution
+
+#print eqv.backward_map
+
+#eval sol -- (![1.664863, 0.031932], 1.159033)
 
 end FittingSphere