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

``````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 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

Be a part of the DaniWeb community

We're a friendly, industry-focused community of 1.18 million developers, IT pros, digital marketers, and technology enthusiasts learning and sharing knowledge.