(*Example of Thm 2*)
n=5;
A=Table[RandomInteger[{-10,10}],{i,1,n},{j,1,n}];
A//MatrixForm
(*Characteristic polynomial*)
Det[x IdentityMatrix[n]-A]
Sum[(-1)^k Tr[Minors[A,k]]x^(n-k),{k,0,n}]

(*Examples of Thm 1 *)

(* Auxiliary function: Generates random orthogonal matrix with QR decomposition*)RandomOrthogonalMatrix[n_Integer]:=Module[{A,q,r,d},(*Generate a random matrix with normally distributed entries*)A=RandomVariate[NormalDistribution[],{n,n}];
(*QR decomposition*)
{q,r}=QRDecomposition[A];
(*Adjust signs to ensure uniform distribution (Haar measure)*)
d=Diagonal[r];
q=Transpose[q].DiagonalMatrix[Sign[d]]];

(*Example of diagonalizing a 'nice' symmetric matrix: integer eigenvalues*)
diag=DiagonalMatrix[Sort[Table[RandomInteger[{-10,10}],n]]];
randOrth=RandomOrthogonalMatrix[n];
A=randOrth.diag.Transpose[randOrth];
A//MatrixForm
(*Characteristic polynomial*)
Det[x IdentityMatrix[n]-A]
Sum[(-1)^k Tr[Minors[A,k]]x^(n-k),{k,0,n}]

(*Eigenvalues*)
eigenvalues=Sort[SolveValues[Det[x IdentityMatrix[n]-A]==0,x]]
eigenspaces=Table[Orthogonalize[NullSpace[eigenvalues[[i]] IdentityMatrix[n]-A]],{i,1,n}];
eigenspaces

(*Checking first and second pairs of eigenvalue/eigenvector:*)
A.Transpose[eigenspaces[[1,All]]]
eigenvalues[[1]] Transpose[eigenspaces[[1,All]]]

A.Transpose[eigenspaces[[2,All]]]
eigenvalues[[2]] Transpose[eigenspaces[[2,All]]]

(*Building the matrix P using eigenvectors found above, accounting for possible multiplicities*)
P={};
For[i=1,i<=n,i++,
If[i==1||eigenvalues[[i-1]]!=eigenvalues[[i]],
P=Join[P,eigenspaces[[i]]];
]];
(*Note P is already orthogonal because we already applied Gram-Schmidt to each eigenspace! *)
P=Transpose[P];
P//MatrixForm
OrthogonalMatrixQ[P]

(*Diagonalization produces a decomposition of the input matrix A as P D P^{-1}  *)
P.DiagonalMatrix[eigenvalues].Transpose[P]//MatrixForm
A//MatrixForm

(*Note that the random orthogonal matrix used and the orthogonal matrix P found above need not the same if some eigenvalue is not simple: *)
randOrth//MatrixForm
P//MatrixForm
OrthogonalMatrixQ[randOrth]
OrthogonalMatrixQ[P]

(*Example of diagonalizing a generic symmetric matrix*)

(*Generate a random n×n symmetric matrix*)RandomSymmetricMatrix[n_Integer,dist_:NormalDistribution[]]:=Module[{A},A=RandomVariate[dist,{n,n}];
(A+Transpose[A])/2];
A=RandomSymmetricMatrix[n];
A//MatrixForm
A==Transpose[A]

(*Characteristic polynomial*)
Det[x IdentityMatrix[n]-A]
Sum[(-1)^k Tr[Minors[A,k]]x^(n-k),{k,0,n}]

(*Eigenvalues*)
eigenvalues=Sort[SolveValues[Det[x IdentityMatrix[n]-A]==0,x]]
eigenspaces=Table[Orthogonalize[NullSpace[eigenvalues[[i]] IdentityMatrix[n]-A]],{i,1,n}];
eigenspaces

(*Building the matrix P using eigenvectors found above, accounting for possible multiplicities*)
P={};
For[i=1,i<=n,i++,
If[i==1||eigenvalues[[i-1]]!=eigenvalues[[i]],
P=Join[P,eigenspaces[[i]]];
]];
(*Note P is already orthogonal because we already applied Gram-Schmidt to each eigenspace! *)
P=Transpose[P];
P//MatrixForm
OrthogonalMatrixQ[P]

(*Diagonalization produces a decomposition of the input matrix A as P D P^{-1}  *)
P.DiagonalMatrix[eigenvalues].Transpose[P]//MatrixForm
A//MatrixForm

eigenvalues^2

(*Comparing efficiency in MatrixPower - depends on version of Mathematica and machine running it!*)
p=1000000;
Timing[MatrixPower[A,p]]
Timing[P.DiagonalMatrix[eigenvalues^p].Transpose[P]]

(*Generically, eigenvalues are simple*)
Clear[A,eigenvalues]
n=15;

(*Two random symmetric endpoints*)
RandomSymmetricMatrix[dim_]:=With[{A=RandomVariate[NormalDistribution[],{dim,dim}]},(A+Transpose[A])/2]

A0=RandomSymmetricMatrix[n];
A1=RandomSymmetricMatrix[n];

(*Matrix family:linear path through matrix space*)
A[t_]:=(1-t)*A0+t*A1

(*Compute sorted eigenvalues as a function of t*)
eigenvalues[t_]:=Sort[Eigenvalues[A[t]]]

(*Sample over t and plot*)
tValues=Range[0,1,0.005];
eigData=Table[eigenvalues[t],{t,tValues}];

(*Reshape for plotting: each eigenvalue as its own series*)
eigSeries=Table[Transpose[{tValues,eigData[[All,k]]}],{k,1,n}];

ListLinePlot[eigSeries,PlotRange->All,PlotStyle->Thick,PlotLegends->Table["λ"<>ToString[k],{k,1,n}],AxesLabel->{"t","Eigenvalue"},PlotLabel->"Eigenvalues of A(t) = (1-t)A0 + t A1",GridLines->Automatic,GridLinesStyle->LightGray,ImageSize->Large]

(*This phenomenon inspires the cover of Peter Lax's book "Linear Algebra and its applications" by Wiley.*)

(*Explanation: The space of real symmetric n×n matrices has dimension n(n+1)/2. Within this space,the subset of matrices with a repeated eigenvalue has codimension 2, so a 1-parameter family A(t)
traces a curve in this space, and curves generically miss subsets of codimension 2, so crossings do not occur.
This can be made rigorous using perturbation theory.
*)

(*Examples of Sylvester's criterion: done in class*)