feat(linear_algebra/matrix): LDL Decomposition (#15220)

Any positive definite matrix `S` can be decomposed as `S = LDLᴴ` where `L` is a lower-triangular matrix and `D` is a diagonal matrix.

src/data/matrix/basic.lean

@@ -1244,6 +1244,10 @@ lemma vec_mul_conj_transpose [fintype n] [star_ring α] (A : matrix m n α) (x :
   vec_mul x Aᴴ = star (mul_vec A (star x)) :=
 funext $ λ i, dot_product_star _ _
 
+lemma mul_mul_apply [fintype n] (A B C : matrix n n α) (i j : n) :
+  (A ⬝ B ⬝ C) i j = A i ⬝ᵥ (B.mul_vec (Cᵀ j)) :=
+by { rw matrix.mul_assoc, simpa only [mul_apply, dot_product, mul_vec] }
+
 end non_unital_semiring
 
 section non_assoc_semiring

src/linear_algebra/matrix/basis.lean

@@ -3,6 +3,7 @@ Copyright (c) 2019 Johannes Hölzl. All rights reserved.
 Released under Apache 2.0 license as described in the file LICENSE.
 Authors: Johannes Hölzl, Patrick Massot, Casper Putz, Anne Baanen
 -/
+import linear_algebra.matrix.nonsingular_inverse
 import linear_algebra.matrix.reindex
 import linear_algebra.matrix.to_lin
 
@@ -234,6 +235,11 @@ lemma basis.to_matrix_mul_to_matrix_flip [decidable_eq ι] [fintype ι'] :
   b.to_matrix b' ⬝ b'.to_matrix b = 1 :=
 by rw [basis.to_matrix_mul_to_matrix, basis.to_matrix_self]
 
+/-- A matrix whose columns form a basis `b'`, expressed w.r.t. a basis `b`, is invertible. -/
+def basis.invertible_to_matrix [decidable_eq ι] [fintype ι] (b b' : basis ι R₂ M₂) :
+  invertible (b.to_matrix b') :=
+matrix.invertible_of_left_inverse _ _ (basis.to_matrix_mul_to_matrix_flip _ _)
+
 @[simp]
 lemma basis.to_matrix_reindex
   (b : basis ι R M) (v : ι' → M) (e : ι ≃ ι') :

src/linear_algebra/matrix/ldl.lean

@@ -0,0 +1,113 @@
+/-
+Copyright (c) 2022 Alexander Bentkamp. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Alexander Bentkamp
+-/
+import analysis.inner_product_space.gram_schmidt_ortho
+import linear_algebra.matrix.pos_def
+
+/-! # LDL decomposition
+
+This file proves the LDL-decomposition of matricies: Any positive definite matrix `S` can be
+decomposed as `S = LDLᴴ` where `L` is a lower-triangular matrix and `D` is a diagonal matrix.
+
+## Main definitions
+
+ * `LDL.lower` is the lower triangular matrix `L`.
+ * `LDL.lower_inv` is the inverse of the lower triangular matrix `L`.
+ * `LDL.diag` is the diagonal matrix `D`.
+
+## Main result
+
+* `ldl_decomposition` states that any positive definite matrix can be decomposed as `LDLᴴ`.
+
+## TODO
+
+* Prove that `LDL.lower` is lower triangular from `LDL.lower_inv_triangular`.
+
+-/
+
+variables {𝕜 : Type*} [is_R_or_C 𝕜]
+variables {n : Type*} [linear_order n] [is_well_order n (<)] [locally_finite_order_bot n]
+
+local notation `⟪`x`, `y`⟫` :=
+@inner 𝕜 (n → 𝕜) (pi_Lp.inner_product_space (λ _, 𝕜)).to_has_inner x y
+
+open matrix
+open_locale matrix
+variables {S : matrix n n 𝕜} [fintype n] (hS : S.pos_def)
+
+/-- The inverse of the lower triangular matrix `L` of the LDL-decomposition. It is obtained by
+applying Gram-Schmidt-Orthogonalization w.r.t. the inner product induced by `Sᵀ` on the standard
+basis vectors `pi.basis_fun`. -/
+noncomputable def LDL.lower_inv : matrix n n 𝕜 :=
+@gram_schmidt
+  𝕜 (n → 𝕜) _ (inner_product_space.of_matrix hS.transpose) n _ _ _ (λ k, pi.basis_fun 𝕜 n k)
+
+lemma LDL.lower_inv_eq_gram_schmidt_basis :
+  LDL.lower_inv hS = ((pi.basis_fun 𝕜 n).to_matrix
+    (@gram_schmidt_basis 𝕜 (n → 𝕜) _
+    (inner_product_space.of_matrix hS.transpose) n _ _ _ (pi.basis_fun 𝕜 n)))ᵀ :=
+begin
+  ext i j,
+  rw [LDL.lower_inv, basis.coe_pi_basis_fun.to_matrix_eq_transpose, coe_gram_schmidt_basis],
+  refl
+end
+
+noncomputable instance LDL.invertible_lower_inv : invertible (LDL.lower_inv hS) :=
+begin
+  rw [LDL.lower_inv_eq_gram_schmidt_basis],
+  haveI := basis.invertible_to_matrix (pi.basis_fun 𝕜 n)
+    (@gram_schmidt_basis 𝕜 (n → 𝕜) _ (inner_product_space.of_matrix hS.transpose)
+      n _ _ _ (pi.basis_fun 𝕜 n)),
+  apply_instance
+end
+
+lemma LDL.lower_inv_orthogonal {i j : n} (h₀ : i ≠ j) :
+  ⟪(LDL.lower_inv hS i), Sᵀ.mul_vec (LDL.lower_inv hS j)⟫ = 0 :=
+show @inner 𝕜 (n → 𝕜) (inner_product_space.of_matrix hS.transpose).to_has_inner
+    (LDL.lower_inv hS i)
+    (LDL.lower_inv hS j) = 0,
+by apply gram_schmidt_orthogonal _ _ h₀
+
+/-- The entries of the diagonal matrix `D` of the LDL decomposition. -/
+noncomputable def LDL.diag_entries : n → 𝕜 :=
+  λ i, ⟪star (LDL.lower_inv hS i), S.mul_vec (star (LDL.lower_inv hS i))⟫
+
+/-- The diagonal matrix `D` of the LDL decomposition. -/
+noncomputable def LDL.diag : matrix n n 𝕜 := matrix.diagonal (LDL.diag_entries hS)
+
+lemma LDL.lower_inv_triangular {i j : n} (hij : i < j) :
+  LDL.lower_inv hS i j = 0 :=
+by rw [← @gram_schmidt_triangular 𝕜 (n → 𝕜) _ (inner_product_space.of_matrix hS.transpose) n _ _ _
+    i j hij (pi.basis_fun 𝕜 n), pi.basis_fun_repr, LDL.lower_inv]
+
+/-- Inverse statement of **LDL decomposition**: we can conjugate a positive definite matrix
+by some lower triangular matrix and get a diagonal matrix. -/
+lemma LDL.diag_eq_lower_inv_conj : LDL.diag hS = LDL.lower_inv hS ⬝ S ⬝ (LDL.lower_inv hS)ᴴ :=
+begin
+  ext i j,
+  by_cases hij : i = j,
+  { simpa only [hij, LDL.diag, diagonal_apply_eq, LDL.diag_entries, matrix.mul_assoc, inner,
+      pi.star_apply, is_R_or_C.star_def, star_ring_end_self_apply] },
+  { simp only [LDL.diag, hij, diagonal_apply_ne, ne.def, not_false_iff, mul_mul_apply],
+    rw [conj_transpose, transpose_map, transpose_transpose, dot_product_mul_vec,
+      (LDL.lower_inv_orthogonal hS (λ h : j = i, hij h.symm)).symm,
+      ← inner_conj_sym, mul_vec_transpose, euclidean_space.inner_eq_star_dot_product,
+      ← is_R_or_C.star_def, ← star_dot_product_star, dot_product_comm, star_star],
+    refl }
+end
+
+/-- The lower triangular matrix `L` of the LDL decomposition. -/
+noncomputable def LDL.lower := (LDL.lower_inv hS)⁻¹
+
+/-- **LDL decomposition**: any positive definite matrix `S` can be
+decomposed as `S = LDLᴴ` where `L` is a lower-triangular matrix and `D` is a diagonal matrix.  -/
+theorem LDL.lower_conj_diag :
+  LDL.lower hS ⬝ LDL.diag hS ⬝ (LDL.lower hS)ᴴ = S :=
+begin
+  rw [LDL.lower, conj_transpose_nonsing_inv, matrix.mul_assoc,
+    matrix.inv_mul_eq_iff_eq_mul_of_invertible (LDL.lower_inv hS),
+    matrix.mul_inv_eq_iff_eq_mul_of_invertible],
+  exact LDL.diag_eq_lower_inv_conj hS,
+end

src/linear_algebra/matrix/nonsingular_inverse.lean

@@ -149,6 +149,29 @@ def invertible_of_left_inverse (h : B ⬝ A = 1) : invertible A :=
 def invertible_of_right_inverse (h : A ⬝ B = 1) : invertible A :=
 ⟨B, mul_eq_one_comm.mp h, h⟩
 
+/-- The transpose of an invertible matrix is invertible. -/
+instance invertible_transpose [invertible A] : invertible Aᵀ :=
+begin
+  haveI : invertible Aᵀ.det,
+    by simpa using det_invertible_of_invertible A,
+  exact invertible_of_det_invertible Aᵀ
+end
+
+/-- A matrix is invertible if the transpose is invertible. -/
+def invertible__of_invertible_transpose [invertible Aᵀ] : invertible A :=
+begin
+  rw ←transpose_transpose A,
+  apply_instance
+end
+
+/-- A matrix is invertible if the conjugate transpose is invertible. -/
+def invertible_of_invertible_conj_transpose [star_ring α] [invertible Aᴴ] :
+  invertible A :=
+begin
+  rw ←conj_transpose_conj_transpose A,
+  apply_instance
+end
+
 /-- Given a proof that `A.det` has a constructive inverse, lift `A` to `(matrix n n α)ˣ`-/
 def unit_of_det_invertible [invertible A.det] : (matrix n n α)ˣ :=
 @unit_of_invertible _ _ A (invertible_of_det_invertible A)
@@ -164,6 +187,12 @@ lemma is_unit_det_of_invertible [invertible A] : is_unit A.det :=
 
 variables {A B}
 
+lemma is_unit_of_left_inverse (h : B ⬝ A = 1) : is_unit A :=
+⟨⟨A, B, mul_eq_one_comm.mp h, h⟩, rfl⟩
+
+lemma is_unit_of_right_inverse (h : A ⬝ B = 1) : is_unit A :=
+⟨⟨A, B, h, mul_eq_one_comm.mp h⟩, rfl⟩
+
 lemma is_unit_det_of_left_inverse (h : B ⬝ A = 1) : is_unit A.det :=
 @is_unit_of_invertible _ _ _ (det_invertible_of_left_inverse _ _ h)
 
@@ -284,6 +313,16 @@ nonsing_inv_mul_cancel_right A B (is_unit_det_of_invertible A)
   A⁻¹ ⬝ (A ⬝ B) = B :=
 nonsing_inv_mul_cancel_left A B (is_unit_det_of_invertible A)
 
+lemma inv_mul_eq_iff_eq_mul_of_invertible (A B C : matrix n n α) [invertible A] :
+  A⁻¹ ⬝ B = C ↔ B = A ⬝ C :=
+⟨λ h, by rw [←h, mul_inv_cancel_left_of_invertible],
+ λ h, by rw [h, inv_mul_cancel_left_of_invertible]⟩
+
+lemma mul_inv_eq_iff_eq_mul_of_invertible (A B C : matrix n n α) [invertible A] :
+  B ⬝ A⁻¹ = C ↔ B = C ⬝ A :=
+⟨λ h, by rw [←h, inv_mul_cancel_right_of_invertible],
+ λ h, by rw [h, mul_inv_cancel_right_of_invertible]⟩
+
 lemma nonsing_inv_cancel_or_zero :
   (A⁻¹ ⬝ A = 1 ∧ A ⬝ A⁻¹ = 1) ∨ A⁻¹ = 0 :=
 begin

src/linear_algebra/matrix/pos_def.lean

@@ -31,6 +31,13 @@ M.is_hermitian ∧ ∀ x : n → 𝕜, x ≠ 0 → 0 < is_R_or_C.re (dot_product
 
 lemma pos_def.is_hermitian {M : matrix n n 𝕜} (hM : M.pos_def) : M.is_hermitian := hM.1
 
+lemma pos_def.transpose {M : matrix n n 𝕜} (hM : M.pos_def) : Mᵀ.pos_def :=
+begin
+  refine ⟨is_hermitian.transpose hM.1, λ x hx, _⟩,
+  convert hM.2 (star x) (star_ne_zero.2 hx) using 2,
+  rw [mul_vec_transpose, matrix.dot_product_mul_vec, star_star, dot_product_comm]
+end
+
 lemma pos_def_of_to_quadratic_form' [decidable_eq n] {M : matrix n n ℝ}
   (hM : M.is_symm) (hMq : M.to_quadratic_form'.pos_def) :
   M.pos_def :=