Theory SquareMatrix

header {* Square Matrices *}

theory SquareMatrix
imports VectorType
begin

typedef (open) ('a,'n) smatrix = "UNIV :: ('n => 'n => 'a) set" ..

declare Abs_smatrix_inverse [simplified, iff]

lemma expand_smatrix_eq:
  "A = B <-> (∀i j. Rep_smatrix A i j = Rep_smatrix B i j)"
by (simp add: Rep_smatrix_inject [symmetric] expand_fun_eq)

lemma smatrix_ext:
  "(!!i j. Rep_smatrix A i j = Rep_smatrix B i j) ==> A = B"
by (simp add: expand_smatrix_eq)


subsection {* Matrix arithmetic *}

instantiation smatrix :: (zero, type) zero
begin

definition
  zero_smatrix_def: "0 = Abs_smatrix (λi j. 0)"

instance ..
end

instantiation smatrix :: (plus, type) plus
begin

definition
  plus_smatrix_def:
    "A + B = Abs_smatrix (λi j. Rep_smatrix A i j + Rep_smatrix B i j)"

instance ..
end

instantiation smatrix :: (minus, type) minus
begin

definition
  minus_smatrix_def:
    "A - B = Abs_smatrix (λi j. Rep_smatrix A i j - Rep_smatrix B i j)"

instance ..
end

instantiation smatrix :: (uminus, type) uminus
begin

definition
  uminus_smatrix_def:
    "- A = Abs_smatrix (λi j. - Rep_smatrix A i j)"

instance ..
end

lemma Rep_smatrix_zero [simp]: "Rep_smatrix 0 i j = 0"
  unfolding zero_smatrix_def by simp

lemma Rep_smatrix_add [simp]:
  "Rep_smatrix (A + B) i j = Rep_smatrix A i j + Rep_smatrix B i j"
  unfolding plus_smatrix_def by simp

lemma Rep_smatrix_diff [simp]:
  "Rep_smatrix (A - B) i j = Rep_smatrix A i j - Rep_smatrix B i j"
  unfolding minus_smatrix_def by simp

lemma Rep_smatrix_uminus [simp]:
  "Rep_smatrix (- A) i j = - Rep_smatrix A i j"
  unfolding uminus_smatrix_def by simp

instance smatrix :: (semigroup_add, type) semigroup_add
by (intro_classes, simp add: expand_smatrix_eq add_assoc)

instance smatrix :: (ab_semigroup_add, type) ab_semigroup_add
by (intro_classes, simp add: expand_smatrix_eq add_commute)

instance smatrix :: (comm_monoid_add, type) comm_monoid_add
by (intro_classes, simp add: expand_smatrix_eq)

instance smatrix :: (cancel_semigroup_add, type) cancel_semigroup_add
by (intro_classes, simp_all add: expand_smatrix_eq)

instance smatrix :: (cancel_ab_semigroup_add, type) cancel_ab_semigroup_add
by (intro_classes, simp_all add: expand_smatrix_eq)

instance smatrix :: (ab_group_add, type) ab_group_add
by (intro_classes, simp_all add: expand_smatrix_eq)

subsection {* Smatrix multiplication *}

instantiation smatrix :: (semiring_0, finite) semiring_0
begin

definition
  times_smatrix_def:
    "A * B = Abs_smatrix
      (λi j. ∑k∈UNIV. Rep_smatrix A i k * Rep_smatrix B k j)"

lemma Rep_smatrix_mult [simp]:
  "Rep_smatrix (A * B) i j =
    (∑k∈UNIV. Rep_smatrix A i k * Rep_smatrix B k j)"
  unfolding times_smatrix_def by simp

instance proof
  fix A B C :: "('a, 'b) smatrix"
  show "0 * A = 0"
    by (rule smatrix_ext) simp
  show "A * 0 = 0"
    by (rule smatrix_ext) simp
  show "(A * B) * C = A * (B * C)"
    apply (rule smatrix_ext, simp)
    apply (subst setsum_left_distrib)
    apply (subst setsum_right_distrib)
    apply (subst mult_assoc)
    apply (rule setsum_commute)
    done
  show "(A + B) * C = A * C + B * C"
    by (rule smatrix_ext)
       (simp add: left_distrib setsum_addf)
  show "A * (B + C) = A * B + A * C"
    by (rule smatrix_ext)
       (simp add: right_distrib setsum_addf)
qed

