SOR and GS Oh my

m=4;
A= RandomReal[{-1,1},{m,m}]+DiagonalMatrix[RandomReal[{m/2,m},{m}]];
A=A+Transpose[A];
b = RandomReal[{-1,1},m];
x=LinearSolve[A,b]
A.x-b

Clear[JacobiStep1]
JacobiStep1[ xIn_,b_]:=Module[
{i,m=Length[b]},
Table[(1/(A[[i,i]]))(b[[i]]+A[[i,i]] xIn[[i]]-Sum[A[[i,j]] xIn[[j]],{j,1,m}]),
{i,1,m}]
]
Clear[JacobiStep2]
SetAttributes[JacobiStep2,HoldFirst]
JacobiStep2[xOut_, xIn_,b_]:=Module[
{i,m=Length[b]},
Do[
xOut[[i]]=(1/(A[[i,i]]))(b[[i]]+A[[i,i]] xIn[[i]]-Sum[A[[i,j]] xIn[[j]],{j,1,m}]),
{i,1,m}]
]
SetAttributes[GSStep,HoldRest]
GSStep[b_,xOut_, xIn_]:=Module[
{i,m=Length[b]},
Do[
xOut[[i]]=(1/(A[[i,i]]))(b[[i]]+A[[i,i]] xIn[[i]]-Sum[A[[i,j]] xIn[[j]],{j,1,m}]),
{i,1,m}]
]
SetAttributes[SORStep,HoldFirst]
SORStep[x_,{A_,b_},ω_]:=Module[
{i,m=Length[b]},
Do[
x[[i]]=(1-ω)x[[i]] + (ω/(A[[i,i]]))(b[[i]]+A[[i,i]] x[[i]]-Sum[A[[i,j]] x[[j]],{j,1,m}]),
{i,1,m}]
]

MaxIter=56;
xOut = 0.0*x;
JData=Table[
xOut=JacobiStep1[xOut,b];
Norm[xOut-x],
{MaxIter}];
xOut = 0.0*x;
AltJData=Table[
JacobiStep2[xOut,xOut,b];
Norm[xOut-x],
{MaxIter}];
xOut = 0.0*x;
GSData=Table[
GSStep[b,xOut,xOut];
Norm[xOut-x],
{MaxIter}];
xOut = 0.0*x;
SORData=Table[
SORStep[xOut,{A,b},1.2];
Norm[xOut-x],
{MaxIter}];
ListLogPlot[{JData,AltJData, GSData,SORData},
PlotStyle->{{PointSize[0.02],Red},Blue, Green, Cyan}]

As it should GS converges twice as fast as Jacobi.   The SOR cnverges but I am not beating GS with random choices of ω.

The optimal SOR relaxation parameter for the Poisson equation in 2D is 
	ω=(2/(1+Sin[π h]))
where h=Δx=Δy.   So for our internal solver h=(κ/2) ((Δt)/((Δx)^2))=(ν/2) with the half comming from the Crank Nickolson symmetry trick.  
Yang and Gobbert Applied Mathematics Letters, 2007, "The Optimal Relaxation Parameter for the SOR Method Applied to
the Poisson Equation in Any Space Dimensions" provide the analogous formula

The general convergnce theory gives 
	(Subscript[ω, opt])=(2/(1+(Sqrt[1-(ρ^2)])))
for the SOR scheme wher ρ is the spectral radius of the Jacobi iteration matrix (Subscript[G, j]). with A=D-L-U  we have
(Subscript[G, j])=I-(D^(-1))A 

See James Demmel, Applied Numerical Linear Algebra, SIAM, 1997.

Clear[AFD, ν,m]
AFD[ν_,m_]:=SparseArray[{
Band[{1,1},{m*m,m*m}]->1-4ν,
Band[{2,1},{m*m,m*m-1}]->ν,
Band[{1,2},{m*m,m*m-1}]->ν,
Band[{1,m+1},{m*(m-1),m*m}]->ν,
Band[{m+1,1},{m*m,m*(m-1)}]->ν
}]

m=4; ν=0.1;
A=AFD[ν,m];
b= RandomReal[{-1,1},{m*m}];
x=LinearSolve[A,b];

ωOpt=(2/(1+Sin[π ν] ))

Clear[SORData,ω,ν,m]
SORData[ω_,ν_,m_]:=Module[
{
xOut=ConstantArray[0.0,{m*m}],
ωOpt=(2/(1+Sin[π ν])),
A=SparseArray[{
Band[{1,1},{m*m,m*m}]->1-4ν,
Band[{2,1},{m*m,m*m-1}]->ν,
Band[{1,2},{m*m,m*m-1}]->ν,
Band[{1,m+1},{m*(m-1),m*m}]->ν,
Band[{m+1,1},{m*m,m*(m-1)}]->ν
},
{m*m,m*m}],
b=ConstantArray[1.0,{m*m}],
x=LinearSolve[A,b]
},
Print[x];
Table[
SORStep[xOut,{A,b},1.2];
Norm[xOut-x];
xOut,
{7}] 
]

ListPlot[Transpose[
data=SORData[1.2,0.1, 4]]]

data[[-1]]

2/(1 + Sin[π 0.03])

Break

