Theory Expander_Graphs.Expander_Graphs_Algebra

section ‹Algebra-only Theorems›

text ‹This section verifies the linear algebraic counter-parts of the graph-theoretic theorems
about Random walks. The graph-theoretic results are then derived in Section~\ref{sec:random_walks}.›

theory Expander_Graphs_Algebra
  imports
    "HOL-Library.Monad_Syntax"
    Expander_Graphs_TTS
begin

lemma pythagoras:
  fixes v w :: "'a::real_inner"
  assumes "v  w  = 0"
  shows "norm (v+w)^2 = norm v^2 + norm w^2"
  using assms by (simp add:power2_norm_eq_inner algebra_simps inner_commute)

definition diag :: "('a :: zero)^'n  'a^'n^'n"
  where "diag v = (χ i j. if i = j then (v $ i) else 0)"

definition ind_vec :: "'n set  real^'n"
  where "ind_vec S = (χ i. of_bool( i  S))"

lemma diag_mult_eq: "diag x ** diag y = diag (x * y)"
  unfolding diag_def
  by (vector matrix_matrix_mult_def)
   (auto simp add:if_distrib if_distribR sum.If_cases)

lemma diag_vec_mult_eq: "diag x *v y = x * y"
  unfolding diag_def matrix_vector_mult_def
  by (simp add:if_distrib if_distribR sum.If_cases times_vec_def)

definition matrix_norm_bound :: "real^'n^'m  real  bool"
  where "matrix_norm_bound A l = (x. norm (A *v x)  l * norm x)"

lemma  matrix_norm_boundI:
  assumes "x. norm (A *v x)  l * norm x"
  shows "matrix_norm_bound A l"
  using assms unfolding matrix_norm_bound_def by simp

lemma matrix_norm_boundD:
  assumes "matrix_norm_bound A l"
  shows "norm (A *v x)  l * norm x"
  using assms unfolding matrix_norm_bound_def by simp

lemma matrix_norm_bound_nonneg:
  fixes A :: "real^'n^'m"
  assumes "matrix_norm_bound A l"
  shows "l  0"
proof -
  have "0  norm (A *v 1)" by simp
  also have "...  l * norm (1::real^'n)"
    using assms(1) unfolding matrix_norm_bound_def by simp
  finally have "0  l  * norm (1::real^'n)"
    by simp
  moreover have "norm (1::real^'n) > 0"
    by simp
  ultimately show ?thesis
    by (simp add: zero_le_mult_iff)
qed

lemma  matrix_norm_bound_0:
  assumes "matrix_norm_bound A 0"
  shows "A = (0::real^'n^'m)"
proof (intro iffD2[OF matrix_eq] allI)
  fix x :: "real^'n"
  have "norm (A *v x) = 0"
    using assms unfolding matrix_norm_bound_def by simp
  thus "A *v x = 0 *v x"
    by simp
qed

lemma matrix_norm_bound_diag:
  fixes x :: "real^'n"
  assumes "i. ¦x $ i¦  l"
  shows "matrix_norm_bound (diag x) l"
proof (rule matrix_norm_boundI)
  fix y :: "real^'n"

  have l_ge_0: "l  0" using assms by fastforce

  have a: "¦x $ i * v¦  ¦l * v¦" for v i
    using l_ge_0 assms by (simp add:abs_mult mult_right_mono)

  have "norm (diag x *v y) = sqrt (i  UNIV. (x $ i * y $ i)^2)"
    unfolding matrix_vector_mult_def diag_def norm_vec_def L2_set_def
    by (auto simp add:if_distrib if_distribR sum.If_cases)
  also have "...  sqrt (i  UNIV. (l * y $ i)^2)"
    by (intro real_sqrt_le_mono sum_mono iffD1[OF abs_le_square_iff] a)
  also have "... = l * norm y"
    using l_ge_0 by (simp add:norm_vec_def L2_set_def algebra_simps
        sum_distrib_left[symmetric] real_sqrt_mult)
  finally show "norm (diag x *v y)  l * norm y" by simp
qed

lemma vector_scaleR_matrix_ac_2: "b *R (A::real^'n^'m) *v x = b *R (A *v x)"
  unfolding vector_transpose_matrix[symmetric]  transpose_scalar
  by (intro vector_scaleR_matrix_ac)

lemma  matrix_norm_bound_scale:
  assumes "matrix_norm_bound A l"
  shows "matrix_norm_bound (b *R A) (¦b¦ * l)"
proof (intro matrix_norm_boundI)
  fix x
  have "norm (b *R A *v x) = norm (b *R (A *v x))"
    by (metis transpose_scalar vector_scaleR_matrix_ac vector_transpose_matrix)
  also have "... = ¦b¦ * norm (A *v x)"
    by simp
  also have "...  ¦b¦ * (l * norm x)"
    using assms matrix_norm_bound_def by (intro mult_left_mono) auto
  also have "...  (¦b¦ * l) * norm x" by simp
  finally show "norm (b *R A *v x)  (¦b¦ * l) * norm x" by simp
qed

definition nonneg_mat :: "real^'n^'m  bool"
  where "nonneg_mat A = (i j. A $ i $ j  0)"

lemma nonneg_mat_1:
  shows "nonneg_mat (mat 1)"
  unfolding nonneg_mat_def mat_def by auto

lemma nonneg_mat_prod:
  assumes "nonneg_mat A" "nonneg_mat B"
  shows "nonneg_mat (A ** B)"
  using assms unfolding nonneg_mat_def matrix_matrix_mult_def
  by (auto intro:sum_nonneg)

lemma nonneg_mat_transpose:
  "nonneg_mat (transpose A) = nonneg_mat A"
  unfolding nonneg_mat_def transpose_def
  by auto

