Theory Spectral_Radius_Theory

theory Spectral_Radius_Theory
imports Square_Free_Factorization Perron_Frobenius Field_as_Ring
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

text ‹Note that @{const yun_factorization} has an offset by 1, 
  so the pair @{term "(fi,i) ∈ set fis"} encodes @{term "fi^(Suc i)"}.›
lemma perron_frobenius_spectral_radius_yun: 
  assumes bnd: "⋀ fi i. (fi,i) ∈ set fis 
    ⟹ (∃ x :: complex. poly (map_poly of_real fi) x = 0 ∧ norm x = 1) 
    ⟹ Suc 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,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 fi) 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 fi},
    since @{term fi} 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: "⋀ fi i. (fi,i) ∈ set fis 
    ⟹ card { x. poly fi x = 0} ≠ degree fi ∨ poly fi 1 = 0 ∨ poly fi (-1) = 0 
    ⟹ Suc 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 "Suc 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

thm perron_frobenius_spectral_radius_yun_real_roots
end