Library mathcomp.algebra.spectral

From HB Require Import structures.
From mathcomp Require Import ssreflect ssrfun ssrbool eqtype choice ssrnat.
From mathcomp Require Import seq div fintype bigop ssralg finset fingroup zmodp.
From mathcomp Require Import poly polydiv order ssrnum matrix mxalgebra vector.
From mathcomp Require Import mxpoly mxred sesquilinear.

Spectral theory This file provides a formalization of Gram-Schmidt orthonormalization, Schur decomposition, etc. M ^t* := M ^t conjC M \is unitarymx == M is a unitary matrix M : 'M[C]_(m, n) with C : numClosedFieldType M \is normalmx == M is a normal matrix M : 'M[C]_n with C : numClosedFieldType dotmx u v == dot product u and v are row vectors over a numClosedFieldType Local notations: ' [u, v] := dotmx u v, ' [u] := ' [u, u] proj_ortho Y := proj_mx <<U>>%MS U^!%MS where U^! is a 1-orthogonal completement of U schmidt A == Gram-Schmidt basis A : 'M[C]_(m, n) schmidt_complete V := col_mx (schmidt (row_base V)) (schmidt (row_base V^!%MS)) spectralmx A, spectral_diag A == (M,X) s.t. A = M^-1 *m diag_mx X *m M A : 'M[C]_n

Set Implicit Arguments.

Import GRing.Theory Order.Theory Num.Theory.
Local Open Scope ring_scope.

