From 0fc5496584478976b86f17941b3ff25ad7c4ebd6 Mon Sep 17 00:00:00 2001 From: Alexander Bentkamp Date: Fri, 19 Aug 2022 19:30:31 +0000 Subject: [PATCH] feat(linear_algebra/matrix): LDL Decomposition (#15220) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 | 4 + src/linear_algebra/matrix/basis.lean | 6 + src/linear_algebra/matrix/ldl.lean | 113 ++++++++++++++++++ .../matrix/nonsingular_inverse.lean | 39 ++++++ src/linear_algebra/matrix/pos_def.lean | 7 ++ 5 files changed, 169 insertions(+) create mode 100644 src/linear_algebra/matrix/ldl.lean diff --git a/src/data/matrix/basic.lean b/src/data/matrix/basic.lean index 56995376e726e..b9cd3f40a67b2 100644 --- a/src/data/matrix/basic.lean +++ b/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 diff --git a/src/linear_algebra/matrix/basis.lean b/src/linear_algebra/matrix/basis.lean index 195d49b6c7cc6..dd1d97858ab9c 100644 --- a/src/linear_algebra/matrix/basis.lean +++ b/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 : ι ≃ ι') : diff --git a/src/linear_algebra/matrix/ldl.lean b/src/linear_algebra/matrix/ldl.lean new file mode 100644 index 0000000000000..4c5ba1efe5525 --- /dev/null +++ b/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 diff --git a/src/linear_algebra/matrix/nonsingular_inverse.lean b/src/linear_algebra/matrix/nonsingular_inverse.lean index 0d9c48d78bd74..48634f7d65856 100644 --- a/src/linear_algebra/matrix/nonsingular_inverse.lean +++ b/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 diff --git a/src/linear_algebra/matrix/pos_def.lean b/src/linear_algebra/matrix/pos_def.lean index 036ba42ff1c3c..c584f1d8e917c 100644 --- a/src/linear_algebra/matrix/pos_def.lean +++ b/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 :=