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
0001 function [Q, R] = qr_unique(A) 0002 % Thin QR factorization ensuring diagonal of R is real, positive if possible. 0003 % 0004 % function [Q, R] = qr_unique(A) 0005 % 0006 % If A is a matrix, then Q, R are matrices such that A = QR where Q'*Q = I 0007 % and R is upper triangular. If A is real, then so are Q and R. 0008 % This is a thin QR factorization in the sense that if A has more rows than 0009 % columns, than Q has the same size as A. 0010 % 0011 % If A has full column rank, then R has positive reals on its diagonal. 0012 % Otherwise, it may have zeros on its diagonal. 0013 % 0014 % This is equivalent to a call to Matlab's qr(A, 0), up to possible 0015 % sign/phase changes of the columns of Q and of the rows of R to ensure 0016 % the stated properties of the diagonal of R. If A has full column rank, 0017 % this decomposition is unique, hence the name of the function. 0018 % 0019 % If A is a 3D array, then Q, R are also 3D arrays and this function is 0020 % applied to each slice separately. 0021 % 0022 % See also: qr randrot randunitary 0023 0024 % This file is part of Manopt: www.manopt.org. 0025 % Original author: Nicolas Boumal, June 18, 2019. 0026 % Contributors: 0027 % Change log: 0028 0029 [m, n, N] = size(A); 0030 if m >= n % A (or its slices) has more rows than columns 0031 Q = zeros(m, n, N, class(A)); 0032 R = zeros(n, n, N, class(A)); 0033 else 0034 Q = zeros(m, m, N, class(A)); 0035 R = zeros(m, n, N, class(A)); 0036 end 0037 0038 for k = 1 : N 0039 0040 [q, r] = qr(A(:, :, k), 0); 0041 0042 % In the real case, s holds the signs of the diagonal entries of R. 0043 % In the complex case, s holds the unit-modulus phases of these 0044 % entries. In both cases, d = diag(s) is a unitary matrix, and 0045 % its inverse is d* = diag(conj(s)). 0046 s = sign(diag(r)); 0047 0048 % Since a = qr (with 'a' the slice of A currently being processed), 0049 % it is also true that a = (qd)(d*r). By construction, qd still has 0050 % orthonormal columns, and d*r has positive real entries on its 0051 % diagonal, /unless/ s contains zeros. The latter can only occur if 0052 % slice a does not have full column rank, so that the decomposition 0053 % is not unique: we make an arbitrary choice in that scenario. 0054 % While exact zeros are unlikely, they may occur if, for example, 0055 % the slice a contains repeated columns, or columns that are equal 0056 % to zero. If an entry should be mathematically zero but is only 0057 % close to zero numerically, then it is attributed an arbitrary 0058 % sign dictated by the numerical noise: this is also fine. 0059 s(s == 0) = 1; 0060 0061 Q(:, :, k) = bsxfun(@times, q, s.'); 0062 R(:, :, k) = bsxfun(@times, r, conj(s)); 0063 0064 end 0065 0066 end