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