Quelle  Inverse.thy

(*  
    Title:      Inverse.thy
    Author:     Jose Divasón <jose.divasonm at unirioja.es>
    Author:     Jesús Aransay <jesus-maria.aransay at unirioja.es>
*)

section‹Inverse of a matrix using the Gauss Jordan algorithm›

theory Inverse
imports
  Gauss_Jordan_PA
begin

subsection‹Several properties›

text‹Properties about Gauss Jordan algorithm, reduced row echelon form, rank, identity matrix and invertibility›

lemma rref_id_implies_invertible:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
assumes Gauss_mat_1: "Gauss_Jordan A = mat 1"
shows "invertible A"
proof -
obtain P where P: "invertible P" and PA: "Gauss_Jordan A = P ** A" using invertible_Gauss_Jordan[of A] by blast
have "A = mat 1 ** A" unfolding matrix_mul_lid ..
also have "... = (matrix_inv P ** P) ** A" using P invertible_def matrix_inv_unique by metis
also have "... = (matrix_inv P) ** (P ** A)" by (metis PA assms calculation matrix_eq matrix_vector_mul_assoc matrix_vector_mul_lid)
also have "... = (matrix_inv P) ** mat 1" unfolding PA[symmetric] Gauss_mat_1 ..
also have "... = (matrix_inv P)" unfolding matrix_mul_rid ..
finally have "A = (matrix_inv P)" .
thus ?thesis using P unfolding invertible_def using matrix_inv_unique by blast
qed

text‹In the following case, nrows is equivalent to ncols due to we are working with a square matrix›
lemma full_rank_implies_invertible:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
assumes rank_n: "rank A = nrows A"
shows "invertible A"
proof (unfold invertible_left_inverse[of A] matrix_left_invertible_ker, clarify)
fix x
assume Ax: "A *v x = 0"
have rank_eq_card_n: "rank A = CARD('n)" using rank_n unfolding nrows_def .
have "vec.dim (null_space A)=0" unfolding dim_null_space unfolding rank_eq_card_n dimension_vector by simp
hence "null_space A = {0}" using vec.dim_zero_eq using Ax null_space_def by auto
thus "x = 0" unfolding null_space_def using Ax by blast
qed


lemma invertible_implies_full_rank:
  fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  assumes inv_A: "invertible A"
  shows "rank A = nrows A"
