(* May 28, 2013,   Mathematica  Cauchy-Binet  and Pythagoras experiments                             *)
(* Paper: Cauchy_Binet for Pseudo Determinants     http://arxiv.org/abs/1306.0062                    *)
(* June 10, This code has now been added to the ArXiv paper                                          *)

(* Pythagoras:   Det^2(A) = sum_P det^2(A_P)      if A is selfadjoint of rank k and P are kxk masks  *)
(* Cauchy-Binet:   p_k( F^T G)  = (-1)^k sum_P det(F_P) det(G_P)  for nxm matrices F,G               *)

(* The Pseudo determinant is up to a sign the first nonzero element of the Characteristic polynomial *)
(* except if p(x) = x^n, the nilpotent situation in which case PDet is equal to 0                    *)

f[s_,n_]:=If[n==Length[s]-1,0,s[[1+n]](-1)^n];
TopNoZero[s_]:=f[s,Last[Position[s,0]][[1]]];
PDet[A_]:=Module[{s=CoefficientList[CharacteristicPolynomial[A,x],x],n},
   n=Sum[Abs[Sign[s[[k]]]],{k,Length[s]}]; d=Det[A];
   If[n==1 || d !=0,res=d,res=TopNoZero[s]];res];

PTrace[A_,B_,k_]:=Module[{U,V,m=Length[A],n=Length[A[[1]]],s={},t={}},
 U=Partition[Flatten[Subsets[Range[m],{k,k}]],k]; u=Length[U];
 V=Partition[Flatten[Subsets[Range[n],{k,k}]],k]; v=Length[V];
 Do[s=Append[s,Table[A[[U[[i,u]],V[[j,v]]]],{u,k},{v,k}]],{i,u},{j,v}];
 Do[t=Append[t,Table[B[[U[[i,u]],V[[j,v]]]],{u,k},{v,k}]],{i,u},{j,v}];
 Sum[Det[s[[l]]]*Det[t[[l]]],{l,Length[s]}]];
PTraces[A_,B_]:=(-x)^m+Sum[(-x)^(m-k) PTrace[A,B,k],{k,Min[n,m]}];
n=5; m=7; r=4; RandomMatrix[n_,m_]:=Table[RandomInteger[2r]-r,{n},{m}];
F=RandomMatrix[n,m]; G=RandomMatrix[n,m];
CharacteristicPolynomial[Transpose[F].G,x]
PTraces[F,G]

(* Cauchy-Binet experiments, We make sure to get examples reaching new territory                    *)

R:=RandomMatrix[n,m]; 
Shuffle:=Module[{},F=R;G=R;A=Transpose[F].G;B=F.Transpose[G]];
Shuffle; While[Abs[Det[A]]+Abs[Det[B]]>0,Shuffle]; 
CharacteristicPolynomial[Transpose[F].G,x]
PTraces[F,G]

(* Pythagoras experiments, again we shuffle until we get examples outside of classical Cauchy-Binet *)

n=4; m=5; R:=RandomMatrix[n,m]; 
Shuffle:=Module[{},F=R;A=Transpose[F].F;B=F.Transpose[F]];
Shuffle; While[Abs[Det[A]]+Abs[Det[B]]>0,Shuffle]; 
{k,l}=Map[MatrixRank,{A,B}];{n,m,k,l}
PDet[A]
PTrace[F,F,k]


(* Some extreme test cases for Pseudo determinant  *)


A={{1,2},{3,1}}; {PDet[A],Det[A]}
A={{0,1},{0,0}}; {Eigenvalues[A],PDet[A],Det[A]}
A={{0}}; {Eigenvalues[A],PDet[A],Det[A]}
A={{0,0},{0,0}}; {Eigenvalues[A],PDet[A],Det[A]}
A={{0,0,0},{0,0,1},{0,0,0}}; {Eigenvalues[A],PDet[A],Det[A]}
A={{1,1,1},{1,1,1},{1,1,1}}; {Eigenvalues[A],PDet[A],Det[A]}

