1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
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;
    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
    Code (Text):
    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
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you help with the solution or looking for help too?
Draft saved Draft deleted



Similar Discussions: Matlab: Solving linear system with QR/Householder
Loading...