definition spec_bound :: "real^'n^'n  real  bool"
  where "spec_bound M l = (l  0  (v. v  1 = 0  norm (M *v v)  l * norm v))"

lemma spec_boundD1:
  assumes "spec_bound M l"
  shows "0  l"
  using assms unfolding spec_bound_def by simp

lemma spec_boundD2:
  assumes "spec_bound M l"
  assumes "v  1 = 0 "
  shows "norm (M *v v)  l * norm v"
  using assms unfolding spec_bound_def by simp

lemma spec_bound_mono:
  assumes "spec_bound M α" "α  β"
  shows "spec_bound M β"
proof -
  have "norm (M *v v)  β * norm v" if "inner v 1 = 0"  for v
  proof -
    have "norm (M *v v)  α * norm v"
      by (intro spec_boundD2[OF assms(1)] that)
    also have "...  β * norm v"
      by (intro mult_right_mono assms(2)) auto
    finally show ?thesis by simp
  qed
  moreover have "β  0"
    using assms(2) spec_boundD1[OF assms(1)] by simp
  ultimately show ?thesis
    unfolding spec_bound_def by simp
qed

definition markov :: "real^'n^'n  bool"
  where "markov M = (nonneg_mat M  M *v 1  = 1  1 v* M = 1)"

lemma markov_symI:
  assumes "nonneg_mat A" "transpose A = A" "A *v 1 = 1"
  shows "markov A"
proof -
  have "1 v* A = transpose A *v 1"
    unfolding vector_transpose_matrix[symmetric] by simp
  also have "... = 1" unfolding assms(2,3) by simp
  finally have "1 v* A = 1" by simp
  thus ?thesis
    unfolding markov_def using assms by auto
qed

lemma markov_apply:
  assumes "markov M"
  shows "M *v 1 = 1" "1 v* M = 1"
  using assms unfolding markov_def by auto

lemma markov_transpose:
  "markov A = markov (transpose A)"
  unfolding markov_def nonneg_mat_transpose by auto
fun matrix_pow where
  "matrix_pow M 0 = mat 1" |
  "matrix_pow M (Suc n) = M ** (matrix_pow M n)"

lemma markov_orth_inv:
  assumes "markov A"
  shows "inner (A *v x) 1 = inner x 1"
proof -
  have "inner (A *v x) 1 = inner x (1 v* A)"
    using dot_lmul_matrix inner_commute by metis
  also have "... = inner x 1"
    using markov_apply[OF assms(1)] by simp
  finally show ?thesis by simp
qed

lemma markov_id:
  "markov (mat 1)"
  unfolding markov_def using nonneg_mat_1 by simp

lemma markov_mult:
  assumes "markov A" "markov B"
  shows "markov (A ** B)"
proof -
  have "nonneg_mat (A ** B)"
    using assms unfolding markov_def by (intro nonneg_mat_prod) auto
  moreover have "(A ** B) *v 1 = 1"
    using assms unfolding markov_def
    unfolding matrix_vector_mul_assoc[symmetric] by simp
  moreover have "1 v* (A ** B) = 1"
    using assms unfolding markov_def
    unfolding vector_matrix_mul_assoc[symmetric] by simp
  ultimately show ?thesis
    unfolding markov_def by simp
qed

lemma markov_matrix_pow:
  assumes "markov A"
  shows "markov (matrix_pow A k)"
  using markov_id assms markov_mult
  by (induction k, auto)

lemma spec_bound_prod:
  assumes "markov A" "markov B"
  assumes "spec_bound A la" "spec_bound B lb"
  shows "spec_bound (A ** B) (la*lb)"
proof -
  have la_ge_0: "la  0" using spec_boundD1[OF assms(3)] by simp

  have "norm ((A ** B) *v x)  (la * lb) * norm x" if "inner x 1 = 0" for x
  proof -
    have "norm ((A ** B) *v x) = norm (A *v (B *v x))"
      by (simp add:matrix_vector_mul_assoc)
    also have "...  la * norm (B *v x)"
      by (intro spec_boundD2[OF assms(3)]) (simp add:markov_orth_inv that assms(2))
    also have "...  la * (lb * norm x)"
      by (intro spec_boundD2[OF assms(4)] mult_left_mono that la_ge_0)
    finally show ?thesis by simp
  qed
  moreover have "la * lb  0"
    using la_ge_0 spec_boundD1[OF assms(4)] by simp
  ultimately show ?thesis
    using spec_bound_def by auto
qed

lemma spec_bound_pow:
  assumes "markov A"
  assumes "spec_bound A l"
  shows "spec_bound (matrix_pow A k) (l^k)"
proof (induction k)
  case 0
  then show ?case unfolding spec_bound_def by simp
next
  case (Suc k)
  have "spec_bound (A ** matrix_pow A k) (l * l ^ k)"
    by (intro spec_bound_prod assms Suc markov_matrix_pow)
  thus ?case by simp
qed

fun intersperse :: "'a  'a list  'a list"
  where
    "intersperse x [] = []" |
    "intersperse x (y#[]) = y#[]" |
    "intersperse x (y#z#zs) = y#x#intersperse x (z#zs)"

lemma intersperse_snoc:
  assumes "xs  []"
  shows "intersperse z (xs@[y]) = intersperse z xs@[z,y]"
  using assms
proof (induction xs rule:list_nonempty_induct)
  case (single x)
  then show ?case by simp
next
  case (cons x xs)
  then obtain xsh xst where t:"xs = xsh#xst"
    by (metis neq_Nil_conv)
  have "intersperse z ((x # xs) @ [y]) = x#z#intersperse z (xs@[y])"
    unfolding t by simp
  also have "... = x#z#intersperse z xs@[z,y]"
    using cons by simp
  also have "... = intersperse z (x#xs)@[z,y]"
    unfolding t by simp
  finally show ?case by simp
qed

