Theory Perron_Frobenius

(* Author: R. Thiemann *)

subsection ‹Perron-Frobenius theorem via Brouwer's fixpoint theorem.›

theory Perron_Frobenius
imports
  "HOL-Analysis.Brouwer_Fixpoint"
  Perron_Frobenius_Aux
begin

text ‹We follow the textbook proof of Serre cite‹Theorem 5.2.1› in "SerreMatrices".›

context
  fixes A :: "complex ^ 'n ^ 'n :: finite"
  assumes rnnA: "real_non_neg_mat A"
begin

private abbreviation(input) sr where "sr  spectral_radius A"

private definition max_v_ev :: "(complex^'n) × complex" where
  "max_v_ev = (SOME v_ev. eigen_vector A (fst v_ev) (snd v_ev)
   norm (snd v_ev) = sr)"

private definition "max_v = (1 / norm1 (fst max_v_ev)) *R fst max_v_ev"
private definition "max_ev = snd max_v_ev"

private lemma max_v_ev:
  "eigen_vector A max_v max_ev"
  "norm max_ev = sr"
  "norm1 max_v = 1"
proof -
  obtain v ev where id: "max_v_ev = (v,ev)" by force
  from spectral_radius_ev[of A] someI_ex[of "λ v_ev. eigen_vector A (fst v_ev) (snd v_ev)
   norm (snd v_ev) = sr", folded max_v_ev_def, unfolded id]
  have v: "eigen_vector A v ev" and ev: "norm ev = sr" by auto
  from normalize_eigen_vector[OF v] ev
  show "eigen_vector A max_v max_ev" "norm max_ev = sr" "norm1 max_v = 1"
    unfolding max_v_def max_ev_def id by auto
qed

text ‹In the definition of S, we use the linear norm instead of the
  default euclidean norm which is defined via the type-class.
  The reason is that S is not convex if one uses the euclidean norm.›

private definition B :: "real ^ 'n ^ 'n" where "B  χ i j. Re (A $ i $ j)"
private definition S where "S = {v :: real ^ 'n . norm1 v = 1  ( i. v $ i  0) 
  ( i. (B *v v) $ i  sr * (v $ i))}"
private definition f :: "real ^ 'n  real ^ 'n" where
  "f v = (1 / norm1 (B *v v)) *R (B *v v)"

private lemma closedS: "closed S"
  unfolding S_def matrix_vector_mult_def[abs_def]
proof (intro closed_Collect_conj closed_Collect_all closed_Collect_le closed_Collect_eq)
  show "continuous_on UNIV norm1"
    by (simp add: continuous_at_imp_continuous_on)
qed (auto intro!: continuous_intros continuous_on_component)

private lemma boundedS: "bounded S"
proof -
  {
    fix v :: "real ^ 'n"
    from norm1_ge_norm[of v] have "norm1 v = 1  norm v  1" by auto
  }
  thus ?thesis
  unfolding S_def bounded_iff
  by (auto intro!: exI[of _ 1])
qed

private lemma compactS: "compact S"
  using boundedS closedS
  by (simp add: compact_eq_bounded_closed)

private lemmas rnn = real_non_neg_matD[OF rnnA]

lemma B_norm: "B $ i $ j = norm (A $ i $ j)"
  using rnn[of i j]
  by (cases "A $ i $ j", auto simp: B_def)

lemma mult_B_mono: assumes " i. v $ i  w $ i"
  shows "(B *v v) $ i  (B *v w) $ i" unfolding matrix_vector_mult_def vec_lambda_beta
  by (rule sum_mono, rule mult_left_mono[OF assms], unfold B_norm, auto)


private lemma non_emptyS: "S  {}"
proof -
  let ?v = "(χ i. norm (max_v $ i)) :: real ^ 'n"
  have "norm1 max_v = 1" by (rule max_v_ev(3))
  hence nv: "norm1 ?v = 1" unfolding norm1_def by auto
  {
    fix i
    have "sr * (?v $ i) = sr * norm (max_v $ i)" by auto
    also have " = (norm max_ev) * norm (max_v $ i)" using max_v_ev by auto
    also have " = norm ((max_ev *s max_v) $ i)" by (auto simp: norm_mult)
    also have "max_ev *s max_v = A *v max_v" using max_v_ev(1)[unfolded eigen_vector_def] by auto
    also have "norm ((A *v max_v) $ i)  (B *v ?v) $ i"
      unfolding matrix_vector_mult_def vec_lambda_beta
      by (rule sum_norm_le, auto simp: norm_mult B_norm)
    finally have "sr * (?v $ i)  (B *v ?v) $ i" .
  } note le = this
  have "?v  S" unfolding S_def using nv le by auto
  thus ?thesis by blast
