QR decomposition with column pivoting, errror

Question

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

http://info.mcs.anl.gov/pub/tech_reports/reports/P551.pdf ,

see page 2, please

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

Thanks for your help...

Answer 1

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.