qr_unique.m

function [Q, R] = qr_unique(A)
% Thin QR factorization ensuring diagonal of R is real, positive if possible.
%
% function [Q, R] = qr_unique(A)
%
% If A is a matrix, then Q, R are matrices such that A = QR where Q'*Q = I
% and R is upper triangular. If A is real, then so are Q and R.
% This is a thin QR factorization in the sense that if A has more rows than
% columns, than Q has the same size as A.
% 
% If A has full column rank, then R has positive reals on its diagonal.
% Otherwise, it may have zeros on its diagonal.
%
% This is equivalent to a call to Matlab's qr(A, 0), up to possible
% sign/phase changes of the columns of Q and of the rows of R to ensure
% the stated properties of the diagonal of R. If A has full column rank,
% this decomposition is unique, hence the name of the function.
%
% If A is a 3D array, then Q, R are also 3D arrays and this function is
% applied to each slice separately.
%
% See also: qr randrot randunitary

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 18, 2019.
% Contributors: 
% Change log: 

    [m, n, N] = size(A);
    if m >= n % A (or its slices) has more rows than columns
        Q = zeros(m, n, N, class(A));
        R = zeros(n, n, N, class(A));
    else
        Q = zeros(m, m, N, class(A));
        R = zeros(m, n, N, class(A));
    end
    
    for k = 1 : N
        
        [q, r] = qr(A(:, :, k), 0);
        
        % In the real case, s holds the signs of the diagonal entries of R.
        % In the complex case, s holds the unit-modulus phases of these
        % entries. In both cases, d = diag(s) is a unitary matrix, and
        % its inverse is d* = diag(conj(s)).
        s = sign(diag(r));
        
        % Since a = qr (with 'a' the slice of A currently being processed),
        % it is also true that a = (qd)(d*r). By construction, qd still has
        % orthonormal columns, and d*r has positive real entries on its
        % diagonal, /unless/ s contains zeros. The latter can only occur if
        % slice a does not have full column rank, so that the decomposition
        % is not unique: we make an arbitrary choice in that scenario.
        % While exact zeros are unlikely, they may occur if, for example,
        % the slice a contains repeated columns, or columns that are equal
        % to zero. If an entry should be mathematically zero but is only
        % close to zero numerically, then it is attributed an arbitrary
        % sign dictated by the numerical noise: this is also fine.
        s(s == 0) = 1;
        
        Q(:, :, k) = bsxfun(@times, q, s.');
        R(:, :, k) = bsxfun(@times, r, conj(s));
        
    end

end