m = 230;
A= RandomReal[{-1,1},{m,m}];
{Q,R}=QRDecomposition[A];
Dia=DiagonalMatrix[ RandomReal[{0.1, 12},{m}]];
A= Q.Dia.Transpose[Q];
ListPlot[Eigenvalues[A]]

Dia

A= SparseArray[ {
Band[{1,1}]:>RandomReal[{-1,1}],
Band[{1,2}]:>2.2,
Band[{2,1}]:>1.2
},{m,m}]

2.0^(-23)

2 Pi/100.

Jen

Laregst entry in matrix.

m=12;
A= RandomReal[{-1,1}, {m,m}];
A=A+Transpose[A]
TestMat=LowerTriangularize[Abs[A],-1]
MaxVal=Max[TestMat]
Position[TestMat,MaxVal]

Largest entry in a specified row!

m=12;
i=6;
A= RandomReal[{-1,1}, {m,m}]
OurRow=Abs[A[[i]]]
MaxVal=Max[OurRow]
Position[OurRow,MaxVal]

Largest k entries in a specified row first.

m=12;
i=6;
k=3;
A= RandomReal[{-1,1}, {m,m}]
OurRow=Abs[A[[i]]]
OurBigs=Ordering[OurRow,-k]
OurRow[[OurBigs]]

Largest k entries in a matrix

