Symmetric Tridiagonal QR eigenvalues

We have the symmetric tridiagonal m×m matrix A∈(ℂ^(m×m)) with diagonal d∈(ℂ^m) and upper diagonal u∈(ℂ^(m-1)).  

	A=(({{(Subscript[d, 1]), (Subscript[u, 1]), 0, …, …, 0}, {(Subscript[u, 1]), (Subscript[d, 2]), (Subscript[u, 2]), ⋱, …, ⋮}, {0, (Subscript[u, 2]), (Subscript[d, 3]), ⋱, 0, ⋮}, {⋮, ⋱, ⋱, ⋱, (Subscript[u, m-2]), 0}, {⋮,  , 0, (Subscript[u, m-2]), (Subscript[d, m-1]), (Subscript[u, m-1])}, {0, …, …, 0, (Subscript[u, m-1]), (Subscript[d, m])}}))
The task is to mininize the work involved in a single step of a QR eigen-value iteration.

Reminder

The effect of pre and post multiplying by the outer product 
	A=a (b^T)
can be computed efficiently.
I)	A.u=a (b.u) 
	In other words A.u is a scaled copy of the vector a.
II)	A.[(Subscript[u, 1])⋮ (Subscript[u, 2])⋮ … ⋮ (Subscript[u, n])] = [A.(Subscript[u, 1])⋮ A.(Subscript[u, 2])⋮ … ⋮ A.(Subscript[u, n])]= [a (b.(Subscript[u, 1]))⋮ a (b.(Subscript[u, 2]))⋮ … ⋮ a (b.(Subscript[u, n]))]
	In other words A.U is a bunch of copies of a. scaled by inner products of the columns with b. 
III) 	U.A works similarly except now it is the rows of U that turn up and in inner products with a and multiply b. 

In a QR eigenvalue iteration that does not remember the Q we can use this to minimize the arithmetic.  You should ntoe that Algorithm 10.1 implements this updateing taking into account the sweep throught the matrix. The reverse updating is contained in the formulae 10.2 an 10.3.  Note:  if you make a complete sweep you need to store all the vectors that describe the rotations.  This is not as bad as it looks because they are pretty sparse.   You need to remember you perform the forward Householder rotations in a specific order using modified values as you go.  This means you need to do the rotations in the correct order when you reverse the order of the factors and post multiply by the product of all the householder rotations.

Classical Iterative Solvers: Jacobi

Split A=L+D+U (with L strictly lower triangular, D diagonal, and U strictly upper triangular) and iterate
	D (Subscript[x, n+1])=b-(L+U)(Subscript[x, n])  
remember "inverting D" is not really a matrix operation.

m=12;
A= RandomReal[{-1,1},{m,m}];
A=A+DiagonalMatrix[RandomReal[{3,12},m]];
Dia=Diagonal[A];
LU=A-DiagonalMatrix[Dia];
b=RandomReal[{-1,1},{m}];
x=RandomReal[{-1,1},{m}];
MatrixForm[
Table[
x= (b-LU.x)/Dia,
{5}
]
]
LinearSolve[A,b]

Classical Iterative Solvers: Gauss-Seidel

Split A=L+D+U (with L strictly lower triangular, D diagonal, and U strictly upper triangular) and iterate
	(L+D) (Subscript[x, n+1])=b-U (Subscript[x, n])  
remember "solving the L+D" equation is not really a matrix operation.It is just a back solve cause L+D is triangular!

m=12;
A= RandomReal[{-1,1},{m,m}];
A=A+DiagonalMatrix[RandomReal[{3,12},m]];
U=UpperTriangularize[A,1];
LD=A-U;
b=RandomReal[{-1,1},{m}];
x=RandomReal[{-1,1},{m}];
MatrixForm[
Table[
x= LinearSolve[LD,b-U.x],
{5}
]
]
LinearSolve[A,b]

Classical Iterative Solvers: SOR

The SOR (Successive Over Relaxation ) idea is just that if a little bit of Gauss Seidel is good then more should be better.  The more is described by a tuneable parameter ω. 
Split A=L+D+U (with L strictly lower triangular, D diagonal, and U strictly upper triangular), write the equation A x = b as
	 (ωL+D+(ω-1)D+ωU)x=ωb
rewrite as
	(D + ωL) x = ωb-[ωU+  (ω - 1) D ]x
and set up the iteration
	(D + ωL) (Subscript[x, n+1]) = ωb-[ωU +  (ω - 1) D ](Subscript[x, n])
remember "solving the D + ωU" equation is not really a matrix operation. It is just a back solve cause D+ω U is upper triangular!

Note with ω=1 this is Gauss-Seidel.  The usefull values of are between 1 and 2.

m=12;
A= RandomReal[{-1,1},{m,m}];
A=A+DiagonalMatrix[RandomReal[{3,12},m]];
ω=1.3;
L=LowerTriangularize[A,-1];
U=UpperTriangularize[A,1];
Dia=DiagonalMatrix[Diagonal[A]];
b=RandomReal[{-1,1},{m}];
x=RandomReal[{-1,1},{m}];
MatrixForm[
Table[
x= LinearSolve[(Dia+ω L),ω b-(ω U +(ω-1)Dia).x],
{12}
]
]
LinearSolve[A,b]

Krylov Subspaces and Arnoldi Iteration

Pick a random vector (OverVector[b]).  Hit it with a matrix A the chances are that you have caught some information about the largest eigenvalue. Hit it again and you have amplified whatever information you got about the largest eigenvalue and picked up some information about the second alrgest eigenvalue, …

The Arnoldi iteration filters and cleans this information at each n for the Krylov subspaces
	(Subscript[𝒦, n])=⟨b,Ab,(A^2)b,…,(A^(n-1))b⟩
by maintaining a list of orthogonal vectors (generated using modified gram-schmidt) that span (Subscript[𝒦, n]) i.e.
	(Subscript[𝒦, n])=⟨(Subscript[q, 1]),(Subscript[q, 2]),…,(Subscript[q, n])⟩.
The code is fairly easy!

m=12;
nMax=4;
A = RandomReal[{-1,1},{m,m}];
b = RandomReal[{-1,1},{m}];
Q=0*A;
Q[[All,1]]=b/Norm[b];
Do[
v=A.Q[[All,n]];
Do[
h = Q[[All,j]].v;
v=v-h Q[[All,j]],
{j,1,n}
];
h=Norm[v];
Q[[All,n+1]]=v/h,
{n,1, nMax-1}
]
MatrixForm[Q]
MatrixForm[Transpose[Q].Q]