2
$\begingroup$

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

1 Answers 1

2

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.