how to avoid floating point error and figure out a way to break out of this loop such that it generates correct answer

function [] = lanczos(A, m)

this is a function whick takes A as input which is a matrix and reduces it to smaller size. And i have got a for loop

for j=2:inf
    w = A*V(:,j) - beta(j)*V(:,j-1);
    alpha(j) = w'*V(:,j);
    w = w - alpha(j)*V(:,j); 
    beta(j+1) = norm(w,2);
    loopcnt = loopcnt + 1;
    if abs(beta(j+1)) = 0
        break
    end
end

now i need the control to come out of loop when beta(j+1) is equal to zero ,but its looping continuously may be because of some floating point error.so i tried something like this

if abs(beta(j+1)) < 0.0001 should come out of loop which worked fine for smaller matrix sizes.

but when matrix size got bigger even its not working and loop is running continuosly can some one suggest me a way to avoid this problem and run the loop properly and break it when beta(j+1) become zero

[David Hill]

If you want help, you need to provide more information on what all your variables are. You should use the code block so it is easy to cut and paste. I also recommend a while loop instead of the for loop you are using.

[Nitish Reddy Kotkur]

function [] = lanczos(A, m)
A = readmatrix('output1.txt','Whitespace',' []'); %output1.txt contains matrix
[n,k] = size(A);
V = zeros(k,k);
V(:,2) = 0;
V(1,2)=1;
V(:,2)=V(:,2)';
V(:,2)=V(:,2)/norm(V(:,2),2);
beta(2)=0;
loopcnt = 0;

%%highlight starts
for j=2:inf
    w = A*V(:,j) - beta(j)*V(:,j-1);
    alpha(j) = w'*V(:,j);
    w = w - alpha(j)*V(:,j); 
    beta(j+1) = norm(w,2);
    loopcnt = loopcnt + 1;
    if abs(beta(j+1)) < 0.0001
    break
end
%%highlight ends
    V(:,j+1) = w/beta(j+1);
    
end
m=loopcnt-1;
T=sparse(m+1,m+1);
for i=2:m+1
    T(i-1,i-1)=alpha(i);
    T(i-1,i)=beta(i+1);
    T(i,i-1)=beta(i+1);
end 
disp(['number of loops: ' num2str(loopcnt)]);
T(m+1,m+1)=alpha(m+2);
V = V(:,2:loopcnt+1)

disp(['approximating eigenvalues are: ', num2str(eigs(T)')]);
T
j
end

This is my entire code,and i have highlighted the loop part. could you help me with it

[David Hill]

j=2;
while  abs(beta(j))>=0.0001
    w = A*V(:,j) - beta(j)*V(:,j-1);
    alpha(j) = w'*V(:,j);
    w = w - alpha(j)*V(:,j); 
    beta(j+1) = norm(w,2);
    j=j+1;
end
    

If you are sure abs(beta(j)) will approach zero, you could troubleshoot by printing beta(j+1) during each loop by removing the semicolon. You should also preallocate alpha and beta if you have some idea how many cycles the loop will make (it will dramatically increase the speed).

[Nitish Reddy Kotkur]

well its not working the problem is figuring out perfect suitable tolerance,which i am unable to find

Answer 1

[Christine Tobler]

There has been a lot of research on how to deal with floating-point error in the Lanczos algorithm, see the Wikipedia page for some discussion.

However, even in exact arithmetic, you should only expect beta(j) to become zero reliably after having performed n iterations (where the matrix A is of size n-by-n). Now usually, the goal of the Lanczos algorithm is to converge in significantly fewer iterations than that, so this is not usually a good convergence criterion.

It's still a good idea to check if beta(j) has become zero, because this would be a breakdown condition - if beta(j) becomes zero, it means that a complete eigenspace of A has been constructed. Either all information you need is already contained in this space, or you need to restart another Lanczos iteration with a new random vector.