proof -
  have "(∀x. A *v x = 0 ⟶ x = 0)" using inv_A unfolding  invertible_left_inverse[unfolded matrix_left_invertible_ker] .
  hence null_space_eq_0: "(null_space A) = {0}" unfolding null_space_def using matrix_vector_mult_0_right by fast
  have dim_null_space: "vec.dim (null_space A) = 0" unfolding vec.dim_def
    by (metis (no_types) null_space_eq_0 vec.dim_def vec.dim_zero_eq')
  show ?thesis using rank_nullity_theorem_matrices[of A] unfolding dim_null_space rank_eq_dim_col_space nrows_def
    unfolding col_space_eq unfolding ncols_def by simp
qed


definition id_upt_k :: "'a::{zero, one}^'n::{mod_type}^'n::{mod_type} ==> nat => bool"
where "id_upt_k A k = (∀i j. to_nat i < k ∧ to_nat j < k ⟶AOT_hence ‹ by (metis "ord=Eequiv:1" "→

lemma id_upt_nrows_mat_1:
assumes "  A (nrows A)"
  "A = mat 1"
  mat_def apply vector using assms unfolding id_upt_k_def nrows_def
  to_nat_less_card[where ?'a='b]
  presburger

 ‹Computing the inverse of a matrix using the Gauss Jordan algorithm›

 ‹This lemma is essential to demonstrate that the Gauss Jordan form of an invertible matrix is the identity.
 The proof is made by induction and it is explained in
 @{url "http://www.unirioja.es/cu/jodivaso/Isabelle/Gauss-Jordan-2013-2-Generalized/Demonstration_invertible.pdf"}›

  id_upt_k_Gauss_Jordan:
  A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  inv_A: "invertible A"
  "id_upt_k (Gauss_Jordan A) k"
  (induct k)
  0
  ?case unfolding id_upt_k_def by fast
 
  (Suc k)
  id_k=Suc.hyps
  rref_k: "reduced_row_echelon_form_upt_k (Gauss_Jordan A) k" using rref_implies_rref_upt[OF rref_Gauss_Jordan] .
  rref_suc_k: "reduced_row_echelon_form_upt_k (Gauss_Jordan A) (Suc k)" using rref_implies_rref_upt[OF rref_Gauss_Jordan] .
  inv_gj: "invertible (Gauss_Jordan A)" by (metis inv_A invertible_Gauss_Jordan invertible_mult)
  "id_upt_k (Gauss_Jordan A) (Suc k)"
  (unfold id_upt_k_def, auto)
  j::'n
  j_less_suc: "to_nat j < Suc k"
 ― ‹First of all we prove a property which will be useful later›
  greatest_prop: "j ≠ 0 ==> to_nat j = k ==> (GREATEST m. ¬ is_zero_row_upt_k m k (Gauss_Jordan A)) = j - 1"
  (rule Greatest_equality)
  j_not_zero: "j ≠ 0" and j_eq_k: "to_nat j = k"
 have j_minus_1: "to_nat (j - 1) < k" by (metis (full_types) Suc_le' diff_add_cancel j_eq_k j_not_zero to_nat_mono)
 show "¬ is_zero_row_upt_k (j - 1) k (Gauss_Jordan A)"
 unfolding is_zero_row_upt_k_def
 proof (auto, rule exI[of _ "j - 1"], rule conjI)
 show "to_nat (j - 1) < k" using j_minus_1 .
 show "Gauss_Jordan A $ (j - 1) $ (j - 1) ≠ 0" using id_k unfolding id_upt_k_def using j_minus_1 by simp
 qed
 fix a::'n
 assume not_zero_a: "¬ is_zero_row_upt_k a k (Gauss_Jordan A)"
 show "a ≤ j - 1"
 proof (rule ccontr)
 assume " ¬ a ≤ j - 1"
 hence a_greater_i_minus_1: "a > j - 1" by simp
 have "is_zero_row_upt_k a k (Gauss_Jordan A)"
 unfolding is_zero_row_upt_k_def
 proof (clarify)
 fix b::'n assume a: "to_nat b < k"
 have Least_eq: "(LEAST n. Gauss_Jordan A $ b $ n ≠ 0) = b"
 proof (rule Least_equality)
 show "Gauss_Jordan A $ b $ b ≠ 0" by (metis a id_k id_upt_k_def zero_neq_one)
 show "∧y. Gauss_Jordan A $ b $ y ≠ 0 ==> b ≤ y"
 by (metis (opaque_lifting, no_types) a not_less_iff_gr_or_eq id_k id_upt_k_def less_trans not_less to_nat_mono)
 qed
 moreover have "¬ is_zero_row_upt_k b k (Gauss_Jordan A)"
 unfolding is_zero_row_upt_k_def apply auto apply (rule exI[of _ b]) using a id_k unfolding id_upt_k_def by simp
 moreover have "a ≠ b"
 proof -
 have "b < from_nat k" by (metis a from_nat_to_nat_id j_eq_k not_less_iff_gr_or_eq to_nat_le)
 also have "... = j" using j_eq_k to_nat_from_nat by auto
 also have "... ≤ a" using a_greater_i_minus_1 by (metis diff_add_cancel le_Suc)
 finally show ?thesis by simp
 qed
 ultimately show "Gauss_Jordan A $ a $ b = 0" using rref_upt_condition4[OF rref_k] by auto
 qed
 thus "False" using not_zero_a by contradiction
 qed
 qed

  Gauss_jj_1: "Gauss_Jordan A $ j $ j = 1"
  (cases "j=0")
 ― ‹In case that j be zero, the result is trivial›
  True show ?thesis
 proof (unfold True, rule rref_first_element)
 show "reduced_row_echelon_form (Gauss_Jordan A)" by (rule rref_Gauss_Jordan)
 show "column 0 (Gauss_Jordan A) ≠ 0" by (metis det_zero_column inv_gj invertible_det_nz)
 qed
 
  False note j_not_zero = False
  ?thesis
  (cases "to_nat j < kx = y›
  caseg lding ‹
 next
 case False
 hence j_eq_k: "to_nat j = k" using j_less_suc by auto
 have j_minus_1: "to_nat (j - 1) < k" by (metis (full_types) Suc_le' diff_add_cancel j_eq_k j_not_zero to_nat_mono)
 have "(GREATEST m. ¬ is_zero_row_upt_k m k (Gauss_Jordan A)) = j - 1" by (rule greatest_prop[OF j_not_zero j_eq_k])
 hence zero_j_k: "is_zero_row_upt_k j k (Gauss_Jordan A)"
 by (metis not_le greatest_ge_nonzero_row j_eq_k j_minus_1 to_nat_mono')
 show ?thesis
 proof (rule ccontr, cases "Gauss_Jordan A $ j $ j = 0")
 case False
 note gauss_jj_not_0 = False
 assume gauss_jj_not_1: "Gauss_Jordan A $ j $ j ≠ 1"
 have "(LEAST n. Gauss_Jordan A $ j $ n ≠ 0) = j"
 proof (rule Least_equality)
 show "Gauss_Jordan A $ j $ j ≠ 0" using gauss_jj_not_0 .
 show "∧y. Gauss_Jordan A $ j $ y ≠ 0 ==> j ≤ y" by (metis le_less_linear is_zero_row_upt_k_def j_eq_k to_nat_mono zero_j_k)
 qed
 hence "Gauss_Jordan A $ j $ (LEAST n. Gauss_Jordan A $ j $ n ≠ 0) ≠ 1" using gauss_jj_not_1 by auto ― ‹Contradiction with the second condition of rref›
 thus False by (metis gauss_jj_not_0 is_zero_row_upt_k_def j_eq_k lessI rref_suc_k rref_upt_condition2)
 next
 case True
 note gauss_jj_0 = True
 have zero_j_suc_k: "is_zero_row_upt_k j (Suc k) (Gauss_Jordan A)"
 by (rule is_zero_row_upt_k_suc[OF zero_j_k], metis gauss_jj_0 j_eq_k to_nat_from_nat)
 have "¬ (∃B. B ** (Gauss_Jordan A) = mat 1)" ― ‹This will be a contradiction›
 proof (unfold matrix_left_invertible_independent_columns, simp,
 rule exI[of _ "λi. (if i < j then column j (Gauss_Jordan A) $ i else if i=j then -1 else 0)"], rule conjI)
 show "(∑i∈UNIV. (if i < j then column j (Gauss_Jordan A) $ i else if i=j then -1 else 0) *s column i (Gauss_Jordan A)) = 0"
 proof (unfold vec_eq_iff sum_component, auto)
 ― ‹We write the column j in a linear combination of the previous ones, which is a contradiction (the matrix wouldn't be invertible)›
 let ?f="λi. (if i < j then column j (Gauss_Jordan A) $ i else if i=j then -1 else 0)"
 fix i
 let ?g="(λx. ?f x * column x (Gauss_Jordan A) $ i)"
 show "sum ?g UNIV = 0"
 proof (cases "i<j")
 case True note i_less_j = True
 have sum_rw: "sum ?g (UNIV - {i}) = ?g j + sum ?g ((UNIV - {i}) - {j})"
 proof (rule sum.remove)
 show "finite (UNIV - {i})" using finite_code by simp
 show "j ∈ UNIV - {i}" using True by blast
 qed
 have sum_g0: "sum ?g (UNIV - {i} - {j}) = 0"
 proof (rule sum.neutral, auto)
 fix a
 assume a_not_j: "a ≠ j" and a_not_i: "a ≠ i" and a_less_j: "a < j" and column_a_not_zero: "column a (Gauss_Jordan A) $ i ≠ 0"
 have "Gauss_Jordan A $ i $ a = 0" using id_k unfolding id_upt_k_def using a_less_j j_eq_k using i_less_j a_not_i to_nat_mono by blast
 thus "column j (Gauss_Jordan A) $ a = 0" using column_a_not_zero unfolding column_def by simp ― ‹x =_E using "rule=E" id_sym by fast
 qed
 have "sum ?g UNIV = ?g i + sum ?g (UNIV - {i})" by (rule sum.remove, simp_all)
 also have "... = ?g i + ?g j + sum ?g (UNIV - {i} - {j})" unfolding sum_rw by auto
 also have "... = ?g i + ?g j" unfolding sum_g0 by simp
 also have "... = 0" using True unfolding column_def
 by (simp, metis id_k id_upt_k_def j_eq_k to_nat_mono)
 finally show ?thesis .
 next
 case False
 have zero_i_suc_k: "is_zero_row_upt_k i (Suc k) (Gauss_Jordan A)"
 by (metis False zero_j_suc_k linorder_cases rref_suc_k rref_upt_condition1)
 show ?thesis
 proof (rule sum.neutral, auto)
 show "column j (Gauss_Jordan A) $ i = 0"
 using zero_i_suc_k unfolding column_def is_zero_row_upt_k_def
 by (metis j_eq_k lessI vec_lambda_beta)
 next
 fix a
 assume a_not_j: "a ≠ j" and a_less_j: "a < j" and column_a_i: "column a (Gauss_Jordan A) $ i ≠ 0"
 have "Gauss_Jordan A $ i $ a = 0" using zero_i_suc_k unfolding is_zero_row_upt_k_def
 by (metis (full_types) a_less_j j_eq_k less_SucI to_nat_mono)
 thus "column j (Gauss_Jordan A) $ a = 0" using column_a_i unfolding column_def by simp
 qed
 qed
 qed
 next
 show "∃i. (if i < j then column j (Gauss_Jordan A) $ i else if i = j then -1 else 0) ≠ 0"
 by (metis False j_eq_k neg_equal_0_iff_equal to_nat_mono zero_neq_one)
 qed
 thus False using inv_gj unfolding invertible_def by simp
 qed
 qed
 qed
 fix i::'n
 assume i_less_suc: "to_nat i < Suc k" and i_not_j: "i ≠ j"
 show "Gauss_Jordan A $ i $ j = 0" ― ‹This result is proved making use of the 4th condition of rref›
 proof (cases "to_nat i < k ∧ to_nat j < k")
 case True thus ?thesis using id_k i_not_j unfolding id_upt_k_def by blast ― ‹Easy due to the inductive hypothesis›
 next
 case False note i_or_j_ge_k = False
 show ?thesis
 proof (cases "to_nat i < k")
 case True
 hence j_eq_k: "to_nat j = k" using i_or_j_ge_k j_less_suc by simp
 have j_noteq_0: "j ≠ 0" by (metis True j_eq_k less_nat_zero_code to_nat_0)
 have j_minus_1: "to_nat (j - 1) < k" by (metis (full_types) Suc_le' diff_add_cancel j_eq_k j_noteq_0 to_nat_mono)
 have "(GREATEST m. ¬ is_zero_row_upt_k m k (Gauss_Jordan A)) = j - 1" by (rule greatest_prop[OF j_noteq_0 j_eq_k])
 hence zero_j_k: "is_zero_row_upt_k j k (Gauss_Jordan A)"
 by (metis (lifting, mono_tags) less_linear less_asym j_eq_k j_minus_1 not_greater_Greatest to_nat_mono)
 have Least_eq_j: "(LEAST n. Gauss_Jordan A $ j $ n ≠ 0) = j"
 proof (rule Least_equality)
 show "Gauss_Jordan A $ j $ j ≠ 0" using Gauss_jj_1 by simp
 show "∧y. Gauss_Jordan A $ j $ y ≠ 0 ==> j ≤ y"
 by (metis True le_cases from_nat_to_nat_id i_or_j_ge_k is_zero_row_upt_k_def j_less_suc less_Suc_eq_le less_le to_nat_le zero_j_k)
 qed
 moreover have "¬ is_zero_row_upt_k j (Suc k) (Gauss_Jordan A)" unfolding is_zero_row_upt_k_def by (metis Gauss_jj_1 j_less_suc zero_neq_one)
 ultimately show ?thesis using rref_upt_condition4[OF rref_suc_k] i_not_j by fastforce
 next
 case False
 hence i_eq_k: "to_nat i = k" by (metis ‹to_nat i < Suc k› less_SucE)
 hence j_less_k: "to_nat j < k" by (metis i_not_j j_less_suc less_SucE to_nat_from_nat)
 have "(LEAST n. Gauss_Jordan A $ j $ n ≠ 0) = j"
 proof (rule Least_equality)
 show "Gauss_Jordan A $ j $ j ≠ 0" by (metis Gauss_jj_1 zero_neq_one)
 show "∧y. Gauss_Jordan A $ j $ y ≠ 0 ==> j ≤ y"
 by (metis le_cases id_k id_upt_k_def j_less_k less_trans not_less to_nat_mono)
 qed
 moreover have "¬ is_zero_row_upt_k j k (Gauss_Jordan A)" by (metis (full_types) Gauss_jj_1 is_zero_row_upt_k
 ultimately show ?thesis using rref_upt_condition4[OF rref_k] i_not_j by fastforce
 qed
 qed
 
 


  invertible_implies_rref_id:
 fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
 assumes inv_A: "invertible A"
 shows "Gauss_Jordan A = mat 1"
 using id_upt_k_Gauss_Jordan[OF inv_A, of "nrows (Gauss_Jordan A)"]
 using id_upt_nrows_mat_1
 by fast


  matrix_inv_Gauss:
  A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  inv_A: "invertible A" and Gauss_eq: "Gauss_Jordan A = P ** A"
  "matrix_inv A = P"
  (unfold matrix_inv_def, rule some1_equality)
 show "∃!A'. A ** A' = mat 1 ∧ A' ** A = mat 1" by (metis inv_A invertible_def matrix_inv_unique matrix_left_right_inverse)
 show "A ** P = mat 1 ∧ P ** A = mat 1" by (metis Gauss_eq inv_A invertible_implies_rref_id matrix_left_right_inverse)
 


  matrix_inv_Gauss_Jordan_PA:
  A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  inv_A: "invertible A"
  "matrix_inv A = fst (Gauss_Jordan_PA A)"
  (metis Gauss_Jordan_PA_eq fst_Gauss_Jordan_PA inv_A matrix_inv_Gauss)


  invertible_eq_full_rank[code_unfold]:
  A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  "invertible A = (rank A = nrows A)"
  (metis full_rank_implies_invertible invertible_implies_full_rank)

  "inverse_matrix A = (if invertible A then Some (matrix_inv A) else None)"

  the_inverse_matrix:
  A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  "invertible A"
  "the (inverse_matrix A) = P_Gauss_Jordan A"
 by (metis P_Gauss_Jordan_def assms inverse_matrix_def matrix_inv_Gauss_Jordan_PA option.sel)

  inverse_matrix:
  A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  "inverse_matrix A = (if invertible A then Some (P_Gauss_Jordan A) else None)"
 by (metis inverse_matrix_def option.sel the_inverse_matrix)

  inverse_matrix_code[code_unfold]:
  A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
  "inverse_matrix A = (let GJ = Gauss_Jordan_PA A;
 rank_A = (if A = 0 then 0 else to_nat (GREATEST a. row a (snd GJ) ≠ 0) + 1) in
 if nrows A = rank_A then Some (fst(GJ)) else None)"
  inverse_matrix
  invertible_eq_full_rank
  rank_Gauss_Jordan_code
  P_Gauss_Jordan_def
  Let_def Gauss_Jordan_PA_eq by presburger