Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I simply rewrote this algorithm

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

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)

%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);

   if ( colnorm(perms) == 0 )

   %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 


% 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'


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));

Thanks for your help...

share|cite|improve this question
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));


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.

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.