In [1]:
##
using Plots
In [2]:
##
using LinearAlgebra
using SparseArrays
mat(n) = sparse(Diagonal(exp2.(0:-1:-n+1)))
mat(5)
Out[2]:
5×5 SparseMatrixCSC{Float64,Int64} with 5 stored entries:
  [1, 1]  =  1.0
  [2, 2]  =  0.5
  [3, 3]  =  0.25
  [4, 4]  =  0.125
  [5, 5]  =  0.0625
In [3]:
## Just to be clear, this is what NOT to do :)
function krylov_simple(A,b,k)
  X = zeros(size(A,1),k+1)
  X[:,1] = b
  for i=1:k
    X[:,i+1] = A*X[:,i]
  end
  return X
end
n = 20
@show cond(Array(mat(n)))
X = krylov_simple(mat(n), ones(n), 12)
@show cond(X)

cond(Array(mat(n))) = 524288.0
cond(X) = 6.879097053164631e17
Out[3]:
6.879097053164631e17
In [4]:
## Let's plot it
n = 20
k = 12
A = mat(n)
X = krylov_simple(A, ones(n), k)
plot(map(i -> cond(X[:,1:i]), 1:k+1), yscale=:log10,
  label="condition number of basis")
title!("matrix condition number $(cond(Array(A)))")
gui()
In [5]:
## Does this happen for a ranodm b?

X = krylov_simple(A, randn(n), k)
Out[5]:
20×13 Array{Float64,2}:
  0.543546    0.543546      0.543546     …   0.543546      0.543546   
  1.00571     0.502856      0.251428         0.000491071   0.000245535
 -0.929344   -0.232336     -0.058084        -2.21573e-7   -5.53932e-8 
  0.822108    0.102764      0.0128454        9.5706e-11    1.19632e-11
  3.07387     0.192117      0.0120073        1.74729e-13   1.09206e-14
 -1.43031    -0.0446972    -0.00139679   …  -3.96991e-17  -1.2406e-18 
  0.0867494   0.00135546    2.11791e-5       1.17567e-21   1.83699e-23
  0.593737    0.00463857    3.62388e-5       3.92902e-24   3.06955e-26
 -0.555166   -0.00216862   -8.47117e-6      -1.79384e-27  -7.00719e-30
 -0.583033   -0.00113874   -2.22409e-6      -9.19864e-31  -1.79661e-33
  0.323987    0.000316394   3.08978e-7   …   2.49591e-34   2.43741e-37
  1.02051     0.000498294   2.43308e-7       3.83872e-37   1.87437e-40
  0.592306    0.000144606   3.53042e-8       1.0879e-40    2.65599e-44
  0.15594     1.90356e-5    2.32369e-9       1.39852e-44   1.70717e-48
 -0.218002   -1.33058e-5   -8.12122e-10     -9.54645e-48  -5.82669e-52
 -0.0623602  -1.90308e-6   -5.80775e-11  …  -1.33339e-51  -4.0692e-56 
  0.478054    7.29452e-6    1.11306e-10      4.99111e-54   7.61584e-59
 -0.245054   -1.86961e-6   -1.4264e-11      -1.24926e-57  -9.53108e-63
  0.195294    7.44987e-7    2.8419e-12       4.86127e-61   1.85443e-66
 -0.231567   -4.41679e-7   -8.42436e-13     -2.81454e-64  -5.36831e-70
In [6]:
## What about a random matrix?
X = krylov_simple(randn(20,20), randn(n), 20)
Out[6]:
20×21 Array{Float64,2}:
  0.175076  -7.5399     15.3406   -59.9074   …   2.58551e11  -3.21132e11
  0.292286  -2.02152    24.0167    42.1057      -2.78872e10   1.15499e12
 -1.47749    0.820739  -17.7809    14.8141       2.22297e11  -9.95827e11
 -0.898173   0.85216    -8.07042   36.9296      -1.19939e10   3.91453e11
  0.38674   -7.17357   -16.59     -36.1547      -5.45763e11   1.61208e12
  2.01662    4.4708    -15.6532    89.4172   …   1.87627e11  -5.76516e11
  0.891571   7.52137    35.4555     4.19145     -3.16829e11   1.61105e12
  0.920549  -4.45145    -5.26734  120.305       -1.14389e11   1.17415e12
 -0.279078   2.20062   -10.2742    33.927       -4.93287e11   1.91727e12
  0.501988   4.37104   -11.7039    12.1633       1.4125e11   -1.11809e12
  1.41609   -3.28994    -7.05227   70.2076   …   2.81336e10  -5.26353e11
  0.816535  -5.85534    28.7678     6.50265     -3.7084e10    6.73427e11
 -0.259387  -1.99098    -3.68624   -3.03523     -4.07693e10   5.45946e11
  1.53055   -1.23563    47.412    -66.1656      -1.86042e11   1.50104e12
 -1.11034   -6.02484    13.7028    32.844       -4.54998e11   2.73407e12
 -0.838749  -1.641     -10.4355   -40.5021   …  -7.12188e10   7.28705e11
 -2.01293    5.74646     7.01269   -6.4355       6.02283e11  -1.29575e12
  0.278806   1.59997   -12.7925     6.2115       1.17598e11  -9.0691e11 
 -1.57591   -0.892614   29.0649     7.42665     -3.93857e11   1.89617e12
 -0.901986   1.52264    -3.79264  -32.1163       1.47565e11  -3.21103e11