## The Arnoldi and Lanczos processes
using LinearAlgebra
function arnoldi(A,b,k)
n = size(A,1)
V = zeros(n,k+1)
H = zeros(k+1,k)
V[:,1] = b/norm(b)
for j=1:k
y = A*V[:,j]
for i=1:k
H[i,j] = V[:,i]'*y
y -= H[i,j]*V[:,i]
end
H[j+1,j] = norm(y)
V[:,j+1] = y/H[j+1,j]
end
return V,H
end
A = randn(5,5)
b = randn(5)
V,H = arnoldi(A,b,3)
@show A*V[:,1:3] - V*H
## Try it on a symmetric matrix.
A = randn(10,10)
A = A + A'
b = randn(10)
V,H = arnoldi(A,b,4)
@show norm(A*V[:,1:4] - V*H)
## What happens with fewer iterations
V,H = arnoldi(A,b,8)
## That is expected as we will/have shown.
function lanczos(A,b,k)
n = size(A,1)
V = zeros(n,k+1)
T = Tridiagonal(zeros(k-1), zeros(k), zeros(k-1))
rho = 0.0
V[:,1] = b/norm(b)
for j=1:k
y = A*V[:,j]
for i=max(j-1,1):j
T[i,j] = V[:,i]'*y
y -= T[i,j]*V[:,i]
end
rho = norm(y)
V[:,j+1] = y/rho
if j < k
T[j+1,j] = rho
end
end
return V,T,rho
end
U,T,rho = lanczos(A,b,4)
@show norm(A*U[:,1:4] - U[:,1:4]*T - rho*U[:,end]*(x = zeros(4); x[4]=1; x)')
##
# This way misses that last column :(
@show norm(A*U[:,1:4] - U[:,1:end-1]*T)
##
function lanczos_tri(A,b,k)
n = size(A,1)
V = zeros(n,k+1)
T = Tridiagonal(zeros(k-1), zeros(k), zeros(k-1))
rho = 0.0
vold = 0.0
v = b/norm(b)
for j=1:k
y = A*v
if j>1
T[j-1,j] = vold'*y
y -= T[j-1,j]*vold
end
T[j,j] = v'*y
y -= T[j,j]*v
vold = v
rho = norm(y)
v = y/rho
if j < k
T[j+1,j] = rho
end
end
return T,rho
end
S,rho = lanczos_tri(A,b,4)
@show norm(S-T)
## We'll see this matrix come up a bunch
V[:,1:4]'*A*V[:,1:4] - T
## Let's try a lot of iterations
S,rho = lanczos_tri(A,b,15)
##
diag(S)
##
diag(S,-1)