Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab: Solving linear system with QR/Householder

  1. Sep 18, 2007 #1
    1. The problem statement, all variables and given/known data
    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.

    3. The attempt at a solution
    The main function is

    Code (Text):
    function x = lin_solve(A,b)
    [R,v] = householder(A);
    y = Qt_times_b(v,b);
    x = R\y;
    Here are the individual functions:

    Code (Text):
    function [R,v] = householder(A)

    [m,n] = size(A);

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

    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};
            v{k} = zeros(m-k+1,1);
    Code (Text):
    function [Q,y] = Qt_times_b(v,b)

    NumberOfReflections = length(v);
    y = b;
    mv = length(v{1});
    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;
    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.

  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted