MATLAB Central Newsreader - RQ decomposition

RQ decomposition

David Doria — Sat, 04 Oct 2008 15:40:36 -0400

Is there an easy way to get an RQ decomposition from the included qr() function? This seems like a useful thing to include for the next release! Its very helpful in camera calibration, because to get the intrinsic parameters from a camera projection matrix, you can simply use the RQ factorization of the projection matrix.

Thanks,

Dave

Re: RQ decomposition

John D'Errico — Sat, 04 Oct 2008 19:16:01 -0400

"David Doria" <daviddoria@gmail.com> wrote in message <gc82pk$e87$1@fred.mathworks.com>...
> Is there an easy way to get an RQ decomposition from the included qr() function? This seems like a useful thing to include for the next release! Its very helpful in camera calibration, because to get the intrinsic parameters from a camera projection matrix, you can simply use the RQ factorization of the projection matrix.

Have you looked on the internet?

Google easily found matlab code for an
RQ decomposition.

So then I looked on the file exchange.
An RQ decomposition is on the FEX. It
is mex code, but still useful.

John

Re: RQ decomposition

Bruno Luong — Sat, 04 Oct 2008 19:57:02 -0400

"David Doria" <daviddoria@gmail.com> wrote in message <gc82pk$e87$1@fred.mathworks.com>...
> Is there an easy way to get an RQ decomposition from the included qr() function?

function [R Q]=rq(A)
% function [R Q]=rq(A)
% A [m x n] with m<=n
% return R [m x n] triangular and
% Q [n x n] unitary (i.e., Q'*Q=I)
% such that A=R*Q
% Author: Bruno Luong
% Last Update: 04/Oct/2008

[m n]=size(A);
if m>n
error('RQ: Number of rows must be smaller than column');
end

[Q R]=qr(flipud(A).');
R=flipud(R.');
R(:,1:m)=R(:,m:-1:1);
Q=Q.';
Q(1:m,:)=Q(m:-1:1,:);

end
% Bruno

Re: RQ decomposition

David Doria — Sat, 04 Oct 2008 23:38:01 -0400

Bruno,

Could you explain in a couple of sentences what that is doing in a little bit "mathy" terms? ie. I've never seen a "flip upside down" step in a linear algebra kind of definition- can you explain why that is done?

Thanks,

Dave

Re: RQ decomposition

Bruno Luong — Sun, 05 Oct 2008 02:23:02 -0400

"David Doria" <daviddoria@gmail.com> wrote in message <gc8uop$2nh$1@fred.mathworks.com>...
> Bruno,
>
> Could you explain in a couple of sentences what that is doing in a little bit "mathy" terms? ie. I've never seen a "flip upside down" step in a linear algebra kind of definition- can you explain why that is done?
>

Hi,

Flipping upside down is multiplying a permutation matrix on the left side. Flipping rows (also in the above but only m rows, not n), is multiplying a permutation matrix on the right side.

I use permutation matrices for the purpose to transform a lower triangular matrix to upper one.

Bruno

Re: RQ decomposition

Bruno Luong — Sun, 05 Oct 2008 13:32:01 -0400

I have extended my rq with more options:

function [R Q]=rq(A, varargin)
% function [R Q]=rq(A)
%
% Perform RQ factorization
% A [m x n]
% returns R [m x n] triangular superior and
% Q [n x n] unitary (i.e., Q'*Q=I)
% such that A = R*Q
%
% Note: Standard RQ requires m<=n.
%
% For n ~=m, default output R has non-zero elements populate
% at the m-first rows
% Use [R Q] = rq(A, 'last') to populate non-sero R-rows at other end.
%
% For "non standard" RQ, i.e. for A with m > n, zero are padded in the
% outputs R and Q. However Q is no longer unitary. More precisely
% Q'*Q = eye(n) but Q*Q' is
% diag([zeros(1,m-n) ones(1,n)]) without 'last' option, and
% diag([ones(1,n) zeros(1,m-n)]) with 'last' and R becomes triangular.
%
% Author: Bruno Luong
% last Update: 05/Oct/2008

[m n]=size(A);

[Q R]=qr(flipud(A).');
R=flipud(R.'); % m x n
Q=Q.'; % n x n
if m>n
    warning('RQ:DimensionBizarre',...
            'RQ: number of rows is larger the number of columns');
    % Padding zeros
    R(end,m)=0; % m x m
    Q(m,end)=0; % m x n
end

R(:,1:m)=R(:,m:-1:1);
Q(1:m,:)=Q(m:-1:1,:);

last=strcmpi(getoption('',varargin{:}),'last');
% Populate R at last rows
if xor(m>n,last)
    R=circshift(R,[0 n-m]);
    Q=circshift(Q,[n-m 0]);
end

end

% Get option if provided
function res=getoption(default, option)
  if nargin<2 || isempty(option)
      res=default;
  else
      res=option;
  end
end

% Bruno

Re: RQ decomposition

David Doria — Sun, 05 Oct 2008 13:36:02 -0400

Hmm I guess what I am wondering is the following question:

Q and R in the QR decomposition of A are the same Q and R in the RQ decomposition of which matrix?

For square 3x3 matices, I got it down to this:

ReverseRows = [0 0 1; 0 1 0 ; 1 0 0]; %right multiply by ReverseRows reverses columns
[Q R] = qr((ReverseRows * A)');
R = ReverseRows * R' * ReverseRows;
Q = ReverseRows * Q';

but I don't know how to put the last two lines INSIDE the rq operator? I guess I'm completely missing the intuition?

Dave

Re: RQ decomposition

Bruno Luong — Sun, 05 Oct 2008 13:59:02 -0400

This "translation" might help you:

QR is Gram-Schmidt orthoganilization of columns of A, started from the first.

RQ is Gram-Schmidt orthoganilization of rows of A, started from the last.

Bruno

Re: RQ decomposition

David Doria — Sun, 05 Oct 2008 14:52:02 -0400

ah excellent! exactly what i was looking for.

Dave