QR decomposition with column pivoting, errror

I simply rewrote this algorithm

http://faculty.nps.edu/borges/Teaching/MA3046/Matlab/qrlsdiary.html

for QR ecomposition with columns pivoting to be applicable to a general matrix A (m, n). The solution is based on Golub algorithm

But it gives good results only for m <= 3. Could somebody help me to fix the code?

R = A;
Q = eye(m);

% I am going to use a permutation matrix.
P = eye(n);

% Compute the norms.
for i = 1 : n
colnorm(i) = R(:,i)'*R(:,i)
end

%Swapping procedure
for i=1:n

%Find max col norm
max_colnorm = colnorm(i); perms = i;
for  j = i + 1 : n
if (colnorm(j) > max_colnorm)
perms = j;
max_colnorm = colnorm(j);
end
end

%Break
if ( colnorm(perms) == 0 )
break;
end

%Swap P
temp = P(:, i);
P(:, i) = P (:, perms)
P(:, perms) = temp

%Swap R
temp = R(:, i);
R(:, i) = R (:, perms)
R(:, perms) = temp

%Swap colnorm
colnorm = colnorm*P

% Get the Householder vector from get_house.
v = get_house(R(:,i),i,m)

% Apply the transformation to R from the left.
R = R-v*(v'*R)

% And also apply it to Q from the right.
Q = Q - (Q*v)*v'

%Norm downdate
if i~=n
colnorm(i+1:n) = colnorm(i+1:n) - R(i, i+1 : n).^2
end

end

% Get the Householder vector from get_house.
v = get_house(R(:,n),n,m)

% Apply the transformation to R from the left.
R = R- v  *(v' * R)

% And also apply it to Q from the right.
Q = Q - (Q * v)* v'


where

function [v] = get_house(x,i,j)

% Initialization.
n = length(x);
v = zeros(n,1);

% Copy that part of x to be worked on to the corresponding positions in v.
v(i:j) = x(i:j);

% Compute the proper Householder vector.
v(i) = v(i) + mysign(x(i)) * norm(x(i:j));

% Normalize the result so that H = I - v*v'. Includes an error check for
% the trivial reflection.

if ((v'*v) > 0)
v = v*sqrt(2/(v'*v));
end


-
For me the code crashed on the line lens = lens*P when i = 1, because lens had not been defined yet. – littleO Nov 16 '12 at 9:38
@ littleO. Thanks for your comment, the code is completely updated... – justik Nov 16 '12 at 11:09

After making a few changes to your code, the code seems to work.

First, I change your mysign() function to -1, as I am not sure how your mysign() is defined. So the line

v(i) = v(i) + mysign(x(i)) * norm(x(i:j));


becomes

v(i) = v(i) - norm(x(i:j));


Second, in your main function body, I commented out the factorization procedure outside the for-loop and added one last line:

% Get the Householder vector from get_house.
%v = get_house(R(:,n),n,m);

% Apply the transformation to R from the left.
%R = R- v  *(v' * R);

% And also apply it to Q from the right.
%Q = Q - (Q * v)* v';

R = R*P';  % put the columns back to its original order!


I don't see why factorization is needed outside the loop. The lower part of R should have been zeroed out by the loop, no matter the loop terminates prematurely or not.

Finally, you forgot to put the columns of R back to its original order.

-