Theory Gauss_Jordan.Inverse
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 ⟶ ((i = j ⟶ A $ i $ j = 1) ∧ (i ≠ j ⟶ A $ i $ j = 0)))"
lemma id_upt_nrows_mat_1:
assumes "id_upt_k A (nrows A)"
shows "A = mat 1"
unfolding mat_def apply vector using assms unfolding id_upt_k_def nrows_def
using to_nat_less_card[where ?'a='b]
by presburger
subsection‹Computing the inverse of a matrix using the Gauss Jordan algorithm›
text‹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"}›
lemma id_upt_k_Gauss_Jordan:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
assumes inv_A: "invertible A"
shows "id_upt_k (Gauss_Jordan A) k"
proof (induct k)
case 0
show ?case unfolding id_upt_k_def by fast
next
case (Suc k)
note id_k=Suc.hyps
have rref_k: "reduced_row_echelon_form_upt_k (Gauss_Jordan A) k" using rref_implies_rref_upt[OF rref_Gauss_Jordan] .
have rref_suc_k: "reduced_row_echelon_form_upt_k (Gauss_Jordan A) (Suc k)" using rref_implies_rref_upt[OF rref_Gauss_Jordan] .
have inv_gj: "invertible (Gauss_Jordan A)" by (metis inv_A invertible_Gauss_Jordan invertible_mult)
show "id_upt_k (Gauss_Jordan A) (Suc k)"
proof (unfold id_upt_k_def, auto)
fix j::'n
assume j_less_suc: "to_nat j < Suc k"
have greatest_prop: "j ≠ 0 ⟹ to_nat j = k ⟹ (GREATEST m. ¬ is_zero_row_upt_k m k (Gauss_Jordan A)) = j - 1"
proof (rule Greatest_equality)
assume 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
show Gauss_jj_1: "Gauss_Jordan A $ j $ j = 1"
proof (cases "j=0")
case 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
next
case False note j_not_zero = False
show ?thesis
proof (cases "to_nat j < k")
case True thus ?thesis using id_k unfolding id_upt_k_def by presburger
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
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)"
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)
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
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"
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
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_def j_less_k zero_neq_one)
ultimately show ?thesis using rref_upt_condition4[OF rref_k] i_not_j by fastforce
qed
qed
qed
qed
lemma 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
lemma matrix_inv_Gauss:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
assumes inv_A: "invertible A" and Gauss_eq: "Gauss_Jordan A = P ** A"
shows "matrix_inv A = P"
proof (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)
qed
lemma matrix_inv_Gauss_Jordan_PA:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
assumes inv_A: "invertible A"
shows "matrix_inv A = fst (Gauss_Jordan_PA A)"
by (metis Gauss_Jordan_PA_eq fst_Gauss_Jordan_PA inv_A matrix_inv_Gauss)
lemma invertible_eq_full_rank[code_unfold]:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
shows "invertible A = (rank A = nrows A)"
by (metis full_rank_implies_invertible invertible_implies_full_rank)
definition "inverse_matrix A = (if invertible A then Some (matrix_inv A) else None)"
lemma the_inverse_matrix:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
assumes "invertible A"
shows "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)
lemma inverse_matrix:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
shows "inverse_matrix A = (if invertible A then Some (P_Gauss_Jordan A) else None)"
by (metis inverse_matrix_def option.sel the_inverse_matrix)
lemma inverse_matrix_code[code_unfold]:
fixes A::"'a::{field}^'n::{mod_type}^'n::{mod_type}"
shows "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)"
unfolding inverse_matrix
unfolding invertible_eq_full_rank
unfolding rank_Gauss_Jordan_code
unfolding P_Gauss_Jordan_def
unfolding Let_def Gauss_Jordan_PA_eq by presburger
end