end

instance smatrix :: (semiring_0_cancel, finite) semiring_0_cancel ..

instance smatrix :: (ring, finite) ring ..

instantiation smatrix :: (semiring_1, finite) semiring_1
begin

definition
  one_smatrix_def:
    "1 = Abs_smatrix (λi j. if i = j then 1 else 0)"

lemma Rep_smatrix_one [simp]:
  "Rep_smatrix 1 i j = (if i = j then 1 else 0)"
  unfolding one_smatrix_def by simp

instance proof
  fix A :: "('a, 'b) smatrix"
  show "1 * A = A"
    apply (rule smatrix_ext, simp)
    apply (cut_tac x=i in UNIV_I)
    apply (erule setsum_diff1' [OF finite, THEN ssubst])
    apply simp
    done
  show "A * 1 = A"
    apply (rule smatrix_ext, simp)
    apply (cut_tac x=j in UNIV_I)
    apply (erule setsum_diff1' [OF finite, THEN ssubst])
    apply simp
    done
  show "(0::('a, 'b) smatrix) ≠ 1"
    by (simp add: expand_smatrix_eq)
qed

end

instance smatrix :: (semiring_1_cancel, finite) semiring_1_cancel ..

instance smatrix :: (ring_1, finite) ring_1 ..


subsection {* Numeral syntax for matrices *}

instantiation smatrix :: (ring_1, finite) number
begin

definition
  number_of_smatrix_def: "(number_of w::(_,_)smatrix) = of_int w"

instance ..
end

(*instance smatrix :: (ring_1, finite) number_ring*)

text {* Unfortunately class @{text number_ring}
 requires commutativity of multiplication. *}


subsection {* Matrices are a real vector space *}

instantiation smatrix :: (scaleR, finite) scaleR
begin

definition
  scaleR_smatrix_def:
   "scaleR r A = Abs_smatrix (λi j. scaleR r (Rep_smatrix A i j))"

instance ..
end

lemma Rep_smatrix_scaleR [simp]:
  "Rep_smatrix (scaleR r A) i j = scaleR r (Rep_smatrix A i j)"
  unfolding scaleR_smatrix_def by simp

instance smatrix :: (real_vector, finite) real_vector
apply (intro_classes)
apply (simp_all add: expand_smatrix_eq)
apply (simp add: scaleR_right_distrib)
apply (simp add: scaleR_left_distrib)
done

instance smatrix :: (real_algebra, finite) real_algebra
by intro_classes (simp_all add: expand_smatrix_eq scaleR_right.setsum)

instance smatrix :: (real_algebra_1, finite) real_algebra_1 ..


subsection {* Applying a matrix to a vector *}

definition
  mmult :: "('a::semiring_0, 'n::finite) smatrix => 'a ^ 'n => 'a ^ 'n"
where
  "mmult A x = Abs_vec (λi. ∑j∈UNIV. Rep_smatrix A i j * Rep_vec x j)"

lemma Rep_vec_mmult [simp]:
  "Rep_vec (mmult A x) i = (∑j∈UNIV. Rep_smatrix A i j * Rep_vec x j)"
  unfolding mmult_def by simp

lemma mmult_mmult: "mmult A (mmult B x) = mmult (A * B) x"
apply (rule vec_ext)
apply (simp add: setsum_right_distrib setsum_left_distrib mult_assoc)
apply (rule setsum_commute)
done

lemma mmult_right_distrib: "mmult A (x + y) = mmult A x + mmult A y"
by (rule vec_ext, simp add: right_distrib setsum_addf)

lemma mmult_left_distrib: "mmult (A + B) x = mmult A x + mmult B x"
by (rule vec_ext, simp add: left_distrib setsum_addf)

interpretation mmult_right:
  additive ["λx. mmult A x :: 'a::ring ^ 'n::finite"]
by unfold_locales (rule mmult_right_distrib)

interpretation mmult_left:
  additive ["λA. mmult A x :: 'a::ring ^ 'n::finite"]
by unfold_locales (rule mmult_left_distrib)

declare mmult_right.zero [simp]
declare mmult_left.zero [simp]

lemma mmult_vec1: "mmult A (vec1 j a) = Abs_vec (λi. Rep_smatrix A i j * a)"
apply (rule vec_ext, simp)
apply (cut_tac x="j" in UNIV_I)
apply (erule setsum_diff1' [OF finite, THEN ssubst])
apply simp
done

text {* Class @{text semiring_1} is required, since
 class @{text semiring_0} would allow multiplication
 in the ring to always return zero. *}

lemma smatrix_mmult_ext:
  fixes A B :: "('a::semiring_1, 'n::finite) smatrix"
  assumes mmult_eq: "!!x. mmult A x = mmult B x"
  shows "A = B"
proof (rule smatrix_ext)
  fix i j
  have "mmult A (vec1 j 1) = mmult B (vec1 j 1)"
    by (rule mmult_eq)
  thus "Rep_smatrix A i j = Rep_smatrix B i j"
    by (simp add: mmult_vec1, simp add: expand_vec_eq)
qed

lemma expand_smatrix_eq_mmult:
  fixes A B :: "('a::semiring_1, 'n::finite) smatrix"
  shows "A = B <-> (∀x. mmult A x = mmult B x)"
by (auto intro!: smatrix_mmult_ext)

lemma mmult_right_scaleR:
  fixes A :: "('a::real_algebra, 'n::finite) smatrix"
  shows "mmult A (scaleR r x) = scaleR r (mmult A x)"
by (rule vec_ext, simp add: scaleR_right.setsum)

lemma mmult_left_scaleR:
  fixes A :: "('a::real_algebra, 'n::finite) smatrix"
  shows "mmult (scaleR r A) x = scaleR r (mmult A x)"
by (rule vec_ext, simp add: scaleR_right.setsum)

lemma mmult_1_left [simp]:
  fixes x :: "'a::semiring_1 ^ 'n::finite"
  shows "mmult 1 x = x"
apply (rule vec_ext)
apply simp
apply (cut_tac x=i in UNIV_I)
apply (erule setsum_diff1' [OF finite, THEN ssubst])
apply simp
done

lemma bounded_mmult:
  fixes A :: "('a::real_normed_algebra,'n::finite) smatrix"
  shows "∃K. ∀x. norm (mmult A x) ≤ norm x * K"
proof (intro exI allI)
  fix x
  have "norm (mmult A x) = set_l2 (λi. norm (Rep_vec (mmult A x) i)) UNIV"
    unfolding norm_vec_def ..
  also have "… ≤ (∑i∈UNIV. ¦norm (Rep_vec (mmult A x) i)¦)"
    by (rule set_l2_le_setsum_abs)
  also have "… = (∑i∈UNIV. norm (∑j∈UNIV. Rep_smatrix A i j * Rep_vec x j))"
    by simp
  also have "… ≤ (∑i∈UNIV. ∑j∈UNIV. norm (Rep_smatrix A i j) * norm x)"
    by (intro setsum_mono order_trans [OF norm_setsum]
              order_trans [OF norm_mult_ineq] mult_left_mono
              norm_Rep_vec_le norm_ge_zero)
  finally show "norm (mmult A x) ≤
                norm x * (∑i∈UNIV. ∑j∈UNIV. norm (Rep_smatrix A i j))"
    by (simp add: setsum_right_distrib mult_commute)
qed

lemma bounded_linear_mmult_right:
  fixes A :: "('a::real_normed_algebra,'n::finite) smatrix"
  shows "bounded_linear (mmult A)"
apply (unfold_locales)
apply (rule mmult_right_scaleR)
apply (rule bounded_mmult)
done

interpretation mmult_right:
  bounded_linear ["mmult (A::('a::real_normed_algebra,'n::finite) smatrix)"]
by (rule bounded_linear_mmult_right)


subsection {* Matrix norm *}

instantiation smatrix :: (real_normed_algebra, finite) norm
begin

definition
  norm_smatrix_def:
    "norm A = (THE r. isLub UNIV {norm (mmult A x) |x. norm x ≤ 1} r)"

instance ..
end

lemma reals_complete1:
  assumes notempty_S: "∃x. x ∈ S"
    and exists_Ub: "∃y. isUb (UNIV::real set) S y"
  shows "∃!t. isLub (UNIV::real set) S t"
apply (rule ex_ex1I)
apply (rule reals_complete [OF notempty_S exists_Ub])
apply (erule (1) real_isLub_unique)
done

lemma isLub_smatrix_norm:
  "isLub UNIV {norm (mmult A x) |x. norm x ≤ 1} (norm A)"
unfolding norm_smatrix_def
apply (rule theI')
apply (rule reals_complete1)
apply (rule_tac x="0" in exI, simp)
apply (rule_tac x="0" in exI, simp)
apply (cut_tac mmult_right.nonneg_bounded [of A])
apply (erule exE, rule_tac x="K" in exI)
apply (rule isUbI)
apply (rule setleI)
apply clarify
apply (drule spec)
apply (erule order_trans)
apply (rule ord_le_eq_trans)
apply (erule (1) mult_right_mono)
apply simp
apply simp
done

lemma smatrix_norm_geI:
  fixes A :: "('a::real_normed_algebra, 'n::finite) smatrix"
  shows "∃x. y = norm (mmult A x) ∧ norm x ≤ 1 ==> y ≤ norm A"
by (rule isLub_smatrix_norm [THEN isLubD2], simp)

lemma smatrix_norm_leI:
  fixes A :: "('a::real_normed_algebra, 'n::finite) smatrix"
  assumes "!!x. norm x ≤ 1 ==> norm (mmult A x) ≤ y"
  shows "norm A ≤ y"
apply (rule isLub_smatrix_norm [THEN isLub_le_isUb])
apply (rule isUbI [OF setleI])
apply (clarify elim!: prems)
apply simp
done

lemma norm_mmult_ineq: "norm (mmult A x) ≤ norm A * norm x"
apply (case_tac "x = 0", simp)
apply (subgoal_tac "norm (mmult A x) / norm x ≤ norm A")
apply (simp add: pos_divide_le_eq)
apply (rule smatrix_norm_geI)
apply (rule_tac x="scaleR (inverse (norm x)) x" in exI)
apply (simp add: norm_scaleR)
apply (simp add: mmult_right_scaleR norm_scaleR)
apply (simp add: nonzero_divide_eq_eq)
done

lemma smatrix_norm_ge_zero:
  fixes A :: "('a::real_normed_algebra, 'n::finite) smatrix"
  shows "0 ≤ norm A"
apply (rule smatrix_norm_geI)
apply (rule_tac x="0" in exI, simp)
done

lemma smatrix_norm_eq_zero:
  fixes A :: "('a::real_normed_algebra_1, 'n::finite) smatrix"
  shows "norm A = 0 <-> A = 0"
proof (safe)
  show "norm (0::('a, 'n) smatrix) = 0"
    apply (rule order_antisym)
    apply (rule smatrix_norm_leI, simp)
    apply (rule smatrix_norm_ge_zero)
    done
next
  assume A: "norm A = 0"
  have "∀x. norm (mmult A x) ≤ norm A * norm x"
      by (simp add: norm_mmult_ineq)
  with A have "∀x. mmult A x = 0" by simp
  thus "A = 0"
    by (simp add: smatrix_mmult_ext)
qed

lemma smatrix_norm_triangle_ineq:
  fixes A B :: "('a::real_normed_algebra, 'n::finite) smatrix"
  shows "norm (A + B) ≤ norm A + norm B"
apply (rule smatrix_norm_leI)
apply (simp only: mmult_left_distrib)
apply (rule order_trans [OF norm_triangle_ineq])
apply (rule add_mono)
apply (rule smatrix_norm_geI, fast)
apply (rule smatrix_norm_geI, fast)
done

lemma smatrix_norm_scaleR:
  fixes A :: "('a::real_normed_algebra, 'n::finite) smatrix"
  shows "norm (scaleR r A) = ¦r¦ * norm A"
apply (rule order_antisym)
apply (rule smatrix_norm_leI)
apply (simp add: mmult_left_scaleR norm_scaleR)
apply (rule mult_left_mono [OF _ abs_ge_zero])
apply (rule smatrix_norm_geI, fast)
apply (case_tac "r = 0", simp add: smatrix_norm_ge_zero)
apply (subgoal_tac "norm A ≤ norm (scaleR r A) / ¦r¦")
apply (simp add: pos_le_divide_eq mult_commute)
apply (rule smatrix_norm_leI)
apply (simp add: pos_le_divide_eq mult_commute)
apply (rule smatrix_norm_geI)
apply (rule_tac x=x in exI)
apply (simp add: mmult_left_scaleR norm_scaleR)
done

lemma smatrix_norm_mult_ineq:
  fixes A B :: "('a::real_normed_algebra, 'n::finite) smatrix"
  shows "norm (A * B) ≤ norm A * norm B"
apply (rule smatrix_norm_leI)
apply (simp only: mmult_mmult [symmetric])
apply (rule order_trans [OF norm_mmult_ineq])
apply (rule mult_left_mono)
apply (rule smatrix_norm_geI, fast)
apply (rule smatrix_norm_ge_zero)
done

lemma smatrix_norm_1:
  "norm (1::('a::real_normed_algebra_1, 'n::finite) smatrix) = 1"
apply (rule order_antisym)
apply (rule smatrix_norm_leI)
apply simp
apply (rule smatrix_norm_geI)
apply simp
apply (rule_tac x="vec1 i 1" in exI, simp)
done

instantiation smatrix :: (real_normed_algebra_1, finite) real_normed_algebra_1
begin

definition
  sgn_smatrix_def: "sgn (x::(_,_)smatrix) = scaleR (inverse (norm x)) x"

instance
apply intro_classes
apply (rule sgn_smatrix_def)
apply (rule smatrix_norm_ge_zero)
apply (rule smatrix_norm_eq_zero)
apply (rule smatrix_norm_triangle_ineq)
apply (rule smatrix_norm_scaleR)
apply (rule smatrix_norm_mult_ineq)
apply (rule smatrix_norm_1)
done

end

subsection {* Vector outer product *}

definition
  outer :: "'a ^ 'n => 'a ^ 'n => ('a::semiring,'n::finite) smatrix" where
  "outer x y = Abs_smatrix (λi j. Rep_vec x i * Rep_vec y j)"

lemma Rep_smatrix_outer [simp]:
  "Rep_smatrix (outer x y) i j = Rep_vec x i * Rep_vec y j"
unfolding outer_def by simp

lemma outer_right_distrib: "outer x (y + z) = outer x y + outer x z"
by (rule smatrix_ext, simp add: right_distrib)

lemma outer_left_distrib: "outer (x + y) z = outer x z + outer y z"
by (rule smatrix_ext, simp add: left_distrib)

interpretation outer_right:
  additive ["λy. outer x y :: ('a::ring, 'n::finite) smatrix"]
by unfold_locales (rule outer_right_distrib)

interpretation outer_left:
  additive ["λx. outer x y :: ('a::ring, 'n::finite) smatrix"]
by unfold_locales (rule outer_left_distrib)

lemma outer_right_scaleR:
  fixes x y :: "'a::real_algebra ^ 'n::finite"
  shows "outer x (scaleR r y) = scaleR r (outer x y)"
by (rule smatrix_ext, simp)

lemma outer_left_scaleR:
  fixes x y :: "'a::real_algebra ^ 'n::finite"
  shows "outer (scaleR r x) y = scaleR r (outer x y)"
by (rule smatrix_ext, simp)

lemma mmult_outer: "mmult (outer x y) z = vec_map (λa. a * dot y z) x"
unfolding dot_def
apply (rule vec_ext)
apply (simp add: mult_assoc)
apply (simp add: setsum_right_distrib [symmetric])
done

lemma norm_mmult_outer_ineq:
  fixes x y z :: "'a::real_normed_algebra ^ 'n::finite"
  shows "norm (mmult (outer x y) z) ≤ norm x * norm y * norm z"
apply (subst mmult_outer)
apply (subst norm_vec_def)
apply simp
apply (rule order_trans)
apply (rule set_l2_mono)
apply (rule norm_mult_ineq)
apply (rule norm_ge_zero)
apply (subst set_l2_left_distrib [symmetric])
apply (rule norm_ge_zero)
apply (subst norm_vec_def [symmetric])
apply (subst mult_assoc)
apply (rule mult_left_mono)
apply (rule norm_dot_ineq)
apply (rule norm_ge_zero)
done

lemma norm_outer_ineq:
  fixes x y z :: "'a::real_normed_algebra ^ 'n::finite"
  shows "norm (outer x y) ≤ norm x * norm y"
apply (rule smatrix_norm_leI)
apply (rule order_trans)
apply (rule norm_mmult_outer_ineq)
apply (rule order_trans)
apply (erule mult_left_mono)
apply (simp add: mult_nonneg_nonneg)
apply simp
done

interpretation outer: bounded_bilinear ["outer"]
apply (unfold_locales)
apply (rule outer_left_distrib)
apply (rule outer_right_distrib)
apply (rule outer_left_scaleR)
apply (rule outer_right_scaleR)
apply (rule_tac x=1 in exI)
apply (simp add: norm_outer_ineq)
done

end