I am trying to implement Lanczos biorthogonalization algorithm to construct a Krylov subspace basis. Given a matrix $A$, and vectors $v,w:(v,w)=1$, it produces three matrices $V,W,T$ such that $V^T W = W^T V = Id, \quad V^T A V = T$.Authors argue, that the biggest problem can occur when a denominator becomes very small during renormalization. However, my code produces results, that satisfy these conditions, only for small dimensions, like 10. If $\dim A =1000$, the code produces rubbish with $||T-V^T A V||\approx 10^{29}$.What should I consider to implement a working (robust) function?Matlab code:function [V, T, W] = lanczos(A, v, m)% Based on: Iterative Solution of Large Linear Systemsv = v / norm(v);w = rand(length(v), 1);tmp = dot(v, w);w = w / tmp;a = [];b = 0;d = 0;V = [zeros(length(v), 1), v];W = [zeros(length(v), 1), w];for j = 1:m    v = V(:, j+1);    w = W(:, j+1);    a(j) = dot(A*v, w);    v_hat = A * v - a(j) * v - d(j) * V(:, j);    w_hat = A' * w - a(j) * w - d(j) * W(:, j);    tmp = dot(v_hat, w_hat);    d(j+1) = sqrt(abs(tmp));    if abs(d(j+1)) < 10^(-10)        error('Delta near zero');    end    b(j+1) = tmp / d(j+1);    if abs(b(j+1)) < 10^(-10)        error('Beta near zero');    end    v_hat = v_hat / d(j+1);    w_hat = w_hat / b(j+1);    W(:, j+2) = w_hat;    V(:, j+2) = v_hat;endV = V(:, 2:(m + 1));W = W(:, 2:(m + 1));b = b(2:m);d = d(2:m);T = diag(a) + diag(b, 1) + diag(d, -1);endJulia code:function lanczos(A,v,m)    v = v/norm(v,2);    len_v = length(v);    w=rand(len_v,1);    tmp = dot(v,w);    w = w/tmp;    a = zeros(m);    b = zeros(m+1);    d = zeros(m+1);    V = zeros(len_v,m+2);    W = zeros(len_v,m+2);    V[:,2] = v;    W[:,2] = w;    for j=1:m        v = V[:, j+1];        w = W[:, j+1];        a[j] = dot(A*v, w);        v_hat = A*v - a[j]*v        w_hat = A'*w - a[j]*w - d[j]*W[:, j];        tmp = dot(v_hat, w_hat);        d[j+1] = sqrt(abs(tmp));        b[j+1] = tmp / d[j+1];        v_hat = v_hat / d[j+1];        w_hat = w_hat / b[j+1];        W[:, j+2] = w_hat;        V[:, j+2] = v_hat;    end    V = V[:, 2:(m + 1)];    W = W[:, 2:(m + 1)];    b = b[2:m];    d = d[2:m];    T = diagm(a,0) + diagm(b, 1) + diagm(d, -1);    println(d')    return V, T, W;enddim = 100;A=rand(dim,dim);A = 0.5*A /norm(A,2);v = rand(dim);m = convert(Int64,dim/5);V, T, W = lanczos(A,v,m)I = eye(m,m);df = [norm(I - V'*W,2), norm(I - W'*V,2) ];df2 = norm(V'*A*V - T,2);println(df)println(df2);[1]: Saad, Y. (1992). Analysis of Some Krylov Subspace Approximations to the Matrix Exponential Operator. SIAM Journal on Numerical Analysis, 29(1), 209228. https://doi.org/10.1137/0729014[2]: Auzinger, W. (2011). Iterative Solution of Large Linear Systems. Lecture Notes, 178.[3]: OLeary, D. P. (2007). Bi-orthogonalization Methods.
How do I implement a working Lanczos biorthogonalization method for Krylov subspace?
linear algebra;reference request;iterative method;krylov method