Banach-Tarski and SO(3)

Stream: general

[Arthur Tilley]

Hi folks,

I’m nearing the end of a formalization of the Banach-Tarski theorem.

Some project-level definitions:

abbrev R3_raw := (Fin 3) → ℝ
abbrev R3 :=  EuclideanSpace ℝ (Fin 3)
def S2 : Set R3 := Metric.sphere (0: R3) 1
abbrev MAT:= Matrix (Fin 3) (Fin 3) ℝ

-- The standard action given by matrix multiplication.
instance SO3_action_on_R3_raw: MulAction SO3 R3_raw where
  smul g v  := (Matrix.mulVec (g : Matrix (Fin 3) (Fin 3) ℝ) v)
  one_smul v := Matrix.one_mulVec v
  mul_smul x y v := (Matrix.mulVec_mulVec v (x : Matrix (Fin 3) (Fin 3) ℝ) (y : Matrix (Fin 3) (Fin 3) ℝ)).symm

instance SO3_action_on_R3: MulAction SO3 R3 := inferInstance

def f {X : Type*} {G: Type*} [Group G] [MulAction G X] (g : G): X → X := fun x: X ↦ g • x

abbrev SO3 := Matrix.specialOrthogonalGroup (Fin 3) ℝ

I have completed what I think should be most of the hard parts, and all that remains is to represent members of SO(3).

Strictly speaking, I need def rot (ax: S2) (θ:ℝ) : SO3 :=... allowing me to prove some fairly basic geometric lemmas.

My current approach is to use the standard Rodrigues rotation matrices as rotation matrices:

-- Rodrigues' formula for the rotation matrix :  I + (sin θ)K + (1-cosθ)K²


def K_mat (a: R3): MAT := !![
  0, -(a 2), (a 1);
  (a 2), 0, -(a 0);
  -(a 1), (a 0), 0;
]

noncomputable def rot_mat (ax: S2) (θ:ℝ) : MAT := (1:MAT) + (Real.sin θ)•(K_mat ax) + (1 - Real.cos θ)•((K_mat ax) ^ 2)

noncomputable def rot (ax: S2) (θ:ℝ) : SO3 :=
  let M := rot_mat ax θ
  have M_is_special : M ∈ SO3 := by
    rw [Matrix.mem_specialOrthogonalGroup_iff]
    constructor
    · rw [Matrix.mem_orthogonalGroup_iff]
      -- M * Mᵀ = 1
      sorry
    · -- M.det = 1
      sorry
  ⟨M, M_is_special⟩

I honestly find these computations are pretty nasty for a generic matrix, but I am still relatively new to Lean.

I need to be able to prove results like

lemma rot_fixed_back_gen_z_single (v w: R3):
f (rot z_axis (z_ang_diff v w).toReal) v = w := by
  sorry

where

noncomputable def z_ang_diff (s t: R3) : Real.Angle :=
  -- NOTE! Projects down to the  z = 0 plane
  let s2: R2 := !₂[(s 0), (s 1)]
  let t2: R2 := !₂[(t 0), (t 1)]
  o.oangle s2 t2

and

lemma rot_comp_add (ax: S2) (t1 t2 : ℝ) : (rot ax t1) * (rot ax t2) = (rot ax (t1 + t2)) := by

lemma rot_fixed (axis: S2) (v: R3): f (rot axis t) v = v → ∃k:ℤ, t = k * 2 * Real.pi := sorry

And to talk about changes of basis

noncomputable def COB_to_Z (axis: S2) : SO3 :=s
  dite (axis.val≠z_axis.val ∧ axis.val≠(-z_axis.val))
  (fun p: _ ↦ rot (s2_cross axis z_axis p.left p.right) (unsigned_ang axis z_axis))
  (fun p: _ ↦ dite (axis.val=z_axis.val)
    (fun _ : _ ↦ (1: SO3))
    (fun _ : _ ↦ rot (x_axis) Real.pi)
  )

allowing results like

lemma ctza_def (axis: S2):   (COB_to_Z axis) • axis.val = z_axis.val := sorry

I have dug into Mathlib quite a bit, but most of the rotation-focused results seem to be for two-dimensional vector spaces.

I’m sure it’s possible to grind through each of these remaining sorrys, but I want to make sure first that I’m not missing anything, or whether anyone knows of any other libraries that handle allow manipulating members of SO(3) as rotations.

The repo is here:

https://github.com/aetilley/banach-tarski

Thanks,
Arthur

[Eric Wieser]

I haven't looked at your repo, and it's a little tricky to give feedback on your message contents without seeing your definition of S2 etc, but one small bit of feedback; dite P (fun p ↦ x) (fun p ↦ y) is much more idiomatic as if p : P then x else y or even just if P then x else y.

[Mirek Olšák]

Hi, congratulations on almost finishing this. Myself, I wouldn't know some parts of how to do Banach-Tarski (why the rotations form a free group).

I understand you have all in your repo but it would be helpful if you paste here all the auxiliary definitions (MAT, S2, ...), so people can play with it immediatelly.

