Grace/Cumulus/LinearAlgebra/Projection.g17

Open in Workbench Download source

Imports

CauchySchwarz

Declarations

Source

1-- Orthogonal projection onto a line and its least-squares characterization2-- Part of Cumulus, the Grace17 mathematics library.3-- Copyright (c) 2026 Grant Passmore.4-- Released under the Apache License 2.0. See LICENSE in the repository root.5-- Authors: Grant Passmore67module Grace.Cumulus.LinearAlgebra.Projection8import CauchySchwarz910-- All vector statements concern the first n coordinates of real sequences.11-- The inverse is total; the projection laws below require a nonzero norm.12def projectionCoefficient (u v : RealSequence, n : Nat) : Real :=13  dotProduct(u, v, n) * realInvValue(squaredNorm(v, n))1415def lineProjection (u v : RealSequence, n : Nat) : RealSequence :=16  mulSequence(constantSequence(projectionCoefficient(u, v, n)), v)1718def lineDistanceSquared (u v : RealSequence, C : Real, n : Nat) : Real :=19  squaredNorm(subSequence(u, mulSequence(constantSequence(C), v)), n)2021theorem projectionCoefficient.real (u v : Set, n : Nat)22  requires sequenceU : RealSequence(u)23  requires sequenceV : RealSequence(v)24  shows Real(projectionCoefficient(u, v, n))25proof26  unfold projectionCoefficient27  typecheck2829theorem lineProjection.sequence (u v : Set, n : Nat)30  requires sequenceU : RealSequence(u)31  requires sequenceV : RealSequence(v)32  shows RealSequence(lineProjection(u, v, n))33proof34  have constants := constantSequence.sequence at projectionCoefficient(u, v, n) with *35  unfold lineProjection36  exact mulSequence.sequence with *3738theorem lineDistanceSquared.real (u v : Set, C : Real, n : Nat)39  requires sequenceU : RealSequence(u)40  requires sequenceV : RealSequence(v)41  shows Real(lineDistanceSquared(u, v, C, n))42proof43  have constants := constantSequence.sequence with *44  have scaled := mulSequence.sequence at constantSequence(C), v with *45  have difference := subSequence.sequence at u, mulSequence(constantSequence(C), v) with *46  unfold lineDistanceSquared47  exact squaredNorm.real with *4849theorem projectionCoefficient.normal_equation (u v : Set, n : Nat)50  requires sequenceU : RealSequence(u)51  requires sequenceV : RealSequence(v)52  requires nonzero : squaredNorm(v, n) ≠ (0:R)53  shows projectionCoefficient(u, v, n) * squaredNorm(v, n) = dotProduct(u, v, n)54proof55  have inverseLaw := realInvValue.mul_right with *56  unfold projectionCoefficient57  linear_combination dotProduct(u, v, n) * inverseLaw5859theorem lineProjection.orthogonal (u v : Set, n : Nat)60  requires sequenceU : RealSequence(u)61  requires sequenceV : RealSequence(v)62  requires nonzero : squaredNorm(v, n) ≠ (0:R)63  shows Orthogonal(subSequence(u, lineProjection(u, v, n)), v, n)64proof65  have projection := lineProjection.sequence at u, v, n with *66  have normal := projectionCoefficient.normal_equation at u, v, n with *67  unfold Orthogonal68  rw dotProduct.sub_left69  unfold lineProjection70  rw dotProduct.scale_left71  unfold squaredNorm at normal72  linear_combination -normal7374theorem lineDistanceSquared.expansion (u v : Set, C : Real, n : Nat)75  requires sequenceU : RealSequence(u)76  requires sequenceV : RealSequence(v)77  shows lineDistanceSquared(u, v, C, n) =78    squaredNorm(u, n) - (2:R) * C * dotProduct(u, v, n) + C * C * squaredNorm(v, n)79proof80  have constants := constantSequence.sequence with *81  have scaled := mulSequence.sequence at constantSequence(C), v with *82  unfold lineDistanceSquared83  simp [squaredNorm.sub, squaredNorm.scale, dotProduct.scale_right]84  ring8586theorem lineProjection.distance_decomposition (u v : Set, C : Real, n : Nat)87  requires sequenceU : RealSequence(u)88  requires sequenceV : RealSequence(v)89  requires nonzero : squaredNorm(v, n) ≠ (0:R)90  shows lineDistanceSquared(u, v, C, n) =91    lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) +92    (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)) * squaredNorm(v, n)93proof94  have normal := projectionCoefficient.normal_equation at u, v, n with *95  rw [lineDistanceSquared.expansion, lineDistanceSquared.expansion]96  linear_combination ((2:R) * C - (2:R) * projectionCoefficient(u, v, n)) * normal9798theorem lineProjection.minimum (u v : Set, C : Real, n : Nat)99  requires sequenceU : RealSequence(u)100  requires sequenceV : RealSequence(v)101  requires nonzero : squaredNorm(v, n) ≠ (0:R)102  shows lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) ≤103    lineDistanceSquared(u, v, C, n)104proof105  have squareNonnegative := realMulValue.square_nonnegative at C - projectionCoefficient(u, v, n) with *106  have normNonnegative := squaredNorm.nonnegative at v, n with *107  have gapNonnegative := realMulValue.nonnegative at (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)), squaredNorm(v, n) with *108  have translated := realLe.add_left at realZero, (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)) * squaredNorm(v, n), lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) with *109  have decomposition := lineProjection.distance_decomposition at u, v, C, n with *110  rw decomposition111  simpa [realAddValue.zero_right] using translated112113theorem lineProjection.strict_minimum (u v : Set, C : Real, n : Nat)114  requires sequenceU : RealSequence(u)115  requires sequenceV : RealSequence(v)116  requires nonzero : squaredNorm(v, n) ≠ (0:R)117  requires different : C ≠ projectionCoefficient(u, v, n)118  shows lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) <119    lineDistanceSquared(u, v, C, n)120proof121  have differenceNonzero : C - projectionCoefficient(u, v, n) ≠ 0122  . intro zero123    have equal : C = projectionCoefficient(u, v, n)124    . linear_combination zero125    exact different with equal126  have squarePositive := realMulValue.square_positive at C - projectionCoefficient(u, v, n) with *127  have normNonnegative := squaredNorm.nonnegative at v, n with *128  have normPositive := Real.positive_of_nonnegative_of_ne_zero at squaredNorm(v, n) with *129  have gapPositive := realLess.mul_right_of_positive at realZero, (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)), squaredNorm(v, n) with *130  simp [realMulValue.zero_left] at gapPositive131  have translated := realLess.add_left at realZero, (C - projectionCoefficient(u, v, n)) * (C - projectionCoefficient(u, v, n)) * squaredNorm(v, n), lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n) with *132  have decomposition := lineProjection.distance_decomposition at u, v, C, n with *133  rw decomposition134  simpa [realAddValue.zero_right] using translated135136theorem lineProjection.unique_minimizer (u v : Set, C : Real, n : Nat)137  requires sequenceU : RealSequence(u)138  requires sequenceV : RealSequence(v)139  requires nonzero : squaredNorm(v, n) ≠ (0:R)140  requires minimizing : lineDistanceSquared(u, v, C, n) ≤141    lineDistanceSquared(u, v, projectionCoefficient(u, v, n), n)142  shows C = projectionCoefficient(u, v, n)143proof144  by_contra different145  have strict := lineProjection.strict_minimum with *146  have candidate := lineDistanceSquared.real at u, v, C, n with *147  have cycle := realLe.transitive_less with *148  have member : lineDistanceSquared(u, v, C, n) ∈ ℝ149  . settle with Real150  have irreflexive := realLess.irreflexive_on151  exact irreflexive with lineDistanceSquared(u, v, C, n), cycle152153-- One orthogonalization step: retain the component perpendicular to v.154def orthogonalResidual (u v : RealSequence, n : Nat) : RealSequence :=155  subSequence(u, lineProjection(u, v, n))156157theorem orthogonalResidual.sequence (u v : Set, n : Nat)158  requires sequenceU : RealSequence(u)159  requires sequenceV : RealSequence(v)160  shows RealSequence(orthogonalResidual(u, v, n))161proof162  have projection := lineProjection.sequence at u, v, n with *163  unfold orthogonalResidual164  exact subSequence.sequence with *165166theorem orthogonalResidual.orthogonal (u v : Set, n : Nat)167  requires sequenceU : RealSequence(u)168  requires sequenceV : RealSequence(v)169  requires nonzero : squaredNorm(v, n) ≠ (0:R)170  shows Orthogonal(orthogonalResidual(u, v, n), v, n)171proof172  unfold orthogonalResidual173  exact lineProjection.orthogonal with *174175theorem orthogonalResidual.reconstruction (u v : Set, n i : Nat)176  requires sequenceU : RealSequence(u)177  requires sequenceV : RealSequence(v)178  shows u[i] = (lineProjection(u, v, n)[i]:R) + (orthogonalResidual(u, v, n)[i]:R)179proof180  have projection := lineProjection.sequence at u, v, n with *181  have ui := RealSequence.value at u, i with *182  have pi := RealSequence.value at lineProjection(u, v, n), i with *183  unfold orthogonalResidual184  rw subSequence.value185  ring186187theorem projectionCoefficient.on_line (v : Set, C : Real, n : Nat)188  requires sequence : RealSequence(v)189  requires nonzero : squaredNorm(v, n) ≠ (0:R)190  shows projectionCoefficient(mulSequence(constantSequence(C), v), v, n) = C191proof192  have constants := constantSequence.sequence with *193  have scaled := mulSequence.sequence at constantSequence(C), v with *194  have inverseLaw := realInvValue.mul_right at squaredNorm(v, n) with *195  unfold projectionCoefficient196  rw dotProduct.scale_left197  unfold squaredNorm at inverseLaw198  unfold squaredNorm199  linear_combination C * inverseLaw200201theorem lineProjection.idempotent (u v : Set, n : Nat)202  requires sequenceU : RealSequence(u)203  requires sequenceV : RealSequence(v)204  requires nonzero : squaredNorm(v, n) ≠ (0:R)205  shows lineProjection(lineProjection(u, v, n), v, n) = lineProjection(u, v, n)206proof207  unfold lineProjection208  rw projectionCoefficient.on_line209  rfl210211theorem projectionCoefficient.of_orthogonal (u v : Set, n : Nat)212  requires sequenceU : RealSequence(u)213  requires sequenceV : RealSequence(v)214  requires perpendicular : Orthogonal(u, v, n)215  shows projectionCoefficient(u, v, n) = (0:R)216proof217  unfold Orthogonal at perpendicular218  simp [projectionCoefficient]219  ring220221theorem orthogonalResidual.projection_coefficient (u v : Set, n : Nat)222  requires sequenceU : RealSequence(u)223  requires sequenceV : RealSequence(v)224  requires nonzero : squaredNorm(v, n) ≠ (0:R)225  shows projectionCoefficient(orthogonalResidual(u, v, n), v, n) = (0:R)226proof227  have residual := orthogonalResidual.sequence at u, v, n with *228  have perpendicular := orthogonalResidual.orthogonal at u, v, n with *229  exact projectionCoefficient.of_orthogonal with *230231theorem orthogonalResidual.pythagorean (u v : Set, n : Nat)232  requires sequenceU : RealSequence(u)233  requires sequenceV : RealSequence(v)234  requires nonzero : squaredNorm(v, n) ≠ (0:R)235  shows squaredNorm(u, n) = squaredNorm(lineProjection(u, v, n), n) +236    squaredNorm(orthogonalResidual(u, v, n), n)237proof238  have projection := lineProjection.sequence at u, v, n with *239  have constants := constantSequence.sequence at projectionCoefficient(u, v, n) with *240  have scaled := mulSequence.sequence at constantSequence(projectionCoefficient(u, v, n)), v with *241  have normal := projectionCoefficient.normal_equation at u, v, n with *242  unfold orthogonalResidual243  rw squaredNorm.sub244  unfold lineProjection245  simp [squaredNorm.scale, dotProduct.scale_right]246  linear_combination -(2:R) * projectionCoefficient(u, v, n) * normal247248theorem projectionCoefficient.add (u w v : Set, n : Nat)249  requires sequenceU : RealSequence(u)250  requires sequenceW : RealSequence(w)251  requires sequenceV : RealSequence(v)252  shows projectionCoefficient(addSequence(u, w), v, n) =253    projectionCoefficient(u, v, n) + projectionCoefficient(w, v, n)254proof255  have sum := addSequence.sequence at u, w with *256  unfold projectionCoefficient257  rw dotProduct.add_left258  ring259260theorem projectionCoefficient.scale (u v : Set, C : Real, n : Nat)261  requires sequenceU : RealSequence(u)262  requires sequenceV : RealSequence(v)263  shows projectionCoefficient(mulSequence(constantSequence(C), u), v, n) =264    C * projectionCoefficient(u, v, n)265proof266  have constants := constantSequence.sequence with *267  have scaled := mulSequence.sequence at constantSequence(C), u with *268  unfold projectionCoefficient269  rw dotProduct.scale_left270  ring271