How to get my QR factorization code to give Q and R and make them stored in the Q and R variables

Hello folks,

I'm new to Matlab and tackling numerical linear algebra with it.

I am trying to get some practice ahead and right now, I'm having issues with my basic QR factorization.

Basically, I'm writing my own general Gram Schmidt code that would take in a matrix A of size mxn and spit out Q(mxn) and R(nxn)

to give you the context of the calculation:

1_ given a matrix A =[a_1, a_2,..... a_n] where a_i are vectors of size (m, 1)

2_ I calculate the elements of Q=[q_1, q_2,... q_n] where q_i=h_i/||h_i||

3_Where h_i is a vector built from the following expression : h_i=a_i-sum(j=0 to i, q_j*(q_j)'*a_i

4_ I calculate the elements of R=[r_1, r_2,..... r_n] as follows : r_11=||a_1||=a_1*(a_1)', for i<j r_ij=(q_i)'*(a_j), for i=j r_ij=||h_i|| and else r_ij=0

I finally got my code running but it spits out intermediary R that don't match the final R I am supposed to get and there is no sign of Q.

Any insight on what is wrong is welcome, anything to improve on it as well.

attached is my code. thank you!

Answer 1

[Shantanu Dixit]

Hi Douji,

It looks like you're facing some issues with the correctness of R and the absence of Q while implementing the Gram-Schmidt process.

Briefly the QR factorization can be broken down into following steps:

1. Orthogonalization of columns: For each column of matrix A, subtract the projection of onto all previously computed using:

v = v - (Q(:, j)' * A(:, i)) * Q(:, j);

2. Normalization for Orthonormal Vectors: After making v orthogonal to all prior vectors, normalize v to get the next column of Q:

Q(:, i) = v / norm(v);

3. Constructing R: Once all columns of Q are computed, R can be computed at once as (follows from the fact that :

R = Q' * A;

Combining the above steps QR factorization can be computed as:

function [Q, R] = qr_factorization(A)

    [m, n] = size(A);
    Q = zeros(m, n);
    for i = 1:n
        % Compute the ith column of Q
        v = A(:, i);
        for j = 1:i-1
            % Subtract the projection of A(:,i) onto the previous Q vectors
            v = v - (Q(:, j)' * A(:, i)) * Q(:, j);
        end
        % Normalize to get q_i
        Q(:, i) = v / norm(v);
    end
    
    % Compute R as Q' * A
    R = Q' * A;
end

% Example Matrix
A = [1, 1, 1; 1, 2, 3; 1, 3, 6];
[Q, R] = qr_factorization(A);

disp('Computed Q:');
disp(Q);

disp('Computed R:');
disp(R);