imports Perron_Frobenius_Irreducible

(* Author: Thiemann *) subsection ‹Handling Non-Irreducible Matrices as Well› theory Perron_Frobenius_General imports Perron_Frobenius_Irreducible begin text ‹We will need to take sub-matrices and permutations of matrices where the former can best be done via JNF-matrices. So, we first need the Perron-Frobenius theorem in the JNF-world. So, we first define irreducibility of a JNF-matrix.› definition graph_of_mat where "graph_of_mat A = (let n = dim_row A; U = {..<n} in { ij. A $$ ij ≠ 0} ∩ U × U)" definition irreducible_mat where "irreducible_mat A = (let n = dim_row A in (∀ i j. i < n ⟶ j < n ⟶ (i,j) ∈ (graph_of_mat A)^+))" definition "nonneg_irreducible_mat A = (nonneg_mat A ∧ irreducible_mat A)" text ‹Next, we have to install transfer rules› context includes lifting_syntax begin lemma HMA_irreducible[transfer_rule]: "((HMA_M :: _ ⇒ _ ^ 'n ^ 'n ⇒ _) ===> (=)) irreducible_mat fixed_mat.irreducible" proof (intro rel_funI, goal_cases) case (1 a A) interpret fixed_mat A . let ?t = "Bij_Nat.to_nat :: 'n ⇒ nat" let ?f = "Bij_Nat.from_nat :: nat ⇒ 'n" from 1[unfolded HMA_M_def] have a: "a = from_hma⇩_{m}A" (is "_ = ?A") by auto let ?n = "CARD('n)" have dim: "dim_row a = ?n" unfolding a by simp have id: "{..<?n} = {0..<?n}" by auto have Aij: "A $ i $ j = ?A $$ (?t i, ?t j)" for i j by (metis (no_types, lifting) to_hma⇩_{m}_def to_hma_from_hma⇩_{m}vec_lambda_beta) have graph: "graph_of_mat a = {(?t i,?t j) | i j. A $ i $ j ≠ 0}" (is "?G = _") unfolding graph_of_mat_def dim Let_def id range_to_nat[symmetric] unfolding a Aij by auto have "irreducible_mat a = (∀i j. i ∈ range ?t ⟶ j ∈ range ?t ⟶ (i,j) ∈ ?G^+)" unfolding irreducible_mat_def dim Let_def range_to_nat by auto also have "… = (∀ i j. (?t i, ?t j) ∈ ?G^+)" by auto also note part1 = calculation have G: "?G = map_prod ?t ?t ` G" unfolding graph G_def by auto have part2: "(?t i, ?t j) ∈ ?G^+ ⟷ (i,j) ∈ G^+" for i j unfolding G by (rule inj_trancl_image, simp add: inj_on_def) show ?case unfolding part1 part2 irreducible_def by auto qed lemma HMA_nonneg_irreducible_mat[transfer_rule]: "(HMA_M ===> (=)) nonneg_irreducible_mat perron_frobenius" unfolding perron_frobenius_def pf_nonneg_mat_def perron_frobenius_axioms_def nonneg_irreducible_mat_def by transfer_prover end text ‹The main statements of Perron-Frobenius can now be transferred to JNF-matrices› lemma perron_frobenius_irreducible: fixes A :: "real Matrix.mat" and cA :: "complex Matrix.mat" assumes A: "A ∈ carrier_mat n n" and n: "n ≠ 0" and nonneg: "nonneg_mat A" and irr: "irreducible_mat A" and cA: "cA = map_mat of_real A" and sr: "sr = Spectral_Radius.spectral_radius cA" shows "eigenvalue A sr" "order sr (char_poly A) = 1" "0 < sr" "eigenvalue cA α ⟹ cmod α ≤ sr" "eigenvalue cA α ⟹ cmod α = sr ⟹ order α (char_poly cA) = 1" "∃ k f. k ≠ 0 ∧ k ≤ n ∧ char_poly A = (monom 1 k - [:sr ^ k:]) * f ∧ (∀x. poly (map_poly complex_of_real f) x = 0 ⟶ cmod x < sr)" proof (atomize (full), goal_cases) case 1 from nonneg irr have irr: "nonneg_irreducible_mat A" unfolding nonneg_irreducible_mat_def by auto note main = perron_frobenius.pf_main_connect[untransferred, cancel_card_constraint, OF A irr, folded sr cA] from main(5,6,7)[OF refl refl n] have "∃ k f. k ≠ 0 ∧ k ≤ n ∧ char_poly A = (monom 1 k - [:sr ^ k:]) * f ∧ (∀x. poly (map_poly complex_of_real f) x = 0 ⟶ cmod x < sr)" by blast with main(1,3,8)[OF n] main(2)[OF _ n] main(4)[OF _ _ n] show ?case by auto qed text ‹We now need permutations on matrices to show that a matrix if a matrix is not irreducible, then it can be turned into a four-block-matrix by a permutation, where the lower left block is 0.› definition permutation_mat :: "nat ⇒ (nat ⇒ nat) ⇒ 'a :: semiring_1 mat" where "permutation_mat n p = Matrix.mat n n (λ (i,j). (if i = p j then 1 else 0))" no_notation m_inv ("invı _" [81] 80) lemma permutation_mat_dim[simp]: "permutation_mat n p ∈ carrier_mat n n" "dim_row (permutation_mat n p) = n" "dim_col (permutation_mat n p) = n" unfolding permutation_mat_def by auto lemma permutation_mat_row[simp]: "p permutes {..<n} ⟹ i < n ⟹ Matrix.row (permutation_mat n p) i = unit_vec n (inv p i)" unfolding permutation_mat_def unit_vec_def by (intro eq_vecI, auto simp: permutes_inverses) lemma permutation_mat_col[simp]: "p permutes {..<n} ⟹ i < n ⟹ Matrix.col (permutation_mat n p) i = unit_vec n (p i)" unfolding permutation_mat_def unit_vec_def by (intro eq_vecI, auto simp: permutes_inverses) lemma permutation_mat_left: assumes A: "A ∈ carrier_mat n nc" and p: "p permutes {..<n}" shows "permutation_mat n p * A = Matrix.mat n nc (λ (i,j). A $$ (inv p i, j))" proof - { fix i j assume ij: "i < n" "j < nc" from p ij(1) have i: "inv p i < n" by (simp add: permutes_def) have "(permutation_mat n p * A) $$ (i,j) = scalar_prod (unit_vec n (inv p i)) (col A j)" by (subst index_mult_mat, insert ij A p, auto) also have "… = A $$ (inv p i, j)" by (subst scalar_prod_left_unit, insert A ij i, auto) also note calculation } thus ?thesis using A by (intro eq_matI, auto) qed lemma permutation_mat_right: assumes A: "A ∈ carrier_mat nr n" and p: "p permutes {..<n}" shows "A * permutation_mat n p = Matrix.mat nr n (λ (i,j). A $$ (i, p j))" proof - { fix i j assume ij: "i < nr" "j < n" from p ij(2) have j: "p j < n" by (simp add: permutes_def) have "(A * permutation_mat n p) $$ (i,j) = scalar_prod (Matrix.row A i) (unit_vec n (p j))" by (subst index_mult_mat, insert ij A p, auto) also have "… = A $$ (i, p j)" by (subst scalar_prod_right_unit, insert A ij j, auto) also note calculation } thus ?thesis using A by (intro eq_matI, auto) qed lemma permutes_lt: "p permutes {..<n} ⟹ i < n ⟹ p i < n" by (meson lessThan_iff permutes_in_image) lemma permutes_iff: "p permutes {..<n} ⟹ i < n ⟹ j < n ⟹ p i = p j ⟷ i = j" by (metis permutes_inverses(2)) lemma permutation_mat_id_1: assumes p: "p permutes {..<n}" shows "permutation_mat n p * permutation_mat n (inv p) = 1⇩_{m}n" by (subst permutation_mat_left[OF _ p, of _ n], force, unfold permutation_mat_def, rule eq_matI, auto simp: permutes_lt[OF permutes_inv[OF p]] permutes_iff[OF permutes_inv[OF p]]) lemma permutation_mat_id_2: assumes p: "p permutes {..<n}" shows "permutation_mat n (inv p) * permutation_mat n p = 1⇩_{m}n" by (subst permutation_mat_right[OF _ p, of _ n], force, unfold permutation_mat_def, rule eq_matI, insert p, auto simp: permutes_lt[OF p] permutes_inverses) lemma permutation_mat_both: assumes A: "A ∈ carrier_mat n n" and p: "p permutes {..<n}" shows "permutation_mat n p * Matrix.mat n n (λ (i,j). A $$ (p i, p j)) * permutation_mat n (inv p) = A" unfolding permutation_mat_left[OF mat_carrier p] by (subst permutation_mat_right[OF _ permutes_inv[OF p], of _ n], force, insert A p, auto intro!: eq_matI simp: permutes_inverses permutes_lt[OF permutes_inv[OF p]]) lemma permutation_similar_mat: assumes A: "A ∈ carrier_mat n n" and p: "p permutes {..<n}" shows "similar_mat A (Matrix.mat n n (λ (i,j). A $$ (p i, p j)))" by (rule similar_matI[OF _ permutation_mat_id_1[OF p] permutation_mat_id_2[OF p] permutation_mat_both[symmetric, OF A p]], insert A, auto) lemma det_four_block_mat_lower_left_zero: fixes A1 :: "'a :: idom mat" assumes A1: "A1 ∈ carrier_mat n n" and A2: "A2 ∈ carrier_mat n m" and A30: "A3 = 0⇩_{m}m n" and A4: "A4 ∈ carrier_mat m m" shows "Determinant.det (four_block_mat A1 A2 A3 A4) = Determinant.det A1 * Determinant.det A4" proof - let ?det = Determinant.det let ?t = "transpose_mat" let ?A = "four_block_mat A1 A2 A3 A4" let ?k = "n + m" have A3: "A3 ∈ carrier_mat m n" unfolding A30 by auto have A: "?A ∈ carrier_mat ?k ?k" by (rule four_block_carrier_mat[OF A1 A4]) have "?det ?A = ?det (?t ?A)" by (rule sym, rule Determinant.det_transpose[OF A]) also have "?t ?A = four_block_mat (?t A1) (?t A3) (?t A2) (?t A4)" by (rule transpose_four_block_mat[OF A1 A2 A3 A4]) also have "?det … = ?det (?t A1) * ?det (?t A4)" by (rule det_four_block_mat_upper_right_zero[of _ n _ m], insert A1 A2 A30 A4, auto) also have "?det (?t A1) = ?det A1" by (rule Determinant.det_transpose[OF A1]) also have "?det (?t A4) = ?det A4" by (rule Determinant.det_transpose[OF A4]) finally show ?thesis . qed lemma char_poly_matrix_four_block_mat: assumes A1: "A1 ∈ carrier_mat n n" and A2: "A2 ∈ carrier_mat n m" and A3: "A3 ∈ carrier_mat m n" and A4: "A4 ∈ carrier_mat m m" shows "char_poly_matrix (four_block_mat A1 A2 A3 A4) = four_block_mat (char_poly_matrix A1) (map_mat (λ x. [:-x:]) A2) (map_mat (λ x. [:-x:]) A3) (char_poly_matrix A4)" proof - from A1 A4 have dim[simp]: "dim_row A1 = n" "dim_col A1 = n" "dim_row A4 = m" "dim_col A4 = m" by auto show ?thesis unfolding char_poly_matrix_def four_block_mat_def Let_def dim by (rule eq_matI, insert A2 A3, auto) qed lemma char_poly_four_block_mat_lower_left_zero: fixes A :: "'a :: idom mat" assumes A: "A = four_block_mat B C (0⇩_{m}m n) D" and B: "B ∈ carrier_mat n n" and C: "C ∈ carrier_mat n m" and D: "D ∈ carrier_mat m m" shows "char_poly A = char_poly B * char_poly D" unfolding A char_poly_def by (subst char_poly_matrix_four_block_mat[OF B C _ D], force, rule det_four_block_mat_lower_left_zero[of _ n _ m], insert B C D, auto) lemma elements_mat_permutes: assumes p: "p permutes {..< n}" and A: "A ∈ carrier_mat n n" and B: "B = Matrix.mat n n (λ (i,j). A $$ (p i, p j))" shows "elements_mat A = elements_mat B" proof - from A B have [simp]: "dim_row A = n" "dim_col A = n" "dim_row B = n" "dim_col B = n" by auto { fix i j assume ij: "i < n" "j < n" let ?p = "inv p" from permutes_lt[OF p] ij have *: "p i < n" "p j < n" by auto from permutes_lt[OF permutes_inv[OF p]] ij have **: "?p i < n" "?p j < n" by auto have "∃ i' j'. B $$ (i,j) = A $$ (i',j') ∧ i' < n ∧ j' < n" "∃ i' j'. A $$ (i,j) = B $$ (i',j') ∧ i' < n ∧ j' < n" by (rule exI[of _ "p i"], rule exI[of _ "p j"], insert ij *, simp add: B, rule exI[of _ "?p i"], rule exI[of _ "?p j"], insert ** p, simp add: B permutes_inverses) } thus ?thesis unfolding elements_mat by auto qed lemma elements_mat_four_block_mat_supseteq: assumes A1: "A1 ∈ carrier_mat n n" and A2: "A2 ∈ carrier_mat n m" and A3: "A3 ∈ carrier_mat m n" and A4: "A4 ∈ carrier_mat m m" shows "elements_mat (four_block_mat A1 A2 A3 A4) ⊇ (elements_mat A1 ∪ elements_mat A2 ∪ elements_mat A3 ∪ elements_mat A4)" proof let ?A = "four_block_mat A1 A2 A3 A4" have A: "?A ∈ carrier_mat (n + m) (n + m)" using A1 A2 A3 A4 by simp from A1 A4 have dim[simp]: "dim_row A1 = n" "dim_col A1 = n" "dim_row A4 = m" "dim_col A4 = m" by auto fix x assume x: "x ∈ elements_mat A1 ∪ elements_mat A2 ∪ elements_mat A3 ∪ elements_mat A4" { assume "x ∈ elements_mat A1" from this[unfolded elements_mat] A1 obtain i j where x: "x = A1 $$ (i,j)" and ij: "i < n" "j < n" by auto have "x = ?A $$ (i,j)" using ij unfolding x four_block_mat_def Let_def by simp from elements_matI[OF A _ _ this] ij have "x ∈ elements_mat ?A" by auto } moreover { assume "x ∈ elements_mat A2" from this[unfolded elements_mat] A2 obtain i j where x: "x = A2 $$ (i,j)" and ij: "i < n" "j < m" by auto have "x = ?A $$ (i,j + n)" using ij unfolding x four_block_mat_def Let_def by simp from elements_matI[OF A _ _ this] ij have "x ∈ elements_mat ?A" by auto } moreover { assume "x ∈ elements_mat A3" from this[unfolded elements_mat] A3 obtain i j where x: "x = A3 $$ (i,j)" and ij: "i < m" "j < n" by auto have "x = ?A $$ (i+n,j)" using ij unfolding x four_block_mat_def Let_def by simp from elements_matI[OF A _ _ this] ij have "x ∈ elements_mat ?A" by auto } moreover { assume "x ∈ elements_mat A4" from this[unfolded elements_mat] A4 obtain i j where x: "x = A4 $$ (i,j)" and ij: "i < m" "j < m" by auto have "x = ?A $$ (i+n,j + n)" using ij unfolding x four_block_mat_def Let_def by simp from elements_matI[OF A _ _ this] ij have "x ∈ elements_mat ?A" by auto } ultimately show "x ∈ elements_mat ?A" using x by blast qed lemma non_irreducible_mat_split: fixes A :: "'a :: idom mat" assumes A: "A ∈ carrier_mat n n" and not: "¬ irreducible_mat A" and n: "n > 1" shows "∃ n1 n2 B B1 B2 B4. similar_mat A B ∧ elements_mat A = elements_mat B ∧ B = four_block_mat B1 B2 (0⇩_{m}n2 n1) B4 ∧ B1 ∈ carrier_mat n1 n1 ∧ B2 ∈ carrier_mat n1 n2 ∧ B4 ∈ carrier_mat n2 n2 ∧ 0 < n1 ∧ n1 < n ∧ 0 < n2 ∧ n2 < n ∧ n1 + n2 = n" proof - from A have [simp]: "dim_row A = n" by auto let ?G = "graph_of_mat A" let ?reachp = "λ i j. (i,j) ∈ ?G^+" let ?reach = "λ i j. (i,j) ∈ ?G^*" have "∃ i j. i < n ∧ j < n ∧ ¬ ?reach i j" proof (rule ccontr) assume "¬ ?thesis" hence reach: "⋀ i j. i < n ⟹ j < n ⟹ ?reach i j" by auto from not[unfolded irreducible_mat_def Let_def] obtain i j where i: "i < n" and j: "j < n" and nreach: "¬ ?reachp i j" by auto from reach[OF i j] nreach have ij: "i = j" by (simp add: rtrancl_eq_or_trancl) from n j obtain k where k: "k < n" and diff: "j ≠ k" by auto from reach[OF j k] diff reach[OF k j] have "?reachp j j" by (simp add: rtrancl_eq_or_trancl) with nreach ij show False by auto qed then obtain i j where i: "i < n" and j: "j < n" and nreach: "¬ ?reach i j" by auto define I where "I = {k. k < n ∧ ?reach i k}" have iI: "i ∈ I" unfolding I_def using nreach i by auto have jI: "j ∉ I" unfolding I_def using nreach j by auto define f where "f = (λ i. if i ∈ I then 1 else 0 :: nat)" let ?xs = "[0 ..< n]" from mset_eq_permutation[OF mset_sort, of ?xs f] obtain p where p: "p permutes {..< n}" and perm: "permute_list p ?xs = sort_key f ?xs" by auto from p have lt[simp]: "i < n ⟹ p i < n" for i by (rule permutes_lt) let ?p = "inv p" have ip: "?p permutes {..< n}" using permutes_inv[OF p] . from ip have ilt[simp]: "i < n ⟹ ?p i < n" for i by (rule permutes_lt) let ?B = "Matrix.mat n n (λ (i,j). A $$ (p i, p j))" define B where "B = ?B" from permutation_similar_mat[OF A p] have sim: "similar_mat A B" unfolding B_def . let ?ys = "permute_list p ?xs" define ys where "ys = ?ys" have len_ys: "length ys = n" unfolding ys_def by simp let ?k = "length (filter (λ i. f i = 0) ys)" define k where "k = ?k" have kn: "k ≤ n" unfolding k_def using len_ys using length_filter_le[of _ ys] by auto have ys_p: "i < n ⟹ ys ! i = p i" for i unfolding ys_def permute_list_def by simp have ys: "ys = map (λ i. ys ! i) [0 ..< n]" unfolding len_ys[symmetric] by (simp add: map_nth) also have "… = map p [0 ..< n]" by (rule map_cong, insert ys_p, auto) also have "[0 ..< n] = [0 ..< k] @ [k ..< n]" using kn using le_Suc_ex upt_add_eq_append by blast finally have ys: "ys = map p [0 ..< k] @ map p [k ..< n]" by simp { fix i assume i: "i < n" let ?g = "(λ i. f i = 0)" let ?f = "filter ?g" from i have pi: "p i < n" using p by simp have "k = length (?f ys)" by fact also have "?f ys = ?f (map p [0 ..< k]) @ ?f (map p [k ..< n])" unfolding ys by simp also note k = calculation finally have True by blast from perm[symmetric, folded ys_def] have "sorted (map f ys)" using sorted_sort_key by metis from this[unfolded ys map_append sorted_append set_map] have sorted: "⋀ x y. x < k ⟹ y ∈ {k..<n} ⟹ f (p x) ≤ f (p y)" by auto have 0: "∀ i < k. f (p i) = 0" proof (rule ccontr) assume "¬ ?thesis" then obtain i where i: "i < k" and zero: "f (p i) ≠ 0" by auto hence "f (p i) = 1" unfolding f_def by (auto split: if_splits) from sorted[OF i, unfolded this] have 1: "j ∈ {k..<n} ⟹ f (p j) ≥ 1" for j by auto have le: "j ∈ {k ..< n} ⟹ f (p j) = 1" for j using 1[of j] unfolding f_def by (auto split: if_splits) also have "?f (map p [k ..< n]) = []" using le by auto from k[unfolded this] have "length (?f (map p [0..<k])) = k" by simp from length_filter_less[of "p i" "map p [0 ..< k]" ?g, unfolded this] i zero show False by auto qed hence "?f (map p [0..<k]) = map p [0..<k]" by auto from arg_cong[OF k[unfolded this, simplified], of set] have 1: "⋀ i. i ∈ {k ..< n} ⟹ f (p i) ≠ 0" by auto have 1: "i < n ⟹ ¬ i < k ⟹ f (p i) ≠ 0" for i using 1[of i] by auto have 0: "i < n ⟹ (f (p i) = 0) = (i < k)" for i using 1[of i] 0[rule_format, of i] by blast have main: "(f i = 0) = (?p i < k)" using 0[of "?p i"] i p by (auto simp: permutes_inverses) have "i ∈ I ⟷ f i ≠ 0" unfolding f_def by simp also have "(f i = 0) ⟷ ?p i < k" using main by auto finally have "i ∈ I ⟷ ?p i ≥ k" by auto } note main = this from main[OF j] jI have k0: "k ≠ 0" by auto from iI main[OF i] have "?p i ≥ k" by auto with ilt[OF i] have kn: "k < n" by auto { fix i j assume i: "i < n" and ik: "k ≤ i" and jk: "j < k" with kn have j: "j < n" by auto have jI: "p j ∉ I" by (subst main, insert jk j p, auto simp: permutes_inverses) have iI: "p i ∈ I" by (subst main, insert i ik p, auto simp: permutes_inverses) from i j have "B $$ (i,j) = A $$ (p i, p j)" unfolding B_def by auto also have "… = 0" proof (rule ccontr) assume "A $$ (p i, p j) ≠ 0" hence "(p i, p j) ∈ ?G" unfolding graph_of_mat_def Let_def using i j p by auto with iI j have "p j ∈ I" unfolding I_def by auto with jI show False by simp qed finally have "B $$ (i,j) = 0" . } note zero = this have dimB[simp]: "dim_row B = n" "dim_col B = n" unfolding B_def by auto have dim: "dim_row B = k + (n - k)" "dim_col B = k + (n - k)" using kn by auto obtain B1 B2 B3 B4 where spl: "split_block B k k = (B1,B2,B3,B4)" (is "?tmp = _") by (cases ?tmp, auto) from split_block[OF this dim] have Bs: "B1 ∈ carrier_mat k k" "B2 ∈ carrier_mat k (n - k)" "B3 ∈ carrier_mat (n - k) k" "B4 ∈ carrier_mat (n - k) (n - k)" and B: "B = four_block_mat B1 B2 B3 B4" by auto have B3: "B3 = 0⇩_{m}(n - k) k" unfolding arg_cong[OF spl[symmetric], of "λ (_,_,B,_). B", unfolded split] unfolding split_block_def Let_def split by (rule eq_matI, auto simp: kn zero) from elements_mat_permutes[OF p A B_def] have elem: "elements_mat A = elements_mat B" . show ?thesis by (intro exI conjI, rule sim, rule elem, rule B[unfolded B3], insert Bs k0 kn, auto) qed lemma non_irreducible_nonneg_mat_split: fixes A :: "'a :: linordered_idom mat" assumes A: "A ∈ carrier_mat n n" and nonneg: "nonneg_mat A" and not: "¬ irreducible_mat A" and n: "n > 1" shows "∃ n1 n2 A1 A2. char_poly A = char_poly A1 * char_poly A2 ∧ nonneg_mat A1 ∧ nonneg_mat A2 ∧ A1 ∈ carrier_mat n1 n1 ∧ A2 ∈ carrier_mat n2 n2 ∧ 0 < n1 ∧ n1 < n ∧ 0 < n2 ∧ n2 < n ∧ n1 + n2 = n" proof - from non_irreducible_mat_split[OF A not n] obtain n1 n2 B B1 B2 B4 where sim: "similar_mat A B" and elem: "elements_mat A = elements_mat B" and B: "B = four_block_mat B1 B2 (0⇩_{m}n2 n1) B4" and Bs: "B1 ∈ carrier_mat n1 n1" "B2 ∈ carrier_mat n1 n2" "B4 ∈ carrier_mat n2 n2" and n: "0 < n1" "n1 < n" "0 < n2" "n2 < n" "n1 + n2 = n" by auto from char_poly_similar[OF sim] have AB: "char_poly A = char_poly B" . from nonneg have nonneg: "nonneg_mat B" unfolding nonneg_mat_def elem by auto have cB: "char_poly B = char_poly B1 * char_poly B4" by (rule char_poly_four_block_mat_lower_left_zero[OF B Bs]) from nonneg have B1_B4: "nonneg_mat B1" "nonneg_mat B4" unfolding B nonneg_mat_def using elements_mat_four_block_mat_supseteq[OF Bs(1-2) _ Bs(3), of "0⇩_{m}n2 n1"] by auto show ?thesis by (intro exI conjI, rule AB[unfolded cB], rule B1_B4, rule B1_B4, rule Bs, rule Bs, insert n, auto) qed text ‹The main generalized theorem. The characteristic polynomial of a non-negative real matrix can be represented as a product of roots of unitys (scaled by the the spectral radius sr) and a polynomial where all roots are smaller than the spectral radius.› theorem perron_frobenius_nonneg: fixes A :: "real Matrix.mat" assumes A: "A ∈ carrier_mat n n" and pos: "nonneg_mat A" and n: "n ≠ 0" shows "∃ sr ks f. sr ≥ 0 ∧ 0 ∉ set ks ∧ ks ≠ [] ∧ char_poly A = prod_list (map (λ k. monom 1 k - [:sr ^ k:]) ks) * f ∧ (∀ x. poly (map_poly complex_of_real f) x = 0 ⟶ cmod x < sr)" proof - define p where "p = (λ sr k. monom 1 k - [: (sr :: real) ^ k:])" let ?small = "λ f sr. (∀ x. poly (map_poly complex_of_real f) x = 0 ⟶ cmod x < sr)" let ?wit = "λ A sr ks f. sr ≥ 0 ∧ 0 ∉ set ks ∧ ks ≠ [] ∧ char_poly A = prod_list (map (p sr) ks) * f ∧ ?small f sr" let ?c = "complex_of_real" interpret c: field_hom ?c .. interpret p: map_poly_inj_idom_divide_hom ?c .. have map_p: "map_poly ?c (p sr k) = (monom 1 k - [:?c sr^k:])" for sr k unfolding p_def by (simp add: hom_distribs) { fix k x sr assume 0: "poly (map_poly ?c (p sr k)) x = 0" and k: "k ≠ 0" and sr: "sr ≥ 0" note 0 also note map_p finally have "x^k = (?c sr)^k" by (simp add: poly_monom) from arg_cong[OF this, of "λ c. root k (cmod c)", unfolded norm_power] k have "cmod x = cmod (?c sr)" using real_root_pos2 by auto also have "… = sr" using sr by auto finally have "cmod x = sr" . } note p_conv = this have "∃ sr ks f. ?wit A sr ks f" using A pos n proof (induct n arbitrary: A rule: less_induct) case (less n A) note pos = less(3) note A = less(2) note IH = less(1) note n = less(4) from n consider (1) "n = 1" | (irr) "irreducible_mat A" | (red) "¬ irreducible_mat A" "n > 1" by force thus "∃ sr ks f. ?wit A sr ks f" proof cases case irr from perron_frobenius_irreducible(3,6)[OF A n pos irr refl refl] obtain sr k f where *: "sr > 0" "k ≠ 0" "char_poly A = p sr k * f" "?small f sr" unfolding p_def by auto hence "?wit A sr [k] f" by auto thus ?thesis by blast next case red from non_irreducible_nonneg_mat_split[OF A pos red] obtain n1 n2 A1 A2 where char: "char_poly A = char_poly A1 * char_poly A2" and pos: "nonneg_mat A1" "nonneg_mat A2" and A: "A1 ∈ carrier_mat n1 n1" "A2 ∈ carrier_mat n2 n2" and n: "n1 < n" "n2 < n" and n0: "n1 ≠ 0" "n2 ≠ 0" by auto from IH[OF n(1) A(1) pos(1) n0(1)] obtain sr1 ks1 f1 where 1: "?wit A1 sr1 ks1 f1" by blast from IH[OF n(2) A(2) pos(2) n0(2)] obtain sr2 ks2 f2 where 2: "?wit A2 sr2 ks2 f2" by blast have "∃ A1 A2 sr1 ks1 f1 sr2 ks2 f2. ?wit A1 sr1 ks1 f1 ∧ ?wit A2 sr2 ks2 f2 ∧ sr1 ≥ sr2 ∧ char_poly A = char_poly A1 * char_poly A2" proof (cases "sr1 ≥ sr2") case True show ?thesis unfolding char by (intro exI, rule conjI[OF 1 conjI[OF 2]], insert True, auto) next case False show ?thesis unfolding char by (intro exI, rule conjI[OF 2 conjI[OF 1]], insert False, auto) qed then obtain A1 A2 sr1 ks1 f1 sr2 ks2 f2 where 1: "?wit A1 sr1 ks1 f1" and 2: "?wit A2 sr2 ks2 f2" and sr: "sr1 ≥ sr2" and char: "char_poly A = char_poly A1 * char_poly A2" by blast show ?thesis proof (cases "sr1 = sr2") case True have "?wit A sr2 (ks1 @ ks2) (f1 * f2)" unfolding char by (insert 1 2 True, auto simp: True p.hom_mult) thus ?thesis by blast next case False with sr have sr1: "sr1 > sr2" by auto have lt: "poly (map_poly ?c (p sr2 k)) x = 0 ⟹ k ∈ set ks2 ⟹ cmod x < sr1" for k x using sr1 p_conv[of sr2 k x] 2 by auto have "?wit A sr1 ks1 (f1 * f2 * prod_list (map (p sr2) ks2))" unfolding char by (insert 1 2 sr1 lt, auto simp: p.hom_mult p.hom_prod_list poly_prod_list prod_list_zero_iff) thus ?thesis by blast qed next case 1 define a where "a = A $$ (0,0)" have A: "A = Matrix.mat 1 1 (λ x. a)" by (rule eq_matI, unfold a_def, insert A 1(1), auto) have char: "char_poly A = [: - a, 1 :]" unfolding A by (auto simp: Determinant.det_def char_poly_def char_poly_matrix_def) from pos A have a: "a ≥ 0" unfolding nonneg_mat_def elements_mat by auto have "?wit A a [1] 1" unfolding char using a by (auto simp: p_def monom_Suc) thus ?thesis by blast qed qed then obtain sr ks f where wit: "?wit A sr ks f" by blast thus ?thesis using wit unfolding p_def by auto qed text ‹And back to HMA world via transfer.› theorem perron_frobenius_non_neg: fixes A :: "real ^ 'n ^ 'n" assumes pos: "non_neg_mat A" shows "∃ sr ks f. sr ≥ 0 ∧ 0 ∉ set ks ∧ ks ≠ [] ∧ charpoly A = prod_list (map (λ k. monom 1 k - [:sr ^ k:]) ks) * f ∧ (∀ x. poly (map_poly complex_of_real f) x = 0 ⟶ cmod x < sr)" using pos proof (transfer, goal_cases) case (1 A) from perron_frobenius_nonneg[OF 1] show ?case by auto qed text ‹We now specialize the theorem for complexity analysis where we are mainly interested in the case where the spectral radius is as most 1. Note that this can be checked by tested that there are no real roots of the characteristic polynomial which exceed 1. Moreover, here the existential quantifier over the factorization is replaced by @{const decompose_prod_root_unity}, an algorithm which computes this factorization in an efficient way.› lemma perron_frobenius_for_complexity: fixes A :: "real ^ 'n ^ 'n" and f :: "real poly" defines "cA ≡ map_matrix complex_of_real A" defines "cf ≡ map_poly complex_of_real f" assumes pos: "non_neg_mat A" and sr: "⋀ x. poly (charpoly A) x = 0 ⟹ x ≤ 1" and decomp: "decompose_prod_root_unity (charpoly A) = (ks, f)" shows "0 ∉ set ks" "charpoly A = prod_root_unity ks * f" "charpoly cA = prod_root_unity ks * cf" "⋀ x. poly (charpoly cA) x = 0 ⟹ cmod x ≤ 1" "⋀ x. poly cf x = 0 ⟹ cmod x < 1" "⋀ x. cmod x = 1 ⟹ order x (charpoly cA) = length [k←ks . x ^ k = 1]" "⋀ x. cmod x = 1 ⟹ poly (charpoly cA) x = 0 ⟹ ∃ k ∈ set ks. x^k = 1" unfolding cf_def cA_def proof (atomize(full), goal_cases) case 1 let ?c = "complex_of_real" let ?cp = "map_poly ?c" let ?A = "map_matrix ?c A" let ?wit = "λ ks f. 0 ∉ set ks ∧ charpoly A = prod_root_unity ks * f ∧ charpoly ?A = prod_root_unity ks * map_poly of_real f ∧ (∀ x. poly (charpoly ?A) x = 0 ⟶ cmod x ≤ 1) ∧ (∀ x. poly (?cp f) x = 0 ⟶ cmod x < 1)" interpret field_hom ?c .. interpret p: map_poly_inj_idom_divide_hom ?c .. { from perron_frobenius_non_neg[OF pos] obtain sr ks f where *: "sr ≥ 0" "0 ∉ set ks" "ks ≠ []" and cp: "charpoly A = prod_list (map (λ k. monom 1 k - [:sr ^ k:]) ks) * f" and small: "⋀ x. poly (?cp f) x = 0 ⟹ cmod x < sr" by blast from arg_cong[OF cp, of "map_poly ?c"] have cpc: "charpoly ?A = prod_list (map (λ k. monom 1 k - [:?c sr ^ k:]) ks) * map_poly ?c f" by (simp add: charpoly_of_real hom_distribs p.prod_list_map_hom[symmetric] o_def) have sr_le_1: "sr ≤ 1" by (rule sr, unfold cp, insert *, cases ks, auto simp: poly_monom) { fix x note [simp] = prod_list_zero_iff o_def poly_monom assume "poly (charpoly ?A) x = 0" from this[unfolded cpc poly_mult poly_prod_list] small[of x] consider (lt) "cmod x < sr" | (mem) k where "k ∈ set ks" "x ^ k = (?c sr) ^ k" by force hence "cmod x ≤ sr" proof (cases) case (mem k) with * have k: "k ≠ 0" by metis with arg_cong[OF mem(2), of "λ x. root k (cmod x)", unfolded norm_power] real_root_pos2[of k] *(1) have "cmod x = sr" by auto thus ?thesis by auto qed simp } note root = this have "∃ ks f. ?wit ks f" proof (cases "sr = 1") case False with sr_le_1 have *: "cmod x ≤ sr ⟹ cmod x < 1" "cmod x ≤ sr ⟹ cmod x ≤ 1" for x by auto show ?thesis by (rule exI[of _ Nil], rule exI[of _ "charpoly A"], insert * root, auto simp: prod_root_unity_def charpoly_of_real) next case sr: True from * cp cpc small root show ?thesis unfolding sr root_unity_def prod_root_unity_def by (auto simp: pCons_one) qed } then obtain Ks F where wit: "?wit Ks F" by auto have cA0: "charpoly ?A ≠ 0" using degree_monic_charpoly[of ?A] by auto from wit have id: "charpoly ?A = prod_root_unity Ks * ?cp F" by auto from of_real_hom.hom_decompose_prod_root_unity[of "charpoly A", unfolded decomp] have decompc: "decompose_prod_root_unity (charpoly ?A) = (ks, ?cp f)" by (auto simp: charpoly_of_real) from wit have small: "cmod x = 1 ⟹ poly (?cp F) x ≠ 0" for x by auto from decompose_prod_root_unity[OF id decompc this cA0] have id: "charpoly ?A = prod_root_unity ks * ?cp F" "F = f" "set Ks = set ks" by auto have "?cp (charpoly A) = ?cp (prod_root_unity ks * f)" unfolding id unfolding charpoly_of_real[symmetric] id p.hom_mult of_real_hom.hom_prod_root_unity .. hence idr: "charpoly A = prod_root_unity ks * f" by auto have wit: "?wit ks f" and idc: "charpoly ?A = prod_root_unity ks * ?cp f" using wit unfolding id idr by auto { fix x assume "cmod x = 1" from small[OF this, unfolded id] have "poly (?cp f) x ≠ 0" by auto from order_0I[OF this] this have ord: "order x (?cp f) = 0" and cf0: "?cp f ≠ 0" by auto have "order x (charpoly ?A) = order x (prod_root_unity ks)" unfolding idc by (subst order_mult, insert cf0 wit ord, auto) also have "… = length [k←ks . x ^ k = 1]" by (subst order_prod_root_unity, insert wit, auto) finally have ord: "order x (charpoly ?A) = length [k←ks . x ^ k = 1]" . { assume "poly (charpoly ?A) x = 0" with cA0 have "order x (charpoly ?A) ≠ 0" unfolding order_root by auto from this[unfolded ord] have "∃ k ∈ set ks. x ^ k = 1" by (cases "[k←ks . x ^ k = 1]", force+) } note this ord } with wit show ?case by blast qed text ‹and convert to JNF-world› lemmas perron_frobenius_for_complexity_jnf = perron_frobenius_for_complexity[unfolded atomize_imp atomize_all, untransferred, cancel_card_constraint, rule_format] end