qed

private lemma convexS: "convex S"
proof (rule convexI)
  fix v w a b
  assume *: "v  S" "w  S" "0  a" "0  b" "a + b = (1 :: real)"
  let ?lin = "a *R v + b *R w"
  from * have 1: "norm1 v = 1" "norm1 w = 1" unfolding S_def by auto
  have "norm1 ?lin = a * norm1 v + b * norm1 w"
    unfolding norm1_def sum_distrib_left sum.distrib[symmetric]
  proof (rule sum.cong)
    fix i :: 'n
    from * have "v $ i  0" "w $ i  0" unfolding S_def by auto
    thus "norm (?lin $ i) = a * norm (v $ i) + b * norm (w $ i)"
      using *(3-4) by auto
  qed simp
  also have " = 1" using *(5) 1 by auto
  finally have norm1: "norm1 ?lin = 1" .
  {
    fix i
    from * have "0  v $ i" "sr * v $ i  (B *v v) $ i" unfolding S_def by auto
    with a  0 have a: "a * (sr * v $ i)  a * (B *v v) $ i" by (intro mult_left_mono)
    from * have "0  w $ i" "sr * w $ i  (B *v w) $ i" unfolding S_def by auto
    with b  0 have b: "b * (sr * w $ i)  b * (B *v w) $ i" by (intro mult_left_mono)
    from a b have "a * (sr * v $ i) + b * (sr * w $ i)  a * (B *v v) $ i + b * (B *v w) $ i" by auto
  } note le = this
  have switch[simp]: " x y. x * a * y = a * x * y"  " x y. x * b * y = b * x * y" by auto
  have [simp]: "x  {v,w}  a * (r * x $h i) = r * (a * x $h i)" for a r i x by auto
  show "a *R v + b *R w  S" using * norm1 le unfolding S_def
    by (auto simp: matrix_vect_scaleR matrix_vector_right_distrib ring_distribs)
qed

private abbreviation (input) r :: "real  complex" where
  "r  of_real"

private abbreviation rv :: "real ^'n  complex ^'n" where
  "rv v  χ i. r (v $ i)"

private lemma rv_0: "(rv v = 0) = (v = 0)"
  by (simp add: of_real_hom.map_vector_0 map_vector_def vec_eq_iff)