I am rather surprised if you need some calculations manually multiplying 3 x 3 matrices. My approach would be to try to have

• Some conversion matrix in SO3 that aligns a variable axis with a specific axis such as z = (0,0,1), a composition of two rotations around z, and around x.
• Rotation around that z, so essentially a 2d-rotation -- there I could compose angles, as well as argue that it doesn't matter which way I moved my arbitrary axis to z, because rotations in 2d commute.
• If I need to obtain an axis of an arbitrary element of SO3, I would hope to get it as the eigen-vector.

On the other hand, I am not an expert in matrices in Lean, so I might be wrong...

[Arthur Tilley]

@Eric Wieser
I have added some auxilliary definitiions to hopefully make it more understandable.
I'm walking a fine line here by both 1) not wanting people to have to look at my code and 2) not wanting my first post on general to be a book.

@Mirek Olšák Don't worry about the free group stuff. That's done. Really all that is remaining is to prove some basic geometry concerning about rotations in SO3. What you mention at the end is part of it. It's just that but proving theorems that involve an arbitrary rotation matrix (for a given axis and angle) turns into a mess. I'm crossing my fingers that there may be a library somewhere. Or maybe some matlib theorems that I'm missing.

[Mirek Olšák]

I wanted to translate any rotation to a 2d rotation. (or composition of them)

[Mirek Olšák]

Actually, I realized, it probably doesn't have to be a composition of rotation along x & z. It should be possible to take a subspace generated with axis and z, rotate axis to z there, and then extend it with LinearIsometry.extend

[Mirek Olšák]

But I would try to stay abstract, ideally without explicitly writing 3x3 matrices.

[Michael Rothgang]

@Heather Macbeth wrote some code for this for the sphere eversion project: it's in https://github.com/leanprover-community/sphere-eversion/blob/master/SphereEversion/ToMathlib/Analysis/InnerProductSpace/Rotation.lean

[Michael Rothgang]

If that code is useful, a PR adding it to mathlib would be welcome.

[Michael Rothgang]

(Perhaps it should be generalised to n dimensions first. I'm unsure about that part.)

[Arthur Tilley]

@Michael Rothgang @Heather Macbeth This looks like it could be very helpful. Thanks so much.

[Moritz Doll]

(deleted)

[Joseph Myers]

For me the natural abstract approach for rotations in arbitrary dimensions starts with docs#Orientation.rotation (2d rotation) which we already have (plus the identity docs#LinearIsometryEquiv.refl on subspaces where you're not rotating) and then combines the isometries on orthogonal subspaces into an isometry of the whole space.

Combining isometries on subspaces (which in general applies to both isometries and isometric equivalences, and with domain and codomain that need not be the same) also decomposes into smaller operations. One operation takes a pair or indexed family of isometries or isometric equivalences and produces an isometry or isometric equivalent of the pairwise or indexed product spaces (the WithLp 2 / PiLp 2 versions thereof, that is). The other operation is an isometric equivalence between an inner product space and the pairwise or indexed product of subspaces that are orthogonal and span the whole space.

(Then if you wish you could show that all orientation-preserving linear isometric equivalences of a finite-dimensional real inner product space can be constructed this way from 2d rotations on orthogonal subspaces, though you probably don't need that for this application.)

[Arthur Tilley]

@Joseph Myers It sounds like you and @Mirek Olšák might be talking about roughly the same approach? That is to define a rotation in the (2d) space orthogonal to the axis then extend that iso LinearIsometry.extend

I might give this a go.

While the code that @Michael Rothgang links is in the same direction, it's a bit hard for my small brain to be able to generalize what is being done there to other results that I need talking about exactly which isometries take some vector v to some vector w.

Also in places I think I need to talk about oriented angle between two vectors, and I think that a definition for that is currently only defined for two dimensional spaces.

[Joseph Myers]

The documentation for docs#LinearIsometry.extend suggests it's an arbitrary extension of an isometry from a subspace to the whole space. In the case of an isometry mapping the subspace to itself I think it will end up extending by the identity on the orthogonal subspace, but you'd be relying on the choice of junk values that isn't currently documented and doesn't currently have any API lemmas to say what that choice is.

[Joseph Myers]

Oriented angles are only meaningful in two dimensions. You can of course put an orientation on a two-dimensional subspace and then take angles between vectors in that subspace.

[Mirek Olšák]

Oriented angles are also meaningful if you have an orientation on the fixed subspace (which you do in 3d, if the axis is represented with an element of the unit sphere).

[Joseph Myers]

I think of that as a matter of API for combining orientations on subspaces. Given orientations of two spaces, you can map them to an orientation of the product space (indexed by the sum of the index types of the two subspace orientations) in the obvious way, and there is also a reverse operation that you can use to map from an orientation of the product space and an orientation of one of the two spaces to an orientation of the other one of the two spaces. In particular, an orientation of $(n+2)$-space and one of $n$-space yield an orientation of the $2$-space in which the rotation takes place. (If those orientations are indexed by Fin (n + 2), Fin n and Fin 2, this also involves a choice of equivalence Fin (n + 2) ≃ (Fin n ⊕ Fin 2).)