In [1]:
## 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

A * V[:, 1:3] - V * H = [1.1102230246251565e-16 -2.220446049250313e-16 -1.6653345369377348e-16; 4.440892098500626e-16 2.220446049250313e-16 1.1102230246251565e-16; 5.551115123125783e-17 0.0 4.440892098500626e-16; 0.0 -2.7755575615628914e-17 -1.1102230246251565e-16; 1.1102230246251565e-16 0.0 0.0]
Out[1]:
5×3 Array{Float64,2}:
 1.11022e-16  -2.22045e-16  -1.66533e-16
 4.44089e-16   2.22045e-16   1.11022e-16
 5.55112e-17   0.0           4.44089e-16
 0.0          -2.77556e-17  -1.11022e-16
 1.11022e-16   0.0           0.0
In [2]:
## 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)

norm(A * V[:, 1:4] - V * H) = 1.4613251254722773e-15
Out[2]:
1.4613251254722773e-15
In [3]:
## What happens with fewer iterations
V,H = arnoldi(A,b,8)
Out[3]:
([0.06863738382997085 -0.21042145929283715 … -0.3414299337473762 -0.3450331087921401; -0.2987848820376377 -0.4087175275000519 … -0.057835021876971295 -0.41222027773848746; … ; -0.16421748776980216 0.35072834280629595 … -0.17700527828916574 0.010619088553484993; 0.006854712294678334 -0.07654914896302442 … 0.4924938542184442 -0.4076980011148628], [-3.450253435088775 5.657668132488314 … 7.75855074630627e-15 1.1521165965699964e-14; 5.657668132488315 3.3490956806035097 … -4.315992008230296e-15 -4.6351811278100286e-15; … ; 0.0 0.0 … 1.2991226740397241 2.1153821852921357; 0.0 0.0 … 0.0 1.7187413869914179])
In [4]:
## 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)')

norm((A * U[:, 1:4] - U[:, 1:4] * T) - rho * U[:, end] * (begin
    x = zeros(4)
    #= In[4]:23 =#
    x[4] = 1
    #= In[4]:23 =#
    x
end)') = 1.5356134357017496e-15
Out[4]:
1.5356134357017496e-15
In [5]:
##
# This way misses that last column :(
@show norm(A*U[:,1:4] - U[:,1:end-1]*T)

norm(A * U[:, 1:4] - U[:, 1:end - 1] * T) = 3.0626091849166674
Out[5]:
3.0626091849166674
In [6]:
##
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)

norm(S - T) = 0.0
Out[6]:
0.0
In [7]:
## We'll see this matrix come up a bunch
V[:,1:4]'*A*V[:,1:4] - T
Out[7]:
4×4 Array{Float64,2}:
  0.0           0.0          1.8348e-15    6.95465e-15
 -8.88178e-16   1.33227e-15  8.88178e-16  -2.72636e-15
  2.20841e-15  -1.77636e-15  0.0           2.22045e-16
  7.27566e-15  -3.22048e-15  3.55271e-15   5.55112e-16
In [8]:
## Let's try a lot of iterations
S,rho = lanczos_tri(A,b,15)
Out[8]:
([-3.450253435088775 5.657668132488314 … 0.0 0.0; 5.657668132488315 3.3490956806035097 … 0.0 0.0; … ; 0.0 0.0 … -0.8411914392828387 3.286310941311094; 0.0 0.0 … 3.2863109413110942 3.8537241463227527], 2.898031105009008)
In [9]:
##
diag(S)
Out[9]:
15-element Array{Float64,1}:
 -3.450253435088775 
  3.3490956806035097
 -2.5300801520141003
 -0.9307081008793854
 -0.3382489384144162
  3.771554559668153 
  0.6817533430301652
  2.1153821852921384
  3.3021281486942047
  2.910235387848049 
 -8.158622997562341 
  8.143875450409851 
 -4.690488219392312 
 -0.8411914392828387
  3.8537241463227527
In [10]:
##
diag(S,-1)
Out[10]:
14-element Array{Float64,1}:
 5.657668132488315    
 5.062750426882765    
 1.8087932613054551   
 3.062609184916668    
 3.3077810972396398   
 3.359179065503639    
 1.299122674039725    
 1.7187413869914197   
 2.064129896278676    
 5.105754512564956e-12
 1.2423633419604647   
 1.7223181524452278   
 0.12645144110000434  
 3.2863109413110942