private lemma rv_mult: "A *v rv v = rv (B *v v)"
proof -
  have "map_matrix r B = A"
    using rnnA unfolding map_matrix_def B_def real_non_neg_mat_def map_vector_def elements_mat_h_def
    by vector
  thus ?thesis
    using of_real_hom.matrix_vector_mult_hom[of B, where 'a = complex]
    unfolding map_vector_def by auto
qed

context
  assumes zero_no_ev: " v. v  S  A *v rv v  0"
begin
private lemma normB_S: assumes v: "v  S"
  shows "norm1 (B *v v)  0"
proof -
  from zero_no_ev[OF v, unfolded rv_mult rv_0]
  show ?thesis by auto
qed

private lemma image_f: "f  S  S"
proof -
  {
    fix v
    assume v: "v  S"
    hence norm: "norm1 v = 1" and ge: " i. v $ i  0" " i. sr * v $ i  (B *v v) $ i" unfolding S_def by auto
    from normB_S[OF v] have normB: "norm1 (B *v v) > 0" using norm1_nonzero by auto
    have fv: "f v = (1 / norm1 (B *v v)) *R (B *v v)" unfolding f_def by auto
    from normB have Bv0: "B *v v  0" unfolding norm1_0_iff[symmetric] by linarith
    have norm: "norm1 (f v) = 1" unfolding fv using normB Bv0 by simp
    define c where "c = (1 / norm1 (B *v v))"
    have c: "c > 0" unfolding c_def using normB by auto
    {
      fix i
      have 1: "f v $ i  0" unfolding fv c_def[symmetric] using c ge
        by (auto simp: matrix_vector_mult_def sum_distrib_left B_norm intro!: sum_nonneg)
      have id1: " i. (B *v f v) $ i = c * ((B *v (B *v v)) $ i)"
        unfolding f_def c_def matrix_vect_scaleR by simp
      have id3: " i. sr * f v $ i = c * ((B *v (sr *R v)) $ i)"
        unfolding f_def c_def[symmetric] matrix_vect_scaleR by auto
      have 2: "sr * f v $ i  (B *v f v) $ i" unfolding id1 id3
        unfolding mult_le_cancel_iff2[OF c > 0]
        by (rule mult_B_mono, insert ge(2), auto)
      note 1 2
    }
    with norm have "f v  S" unfolding S_def by auto
  }
  thus ?thesis by blast
qed

private lemma cont_f: "continuous_on S f"
  unfolding f_def[abs_def] continuous_on using normB_S
  unfolding norm1_def
  by (auto intro!: tendsto_eq_intros)

qualified lemma perron_frobenius_positive_ev:
  " v. eigen_vector A v (r sr)  real_non_neg_vec v"
proof -
  from brouwer[OF compactS convexS non_emptyS cont_f image_f]
  obtain v where v: "v  S" and fv: "f v = v" by auto
  define ev where "ev = norm1 (B *v v)"
  from normB_S[OF v] have "ev  0" unfolding ev_def by auto
  with norm1_ge_0[of "B *v v", folded ev_def] have norm: "ev > 0" by auto
  from arg_cong[OF fv[unfolded f_def], of "λ (w :: real ^ 'n). ev *R w"] norm
  have ev: "B *v v = ev *s v" unfolding ev_def[symmetric] scalar_mult_eq_scaleR by simp
  with v[unfolded S_def] have ge: " i. sr * v $ i  ev * v $ i" by auto
  have "A *v rv v = rv (B *v v)" unfolding rv_mult ..
  also have " = ev *s rv v" unfolding ev vec_eq_iff
    by (simp add: scaleR_conv_of_real scaleR_vec_def)
  finally have ev: "A *v rv v = ev *s rv v" .
  from v have v0: "v  0" unfolding S_def by auto
  hence "rv v  0" unfolding rv_0 .
  with ev have ev: "eigen_vector A (rv v) ev" unfolding eigen_vector_def by auto
  hence "eigen_value A ev" unfolding eigen_value_def by auto
  from spectral_radius_max[OF this] have le: "norm (r ev)  sr" .
  from v0 obtain i where "v $ i  0" unfolding vec_eq_iff by auto
  from v have "v $ i  0" unfolding S_def by auto
  with v $ i  0 have "v $ i > 0" by auto
  with ge[of i] have ge: "sr  ev" by auto
  with le have sr: "r sr = ev" by auto
  from v have *: "real_non_neg_vec (rv v)" unfolding S_def real_non_neg_vec_def vec_elements_h_def by auto
  show ?thesis unfolding sr
    by (rule exI[of _ "rv v"], insert * ev norm, auto)
qed
end

qualified lemma perron_frobenius_both:
  " v. eigen_vector A v (r sr)  real_non_neg_vec v"
proof (cases " v  S. A *v rv v  0")
  case True
  show ?thesis
    by (rule Perron_Frobenius.perron_frobenius_positive_ev[OF rnnA], insert True, auto)
next
  case False
  then obtain v where v: "v  S" and A0: "A *v rv v = 0" by auto
  hence id: "A *v rv v = 0 *s rv v" and v0: "v  0" unfolding S_def by auto
  from v0 have "rv v  0" unfolding rv_0 .
  with id have ev: "eigen_vector A (rv v) 0" unfolding eigen_vector_def by auto
  hence "eigen_value A 0" unfolding eigen_value_def ..
  from spectral_radius_max[OF this] have 0: "0  sr" by auto
  from v[unfolded S_def] have ge: " i. sr * v $ i  (B *v v) $ i" by auto
  from v[unfolded S_def] have rnn: "real_non_neg_vec (rv v)"
    unfolding real_non_neg_vec_def vec_elements_h_def by auto
  from v0 obtain i where "v $ i  0" unfolding vec_eq_iff by auto
  from v have "v $ i  0" unfolding S_def by auto
  with v $ i  0 have vi: "v $ i > 0" by auto
  from rv_mult[of v, unfolded A0] have "rv (B *v v) = 0" by simp
  hence "B *v v = 0" unfolding rv_0 .
  from ge[of i, unfolded this] vi have ge: "sr  0" by (simp add: mult_le_0_iff)
  with 0  sr have "sr = 0" by auto
  show ?thesis unfolding sr = 0 using rnn ev by auto
qed
end

text ‹Perron Frobenius: The largest complex eigenvalue of a real-valued non-negative matrix
  is a real one, and it has a real-valued non-negative eigenvector.›

lemma perron_frobenius:
  assumes "real_non_neg_mat A"
  shows "v. eigen_vector A v (of_real (spectral_radius A))  real_non_neg_vec v"
  by (rule Perron_Frobenius.perron_frobenius_both[OF assms])

text ‹And a version which ignores the eigenvector.›

lemma perron_frobenius_eigen_value:
  assumes "real_non_neg_mat A"
  shows "eigen_value A (of_real (spectral_radius A))"
  using perron_frobenius[OF assms] unfolding eigen_value_def by blast

end