m=12;
A= RandomReal[{-1,1}, {m,m}];
A=A+Transpose[A]
TestMat=LowerTriangularize[Abs[A],-1];
k=4;
Thing=Ordering[TestMat,k, Max[#1]>Max[#2]&]

Clear[FindMaxInVec]
FindMaxInVec[x_]:= Position[x,Max[x]][[1,1]]
FindMaxInVec[{1.1614968234777807`,0.07536543558863018`,0.2951807095580832`,0.821704015037203`,1.5770034045514514`,0.2752890292528325`,1.1597790647062998`,0.24853417990316107`,0.`,0.`,0.`,0.`}]

Map[FindMaxInVec,TestMat[[Thing]]]

Ordering[TestMat,k, Max[#1]>Max[#2]&]

?Ordering

TestMat[[2]]
TestMat[[6]]
TestMat[[5]]

Crud

SetOptions[
EvaluationNotebook[],DockedCells->{
Cell[
BoxData[ToBoxes[
Button["wombat",FrontEnd`FrontEndExecute[FrontEnd`FrontEndToken[ButtonNotebook[],"HyperlinkGoBack"]],Evaluator->None]]],
"DockedCell"],
Cell[
BoxData[ToBoxes[
Button["back",FrontEnd`FrontEndExecute[FrontEnd`FrontEndToken[ButtonNotebook[],"HyperlinkGoBack"]],Evaluator->None]]],
"DockedCell"]}]

?Placeholder

MatrixForm[

?Slider

{Slider[Dynamic[x],{-23,1}],Dynamic[x]}

?Switch

?RandomReal

?UniformDistribution

?RandomInteger

?UniformDistribution

?Band

A= RandomReal[{-1,1},{4,4}]
A-DiagonalMatrix[Diagonal
Diagonal[A]={1,2,3,4}

Manipulate[
{μ,s}={0,1};
Switch[DiagonalDistribution,
"Uniform",A=RandomReal[UniformDistribution[{μ-s,μ+s}],{m,m}],
"Normal",  A=RandomReal[NormalDistribution[μ,s],{m,m}],
"Integer",A= RandomInteger[Evaluate[{μ-s,μ+s}],{m,m}],
"Blank",    A= ConstantArray[Placeholder[],{m,m}]
];
Switch[OffDiagonalDistribution,
"Uniform",A=RandomReal[UniformDistribution[{μ-s,μ+s}],{m,m}],
"Normal",  A=RandomReal[NormalDistribution[μ,s],{m,m}],
"Integer",A= RandomInteger[Evaluate[{μ-s,μ+s}],{m,m}],
"Blank",    A= ConstantArray[Placeholder[],{m,m}]
];
Switch[Type,
"General", B=A,
"Symmetric",  B=A + Transpose[A],
"Symmetric+DD",B= A+Transpose[A] + 234 IdentityMatrix[m],
"Skew Symmetric",B= A-Transpose[A]
];
MatrixForm[B],
{{m,6}, 3, 9,1},
{DiagonalDistribution,{"Uniform", "Normal","Integer"}},
{OffDiagonalDistribution,{"Uniform", "Normal","Integer"}},
{Type,{"General","Symmetric", "Symmetric+DD","Skew Symmetric"}}
]

1+1

Twiddlers

I want to try to make a twiddler that shows linear algebra.
It needs to have an inititialization process and a modification process.

Things that might be usefulll "DockedCells"  appear at the top of a notebook.  
Dumping active code downwards is probably good.   I need a dragger and dropper for values into matrices
Scrolling downwards

m=4
A=Table[Placeholder[],{m},{m}];
MatrixForm[A]

(({{1, 1, \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}, {\[Placeholder], \[Placeholder], \[Placeholder], \[Placeholder]}}))

m=4;
A = RandomReal[{-1,1},{m,m}];
L=IdentityMatrix[m];
L[[1,2]]=-A[[1,2]]/A[[1,1]];
PaletteNotebook[StringForm["(Underscript[(Underscript[``, ⏟]), A]).(Underscript[(Underscript[``, ⏟]), L]) ⟹ (Underscript[(Underscript[``, ⏟]), A . L])",MatrixForm[A],MatrixForm[L],MatrixForm[A.L]]]

CreateWindow[DockedCells->Cell[StringForm["(Underscript[(Underscript[``, ⏟]), A]).(Underscript[(Underscript[``, ⏟]), L]) ⟹ (Underscript[(Underscript[``, ⏟]), A . L])",MatrixForm[A],MatrixForm[L],MatrixForm[A.L]]]]

StringForm["(Underscript[(Underscript[``, ⏟]), A]).(Underscript[(Underscript[``, ⏟]), L]) ⟹ (Underscript[(Underscript[``, ⏟]), A . L])",MatrixForm[A],MatrixForm[L],MatrixForm[A.L]]

CreateWindow[DockedCells->Cell["XXXX","DockedCell"]];
CreateWindow[DockedCells->Cell[A,"DockedCell"]];

Code Complexity

Students find the code very hard and once it gets too complicated for them they shut down.
Students do not get swamped with the Sparse Array structure.  They like band.

OP Counts

OP Counts and Timing are OK.

Linear Algebra

Students do not understand what linear algebra is for.
Students do not understand eigen/SVD stuff
Students do not have any geometric understanding

?Arrow

Graphics3D[

A= {
{-1,0},
{0,1}
};
Graphics[{
Arrowheads[{0.01,0.1}],
Hmm=Table[
Arrow[{{x1,x2},A.{x1,x2}}],
{x1,-1,1,0.5},{x2,-1,1,0.5}]
},
Frame->True,
AspectRatio->Automatic
]

A.A

Hmm

Graphics[Arrow[{{0,0},{1,2}}]]

A.{x1,x2}

?Arrow

Manipulate[
A.{x1,x2},

Visualizing Projections and Rotations

Visualizing stuff is hard and 2D does not make enough distinction between the two spaces for projections.

Clear[ProjectionBuilder,nIn]
ProjectionBuilder[nIn_List]:=Module[
{Q,i},
Q=IdentityMatrix[3];
Q[[1]]=N[nIn];
Q[[2]]=Cross[ nIn,{1.0,0,0}];
Q[[3]]=Cross[nIn,Q[[2]]];
Do[
Q[[i]] = Q[[i]]/Norm[Q[[i]]],
{i,3}];
Q
]

nIn={1,1,1};
{n,t1,t2}= ProjectionBuilder[nIn]

House

v={0.3,-0.3,0.7};
n={1.2,3.4,5.6};
(* n=v+Norm[v] {1,0,0}; *)
n=n/Norm[n];

P = IdentityMatrix[3]-Transpose[{n}].{n};
Q=  IdentityMatrix[3]-2Transpose[{n}].{n};
Show[
ContourPlot3D[ {
(x^2)+(y^2)+(z^2)==Norm[v],
z==0,
n.{x,y,z}==0
},
{x,-1,1},
{y,-1,1},
{z,-1,1},
ContourStyle->{{Pink,Opacity[0.2]},LightGreen, LightBlue},
FaceGrids->None,
Mesh->None
],
Graphics3D[{PointSize[0.02],Thickness[0.01],
{Red,Point[v], Point[P.v],Point[Q.v]},
{Red,Arrow[{{0,0,0},v}]},
{Blue,Arrow[{{0,0,0},P.v}]},
{Green,Arrow[{v, P.v}],Arrow[{P.v,Q.v}]}
}]
]

Cholesky

The Cholesky decomposition can be implemented in a number of different ways.

Outer Product Cholesky

OuterProductCholesky[A_]:= Module[
{m=Length[A],L=IdentityMatrix[m],B=A,a,b},
Do[
a=(Sqrt[B[[i,i]]]);
b=B[[i,i+1;;m]]/a;
B[[i+1;;m,i+1;;m]] -=Outer[Times,b,b];
L[[i,i]]=a;
L[[i+1;;m,i]]=b,
{i,1,m}];
L
]

m = 6;
A = RandomReal[{-1,1},{m,m}] + m IdentityMatrix[m];
A=A+Transpose[A];
L=OuterProductCholesky[A];
MatrixForm[L]
Norm[A-L.Transpose[L]]

Loopy Cholesky

LoopyCholesky[A_]:= Module[
{m=Length[A],L=ConstantArray[0,{m,m}],i,j},
Do[
Do[
L[[i,j]]=(1/(L[[j,j]]))(A[[i,j]]-Sum[L[[i,k]]L[[j,k]],{k,1,i-1}]),
{j,1,i-1}];
L[[i,i]]=(Sqrt[A[[i,i]]-Sum[((L[[i,k]])^2),{k,1,i-1}]]),
{i,1,m}];
L
]

m = 12;
A = RandomReal[{-1,1},{m,m}] + m IdentityMatrix[m];
A=A+Transpose[A];
L=LoopyCholesky[A];
Norm[A-L.Transpose[L]]

Time Trials

m = 12;
A = RandomReal[{-1,1},{m,m}] + m IdentityMatrix[m];
A=A+Transpose[A];
Timing[LoopyL=LoopyCholesky[A];]

Timing[OuterProductL=OuterProductCholesky[A];]
Norm[LoopyL-OuterProductL]

Built-In Commands

Built in commands may do more than you think.  They may also store stuff in compressed formats.

?LUDecomposition

A={
{1.0,2.0,3.0},
{3.,4.0,5.0},
{45.0,7.0,9.0}
};
{lu,p,c}=LUDecomposition[A]

m=3;
L=LowerTriangularize[lu,-1]+IdentityMatrix[m];
U=UpperTriangularize[lu];
L.U-A

?QRDecomposition

{Q,R}=QRDecomposition[A];
Norm[A-Transpose[Q].R]
Norm[Q.A-R]

Lots of other libraries do other things.  You always need to read the documentation!

gsl library web page
14.2 QR Decomposition

A general rectangular M-by-N matrix A has a QR decomposition into the product of an orthogonal M-by-M square matrix Q (where Q^T Q = I) and an M-by-N right-triangular matrix R,

     A = Q R

This decomposition can be used to convert the linear system A x = b into the triangular system R x = Q^T b, which can be solved by back-substitution. Another use of the QR decomposition is to compute an orthonormal basis for a set of vectors. The first N columns of Q form an orthonormal basis for the range of A, ran(A), when A has full column rank.
— Function: int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)

    This function factorizes the M-by-N matrix A into the QR decomposition A = Q R. On output the diagonal and upper triangular part of the input matrix contain the matrix R. The vector tau and the columns of the lower triangular part of the matrix A contain the Householder coefficients and Householder vectors which encode the orthogonal matrix Q. The vector tau must be of length k=\min(M,N). The matrix Q is related to these components by, Q = Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i v_i^T and v_i is the Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same storage scheme as used by lapack.

    The algorithm used to perform the decomposition is Householder QR (Golub & Van Loan, Matrix Computations, Algorithm 5.2.1).

?CholeskyDecomposition

New Jacobi

J[{i_,j_},{p1_,p2_},m_]:= Module[
{n=Norm[{p1,p2}]},
(* (({{ p1, p2}, {-p2, P1}}))*)
 SparseArray[{
{i,i}->p1,
{j,j}->p1,
{i,j}->p2,
{j,i}->-1p2,
Band[{1,1}]->1
},
{m,m}
]
]
J[{i_,j_,k_,l_},{p1_,p2_,p3_,p4_},m_]:= Module[
{n=Norm[{p1,p2,p3,p4}],pp1,pp2,pp3,pp4},
(* LinAlg+App 332(2001) 37-80 Fasbender p65 *)
(* (({{ p1,  p2,  p3,  p4}, {-p2,  p1,  p4, -p3}, {-p3, -p4,  p1,  p2}, {-p4,  p3, -p2,  p1}})) *)
{pp1,pp2,pp3,pp4}={p1,p2,p3,p4}/n;
 SparseArray[{
{i,i}->pp1,  {i,j}->pp2,{i,k}->pp3, {i,l}->pp4,
{j,i}->-pp2,{j,j}->pp1,{j,k}->pp4,{j,l}->-pp3,
{k,i}->-pp3,{k,j}->-pp4,{k,k}->pp1,{k,l}->pp2,
{l,i}->-pp4,{l,j}->pp3,{l,k}->-pp2,{l,l}->pp1,
Band[{1,1}]->1
},
{m,m}
]
]

J[
{i1_,i2_,i3_,i4_,i5_,i6_,i7_,i8_},
{p1_,p2_,p3_,p4_,p5_,p6_,p7_,p8_},m_]:= Module[
{n=Norm[{p1,p2,p3,p4,p5,p6,p7,p8}],
pp1,pp2,pp3,pp4,pp5,pp6,pp7,pp8},
(* LinAlg+App 332(2001) 37-80 Fasbender p66 *)
(* (({{ p1,  p2,  p3,  p4,  p5,  p6,  p7,  p8}, {-p2,  p1,  p4, -p3,  p6, -p5, -p8,  p7}, {-p3, -p4,  p1,  p2,  p7,  p8, -p5, -p6}, {-p4,  p3, -p2,  p1,  p8, -p7,  p6, -p5}, {-p5, -p6, -p7, -p8,  p1,  p2,  p3,  p4}, {-p6,  p5, -p8,  p7, -p2,  p1,  p4, -p3}, {-p7,  p8,  p5, -p6, -p3, -p4,  p1,  p2}, {-p8, -p7,  p6, -p5, -p4,  p3, -p2,  p1}})) *)
{pp1,pp2,pp3,pp4}={p1,p2,p3,p4}/n;
 SparseArray[{
{i,i}->pp1,  {i,j}->pp2,{i,k}->pp3, {i,l}->pp4,
{j,i}->-pp2,{j,j}->pp1,{j,k}->pp4,{j,l}->-pp3,
{k,i}->-pp3,{k,j}->-pp4,{k,k}->pp1,{k,l}->pp2,
{l,i}->-pp4,{l,j}->pp3,{l,k}->-pp2,{l,l}->pp1,
Band[{1,1}]->1
},
{m,m}
]
]

m=12;
A=RandomReal[{-1,1},{m,m}];
A=A+Transpose[A];
{i,j,k,l}={1,2,3,9};
p=A[[{i,j,k,l},i]];
JJ=J[{i,j,k,l}, p,m];
MatrixForm[JJ]
MatrixForm[A]
MatrixForm[JJ.A]

MatrixForm[A]
MatrixForm[Chop[JJ.A.Transpose[JJ]]]

Optimization Jacobi

What I have in the books etc are definitions of convenient and easily computed rotations.  The default versions zero entries in a column in the forward mode.
	A⟶J.A
The ArcTan version for the 2×2 does simultaneous zeroing in a symetric matrix for two entries in the similarity case
	A⟶J.A.(J^T)
No matter what I plug in J is always a rotation.

The questiona are:
	1) Are there simultaneous zeroing condition for the 4×4 an 8×8 cases?
	2) Is it better to make the off diagonal small in the block rather than make some elenments zero?
	
This starts to sound like an optimization problem to maximize the diagonal entries in the block.    The deal here is that we are hoping to have a simple analytic answer to the optimization problem

A= Array[a,{2,2}](*/.a[2,1]->a[1,2];*) 
JJ = J[{1,2},{p1,p2},2];
MatrixForm[JJ]
AA=JJ.A.Transpose[JJ]
F[{p1_,p2_}]=Simplify[((Diagonal[AA])^2)]/.List->Plus
Minimize[{
F[{p1,p2}],
(p1^2)+(p2^2)==1
},
{p1,p2}]

Simplify[
Simplify[
AA[[1,2]]/.{p2->ϵ p1}
]/(p1^2)]==0
Simplify[
Simplify[
AA[[2,1]]/.{p2->ϵ p1}
]/(p1^2)]==0

Solve[
Simplify[
Simplify[
AA[[1,2]]/.{p2->ϵ p1}
]/(p1^2)]==0, ϵ]//Simplify

?Solve

Clear[p1,p2]
Solve[
{
AA[[1,2]]==0,
(p1^2)+(p2^2)==1
},
p2]

Solve[ AA[[1,2]]==0,p2]

Minimize[{
F[{p1,p2}],
(p1^2)+(p2^2)==1
},
{p1,p2}]

D

Clear[c,s]
JJ=J[{1,2},{p1,p2},2]
A={
{p1,p2},
{p2,p1}
}
AA=JJ.A
Simplify[((Norm[Expand[Diagonal[AA]]/.{(s^2)->1-(c^2),s->(Sqrt[1-(c^2)])}])^2).
Assumptions->{0≤c≤1}
]

MatrixForm[Simplify[AA/.a[2,1]->a[1,2]]]

Expand[AA[[1,2]]]/(c s) == 0

Solve[ 
AA[[1,2]] == 0,
c]

QR with Shift

The QRShift step is simple.  Here it is:

QRShiftStep[A_,μ_] := Module[
{Q,R, Id = IdentityMatrix[Length[A]]},
{Q,R}=QRDecomposition[A-μ Id];
(* Note MMa's QR decomp gives the transpose of the Q described in the text. *)
R.Transpose[Q] + μ Id
]

We need to test it.

m=4;
A=RandomReal[ {-1,1},{m,m}];
A=A + Transpose[A];
Eigenvalues[A]

The next question is how to choose the shifts.  Two different questioms: 
1) If we knew good approximations to the eigenvalues we would use one of these a shift! 
As  pointed out on Friday this looks useless untill we realize that the diagonal entries (particularly the last one) are good approximations to eigenvalues.

B=QRShiftStep[A,3.43967];
BB=QRShiftStep[A,0.21838];
BBB=QRShiftStep[A,0.391285];
MatrixForm[Map[Eigenvalues,{A,B,BB,BBB}]]
Map[MatrixForm,{B,BB,BBB}]

QRIterationNC[A_, NIter_] := Module[
{Q, R, AA=A,i=0},
Do[
Print[ 
"Iteration=",i++,
", Lower Triangular F Norm=",Norm[LowerTriangularize[AA,-1],"Frobenius"]];
Print[Diagonal[AA]];
(* Implementing a shift simply by the trailing entry *)
AA=QRShiftStep[AA,AA[[-1,-1]]],
{NIter}];
AA
]

Γ=QRIterationNC[A,18];
MatrixForm[Γ]

Other shifting startegies are possible.

2×2 Jacobi Iteration: Exact Zeros

As we discussed before it is possible to use similarity rotations to introduce pairs of zeros into a matrix

J2[{i_,j_},θ_,m_] := SparseArray[{
{i,i}->Cos[θ],
{j,j}->Cos[θ],
{i,j}->-Sin[θ],
{j,i}->Sin[θ],
Band[{1,1}]->1.0
},
{m,m}]

m= 5;
A=RandomReal[{-1,1},{m,m}];
A= A + Transpose[A];
{i,j}={2,4};
θ= 0.5 ArcTan[ 2 ((A[[i,j]])/(A[[j,j]]-A[[i,i]]))]
JTemp=J2[{i,j},θ,m];
MatrixForm[JTemp.A.Transpose[JTemp]]

A sweep is a pass throught the matrix reducing all the off diagonals.  Lets make a decent sized test matrix and see how we could organizze a sweep.

m= 100;
A=RandomReal[{-1,1},{m,m}];
A= A + Transpose[A];

Repeated sweeps : fixed order.

Here are repeated sweeps through the matrix in a fixed order.

AA=A;
NSweep=12;
Data=Table[
Do[
θ= 0.5 ArcTan[ 2 ((AA[[i,j]])/(AA[[j,j]]-AA[[i,i]]))];
JTemp=J2[{i,j},θ,m];
AA=JTemp.AA.Transpose[JTemp],
{i,1,m,1},
{j,i+1,m}
];
{n*(m-1)*m/2,((Norm[LowerTriangularize[AA,-1],"Frobenius"])/(Norm[AA,"Frobenius"]))},
{n,NSweep}];
ListLogPlot[Data,
FrameLabel->{"J Comps","Frobenius OffDiagona Fraction"},
GridLines->Automatic,
PlotStyle->{PointSize[0.03],Blue},
Frame->True,
PlotLabel->"Jacobi Fixed Order"]
Sort[Diagonal[AA],Abs[#1]>Abs[#2]&]
Eigenvalues[A]

Repeated row comps : Partial Max

Here are repeated passes through the matrix reducing the maximum in each column

AA=A;
NPartialSweep=120;
Data=Table[
Do[
AbsCol =Abs[AA[[All,i]]];
AbsCol[[i]]=0.0;
j=Position[AbsCol,Max[AbsCol],1][[1,1]];
θ= 0.5 ArcTan[ 2 ((AA[[i,j]])/(AA[[j,j]]-AA[[i,i]]))];
JTemp=J2[{i,j},θ,m];
AA=JTemp.AA.Transpose[JTemp],
{i,1,m,1}
];
{n*m,((Norm[LowerTriangularize[AA,-1],"Frobenius"])/(Norm[AA,"Frobenius"]))},
{n,NPartialSweep}];
ListLogPlot[Data,
FrameLabel->{"J Comps","Frobenius OffDiagonal Fraction"},
GridLines->Automatic,
Frame->True,
PlotLabel->"Jacobi Col Max Search"]
Sort[Diagonal[AA],Abs[#1]>Abs[#2]&]
Eigenvalues[A]

Repeated comps : Full Max

Here are repeated runs reducing the maximum of diagonal each time

AA=A;
NOps=12000;
Data=Table[
TempMat=Abs[AA];
Do[TempMat[[k,k]]=0.0,{k,1,m}];
{i,j}=Position[TempMat,Max[TempMat]][[1]];
θ= 0.5 ArcTan[ 2 ((AA[[i,j]])/(AA[[j,j]]-AA[[i,i]]))];
JTemp=J2[{i,j},θ,m];
AA=JTemp.AA.Transpose[JTemp];
{n,((Norm[LowerTriangularize[AA,-1],"Frobenius"])/(Norm[AA,"Frobenius"]))},
{n,NOps}];
ListLogPlot[Data,
FrameLabel->{"J Ops","Frobenius OffDiagonal Fraction"},
GridLines->Automatic,
Frame->True,
PlotLabel->"Jacobi Max Search"]
Sort[Diagonal[AA],Abs[#1]>Abs[#2]&]
Eigenvalues[A]

MatrixPlot[AA]

ListPlot3D[Log[Abs[AA]]]

m

Eigen Stuff for 4×4

m=4;
A = RandomReal[{-1,1},{m,m}];
A=A+Transpose[A];
J4[{i_,j_,k_,l_},{p1_,p2_,p3_,p4_}] := Module[
{pNorm=Norm[{p1,p2,p3,p4}]},
(1/pNorm)(({{ p1,  p2,  p3,  p4}, {-p2,  p1,  p4, -p3}, {-p3, -p4,  p1,  p2}, {-p4,  p3, -p2,  p1}}))
]
J4[{i_,j_,k_,l_},{p1_,p2_,p3_,p4_},noScale] := Module[
{pNorm=Norm[{p1,p2,p3,p4}]},
(({{ p1,  p2,  p3,  p4}, {-p2,  p1,  p4, -p3}, {-p3, -p4,  p1,  p2}, {-p4,  p3, -p2,  p1}}))
]

So I should be able to compute fast the effect of combining a mini sweep in one left and right rotation matrix for a gener

Two competing ideas:  
1) A simple analytic approximation to an eigenvalue ( I can do someting rough and fast either using the analytics of Mathematica or using a low rank one update trick)
2) Just grind the damm thing down with a succesion of Givens or Householder ops.  Accumulate the ops as a small array and then broadcast the bundle of ops up to the entire matrix. 

I may be able to pretend that I have reduced it to completely to a diagonal.  

The analytics do not appear to be tractable in Mathematica with only radicals.

symsubs=Flatten[
Table[
a[i,j]->a[j,i],
{i,1,4},
{j,i+1,4}
]
]

Clear[A,a]
symsubs=Flatten[
Table[
a[j,i]->a[i,j],
{i,1,4},
{j,i+1,4}
]
];
ExpandSub = {a[1,4]->ϵ, a[1,3]->α[1,3] ϵ,a[2,4]-> α[2,4] ϵ,a[2,3]->α[2,3]ϵ}
AltExpandSub ={a[1,1]->ϵ1, a[2,2]->ϵ1 α[2,2], a[3,3]->ϵ1 α[3,3],a[4,4]->ϵ1 α[4,4],a[1,2]->ϵ2,a[3,4]->ϵ3}
A[ϵ1_,ϵ2_,ϵ3_]=Array[a,{4,4}]/.symsubs/.AltExpandSub;
MatrixForm[A[0,0,0]]
F[ϵ1_,ϵ2_,ϵ3_]:=Eigenvalues[A[ϵ1,ϵ2,ϵ3]]

D[F[0,0,ϵ3],ϵ3]/.ϵ3->0

F[0,0,0]

Simplify[F[0]]

A[ϵ]

F'[0]

JJ=J4[{1,2,3,4},A[[All,1]],noScale].A

J2[{i_,j_},θ_,m_] := SparseArray[{
{i,i}->Cos[θ],
{j,j}->Cos[θ],
{i,j}->-Sin[θ],
{j,i}->Sin[θ],
Band[{1,1}]->1.0
},
{m,m}]

Testing Gershorgin and Brauer: Complex Matrices

Gershoring is eay to visualize center is at diagonal radius is the missing stuff

m=4;
A = RandomComplex[{-1-I,1+I},{m,m}]+DiagonalMatrix[RandomComplex[{1+I , 14+3 I},{m}]];
Rs=Table[
Sum[Abs[A[[i,j]]],{j,1,m}]-Abs[A[[i,i]]],
{i,1,m}];
{RegionPlot[
Evaluate[
Flatten[
Table[
Table[
Abs[x+ I y-A[[i,i]]]Abs[x+I y - A[[j,j]]] <= Rs[[i]] Rs[[j]],
{j,i+1,m}],
{i,1,m}]
]
],
{x,-10,20},{y,-10,10},
AxesLabel->{"Re","Im"},
PlotLabel->"Brauer-Cassini",
Epilog->Map[ Point, Eigenvalues[A]/.CompToRealSub]
],
RegionPlot[
Evaluate[
Table[
Abs[x+ I y-A[[i,i]]]<= Rs[[i]],
{i,1,m}]
],
{x,-10,20},{y,-10,10},
AxesLabel->{"Re","Im"},
PlotLabel->"Gershorgin",
Epilog->Map[ Point, Eigenvalues[A]/.{z_Real:>{Re[z],Im[z]},z_Complex:>{Re[z],Im[z]}}]
]
}

MatrixForm[

?RegionPlot

CompToRealSub={z_Real:>{Re[z],Im[z]},z_Complex:>{Re[z],Im[z]}};

m=4;
A = RandomComplex[{-1-I,1+I},{m,m}]+DiagonalMatrix[RandomComplex[{12+10I , 14+3 I},{m}]]
(* A=A+ Transpose[A]; *)
Eval=Eigenvalues[A]/.CompToRealSub
Rads= Table[
Sum[Abs[A[[i,j]]],{j,1,m}]-Abs[A[[i,i]]],
{i,1,m}]
Centers=Table[A[[i,i]],{i,1,m}]/.CompToRealSub
GCircles=Table[
{Hue[i/m],Circle[ Centers[[i]], Rads[[i]]], Point[Eval[[i]]]},
{i,1,m}];

Graphics[
GCircles,
Frame->True,
AspectRatio->Automatic
]

?Circle

Testing Gershorgin and Brauer: Real Non-Symmetric Matrices

Gershoring is eay to visualize center is at diagonal radius is the missing stuff

m=4;
A = RandomReal[{-1,1},{m,m}]+DiagonalMatrix[RandomReal[{1, 2},{m}]];
Rs=Table[
Sum[Abs[A[[i,j]]],{j,1,m}]-Abs[A[[i,i]]],
{i,1,m}];
{RegionPlot[
Evaluate[
Flatten[
Table[
Table[
Abs[x+ I y-A[[i,i]]]Abs[x+I y - A[[j,j]]] <= Rs[[i]] Rs[[j]],
{j,i+1,m}],
{i,1,m}]
]
],
{x,-5,5},{y,-5,5},
AxesLabel->{"Re","Im"},
PlotLabel->"Brauer-Cassini",
Epilog->Map[ Point, Eigenvalues[A]/.CompToRealSub]
],
RegionPlot[
Evaluate[
Table[
Abs[x+ I y-A[[i,i]]]<= Rs[[i]],
{i,1,m}]
],
{x,-5,5},{y,-5,5},
AxesLabel->{"Re","Im"},
PlotLabel->"Gershorgin",
Epilog->Map[ Point, Eigenvalues[A]/.{z_Real:>{Re[z],Im[z]},z_Complex:>{Re[z],Im[z]}}]
]
}

Testing Gershorgin and Brauer: Real 
Symmetric Matrices

Gershoring is eay to visualize center is at diagonal radius is the missing stuff

m=4;
A = RandomReal[{-1,1},{m,m}]+DiagonalMatrix[RandomReal[{1, 2},{m}]];
A= A+Transpose[A];
Rs=Table[
Sum[Abs[A[[i,j]]],{j,1,m}]-Abs[A[[i,i]]],
{i,1,m}];
{RegionPlot[
Evaluate[
Flatten[
Table[
Table[
Abs[x+ I y-A[[i,i]]]Abs[x+I y - A[[j,j]]] <= Rs[[i]] Rs[[j]],
{j,i+1,m}],
{i,1,m}]
]
],
{x,-5,5},{y,-5,5},
AxesLabel->{"Re","Im"},
PlotLabel->"Brauer-Cassini",
Epilog->Map[ Point, Eigenvalues[A]/.CompToRealSub]
],
RegionPlot[
Evaluate[
Table[
Abs[x+ I y-A[[i,i]]]<= Rs[[i]],
{i,1,m}]
],
{x,-5,5},{y,-5,5},
AxesLabel->{"Re","Im"},
PlotLabel->"Gershorgin",
Epilog->Map[ Point, Eigenvalues[A]/.{z_Real:>{Re[z],Im[z]},z_Complex:>{Re[z],Im[z]}}]
]
}

New Idea

{i,j}={1,4};
AA= A[[{i,j},{i,j}]];
Rs=Table[
Sum[ Abs[A[[ii,k]]],{k,1,m}]-Abs[A[[i,i]]]-Abs[A[[j,j]]],
{ii,{i,j}}]
Evals=Eigenvalues[AA]
Eigenvalues[A]
MatrixForm[A]
RegionPlot[
Evaluate[
Table[
Abs[x+ I y-A[[i,i]]]<= Rs[[i]],
{i,1,m}]
],
{x,-5,5},{y,-5,5},
AxesLabel->{"Re","Im"},
PlotLabel->"Gershorgin",
Epilog->Map[ Point, Eigenvalues[A]/.{z_Real:>{Re[z],Im[z]},z_Complex:>{Re[z],Im[z]}}]
]

Break

FullSimplify[Cos[1/2 ArcTan[2 ((d )/(b-a))]]]

?FullSimplify

Cos[ArcTan[θ]/2]

Crude Dominant SV

The following is very crude.  There could be reasons why the special 3×3 in the algorithm work better with some algorihm than you might think.

A= RandomReal[{-1,1},{3,3}];
SampleSize = 100;
Data=Table[
x= RandomReal[{-1,1},{3}];
V=A.Transpose[A].x;
V=V/Norm[V];
U=A.V;
SVEstimate = Norm[U];
U=U/SVEstimate;
SVEstimate,
{SampleSize}];

{u,w,v}=SingularValueDecomposition[A];
SV=w[[1,1]];
ListPlot[Sort[Data],
Epilog->Table[
Line[{{1,w[[i,i]]},{SampleSize,w[[i,i]]}}],
{i,1,3}],
PlotRange->{0, 1.1 SV}
]

?SingularValueDecomposition

4×4 Jacobi with the 3×3 sub problem in Mathematica

J4[ ASub_List] := Module[
{Blah},

?BlockMatrix

Cos[θ] IdentityMatrix[4]

J[8,4,θ_] =SparseArray[
{
Band[{1,1}]->Cos[θ],
Band[{5,1}]->Sin[θ],
Band[{1,5}]->-Sin[θ]
},
{8,8}
]
MatrixForm[J[8,4,θ]]
MatrixForm[J4.Transpose[J4]]

Reduced Info Sparse Vector

Everyone is cool with a low rank approximation to a matrix.  Everyone is cool with a reduced non-zero entry count as an optimality condition for a vector.  Example all a diagonaliing basis is is a minimalist representation of a matrix!  Similarly with an SVD. A partial diagonalization is just something that has not gotten to the end yet.

Bad initial idea

m =1024;
ListPlot[
Sort[
Table[
v = RandomReal[NormalDistribution[1,128],{m}];
p=Ordering[v,54,Abs[#1]>Abs[#2]&];
Norm[v[[p]]]/Norm[v],
{100}]
]
]

The storage requirement for any sane sparse vector storage is one index (absolute or shift) for each non-zero entry along witht he value of the non-zero entry.

?Eigensystem

ListPlot3D[A]

A=Import[ "/home/faculty/struther/Desktop/Link to 4610/TestMatricesETC/cavity03.rua"];
{Evals,EVecs}=Eigensystem[A];
ListPlot[
Table[
Tooltip[{Re[Evals[[i]]],Im[Evals[[i]]]},i],
{i,1,Length[Evals]}],
PlotRange->All
]

Data=Table[
Sort[
Table[
v=EVecs[[eVecNum]];
vC=v[[Ordering[v, CutOff, Abs[#1]>Abs[#2]&]]];
((Norm[vC])/(Norm[v])),
{eVecNum,1,Length[EVecs]}]
],
{CutOff,{1,2,3,4,8,16,32,64,128}}];

Position[EVecs[[200]],1.0]

ListPlot[A[[All,134]]]

MatrixPlot[Abs[EVecs]]
ListDensityPlot[Abs[EVecs]]

ListContourPlot[
Transpose[Data],
Contours->Range[0,1,0.05],
DataRange->{
{1,2,3,4,8,16,32,64,128},
{0,100}
}
]