feat(linear_algebra/matrix): Spectral theorem for matrices (#14231)

Co-authored-by: Eric Wieser <wieser.eric@gmail.com>

Author: abentkamp

src/analysis/inner_product_space/pi_L2.lean

@@ -39,7 +39,7 @@ open_locale big_operators uniformity topological_space nnreal ennreal complex_co
 
 noncomputable theory
 
-variables {ι : Type*}
+variables {ι : Type*} {ι' : Type*}
 variables {𝕜 : Type*} [is_R_or_C 𝕜] {E : Type*} [inner_product_space 𝕜 E]
 variables {E' : Type*} [inner_product_space 𝕜 E']
 variables {F : Type*} [inner_product_space ℝ F]
@@ -112,6 +112,9 @@ instance : inner_product_space 𝕜 (euclidean_space 𝕜 ι) := by apply_instan
 lemma finrank_euclidean_space_fin {n : ℕ} :
   finite_dimensional.finrank 𝕜 (euclidean_space 𝕜 (fin n)) = n := by simp
 
+lemma euclidean_space.inner_eq_star_dot_product (x y : euclidean_space 𝕜 ι) :
+  ⟪x, y⟫ = matrix.dot_product (star $ pi_Lp.equiv _ _ x) (pi_Lp.equiv _ _ y) := rfl
+
 /-- A finite, mutually orthogonal family of subspaces of `E`, which span `E`, induce an isometry
 from `E` to `pi_Lp 2` of the subspaces equipped with the `L2` inner product. -/
 def direct_sum.is_internal.isometry_L2_of_orthogonal_family

src/analysis/normed_space/pi_Lp.lean

@@ -83,8 +83,11 @@ to compare the `L^p` and `L^∞` distances through it. -/
 protected def equiv : pi_Lp p α ≃ Π (i : ι), α i :=
 equiv.refl _
 
-@[simp] lemma equiv_apply (x : pi_Lp p α) (i : ι) : pi_Lp.equiv p α x i = x i := rfl
-@[simp] lemma equiv_symm_apply (x : Π i, α i) (i : ι) : (pi_Lp.equiv p α).symm x i = x i := rfl
+lemma equiv_apply (x : pi_Lp p α) (i : ι) : pi_Lp.equiv p α x i = x i := rfl
+lemma equiv_symm_apply (x : Π i, α i) (i : ι) : (pi_Lp.equiv p α).symm x i = x i := rfl
+
+@[simp] lemma equiv_apply' (x : pi_Lp p α) : pi_Lp.equiv p α x = x := rfl
+@[simp] lemma equiv_symm_apply' (x : Π i, α i) : (pi_Lp.equiv p α).symm x = x := rfl
 
 section
 /-!

src/linear_algebra/matrix/basis.lean

@@ -177,6 +177,32 @@ lemma basis_to_matrix_mul_linear_map_to_matrix_mul_basis_to_matrix
   c.to_matrix c' ⬝ linear_map.to_matrix b' c' f ⬝ b'.to_matrix b = linear_map.to_matrix b c f :=
 by rw [basis_to_matrix_mul_linear_map_to_matrix, linear_map_to_matrix_mul_basis_to_matrix]
 
+lemma basis_to_matrix_mul [decidable_eq κ]
+    (b₁ : basis ι R M) (b₂ : basis ι' R M) (b₃ : basis κ R N) (A : matrix ι' κ R) :
+  b₁.to_matrix b₂ ⬝ A = linear_map.to_matrix b₃ b₁ (to_lin b₃ b₂ A) :=
+begin
+  have := basis_to_matrix_mul_linear_map_to_matrix b₃ b₁ b₂ (matrix.to_lin b₃ b₂ A),
+  rwa [linear_map.to_matrix_to_lin] at this
+end
+
+lemma mul_basis_to_matrix [decidable_eq ι] [decidable_eq ι']
+    (b₁ : basis ι R M) (b₂ : basis ι' R M) (b₃ : basis κ R N) (A : matrix κ ι R) :
+  A ⬝ b₁.to_matrix b₂ = linear_map.to_matrix b₂ b₃ (to_lin b₁ b₃ A) :=
+begin
+  have := linear_map_to_matrix_mul_basis_to_matrix b₂ b₁ b₃ (matrix.to_lin b₁ b₃ A),
+  rwa [linear_map.to_matrix_to_lin] at this
+end
+
+lemma basis_to_matrix_basis_fun_mul (b : basis ι R (ι → R)) (A : matrix ι ι R) :
+  b.to_matrix (pi.basis_fun R ι) ⬝ A = (λ i j, b.repr (Aᵀ j) i) :=
+begin
+  classical,
+  simp only [basis_to_matrix_mul _ _ (pi.basis_fun R ι), matrix.to_lin_eq_to_lin'],
+  ext i j,
+  rw [linear_map.to_matrix_apply, matrix.to_lin'_apply, pi.basis_fun_apply,
+    matrix.mul_vec_std_basis_apply]
+end
+
 /-- A generalization of `linear_map.to_matrix_id`. -/
 @[simp] lemma linear_map.to_matrix_id_eq_basis_to_matrix [decidable_eq ι] :
   linear_map.to_matrix b b' id = b'.to_matrix b :=

src/linear_algebra/matrix/hermitian.lean

@@ -0,0 +1,176 @@
+/-
+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.spectrum
+
+/-! # Hermitian matrices
+
+This file defines hermitian matrices and some basic results about them.
+
+## Main definition
+
+ * `matrix.is_hermitian `: a matrix `A : matrix n n α` is hermitian if `Aᴴ = A`.
+
+## Tags
+
+self-adjoint matrix, hermitian matrix
+
+-/
+
+namespace matrix
+
+variables {α β : Type*} {m n : Type*} {A : matrix n n α}
+
+open_locale matrix
+
+local notation `⟪`x`, `y`⟫` := @inner α (pi_Lp 2 (λ (_ : n), α)) _ x y
+
+section non_unital_semiring
+
+variables [non_unital_semiring α] [star_ring α] [non_unital_semiring β] [star_ring β]
+
+/-- A matrix is hermitian if it is equal to its conjugate transpose. On the reals, this definition
+captures symmetric matrices. -/
+def is_hermitian (A : matrix n n α) : Prop := Aᴴ = A
+
+lemma is_hermitian.eq {A : matrix n n α} (h : A.is_hermitian) : Aᴴ = A := h
+
+@[ext]
+lemma is_hermitian.ext {A : matrix n n α} : (∀ i j, star (A j i) = A i j) → A.is_hermitian :=
+by { intros h, ext i j, exact h i j }
+
+lemma is_hermitian.apply {A : matrix n n α} (h : A.is_hermitian) (i j : n) : star (A j i) = A i j :=
+by { unfold is_hermitian at h, rw [← h, conj_transpose_apply, star_star, h] }
+
+lemma is_hermitian.ext_iff {A : matrix n n α} : A.is_hermitian ↔ ∀ i j, star (A j i) = A i j :=
+⟨is_hermitian.apply, is_hermitian.ext⟩
+
+lemma is_hermitian_mul_conj_transpose_self [fintype n] (A : matrix n n α) :
+  (A ⬝ Aᴴ).is_hermitian :=
+by rw [is_hermitian, conj_transpose_mul, conj_transpose_conj_transpose]
+
+lemma is_hermitian_transpose_mul_self [fintype n] (A : matrix n n α) :
+  (Aᴴ ⬝ A).is_hermitian :=
+by rw [is_hermitian, conj_transpose_mul, conj_transpose_conj_transpose]
+
+lemma is_hermitian_add_transpose_self (A : matrix n n α) :
+  (A + Aᴴ).is_hermitian :=
+by simp [is_hermitian, add_comm]
+
+lemma is_hermitian_transpose_add_self (A : matrix n n α) :
+  (Aᴴ + A).is_hermitian :=
+by simp [is_hermitian, add_comm]
+
+@[simp] lemma is_hermitian_zero :
+  (0 : matrix n n α).is_hermitian :=
+conj_transpose_zero
+
+-- TODO: move
+lemma conj_transpose_map {A : matrix n n α} (f : α → β) (hf : f ∘ star = star ∘ f) :
+  Aᴴ.map f = (A.map f)ᴴ :=
+by rw [conj_transpose, conj_transpose, ←transpose_map, map_map, map_map, hf]
+
+@[simp] lemma is_hermitian.map {A : matrix n n α} (h : A.is_hermitian) (f : α → β)
+    (hf : f ∘ star = star ∘ f) :
+  (A.map f).is_hermitian :=
+by {refine (conj_transpose_map f hf).symm.trans _, rw h.eq }
+
+@[simp] lemma is_hermitian.transpose {A : matrix n n α} (h : A.is_hermitian) :
+  Aᵀ.is_hermitian :=
+by { rw [is_hermitian, conj_transpose, transpose_map], congr, exact h }
+
+@[simp] lemma is_hermitian.conj_transpose {A : matrix n n α} (h : A.is_hermitian) :
+  Aᴴ.is_hermitian :=
+h.transpose.map _ rfl
+
+@[simp] lemma is_hermitian.add {A B : matrix n n α} (hA : A.is_hermitian) (hB : B.is_hermitian) :
+  (A + B).is_hermitian :=
+(conj_transpose_add _ _).trans (hA.symm ▸ hB.symm ▸ rfl)
+
+@[simp] lemma is_hermitian.minor {A : matrix n n α} (h : A.is_hermitian) (f : m → n) :
+  (A.minor f f).is_hermitian :=
+(conj_transpose_minor _ _ _).trans (h.symm ▸ rfl)
+
+/-- The real diagonal matrix `diagonal v` is hermitian. -/
+@[simp] lemma is_hermitian_diagonal [decidable_eq n] (v : n → ℝ) :
+  (diagonal v).is_hermitian :=
+diagonal_conj_transpose _
+
+/-- A block matrix `A.from_blocks B C D` is hermitian,
+    if `A` and `D` are hermitian and `Bᴴ = C`. -/
+lemma is_hermitian.from_blocks
+  {A : matrix m m α} {B : matrix m n α} {C : matrix n m α} {D : matrix n n α}
+  (hA : A.is_hermitian) (hBC : Bᴴ = C) (hD : D.is_hermitian) :
+  (A.from_blocks B C D).is_hermitian :=
+begin
+  have hCB : Cᴴ = B, {rw ← hBC, simp},
+  unfold matrix.is_hermitian,
+  rw from_blocks_conj_transpose,
+  congr;
+  assumption
+end
+
+/-- This is the `iff` version of `matrix.is_hermitian.from_blocks`. -/
+lemma is_hermitian_from_blocks_iff
+  {A : matrix m m α} {B : matrix m n α} {C : matrix n m α} {D : matrix n n α} :
+  (A.from_blocks B C D).is_hermitian ↔ A.is_hermitian ∧ Bᴴ = C ∧ Cᴴ = B ∧ D.is_hermitian :=
+⟨λ h, ⟨congr_arg to_blocks₁₁ h, congr_arg to_blocks₂₁ h,
+       congr_arg to_blocks₁₂ h, congr_arg to_blocks₂₂ h⟩,
+ λ ⟨hA, hBC, hCB, hD⟩, is_hermitian.from_blocks hA hBC hD⟩
+
+end non_unital_semiring
+
+section semiring
+
+variables [semiring α] [star_ring α] [semiring β] [star_ring β]
+
+@[simp] lemma is_hermitian_one [decidable_eq n] :
+  (1 : matrix n n α).is_hermitian :=
+conj_transpose_one
+
+end semiring
+
+section ring
+
+variables [ring α] [star_ring α] [ring β] [star_ring β]
+
+@[simp] lemma is_hermitian.neg {A : matrix n n α} (h : A.is_hermitian) :
+  (-A).is_hermitian :=
+(conj_transpose_neg _).trans (congr_arg _ h)
+
+@[simp] lemma is_hermitian.sub {A B : matrix n n α} (hA : A.is_hermitian) (hB : B.is_hermitian) :
+  (A - B).is_hermitian :=
+(conj_transpose_sub _ _).trans (hA.symm ▸ hB.symm ▸ rfl)
+
+end ring
+
+section is_R_or_C
+
+variables [is_R_or_C α] [is_R_or_C β]
+
+/-- A matrix is hermitian iff the corresponding linear map is self adjoint. -/
+lemma is_hermitian_iff_is_self_adjoint [fintype n] [decidable_eq n] {A : matrix n n α} :
+  is_hermitian A ↔ inner_product_space.is_self_adjoint
+    ((pi_Lp.linear_equiv α (λ _ : n, α)).symm.conj A.to_lin' : module.End α (pi_Lp 2 _)) :=
+begin
+  rw [inner_product_space.is_self_adjoint, (pi_Lp.equiv 2 (λ _ : n, α)).symm.surjective.forall₂],
+  simp only [linear_equiv.conj_apply, linear_map.comp_apply, linear_equiv.coe_coe,
+    pi_Lp.linear_equiv_apply, pi_Lp.linear_equiv_symm_apply, linear_equiv.symm_symm],
+  simp_rw [euclidean_space.inner_eq_star_dot_product, equiv.apply_symm_apply, to_lin'_apply,
+    star_mul_vec, dot_product_mul_vec],
+  split,
+  { rintro (h : Aᴴ = A) x y,
+    rw h },
+  { intro h,
+    ext i j,
+    simpa only [(pi.single_star i 1).symm, ← star_mul_vec, mul_one, dot_product_single,
+      single_vec_mul, star_one, one_mul] using
+        h (@pi.single _ _ _ (λ i, add_zero_class.to_has_zero α) i 1)
+          (@pi.single _ _ _ (λ i, add_zero_class.to_has_zero α) j 1) }
+end
+
+end is_R_or_C
+
+end matrix

src/linear_algebra/matrix/spectrum.lean

@@ -0,0 +1,89 @@
+/-
+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.spectrum
+import linear_algebra.matrix.hermitian
+
+/-! # Spectral theory of hermitian matrices
+
+This file proves the spectral theorem for matrices. The proof of the spectral theorem is based on
+the spectral theorem for linear maps (`diagonalization_basis_apply_self_apply`).
+
+## Tags
+
+spectral theorem, diagonalization theorem
+
+-/
+
+namespace matrix
+
+variables {𝕜 : Type*} [is_R_or_C 𝕜] [decidable_eq 𝕜] {n : Type*} [fintype n] [decidable_eq n]
+variables {A : matrix n n 𝕜}
+
+open_locale matrix
+
+namespace is_hermitian
+
+variables (hA : A.is_hermitian)
+
+/-- The eigenvalues of a hermitian matrix, indexed by `fin (fintype.card n)` where `n` is the index
+type of the matrix. -/
+noncomputable def eigenvalues₀ : fin (fintype.card n) → ℝ :=
+@inner_product_space.is_self_adjoint.eigenvalues 𝕜 _ _ (pi_Lp 2 (λ (_ : n), 𝕜)) _ A.to_lin'
+  (is_hermitian_iff_is_self_adjoint.1 hA) _ (fintype.card n) finrank_euclidean_space
+
+/-- The eigenvalues of a hermitian matrix, reusing the index `n` of the matrix entries. -/
+noncomputable def eigenvalues : n → ℝ :=
+λ i, hA.eigenvalues₀ $ (fintype.equiv_of_card_eq (fintype.card_fin _)).symm i
+
+/-- A choice of an orthonormal basis of eigenvectors of a hermitian matrix. -/
+noncomputable def eigenvector_basis : basis n 𝕜 (n → 𝕜) :=
+(@inner_product_space.is_self_adjoint.eigenvector_basis 𝕜 _ _
+    (pi_Lp 2 (λ (_ : n), 𝕜)) _ A.to_lin' (is_hermitian_iff_is_self_adjoint.1 hA) _
+    (fintype.card n) finrank_euclidean_space).reindex
+  (fintype.equiv_of_card_eq (fintype.card_fin _))
+
+/-- A matrix whose columns are an orthonormal basis of eigenvectors of a hermitian matrix. -/
+noncomputable def eigenvector_matrix : matrix n n 𝕜 :=
+(pi.basis_fun 𝕜 n).to_matrix (eigenvector_basis hA)
+
+/-- The inverse of `eigenvector_matrix` -/
+noncomputable def eigenvector_matrix_inv : matrix n n 𝕜 :=
+(eigenvector_basis hA).to_matrix (pi.basis_fun 𝕜 n)
+
+lemma eigenvector_matrix_mul_inv :
+  hA.eigenvector_matrix ⬝ hA.eigenvector_matrix_inv = 1 :=
+by apply basis.to_matrix_mul_to_matrix_flip
+
+/-- *Diagonalization theorem*, *spectral theorem* for matrices; A hermitian matrix can be
+diagonalized by a change of basis.
+
+For the spectral theorem on linear maps, see `diagonalization_basis_apply_self_apply`. -/
+theorem spectral_theorem :
+  hA.eigenvector_matrix_inv ⬝ A =
+    diagonal (coe ∘ hA.eigenvalues) ⬝ hA.eigenvector_matrix_inv :=
+begin
+  rw [eigenvector_matrix_inv, basis_to_matrix_basis_fun_mul],
+  ext i j,
+  convert @inner_product_space.is_self_adjoint.diagonalization_basis_apply_self_apply 𝕜 _ _
+    (pi_Lp 2 (λ (_ : n), 𝕜)) _ A.to_lin' (is_hermitian_iff_is_self_adjoint.1 hA) _ (fintype.card n)
+    finrank_euclidean_space (euclidean_space.single j 1)
+    ((fintype.equiv_of_card_eq (fintype.card_fin _)).symm i),
+  { rw [eigenvector_basis, inner_product_space.is_self_adjoint.diagonalization_basis,
+      to_lin'_apply],
+    simp only [basis.to_matrix, basis.coe_to_orthonormal_basis_repr, basis.equiv_fun_apply],
+    simp_rw [basis.reindex_repr, euclidean_space.single, pi_Lp.equiv_symm_apply', mul_vec_single,
+      mul_one],
+    refl },
+  { simp only [diagonal_mul, (∘), eigenvalues, eigenvector_basis,
+      inner_product_space.is_self_adjoint.diagonalization_basis],
+    rw [basis.to_matrix_apply, basis.coe_to_orthonormal_basis_repr, basis.reindex_repr,
+      basis.equiv_fun_apply, pi.basis_fun_apply, eigenvalues₀, linear_map.coe_std_basis,
+      euclidean_space.single, pi_Lp.equiv_symm_apply'] }
+end
+
+end is_hermitian
+
+end matrix

src/linear_algebra/matrix/to_lin.lean

@@ -172,10 +172,14 @@ variables {l m n : Type*}
 
 variables [fintype n] [decidable_eq n]
 
-@[simp] lemma matrix.mul_vec_std_basis (M : matrix m n R) (i j) :
+lemma matrix.mul_vec_std_basis (M : matrix m n R) (i j) :
   M.mul_vec (std_basis R (λ _, R) j 1) i = M i j :=
 (congr_fun (matrix.mul_vec_single _ _ (1 : R)) i).trans $ mul_one _
 
+@[simp] lemma matrix.mul_vec_std_basis_apply (M : matrix m n R) (j) :
+  M.mul_vec (std_basis R (λ _, R) j 1) = Mᵀ j :=
+funext $ λ i, matrix.mul_vec_std_basis M i j
+
 /-- Linear maps `(n → R) →ₗ[R] (m → R)` are linearly equivalent to `matrix m n R`. -/
 def linear_map.to_matrix' : ((n → R) →ₗ[R] (m → R)) ≃ₗ[R] matrix m n R :=
 { to_fun := λ f i j, f (std_basis R (λ _, R) j 1) i,