TODO: move?
TODO: move?
Lemma common_eigenvector {C : numClosedFieldType} n (As : seq 'M[C]_n) :
  (n > 0)%N {in As &, A B, comm_mx A B}
  exists2 v : 'rV_n, v != 0 & all (fun Astablemx v A) As.

TODO: move?
Lemma common_eigenvector2 {C : numClosedFieldType}n (A B : 'M[C]_n) :
  (n > 0)%N A ×m B = B ×m A
  exists2 v : 'rV_n, v != 0 & (stablemx v A) && (stablemx v B).

Notation "M ^t*" := (M ^t conjC) (at level 30).
Notation realmx := (mxOver Num.real).

Lemma trmxCK {C : numClosedFieldType} m n (A : 'M[C]_(m, n)) : A ^t× ^t× = A.

Section realmx.
Context {C : numClosedFieldType} {m n : nat}.
Implicit Types A B : 'M[C]_(m, n).

Lemma realmxC A : A \is a realmx A ^ conjC = A.

Lemma realmxD A B : A \is a realmx B \is a realmx A + B \is a realmx.

Lemma Remx_rect : {in realmx &, A B, (A + 'i *: B) ^ (@Re _) = A}.

Lemma Immx_rect : {in realmx &, A B, (A + 'i *: B) ^ (@Im _) = B}.

Lemma eqmx_ReiIm A B A' B' :
  A \is a realmx B \is a realmx A' \is a realmx B' \is a realmx
  (A + 'i *: B) = (A' + 'i *: B') (A, B) = (A', B').

End realmx.

Lemma realsym_hermsym {C : numClosedFieldType} {n} (A : 'M[C]_n) :
  A \is symmetricmx A \is a realmx A \is hermsymmx.

Lemma real_similar {C : numClosedFieldType} {n} (A B : 'M[C]_n) :
  similar_in unitmx A B
  A \is a realmx B \is a realmx similar_in [predI realmx & unitmx] A B.

Section unitarymx.
Context {C : numClosedFieldType}.

Definition unitarymx {m n} := [qualify X : 'M[C]_(m, n) | X ×m X ^t× == 1%:M].
Fact unitarymx_key m n : pred_key (@unitarymx m n).
Canonical unitarymx_keyed m n := KeyedQualifier (unitarymx_key m n).

Lemma unitarymxP m n {M : 'M[C]_(m, n)} :
  reflect (M ×m M^t× = 1%:M) (M \is unitarymx).

Lemma mulmxtVK m1 m2 n (A : 'M[C]_(m1, n)) (B : 'M[C]_(n, m2)) :
  B \is unitarymx A ×m B ×m B^t× = A.

Lemma unitarymx_unit n (M : 'M[C]_n) : M \is unitarymx M \in unitmx.

Lemma invmx_unitary n (M : 'M[C]_n) : M \is unitarymx invmx M = M^t×.

Lemma mulmxKtV m1 m2 n (A : 'M[C]_(m1, n)) (B : 'M[C]_(m2, n)) :
  B \is unitarymx m2 = n A ×m B^t× ×m B = A.

Lemma mxrank_unitary m n (M : 'M[C]_(m, n)) : M \is unitarymx \rank M = m.

Lemma mul_unitarymx m n p (A : 'M[C]_(m, n)) (B : 'M[C]_(n, p)) :
  A \is unitarymx B \is unitarymx A ×m B \is unitarymx.

Lemma pinvmx_unitary n (M : 'M[C]_n) : M \is unitarymx pinvmx M = M^t×.

Lemma conjymx n (P M : 'M[C]_n) : P \is unitarymx conjmx P M = P ×m M ×m P^t×.

Lemma trmx_unitary n (M : 'M[C]_n) : (M ^T \is unitarymx) = (M \is unitarymx).

Lemma conjC_unitary m n (M : 'M[C]_(m, n)) :
  (M ^ conjC \is unitarymx) = (M \is unitarymx).

Lemma trmxC_unitary n (M : 'M[C]_n) : (M ^t× \is unitarymx) = (M \is unitarymx).

End unitarymx.

Section normalmx.
Context {C : numClosedFieldType} {n : nat}.

Definition normalmx := [qualify M : 'M[C]_n | M ×m M ^t× == M ^t× ×m M].
Fact normalmx_key : pred_key normalmx.
Canonical normalmx_keyed := KeyedQualifier normalmx_key.

Lemma normalmxP {M : 'M[C]_n} :
  reflect (M ×m M ^t× = M ^t× ×m M) (M \is normalmx).

Lemma hermitian_normalmx (A : 'M[C]_n) : A \is hermsymmx A \is normalmx.

Lemma symmetric_normalmx (A : 'M[C]_n) : A \is symmetricmx
  A \is a realmx A \is normalmx.

End normalmx.

Section Spectral.
Variable (C : numClosedFieldType).

Local Notation dotmx_def := (form_of_matrix (@conjC _) 1%:M).
Definition dotmx n (u v : 'rV[C]_n) := dotmx_def u%R v%R.

TODO: bug report we were expecting HB.instance Definition _ n := Bilinear.on (@dotmx n). to be sufficient to equip dotmx with the bilinear structure but needed to use .copy in the end as in:

Local Notation "''[' u , v ]" := (dotmx u v) : ring_scope.
Local Notation "''[' u ]" := '[u, u]%R : ring_scope.


Lemma dotmxE n (u v : 'rV[C]_n) : '[u, v] = ( u ×m v ^t× ) 0 0.

Lemma row_unitarymxP m n {M : 'M[C]_(m, n)} :
  reflect ( i j, '[row i M, row j M] = (i == j)%:R) (M \is unitarymx).

Fact dotmx_is_dotmx n (u : 'rV[C]_n) : u != 0 0 < '[u].


Local Notation "B ^!" :=
  (orthomx (@conjC C) (mx_of_hermitian (hermitian1mx _)) B) : matrix_set_scope.
Local Notation "A '_|_ B" := (A%MS B^!)%MS : bool_scope.

Lemma orthomx1E m n p (A : 'M[C]_(m, n)) (B : 'M_(p, n)) :
  (A '_|_ B)%MS = (A ×m B^t× == 0).

Lemma orthomx1P m n p {A : 'M[C]_(m, n)} {B : 'M_(p, n)} :
  reflect (A ×m B^t× = 0) (A '_|_ B).

Lemma orthomx_disj n p q (A : 'M[C]_(p, n)) (B :'M_(q, n)) :
  A '_|_ B (A :&: B = 0)%MS.

Lemma orthomx_ortho_disj n p (A : 'M[C]_(p, n)) : (A :&: A^! = 0)%MS.

Lemma rank_ortho p n (A : 'M[C]_(p, n)) : \rank A^! = (n - \rank A)%N.

Lemma add_rank_ortho p n (A : 'M[C]_(p, n)) : (\rank A + \rank A^!)%N = n.

Lemma addsmx_ortho p n (A : 'M[C]_(p, n)) : (A + A^! :=: 1%:M)%MS.

Lemma ortho_id p n (A : 'M[C]_(p, n)) : (A^!^! :=: A)%MS.

Lemma submx_ortho p m n (U : 'M[C]_(p, n)) (V : 'M_(m, n)) :
  (U^! V^!)%MS = (V U)%MS.

Definition proj_ortho p n (U : 'M[C]_(p, n)) := proj_mx <<U>>%MS U^!%MS.

Lemma sub_adds_genmx_ortho (p m n : nat) (U : 'M[C]_(p, n)) (W : 'M_(m, n)) :
  (W <<U>> + U^!)%MS.
Local Hint Resolve sub_adds_genmx_ortho : core.

Lemma cap_genmx_ortho p n (U : 'M[C]_(p, n)) : (<<U>> :&: U^!)%MS = 0.
Local Hint Resolve cap_genmx_ortho : core.

Lemma proj_ortho_sub p m n (U : 'M_(p, n)) (W : 'M_(m, n)) :
  (W ×m proj_ortho U U)%MS.

Lemma proj_ortho_compl_sub p m n (U : 'M_(p, n)) (W : 'M_(m, n)) :
  (W - W ×m proj_ortho U U^!)%MS.

Lemma proj_ortho_id p m n (U : 'M_(p, n)) (W : 'M_(m, n)) :
  (W U)%MS W ×m proj_ortho U = W.

Lemma proj_ortho_0 p m n (U : 'M_(p, n)) (W : 'M_(m, n)) :
  (W U^!)%MS W ×m proj_ortho U = 0.

Lemma add_proj_ortho p m n (U : 'M_(p, n)) (W : 'M_(m, n)) :
  W ×m proj_ortho U + W ×m proj_ortho U^!%MS = W.

Lemma proj_ortho_proj m n (U : 'M_(m, n)) : let P := proj_ortho U in P ×m P = P.

Lemma proj_orthoE p n (U : 'M_(p, n)) : (proj_ortho U :=: U)%MS.

Lemma orthomx_proj_mx_ortho p p' m m' n
  (A : 'M_(p, n)) (A' : 'M_(p', n))
  (W : 'M_(m, n)) (W' : 'M_(m', n)) :
  A '_|_ A' W ×m proj_ortho A '_|_ W' ×m proj_ortho A'.

Lemma schmidt_subproof m n (A : 'M[C]_(m, n)) : (m n)%N
  exists2 B : 'M_(m, n), B \is unitarymx & [ i : 'I_m,
   (row i A (\sum_(k < m | (k i)%N) <<row k B>>))%MS
   && ('[row i A, row i B] 0) ].

Definition schmidt m n (A : 'M[C]_(m, n)) :=
  if (m n)%N =P true is ReflectT le_mn
  then projT1 (sig2_eqW (schmidt_subproof A le_mn))
  else A.

Lemma schmidt_unitarymx m n (A : 'M[C]_(m, n)) : (m n)%N
  schmidt A \is unitarymx.
Hint Resolve schmidt_unitarymx : core.

Lemma row_schmidt_sub m n (A : 'M[C]_(m, n)) i :
  (row i A (\sum_(k < m | (k i)%N) <<row k (schmidt A)>>))%MS.

Lemma form1_row_schmidt m n (A : 'M[C]_(m, n)) i :
  '[row i A, row i (schmidt A)] 0.

Lemma schmidt_sub m n (A : 'M[C]_(m, n)) : (A schmidt A)%MS.
Hint Resolve schmidt_sub : core.

Lemma eqmx_schmidt_full m n (A : 'M[C]_(m, n)) :
  row_full A (schmidt A :=: A)%MS.

Lemma eqmx_schmidt_free m n (A : 'M[C]_(m, n)) :
  row_free A (schmidt A :=: A)%MS.

Definition schmidt_complete m n (V : 'M[C]_(m, n)) :=
  col_mx (schmidt (row_base V)) (schmidt (row_base V^!%MS)).

Lemma schmidt_complete_unitarymx m n (V : 'M[C]_(m, n)) :
  schmidt_complete V \is unitarymx.

Lemma cotrigonalization n (As : seq 'M[C]_n) :
  {in As &, A B, comm_mx A B}
  cotrigonalizable_in (@unitarymx C n n) As.

Theorem Schur n (A : 'M[C]_n) : (n > 0)%N
  trigonalizable_in (@unitarymx C n n) A.

Lemma cotrigonalization2 n (A B : 'M[C]_n) : A ×m B = B ×m A
  exists2 P : 'M[C]_n, P \is unitarymx &
    similar_trig P A && similar_trig P B.

Theorem orthomx_spectral_subproof n {A : 'M[C]_n} : reflect
  (exists2 sp : 'M_n × 'rV_n,
                sp.1 \is unitarymx &
                A = invmx sp.1 ×m diag_mx sp.2 ×m sp.1)
  (A \is normalmx).

Definition spectralmx n (A : 'M[C]_n) : 'M[C]_n :=
  if @orthomx_spectral_subproof _ A is ReflectT P
  then (projT1 (sig2_eqW P)).1 else 1%:M.

Definition spectral_diag n (A : 'M[C]_n) : 'rV_n :=
  if @orthomx_spectral_subproof _ A is ReflectT P
  then (projT1 (sig2_eqW P)).2 else 0.

Lemma spectral_unitarymx n (A : 'M[C]_n) : spectralmx A \is unitarymx.

Lemma spectral_unit n (A : 'M[C]_n) : spectralmx A \in unitmx.

Theorem orthomx_spectralP {n} {A : 'M[C]_n}
  (P := spectralmx A) (sp := spectral_diag A) :
  reflect (A = invmx P ×m diag_mx sp ×m P) (A \is normalmx).

Lemma hermitian_spectral_diag_real n (A : 'M[C]_n) : A \is hermsymmx
  spectral_diag A \is a realmx.

End Spectral.