lemma foldl_intersperse:
  assumes "xs  []"
  shows "foldl f a ((intersperse x xs)@[x]) = foldl (λy z. f (f y z) x) a xs"
  using assms by (induction xs rule:rev_nonempty_induct) (auto simp add:intersperse_snoc)

lemma foldl_intersperse_2:
  shows "foldl f a (intersperse y (x#xs)) = foldl (λx z. f (f x y) z) (f a x) xs"
proof (induction xs rule:rev_induct)
  case Nil
  then show ?case by simp
next
  case (snoc xst xs)
  have "foldl f a (intersperse y ((x # xs) @ [xst])) = foldl (λx. f (f x y)) (f a x) (xs @ [xst])"
    by (subst intersperse_snoc, auto simp add:snoc)
  then show ?case  by simp
qed


context regular_graph_tts
begin

definition stat :: "real^'n"
  where "stat = (1 / real CARD('n)) *R 1"

definition J :: "('c :: field)^'n^'n"
  where "J = (χ i j. of_nat 1 / of_nat CARD('n))"

lemma inner_1_1: "1  (1::real^'n) = CARD('n)"
  unfolding inner_vec_def by simp

definition proj_unit :: "real^'n  real^'n"
  where "proj_unit v = (1  v) *R stat"

definition proj_rem :: "real^'n  real^'n"
  where "proj_rem v = v - proj_unit v"

lemma proj_rem_orth: "1  (proj_rem v) = 0"
  unfolding proj_rem_def proj_unit_def inner_diff_right stat_def
  by (simp add:inner_1_1)

lemma split_vec: "v = proj_unit v + proj_rem v"
  unfolding proj_rem_def by simp

lemma apply_J: "J *v x = proj_unit x"
proof (intro iffD2[OF vec_eq_iff] allI)
  fix i
  have "(J *v x) $ i = inner (χ j. 1 / real CARD('n)) x"
    unfolding matrix_vector_mul_component J_def by simp
  also have "... = inner stat x"
    unfolding stat_def scaleR_vec_def by auto
  also have "... = (proj_unit x) $ i"
    unfolding proj_unit_def stat_def by simp
  finally show "(J *v x) $ i = (proj_unit x) $ i" by simp
qed

lemma spec_bound_J: "spec_bound (J :: real^'n^'n) 0"
proof -
  have "norm (J *v v) = 0" if "inner v 1 = 0" for v :: "real^'n"
  proof -
    have "inner (proj_unit v + proj_rem v) 1 = 0"
      using that by (subst (asm) split_vec[of "v"], simp)
    hence "inner (proj_unit v) 1 = 0"
      using proj_rem_orth inner_commute unfolding inner_add_left
      by (metis add_cancel_left_right)
    hence "proj_unit v = 0"
      unfolding proj_unit_def stat_def by simp
    hence "J *v v = 0"
      unfolding apply_J by simp
    thus ?thesis by simp
  qed
  thus ?thesis
    unfolding spec_bound_def by simp
qed

lemma matrix_decomposition_lemma_aux:
  fixes A :: "real^'n^'n"
  assumes "markov A"
  shows "spec_bound A l  matrix_norm_bound (A - (1-l) *R J) l" (is "?L  ?R")
