Matlab: Solving linear system with QR/Householder

[wu_weidong]

Homework Statement

Hi all,
I'm trying to implement the QR method for solving the linear system Ax = b. The QR factorization is achieved using Householder method.

The Attempt at a Solution

The main function is

function x = lin_solve(A,b)
[R,v] = householder(A);
y = Qt_times_b(v,b);
x = R\y;

Here are the individual functions:

function [R,v] = householder(A)

[m,n] = size(A);

if m>=n,
    NumberOfReflections = n;
else
    NumberOfReflections = m-1;
end

R = A;

v = cell(NumberOfReflections,1);

for k = 1:NumberOfReflections,
    x = R(k:m,k);

    xnorm = norm(x);

    if xnorm>0,
        % Compute the normal vector of the reflector
        v{k} = -x;
        v{k}(1) = v{k}(1) - sign(x(1))*xnorm;
        v{k} = (sqrt(2)/norm(v{k}))*v{k};

	% Update R
        for j = k:NumberOfReflections,
		R(k:m,j) = R(k:m,j) - (v{k}'*R(k:m,j))*v{k};
	end
    else
        v{k} = zeros(m-k+1,1);
    end
end

function [Q,y] = Qt_times_b(v,b)

NumberOfReflections = length(v);
y = b;
mv = length(v{1});
Q=eye(mv);
for k = 1:NumberOfReflections,
   F = eye(mv-k+1)-2*(v{k}*v{k}')/norm(v{k})^2;
   Qk = [eye(k-1) zeros(k-1,mv-k+1);zeros(mv-k+1,k-1) F];
   y = Qk*y;
end

What I don't understand is the function works for square matrices, i.e. A that are n-by-n, but if m < n, or m > n, then I get the incorrect solution x.

I tried to check my R and Qk against the qr factorization function provided by matlab. Say [QR, RR] = qr(A).

If m < n, then my R is not equal to RR; if m > n, R = RR. In both cases, Q' = QR, where Q' = transpose of matrix Q.

What is wrong with my program?

Thank you.

Regards,
Rayne

 

[mmwave]

Hi Rayne,

It seems like your main issue is with the dimensions of your matrices. In your lin_solve function, you are using the backslash operator to solve for x, which is only valid for square matrices. If m > n, then R is not square and the backslash operator will not work. You can try using the pseudoinverse function, pinv, instead. Also, make sure to check the dimensions of your matrices at each step to ensure they are compatible for the operations you are performing.

Another issue could be with your implementation of the householder function. Make sure to follow the steps correctly and check for any errors in your code. You can also try debugging by using a simpler example and comparing your results to the hand calculations.

Hope this helps! Good luck with your implementation.

 

[tacman]

Your code looks correct, but there may be some issues with how you are calling the functions or how you are setting up your matrices. Make sure that your input matrix A is a full rank matrix, otherwise the QR factorization may not work correctly. Also, check your indexing when updating R in the householder function, as it may be causing some issues when m < n. Another possibility is that you are not transposing the Q matrix correctly in the Qt_times_b function. It may be helpful to compare your code with examples or other implementations to see if there are any differences. Good luck with your implementation!