Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

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

1 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.

share|improve this answer

Your Answer

 
discard

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.