proof
  assume a:"?L"
  hence l_ge_0: "l  0" using spec_boundD1 by auto
  show "?R"
  proof (rule matrix_norm_boundI)
    fix x :: "real^'n"
    have "(A - (1-l) *R J) *v x = A *v x - (1-l) *R (proj_unit x)"
      by (simp add:algebra_simps vector_scaleR_matrix_ac_2 apply_J)
    also have "... = A *v proj_unit x + A *v proj_rem x - (1-l) *R (proj_unit x)"
      by (subst split_vec[of "x"], simp add:algebra_simps)
    also have "... = proj_unit x + A *v proj_rem x - (1-l) *R (proj_unit x)"
      using markov_apply[OF assms(1)]
      unfolding proj_unit_def stat_def by (simp add:algebra_simps)
    also have "... = A *v proj_rem x + l *R proj_unit x" (is "_ = ?R1")
      by (simp add:algebra_simps)
    finally have d:"(A - (1-l) *R J) *v x = ?R1" by simp

    have "inner (l *R proj_unit x) (A *v proj_rem x) =
      inner ((l * inner 1 x / real CARD('n)) *R 1 v* A) (proj_rem x)"
      by (subst dot_lmul_matrix[symmetric]) (simp add:proj_unit_def stat_def)
    also have "... = (l * inner 1 x / real CARD('n)) * inner 1 (proj_rem x)"
      unfolding scaleR_vector_matrix_assoc markov_apply[OF assms] by simp
    also have "... = 0"
      unfolding proj_rem_orth by simp
    finally have b:"inner (l *R proj_unit x) (A *v proj_rem x) = 0" by simp

    have c: "inner (proj_rem x) (proj_unit x) = 0"
      using proj_rem_orth[of "x"]
      unfolding proj_unit_def stat_def by (simp add:inner_commute)

    have "norm (?R1)^2 = norm (A *v proj_rem x)^2 + norm (l *R proj_unit x)^2"
      using b by (intro pythagoras) (simp add:inner_commute)
    also have "...  (l * norm (proj_rem x))^2 + norm (l *R proj_unit x)^2"
      using proj_rem_orth[of "x"]
      by (intro add_mono power_mono spec_boundD2 a) (auto simp add:inner_commute)
    also have "... = l^2 * (norm (proj_rem x)^2 + norm (proj_unit x)^2)"
      by (simp add:algebra_simps)
    also have "... = l^2 * (norm (proj_rem x + proj_unit x)^2)"
      using c by (subst pythagoras) auto
    also have "... = l^2 * norm x^2"
      by (subst (3) split_vec[of "x"]) (simp add:algebra_simps)
    also have "... = (l * norm x)^2"
      by (simp add:algebra_simps)
    finally have "norm (?R1)^2  (l * norm x)^2" by simp
    hence "norm (?R1)  l * norm x"
      using l_ge_0 by (subst (asm) power_mono_iff) auto

    thus "norm ((A - (1-l) *R J) *v x)  l * norm x"
      unfolding d by simp
  qed
next
  assume a:"?R"
  have "norm (A *v x)  l * norm x" if "inner x 1 = 0" for x
  proof -
    have "(1 - l) *R J *v x = (1 - l) *R (proj_unit x)"
      by (simp add:vector_scaleR_matrix_ac_2 apply_J)
    also have "... = 0"
      unfolding proj_unit_def using that by (simp add:inner_commute)
    finally have b: "(1 - l) *R J *v x = 0" by simp

    have "norm (A *v x) = norm ((A - (1-l) *R J) *v x  + ((1-l) *R J) *v x)"
      by (simp add:algebra_simps)
    also have "...  norm ((A - (1-l) *R J) *v x) + norm (((1-l) *R J) *v x)"
      by (intro norm_triangle_ineq)
    also have "...  l * norm x + 0"
      using a b unfolding  matrix_norm_bound_def by (intro add_mono, auto)
    also have "... = l * norm x"
      by simp
    finally show ?thesis by simp
  qed

  moreover have "l  0"
    using a matrix_norm_bound_nonneg by blast

  ultimately show "?L"
    unfolding spec_bound_def by simp
qed

lemma matrix_decomposition_lemma:
  fixes A :: "real^'n^'n"
  assumes "markov A"
  shows "spec_bound A l  (E. A = (1-l) *R J + l *R E  matrix_norm_bound E 1  l  0)"
    (is "?L  ?R")
proof -
  have "?L  matrix_norm_bound (A - (1-l) *R J) l"
    using matrix_decomposition_lemma_aux[OF assms] by simp
  also have "...  ?R"
  proof
    assume a:"matrix_norm_bound (A - (1 - l) *R J) l"
    hence l_ge_0: "l  0" using matrix_norm_bound_nonneg by auto
    define E where "E = (1/l) *R (A - (1-l) *R J)"
    have "A = J" if "l = 0"
    proof -
      have "matrix_norm_bound (A - J) 0"
        using a that by simp
      hence "A - J = 0" using matrix_norm_bound_0 by blast
      thus "A = J" by simp
    qed
    hence "A = (1-l) *R J + l *R E"
      unfolding E_def by simp
    moreover have "matrix_norm_bound E 1"
    proof (cases "l = 0")
      case True
      hence "E = 0" if "l = 0"
        unfolding E_def by simp
      thus "matrix_norm_bound E 1" if "l = 0"
        using that unfolding matrix_norm_bound_def by auto
    next
      case False
      hence "l > 0" using l_ge_0 by simp
      moreover have "matrix_norm_bound E (¦1 / l¦* l)"
        unfolding E_def
        by (intro matrix_norm_bound_scale a)
      ultimately show ?thesis by auto
    qed
    ultimately show ?R using l_ge_0 by auto
  next
    assume a:?R
    then obtain E where E_def: "A = (1 - l) *R J + l *R E"  "matrix_norm_bound E 1" "l  0"
      by auto
    have "matrix_norm_bound (l *R E) (abs l*1)"
      by (intro matrix_norm_bound_scale E_def(2))
    moreover have "l  0" using E_def by simp
    moreover have " l *R E = (A - (1 - l) *R J)"
      using E_def(1) by simp
    ultimately show "matrix_norm_bound (A - (1 - l) *R J) l"
      by simp
  qed
  finally show ?thesis by simp
qed

lemma hitting_property_alg:
  fixes S :: "('n :: finite) set"
  assumes l_range: "l  {0..1}"
  defines "P  diag (ind_vec S)"
  defines "μ  card S / CARD('n)"
  assumes "M. M  set Ms  spec_bound M l  markov M"
  shows "foldl (λx M. P *v (M *v x)) (P *v stat) Ms  1  (μ + l * (1-μ))^(length Ms+1)"
proof -
  define t :: "real^'n" where "t = (χ i. of_bool (i  S))"
  define r where "r = foldl (λx M. P *v (M *v x)) (P *v stat) Ms"
  have P_proj: "P ** P = P"
    unfolding P_def diag_mult_eq ind_vec_def by (intro arg_cong[where f="diag"]) (vector)

  have P_1_left: "1 v* P = t"
    unfolding P_def diag_def ind_vec_def vector_matrix_mult_def t_def by simp

  have P_1_right: "P *v 1 = t"
    unfolding P_def diag_def ind_vec_def matrix_vector_mult_def t_def by simp

  have P_norm :"matrix_norm_bound P 1"
    unfolding P_def ind_vec_def by (intro matrix_norm_bound_diag) simp

  have norm_t: "norm t = sqrt (real (card S))"
    unfolding t_def norm_vec_def L2_set_def of_bool_def
    by (simp add:sum.If_cases if_distrib if_distribR)

  have μ_range: "μ  0" "μ  1"
    unfolding μ_def by (auto simp add:card_mono)

  define condition :: "real^'n  nat  bool"
    where "condition = (λx n. norm x  (μ + l * (1-μ))^n * sqrt (card S)/CARD('n)  P *v x = x)"

  have a:"condition r (length Ms)"
    unfolding r_def using assms(4)
  proof (induction Ms rule:rev_induct)
    case Nil
    have "norm (P *v stat) = (1 / real CARD('n)) * norm t"
      unfolding stat_def matrix_vector_mult_scaleR P_1_right by simp
    also have "...   (1 / real CARD('n)) * sqrt (real (card S))"
      using  norm_t by (intro mult_left_mono) auto
    also have "... = sqrt (card S)/CARD('n)" by simp
    finally have "norm (P *v stat)  sqrt (card S)/CARD('n)" by simp
    moreover have "P *v (P *v stat) = P *v stat"
      unfolding matrix_vector_mul_assoc P_proj by simp
    ultimately show ?case unfolding condition_def by simp
  next
    case (snoc M xs)
    hence "spec_bound M l  markov M"
        using snoc(2) by simp
    then obtain E where E_def: "M = (1-l) *R J + l *R E" "matrix_norm_bound E 1"
      using iffD1[OF matrix_decomposition_lemma] by auto

    define y where "y = foldl (λx M. P *v (M *v x)) (P *v stat) xs"
    have b:"condition y (length xs)"
      using snoc unfolding y_def by simp
    hence a:"P *v y = y" using condition_def by simp

    have "norm (P *v (M *v y)) = norm (P *v ((1-l)*R J *v y) + P *v (l *R E *v y))"
      by (simp add:E_def algebra_simps)
    also have "...  norm (P *v ((1-l)*R J *v y)) + norm (P *v (l *R E *v y)) "
      by (intro norm_triangle_ineq)
    also have "... = (1 - l) * norm (P *v (J *v y)) + l * norm (P *v (E *v y))"
      using l_range
      by (simp add:vector_scaleR_matrix_ac_2 matrix_vector_mult_scaleR)
    also have "... = (1-l) * ¦1  (P *v y)/real CARD('n)¦ * norm t + l * norm (P *v (E *v y))"
      by (subst a[symmetric])
        (simp add:apply_J proj_unit_def stat_def P_1_right matrix_vector_mult_scaleR)
    also have "... = (1-l) * ¦t  y¦/real CARD('n) * norm t + l * norm (P *v (E *v y))"
      by (subst dot_lmul_matrix[symmetric]) (simp add:P_1_left)
    also have "...  (1-l) * (norm t * norm y) / real CARD('n) * norm t + l * (1 * norm (E *v y))"
      using P_norm Cauchy_Schwarz_ineq2 l_range
      by (intro add_mono mult_right_mono mult_left_mono divide_right_mono matrix_norm_boundD) auto
    also have "... = (1-l) * μ * norm y + l * norm (E *v y)"
      unfolding μ_def norm_t by simp
    also have "...  (1-l) * μ * norm y + l * (1 * norm y)"
      using μ_range l_range
      by (intro add_mono matrix_norm_boundD mult_left_mono E_def) auto
    also have "... = (μ + l * (1-μ)) * norm y"
      by (simp add:algebra_simps)
    also have "...  (μ + l * (1-μ)) * ((μ + l * (1-μ))^length xs * sqrt (card S)/CARD('n))"
      using b μ_range l_range unfolding condition_def
      by (intro mult_left_mono) auto
    also have "... = (μ + l * (1-μ))^(length xs +1) * sqrt (card S)/CARD('n)"
      by simp
    finally have "norm (P *v (M *v y))  (μ + l * (1-μ))^(length xs +1) * sqrt (card S)/CARD('n)"
      by simp

    moreover have "P *v (P *v (M *v y)) = P *v (M *v y)"
      unfolding matrix_vector_mul_assoc matrix_mul_assoc P_proj
      by simp

    ultimately have "condition (P *v (M *v y)) (length (xs@[M]))"
      unfolding condition_def by simp

    then show ?case
      unfolding y_def by simp
  qed

  have "inner r 1 = inner (P *v r) 1"
    using a condition_def by simp
  also have "... = inner (1 v* P) r"
    unfolding dot_lmul_matrix by (simp add:inner_commute)
  also have "... = inner t r"
    unfolding P_1_left by simp
  also have "...  norm t * norm r"
    by (intro norm_cauchy_schwarz)
  also have "...  sqrt (card S) * ((μ + l * (1-μ))^(length Ms) * sqrt(card S)/CARD('n))"
    using a unfolding condition_def norm_t
    by (intro mult_mono) auto
  also have "... = (μ + 0) * ((μ + l * (1-μ))^(length Ms))"
    by (simp add:μ_def)
  also have "...  (μ + l * (1-μ)) * (μ + l * (1-μ))^(length Ms)"
    using μ_range l_range
    by (intro mult_right_mono zero_le_power add_mono) auto
  also have "... = (μ + l * (1-μ))^(length Ms+1)" by simp
  finally show ?thesis
    unfolding r_def by simp
qed

lemma upto_append:
  assumes "i  j" "j  k"
  shows  "[i..<j]@[j..<k] = [i..<k]"
  using assms by (metis less_eqE upt_add_eq_append)

definition bool_list_split :: "bool list  (nat list × nat)"
  where "bool_list_split xs = foldl (λ(ys,z) x. (if x then (ys@[z],0) else (ys,z+1))) ([],0) xs"

lemma bool_list_split:
  assumes "bool_list_split xs = (ys,z)"
  shows "xs = concat (map (λk. replicate k False@[True]) ys)@replicate z False"
  using assms
proof (induction xs arbitrary: ys z rule:rev_induct)
  case Nil
  then show ?case unfolding bool_list_split_def by simp
next
  case (snoc x xs)
  obtain u v where uv_def: "bool_list_split xs = (u,v)"
    by (metis surj_pair)

  show ?case
  proof (cases x)
    case True
    have a:"ys = u@[v]" "z = 0"
      using snoc(2) True uv_def unfolding bool_list_split_def by auto
    have "xs@[x] = concat (map (λk. replicate k False@[True]) u)@replicate v False@[True]"
      using snoc(1)[OF uv_def] True by simp
    also have "... = concat (map (λk. replicate k False@[True]) (u@[v]))@replicate 0 False"
      by simp
    also have "... = concat (map (λk. replicate k False@[True]) (ys))@replicate z False"
      using a by simp
    finally show ?thesis by simp
  next
    case False
    have a:"ys = u" "z = v+1"
      using snoc(2) False uv_def unfolding bool_list_split_def by auto
    have "xs@[x] = concat (map (λk. replicate k False@[True]) u)@replicate (v+1) False"
      using snoc(1)[OF uv_def] False unfolding replicate_add by simp
    also have "... = concat (map (λk. replicate k False@[True]) (ys))@replicate z False"
      using a by simp
    finally show ?thesis by simp
  qed
qed

lemma bool_list_split_count:
  assumes "bool_list_split xs = (ys,z)"
  shows "length (filter id xs) = length ys"
  unfolding bool_list_split[OF assms(1)] by (simp add:filter_concat comp_def)

lemma foldl_concat:
  "foldl f a (concat xss) = foldl (λy xs. foldl f y xs) a xss"
  by (induction xss rule:rev_induct, auto)

lemma hitting_property_alg_2:
  fixes S :: "('n :: finite) set" and l :: nat
  fixes M :: "real^'n^'n"
  assumes α_range: "α  {0..1}"
  assumes "I  {..<l}"
  defines "P i  (if i  I then diag (ind_vec S) else mat 1)"
  defines "μ  real (card S) / real (CARD('n))"
  assumes "spec_bound M α" "markov M"
  shows
    "foldl (λx M. M *v x) stat (intersperse M (map P [0..<l]))  1  (μ+α*(1-μ))^card I"
    (is "?L  ?R")
proof (cases "I  {}")
  case True
  define xs where "xs = map (λi. i  I) [0..<l]"
  define Q where "Q = diag (ind_vec S)"
  define P' where "P' = (λx. if x then Q else mat 1)"

  let ?rep = "(λx. replicate x (mat 1))"

  have P_eq: "P i = P' (i  I)" for i
    unfolding P_def P'_def Q_def by simp

  have "l > 0"
    using True assms(2) by auto
  hence xs_ne: "xs  []"
    unfolding xs_def by simp

  obtain ys z where ys_z: "bool_list_split xs = (ys,z)"
    by (metis surj_pair)


  have "length ys = length (filter id xs)"
    using bool_list_split_count[OF ys_z] by simp
  also have "... = card (I  {0..<l})"
    unfolding xs_def filter_map by (simp add:comp_def distinct_length_filter)
  also have "... = card I"
    using Int_absorb2[OF assms(2)] unfolding atLeast0LessThan by simp
  finally have  len_ys: "length ys = card I" by simp

  hence "length ys > 0"
    using True assms(2) by (metis card_gt_0_iff finite_nat_iff_bounded)
  then obtain yh yt where ys_split: "ys = yh#yt"
    by (metis length_greater_0_conv neq_Nil_conv)

  have a:"foldl (λx N. M *v (N *v x)) x (?rep z)  1 = x  1" for x
  proof (induction z)
    case 0
    then show ?case by simp
  next
    case (Suc z)
    have "foldl (λx N. M *v (N *v x)) x (?rep (z+1))  1 = x  1"
      unfolding replicate_add using Suc
      by (simp add:markov_orth_inv[OF assms(6)])
    then show ?case by simp
  qed

  have "M *v stat = stat"
    using assms(6) unfolding stat_def matrix_vector_mult_scaleR markov_def by simp
  hence b: "foldl (λx N. M *v (N *v x)) stat (?rep yh) = stat"
    by (induction yh, auto)

  have "foldl (λx N. N *v (M *v x)) a (?rep x) = matrix_pow M x *v a" for x a
  proof (induction x)
    case 0
    then show ?case by simp
  next
    case (Suc x)
    have "foldl (λx N. N *v (M *v x)) a (?rep (x+1)) =  matrix_pow M (x+1) *v a"
      unfolding replicate_add using Suc by (simp add: matrix_vector_mul_assoc)
    then show ?case by simp
  qed
  hence c: "foldl (λx N. N *v (M *v x)) a (?rep x @ [Q]) = Q *v (matrix_pow M (x+1) *v a)" for x a
    by (simp add:matrix_vector_mul_assoc matrix_mul_assoc)

  have d: "spec_bound N α  markov N" if t1: "N  set (map (λx. matrix_pow M (x + 1)) yt)" for N
  proof -
    obtain y where N_def: "N = matrix_pow M (y+1)"
      using t1 by auto
    hence d1: "spec_bound N (α^(y+1))"
      unfolding N_def using spec_bound_pow assms(5,6) by blast
    have "spec_bound N (α^1)"
      using α_range by (intro spec_bound_mono[OF d1] power_decreasing) auto
    moreover have "markov N"
      unfolding N_def by (intro markov_matrix_pow assms(6))
    ultimately show ?thesis by simp
  qed

  have "?L = foldl (λx M. M *v x) stat (intersperse M (map P' xs))  1"
    unfolding P_eq xs_def map_map by (simp add:comp_def)
  also have "... = foldl (λx M. M *v x) stat (intersperse M (map P' xs)@[M])  1"
    by (simp add:markov_orth_inv[OF assms(6)])
  also have "... = foldl (λx N. M *v (N *v x)) stat (map P' xs)  1"
    using xs_ne by (subst foldl_intersperse) auto
  also have "... = foldl (λx N. M *v (N *v x)) stat ((ys  (λx. ?rep x @ [Q])) @ ?rep z)  1"
    unfolding bool_list_split[OF ys_z] P'_def List.bind_def by (simp add: comp_def map_concat)
  also have "... = foldl (λx N. M *v (N *v x)) stat (ys  (λx. ?rep x @ [Q]))  1"
    by (simp add: a)
  also have "... = foldl (λx N. M *v (N *v x)) stat (?rep yh @[Q]@(yt (λx. ?rep x @ [Q])))  1"
    unfolding ys_split by simp
  also have "... = foldl (λx N. M *v (N *v x)) stat ([Q]@(yt (λx. ?rep x @ [Q])))  1"
    by (simp add:b)
  also have "... = foldl (λx N. N *v x) stat (intersperse M (Q#(yt (λx.?rep x@[Q])))@[M])1"
    by (subst foldl_intersperse, auto)
  also have "... = foldl (λx N. N *v x) stat (intersperse M (Q#(yt (λx.?rep x@[Q]))))  1"
    by (simp add:markov_orth_inv[OF assms(6)])
  also have "... = foldl (λx N. N *v (M *v x)) (Q *v stat) (yt (λx.?rep x@[Q]))  1"
    by (subst foldl_intersperse_2, simp)
  also have "... = foldl (λa x. foldl (λx N. N *v (M *v x)) a (?rep x @ [Q])) (Q *v stat) yt  1"
    unfolding List.bind_def foldl_concat foldl_map by simp
  also have "... = foldl (λa x. Q *v (matrix_pow M (x+1) *v a)) (Q *v stat) yt  1"
    unfolding c by simp
  also have "... = foldl (λa N. Q *v (N *v a)) (Q *v stat) (map (λx. matrix_pow M (x+1)) yt)  1"
    by (simp add:foldl_map)
  also have "...  (μ + α*(1-μ))^(length (map (λx. matrix_pow M (x+1)) yt)+1)"
    unfolding μ_def Q_def by (intro hitting_property_alg α_range d) simp
  also have "... = (μ + α*(1-μ))^(length ys)"
    unfolding ys_split by simp
  also have "... = ?R" unfolding len_ys by simp
  finally show ?thesis by simp
next
  case False
  hence I_empty: "I = {}" by simp

  have "?L = stat  (1 :: real^'n)"
  proof (cases "l > 0")
    case True
    have "?L = foldl (λx M. M *v x) stat ((intersperse M (map P [0..<l]))@[M])  1"
      by (simp add:markov_orth_inv[OF assms(6)])
    also have "... = foldl (λx N. M *v (N *v x)) stat (map P [0..<l])   1"
      using True by (subst foldl_intersperse, auto)
    also have "... = foldl (λx N. M *v (N *v x)) stat (map (λ_. mat 1) [0..<l])   1"
      unfolding  P_def using I_empty by simp
    also have "... = foldl (λx _. M *v x) stat [0..<l]  1"
      unfolding foldl_map by simp
    also have "... = stat  (1 :: real^'n)"
      by (induction l, auto simp add:markov_orth_inv[OF assms(6)])
    finally show ?thesis by simp
  next
    case False
    then show ?thesis by simp
  qed
  also have "... = 1"
    unfolding stat_def by (simp add:inner_vec_def)
  also have "...  ?R" unfolding I_empty by simp
  finally show ?thesis by simp
qed

lemma uniform_property_alg:
  fixes x :: "('n :: finite)" and l :: nat
  assumes "i < l"
  defines "P j  (if j = i then diag (ind_vec {x}) else mat 1)"
  assumes "markov M"
  shows "foldl (λx M. M *v x) stat (intersperse M (map P [0..<l]))  1 = 1 / CARD('n)"
    (is "?L = ?R")
proof -
  have a:"l > 0" using assms(1) by simp

  have 0: "foldl (λx N. M *v (N *v x)) y (xs)  1 = y   1" if "set xs  {mat 1}" for xs y
    using that
  proof (induction xs rule:rev_induct)
    case Nil
    then show ?case by simp
  next
    case (snoc x xs)
    have "x = mat 1"
      using snoc(2) by simp
    hence "foldl (λx N. M *v (N *v x)) y (xs @ [x])  1 = foldl (λx N. M *v (N *v x)) y xs  1"
      by (simp add:markov_orth_inv[OF assms(3)])
    also have "... = y  1"
      using snoc(2) by (intro snoc(1)) auto
    finally show ?case by simp
  qed

  have M_stat: "M *v stat = stat"
    using assms(3) unfolding stat_def matrix_vector_mult_scaleR markov_def by simp

  hence 1: "(foldl (λx N. M *v (N *v x)) stat xs) = stat" if "set xs  {mat 1}" for xs
    using that by (induction xs, auto)

  have "?L = foldl (λx M. M *v x) stat ((intersperse M (map P [0..<l]))@[M])  1"
    by (simp add:markov_orth_inv[OF assms(3)])
  also have "... = foldl (λx N. M *v (N *v x)) stat (map P [0..<l])  1"
    using a by (subst foldl_intersperse) auto
  also have "... = foldl (λx N. M *v (N *v x)) stat (map P ([0..<i+1]@[i+1..<l]))  1"
    using assms(1) by (subst upto_append) auto
  also have "... = foldl (λx N. M *v (N *v x)) stat (map P [0..<i + 1])  1"
    unfolding map_append foldl_append  P_def by (subst 0) auto
  also have "... = foldl (λx N. M *v (N *v x)) stat (map P ([0..<i]@[i]))  1"
    by simp
  also have "... = (M *v (diag (ind_vec {x}) *v stat))  1"
    unfolding map_append foldl_append P_def by (subst 1) auto
  also have "... = (diag (ind_vec {x}) *v stat)  1"
    by (simp add:markov_orth_inv[OF assms(3)])
  also have "... = ((1/CARD('n)) *R ind_vec {x})  1"
    unfolding diag_def ind_vec_def  stat_def matrix_vector_mult_def
    by (intro arg_cong2[where f="(∙)"] refl)
      (vector of_bool_def sum.If_cases if_distrib if_distribR)
  also have "... = (1/CARD('n)) * (ind_vec {x}  1)"
    by simp
  also have "... = (1/CARD('n)) * 1"
    unfolding inner_vec_def ind_vec_def of_bool_def
    by (intro arg_cong2[where f="(*)"] refl) (simp)
  finally show ?thesis by simp
qed

end

lemma foldl_matrix_mult_expand:
  fixes Ms :: "(('r::{semiring_1,comm_monoid_mult})^'a^'a) list"
  shows "(foldl (λx M. M *v x) a Ms) $ k = (x | length x = length Ms+1  x! length Ms = k.
  ( i< length Ms. (Ms ! i) $ (x ! (i+1)) $ (x ! i)) * a $ (x ! 0))"
proof (induction Ms arbitrary: k rule:rev_induct)
  case Nil
  have "length x = Suc 0  x = [x!0]" for x :: "'a list"
    by (cases x, auto)
  hence "{x. length x = Suc 0  x ! 0 = k} = {[k]}"
    by auto
  thus ?case by auto
next
  case (snoc M Ms)
  let ?l = "length Ms"

  have 0: "finite {w. length w = Suc (length Ms)  w ! length Ms = i}" for i :: 'a
    using finite_lists_length_eq[where A="UNIV::'a set" and n="?l +1"] by simp

  have "take (?l+1) x @ [x ! (?l+1)] = x" if "length x = ?l+2" for x :: "'a list"
  proof -
    have "take (?l+1) x @ [x ! (?l+1)] = take (Suc (?l+1)) x"
      using that by (intro take_Suc_conv_app_nth[symmetric], simp)
    also have "... = x"
      using that by simp
    finally show ?thesis by simp
  qed
  hence 1: "bij_betw  (take (?l+1)) {w. length w=?l+2  w!(?l+1) =k} {w. length w = ?l+1}"
    by (intro bij_betwI[where g="λx. x@[k]"]) (auto simp add:nth_append)

  have "foldl (λx M. M *v x) a (Ms @ [M]) $ k = (jUNIV. M$k$j *(foldl (λx M. M *v x) a Ms $ j))"
    by (simp add:matrix_vector_mult_def)
  also have "... =
    (jUNIV. M$k$j * (w|length w=?l+1w!?l=j. (i<?l. Ms!i $ w!(i+1) $ w!i) * a $ w!0))"
    unfolding snoc by simp
  also have "... =
    (jUNIV. (w|length w=?l+1w!?l=j. M$k$w!?l * (i<?l. Ms!i $ w!(i+1) $ w!i) * a $ w!0))"
    by (intro sum.cong refl) (simp add: sum_distrib_left algebra_simps)
  also have "... = (w (j  UNIV. {w. length w=?l+1  w!?l =j}).
    M$k$w!?l*(i<?l. Ms!i $ w!(i+1) $ w!i) * a $ w!0)"
    using 0 by (subst sum.UNION_disjoint, simp, simp) auto
  also have "... = (w | length w=?l+1. M$k$(w!?l)*(i<?l. Ms!i $ w!(i+1) $ w!i) * a $ w!0)"
    by (intro sum.cong arg_cong2[where f="(*)"] refl) auto
  also have "... = (w  take (?l+1) ` {w. length w=?l+2  w!(?l+1) =k}.
    M$k$w!?l*(i<?l. Ms!i $ w!(i+1) $ w!i) * a $ w!0)"
    using 1 unfolding bij_betw_def by (intro sum.cong refl, auto)
  also have "... = (w|length w=?l+2w!(?l+1)=k. M$k$w!?l*(i<?l. Ms!i $ w!(i+1) $ w!i)* a$w!0)"
    using 1 unfolding bij_betw_def by (subst sum.reindex, auto)
  also have "... = (w|length w=?l+2w!(?l+1)=k.
    (Ms@[M])!?l$k$w!?l*(i<?l. (Ms@[M])!i $ w!(i+1) $ w!i)* a$w!0)"
    by (intro sum.cong arg_cong2[where f="(*)"] prod.cong refl) (auto simp add:nth_append)
  also have "... = (w|length w=?l+2w!(?l+1)=k. (i<(?l+1). (Ms@[M])!i $ w!(i+1) $ w!i)* a$w!0)"
    by (intro sum.cong, auto simp add:algebra_simps)
  finally have "foldl (λx M. M *v x) a (Ms @ [M]) $ k =
    ( w | length w = ?l+2  w ! (?l+1) = k. (i<(?l+1). (Ms@[M])!i $ w!(i+1) $ w!i)* a$w!0)"
    by simp
  then show ?case by simp
qed

lemma foldl_matrix_mult_expand_2:
  fixes Ms :: "(real^'a^'a) list"
  shows "(foldl (λx M. M *v x) a Ms)  1 = (x | length x = length Ms+1.
          ( i< length Ms. (Ms ! i) $ (x ! (i+1)) $ (x ! i)) * a $ (x ! 0))"
  (is "?L = ?R")
proof -
  let ?l = "length Ms"
  have "?L = (j  UNIV. (foldl (λx M. M *v x) a Ms) $ j)"
    by (simp add:inner_vec_def)
  also have "... = (jUNIV. x|length x=?l+1  x!?l=j.(i<?l. Ms!i $ x!(i+1) $ x!i) * a $ x!0)"
    unfolding foldl_matrix_mult_expand by simp
  also have "... = (x  (j UNIV.{w. length w = length Ms+1  w ! length Ms = j}).
          ( i< length Ms. (Ms ! i) $ (x ! (i+1)) $ (x ! i)) * a $ (x ! 0))"
    using finite_lists_length_eq[where A="UNIV::'a set" and n="?l +1"]
    by (intro sum.UNION_disjoint[symmetric]) auto
  also have "... = ?R"
    by (intro sum.cong, auto)
  finally show ?thesis by simp
qed

end