Theory Perron_Frobenius_General
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))"
unbundle no m_inv_syntax
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