Theory Spectral_Radius_Theory
section ‹Combining Spectral Radius Theory with Perron Frobenius theorem›
theory Spectral_Radius_Theory
imports
Polynomial_Factorization.Square_Free_Factorization
Jordan_Normal_Form.Spectral_Radius
Jordan_Normal_Form.Char_Poly
Perron_Frobenius
"HOL-Computational_Algebra.Field_as_Ring"
begin
abbreviation spectral_radius where "spectral_radius ≡ Spectral_Radius.spectral_radius"
hide_const (open) Module.smult
text ‹Via JNFs it has been proven that the growth of $A^k$ is polynomially bounded,
if all complex eigenvalues have a norm at most 1, i.e., the spectral radius must be
at most 1. Moreover, the degree of the polynomial growth can be
bounded by the order of those roots which have norm 1, cf. @{thm spectral_radius_poly_bound}.›
text ‹Perron Frobenius theorem tells us that for a real valued non negative matrix,
the largest eigenvalue is a real non-negative one. Hence, we only have to check, that all real
eigenvalues are at most one.›
text ‹We combine both theorems in the following. To be more precise,
the set-based complexity results from JNFs with the type-based
Perron Frobenius theorem in HMA are connected to obtain
a set based complexity criterion for real-valued
non-negative matrices, where one only investigated the real valued eigenvalues for
checking the eigenvalue-at-most-1 condition.
Here, in the precondition of the roots of the polynomial, the type-system ensures
that we only have to look at real-valued eigenvalues, and can ignore the
complex-valued ones.
The linkage between set-and type-based is performed via HMA-connect.›
lemma perron_frobenius_spectral_radius_complex: fixes A :: "complex mat"
assumes A: "A ∈ carrier_mat n n"
and real_nonneg: "real_nonneg_mat A"
and ev_le_1: "⋀ x. poly (char_poly (map_mat Re A)) x = 0 ⟹ x ≤ 1"
and ev_order: "⋀ x. norm x = 1 ⟹ order x (char_poly A) ≤ d"
shows "∃c1 c2. ∀k. norm_bound (A ^⇩m k) (c1 + c2 * real k ^ (d - 1))"
proof (cases "n = 0")
case False
hence n: "n > 0" "n ≠ 0" by auto
define sr where "sr = spectral_radius A"
note sr = spectral_radius_mem_max[OF A n(1), folded sr_def]
show ?thesis
proof (rule spectral_radius_poly_bound[OF A], unfold sr_def[symmetric])
let ?cr = "complex_of_real"
text ‹here is the transition from type-based perron-frobenius to set-based›
from perron_frobenius[untransferred, cancel_card_constraint, OF A real_nonneg n(2)]
obtain v where v: "v ∈ carrier_vec n" and ev: "eigenvector A v (?cr sr)" and
rnn: "real_nonneg_vec v" unfolding sr_def by auto
define B where "B = map_mat Re A"
let ?A = "map_mat ?cr B"
have AB: "A = ?A" unfolding B_def
by (rule eq_matI, insert real_nonneg[unfolded real_nonneg_mat_def elements_mat_def], auto)
define w where "w = map_vec Re v"
let ?v = "map_vec ?cr w"
have vw: "v = ?v" unfolding w_def
by (rule eq_vecI, insert rnn[unfolded real_nonneg_vec_def vec_elements_def], auto)
have B: "B ∈ carrier_mat n n" unfolding B_def using A by auto
from AB vw ev have ev: "eigenvector ?A ?v (?cr sr)" by simp
have "eigenvector B w sr"
by (rule of_real_hom.eigenvector_hom_rev[OF B ev])
hence "eigenvalue B sr" unfolding eigenvalue_def by blast
from ev_le_1[folded B_def, OF this[unfolded eigenvalue_root_char_poly[OF B]]]
show "sr ≤ 1" .
next
fix ev
assume "cmod ev = 1"
thus "order ev (char_poly A) ≤ d" by (rule ev_order)
qed
next
case True
with A show ?thesis
by (intro exI[of _ 0], auto simp: norm_bound_def)
qed
text ‹The following lemma is the same as @{thm perron_frobenius_spectral_radius_complex},
except that now the type @{typ real} is used instead of @{typ complex}.›
lemma perron_frobenius_spectral_radius: fixes A :: "real mat"
assumes A: "A ∈ carrier_mat n n"
and nonneg: "nonneg_mat A"
and ev_le_1: "∀ x. poly (char_poly A) x = 0 ⟶ x ≤ 1"
and ev_order: "∀ x :: complex. norm x = 1 ⟶ order x (map_poly of_real (char_poly A)) ≤ d"
shows "∃c1 c2. ∀k a. a ∈ elements_mat (A ^⇩m k) ⟶ abs a ≤ (c1 + c2 * real k ^ (d - 1))"
proof -
let ?cr = "complex_of_real"
let ?B = "map_mat ?cr A"
have B: "?B ∈ carrier_mat n n" using A by auto
have rnn: "real_nonneg_mat ?B" using nonneg unfolding real_nonneg_mat_def nonneg_mat_def
by (auto simp: elements_mat_def)
have id: "map_mat Re ?B = A"
by (rule eq_matI, auto)
have "∃c1 c2. ∀k. norm_bound (?B ^⇩m k) (c1 + c2 * real k ^ (d - 1))"
by (rule perron_frobenius_spectral_radius_complex[OF B rnn], unfold id,
insert ev_le_1 ev_order, auto simp: of_real_hom.char_poly_hom[OF A])
then obtain c1 c2 where nb: "⋀ k. norm_bound (?B ^⇩m k) (c1 + c2 * real k ^ (d - 1))" by auto
show ?thesis
proof (rule exI[of _ c1], rule exI[of _ c2], intro allI impI)
fix k a
assume "a ∈ elements_mat (A ^⇩m k)"
with pow_carrier_mat[OF A] obtain i j where a: "a = (A ^⇩m k) $$ (i,j)" and ij: "i < n" "j < n"
unfolding elements_mat by force
from ij nb[of k] A have "norm ((?B ^⇩m k) $$ (i,j)) ≤ c1 + c2 * real k ^ (d - 1)"
unfolding norm_bound_def by auto
also have "(?B ^⇩m k) $$ (i,j) = ?cr a"
unfolding of_real_hom.mat_hom_pow[OF A, symmetric] a using ij A by auto
also have "norm (?cr a) = abs a" by auto
finally show "abs a ≤ (c1 + c2 * real k ^ (d - 1))" .
qed
qed
text ‹We can also convert the set-based lemma @{thm perron_frobenius_spectral_radius}
to a type-based version.›
lemma perron_frobenius_spectral_type_based:
assumes "non_neg_mat (A :: real ^ 'n ^ 'n)"
and "∀ x. poly (charpoly A) x = 0 ⟶ x ≤ 1"
and "∀ x :: complex. norm x = 1 ⟶ order x (map_poly of_real (charpoly A)) ≤ d"
shows "∃c1 c2. ∀k a. a ∈ elements_mat_h (matpow A k) ⟶ abs a ≤ (c1 + c2 * real k ^ (d - 1))"
using assms perron_frobenius_spectral_radius
by (transfer, blast)
text ‹And of course, we can also transfer the type-based lemma back to a set-based setting,
only that -- without further case-analysis --
we get the additional assumption @{term "(n :: nat) ≠ 0"}.›
lemma assumes "A ∈ carrier_mat n n"
and "nonneg_mat A"
and "∀ x. poly (char_poly A) x = 0 ⟶ x ≤ 1"
and "∀ x :: complex. norm x = 1 ⟶ order x (map_poly of_real (char_poly A)) ≤ d"
and "n ≠ 0"
shows "∃c1 c2. ∀k a. a ∈ elements_mat (A ^⇩m k) ⟶ abs a ≤ (c1 + c2 * real k ^ (d - 1))"
using perron_frobenius_spectral_type_based[untransferred, cancel_card_constraint, OF assms] .
text ‹Note that the precondition eigenvalue-at-most-1 can easily be formulated as a cardinality
constraints which can be decided by Sturm's theorem.
And in order to obtain a bound on the order, one can
perform a square-free-factorization (via Yun's factorization algorithm)
of the characteristic polynomial into
$f_1^1 \cdot \ldots f_d^d$ where each $f_i$ has precisely the roots of order $i$.›
context
fixes A :: "real mat" and c :: real and fis and n :: nat
assumes A: "A ∈ carrier_mat n n"
and nonneg: "nonneg_mat A"
and yun: "yun_factorization gcd (char_poly A) = (c,fis)"
and ev_le_1: "card {x. poly (char_poly A) x = 0 ∧ x > 1} = 0"
begin
lemma perron_frobenius_spectral_radius_yun:
assumes bnd: "⋀ f⇩i i. (f⇩i,i) ∈ set fis
⟹ (∃ x :: complex. poly (map_poly of_real f⇩i) x = 0 ∧ norm x = 1)
⟹ i ≤ d"
shows "∃c1 c2. ∀k a. a ∈ elements_mat (A ^⇩m k) ⟶ abs a ≤ (c1 + c2 * real k ^ (d - 1))"
proof (rule perron_frobenius_spectral_radius[OF A nonneg]; intro allI impI)
let ?cr = complex_of_real
let ?cp = "map_poly ?cr (char_poly A)"
fix x :: complex
assume x: "norm x = 1"
have A0: "char_poly A ≠ 0" using degree_monic_char_poly[OF A] by auto
interpret field_hom_0' ?cr by (standard, auto)
from A0 have cp0: "?cp ≠ 0" by auto
obtain ox where ox: "order x ?cp = ox" by blast
note sff = square_free_factorization_order_root[OF yun_factorization(1)[OF
yun_factorization_hom[of "char_poly A", unfolded yun map_prod_def split]] cp0, of x ox, unfolded ox]
show "order x ?cp ≤ d" unfolding ox
proof (cases ox)
case (Suc oo)
with sff obtain fi where mem: "(fi,Suc oo) ∈ set fis" and rt: "poly (map_poly ?cr fi) x = 0" by auto
from bnd[OF mem exI[of _ x], OF conjI[OF rt x]]
show "ox ≤ d" unfolding Suc .
qed auto
next
let ?L = "{x. poly (char_poly A) x = 0 ∧ x > 1}"
fix x :: real
assume rt: "poly (char_poly A) x = 0"
have "finite ?L"
by (rule finite_subset[OF _ poly_roots_finite[of "char_poly A"]],
insert degree_monic_char_poly[OF A], auto)
with ev_le_1 have "?L = {}" by simp
with rt show "x ≤ 1" by auto
qed
text ‹Note that the only remaining problem in applying
@{thm perron_frobenius_spectral_radius_yun} is to check the
condition @{term "∃ x :: complex. poly (map_poly of_real f⇩i) x = 0 ∧ norm x = 1"}.
Here, there are at least three possibilities.
First, one can just ignore this precondition and weaken the statement.
Second, one can apply Sturm's theorem to determine whether all roots are real.
This can be done by comparing the number of distinct real roots with the degree of @{term f⇩i},
since @{term f⇩i} is square-free. If all roots are real, then one can decide the criterion
by checking the only two possible real roots with norm equal to 1, namely 1 and -1.
If on the other hand there are complex roots, then we loose precision at this point.
Third, one uses a factorization algorithm (e.g., via complex algebraic numbers) to
precisely determine the complex roots and decide the condition.
The second approach is illustrated in the following theorem. Note that all preconditions --
including the ones from the context --
can easily be checked with the help of Sturm's method.
This method is used as a fast approximative technique in CeTA \<^cite>‹"CeTA"›. Only if the desired degree
cannot be ensured by this method, the more costly complex algebraic number based
factorization is applied.›
lemma perron_frobenius_spectral_radius_yun_real_roots:
assumes bnd: "⋀ f⇩i i. (f⇩i,i) ∈ set fis
⟹ card { x. poly f⇩i x = 0} ≠ degree f⇩i ∨ poly f⇩i 1 = 0 ∨ poly f⇩i (-1) = 0
⟹ i ≤ d"
shows "∃c1 c2. ∀k a. a ∈ elements_mat (A ^⇩m k) ⟶ abs a ≤ (c1 + c2 * real k ^ (d - 1))"
proof (rule perron_frobenius_spectral_radius_yun)
fix fi i
let ?cr = complex_of_real
let ?cp = "map_poly ?cr"
assume fi: "(fi, i) ∈ set fis"
and "∃ x. poly (map_poly ?cr fi) x = 0 ∧ norm x = 1"
then obtain x where rt: "poly (?cp fi) x = 0" and x: "norm x = 1" by auto
show "i ≤ d"
proof (rule bnd[OF fi])
consider (c) "x ∉ ℝ" | (1) "x = 1" | (m1) "x = -1" | (r) "x ∈ ℝ" "x ∉ {1, -1}"
by (cases "x ∈ ℝ"; auto)
thus "card {x. poly fi x = 0} ≠ degree fi ∨ poly fi 1 = 0 ∨ poly fi (- 1) = 0"
proof (cases)
case 1
from rt have "poly fi 1 = 0"
unfolding 1 by simp
thus ?thesis by simp
next
case m1
have id: "-1 = ?cr (-1)" by simp
from rt have "poly fi (-1) = 0"
unfolding m1 id of_real_hom.hom_zero[where 'a=complex,symmetric] of_real_hom.poly_map_poly by simp
thus ?thesis by simp
next
case r
then obtain y where xy: "x = of_real y" unfolding Reals_def by auto
from r(2)[unfolded xy] have y: "y ∉ {1,-1}" by auto
from x[unfolded xy] have "abs y = 1" by auto
with y have False by auto
thus ?thesis ..
next
case c
from yun_factorization(2)[OF yun] fi have "monic fi" by auto
hence fi: "?cp fi ≠ 0" by auto
hence fin: "finite {x. poly (?cp fi) x = 0}" by (rule poly_roots_finite)
have "?cr ` {x. poly (?cp fi) (?cr x) = 0} ⊂ {x. poly (?cp fi) x = 0}" (is "?l ⊂ ?r")
proof (rule, force)
have "x ∈ ?r" using rt by auto
moreover have "x ∉ ?l" using c unfolding Reals_def by auto
ultimately show "?l ≠ ?r" by blast
qed
from psubset_card_mono[OF fin this] have "card ?l < card ?r" .
also have "… ≤ degree (?cp fi)" by (rule poly_roots_degree[OF fi])
also have "… = degree fi" by simp
also have "?l = ?cr ` {x. poly fi x = 0}" by auto
also have "card … = card {x. poly fi x = 0}"
by (rule card_image, auto simp: inj_on_def)
finally have "card {x. poly fi x = 0} ≠ degree fi" by simp
thus ?thesis by auto
qed
qed
qed
end
end