I am implementing the “Lanczos” algorithm which is shown in the Trefethen book (and in the wiki), but shortly here: https://imgur.com/8qE2aR3

The algorithm creates an orthogonal matrix Q and a tridiagonal matrix T such that a matrix A can be written as: A = Q T Q’.

I check my implementation by seeing if my Q matrix is orthogonal: Q Q’ = Identity. However, this is not true as seen in the output full(Q * Q'). Can anyone help me understand what is wrong with my implementation?

MWE:

%% Lanzcos Iteration: following Trefethen
% alpha = main diagonal, beta = super/sub diagonal
% A = Q T_n Q*

n=10;
A = generateSPDmatrix(n);
b = rand(n,1); % vector of length n = 10
[Q,alpha,beta] = Lanczos(A,b,n);
full(Q * Q')

function [Q, alpha, beta] = Lanczos(A,b,iter)
    %% Some initialization
    [row, col] = size(b);
    Q          = zeros(row, iter); % Q = [q_1 | q_2 | ... | q_n]
    Q(:,1)     = b/norm(b);  % q_1 = b/||b||
    alpha = zeros(iter,1);   % Main diagonal
    beta  = zeros(iter-1,1); % beta = (beta_1,..., beta_n). 
    
    %% Perform the Lanczos iteration: three-term recurrence relation
    for i = 1:iter
        %% Construct column vector q_(i+1) as Q = [... | q_(i+1) | ...]
        v        = A*Q(:,i);      % A * q_i 
        alpha(i) = (Q(:,i))' * v; % q_i^T M * A * q_i
    
        v = v - alpha(i)*Q(:,i);  % q_2 = (A q_1 - alpha_1 q_1)
        if i > 1
            v = v - beta(i-1)*Q(:,i-1); % q_i = (A q_i - beta_(i-1) q_(i-1) - alpha_i q_i)
        end
        
        beta(i)  = norm(v); 
        v = v/beta(i); % Normalize 

        Q(:,i+1) = v;
    end
end

%% SPD generation code found on Math.SE
function A = generateSPDmatrix(n)
% Generate a dense n x n symmetric, positive definite matrix

A = rand(n,n); % generate a random n x n matrix

% construct a symmetric matrix using either
A = 0.5*(A+A'); 
%A = A*A';
% The first is significantly faster: O(n^2) compared to O(n^3)

% since A(i,j) < 1 by construction and a symmetric diagonally dominant matrix
%   is symmetric positive definite, which can be ensured by adding nI
A = A + n*eye(n);
end