Matlab: solving the linear system with QR/Householder

Collapse
X
 
  • Time
  • Show
Clear All
new posts
  • lancer6238
    New Member
    • Sep 2007
    • 2

    Matlab: solving the linear system with QR/Householder

    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 main function is

    Code:
    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:
    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 y = Qt_times_b(v,b)
    NumberOfReflections = length(v);
    y = b;
    p = NumberOfReflections;
    for k = 1:NumberOfReflections,
       F = eye(p)-2*(v{k}*v{k}')/norm(v{k})^2;
       p = p -1;
       % Put F into "Q" to get Qk
       y = Qk*y;
    end
    I'm having trouble with the "Put F into Q to get Q_k" part.

    I know Q is implicitly defined by the Householder reflectors stored in v, computed in the householder function, and y can be obtained by Q'b = Q_n * ... * Q_3 * Q_2 * Q_1 * b

    Also, Qk =
    [ I_(k-1) 0 ]
    [ 0 F_k ]
    matrix, where F_k = I - 2 (v * v') / norm(v)^2

    But how do I put F into Q?

    Thank you.

    Regards,
    Rayne
Working...