1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
function [R, Q] = rq (A, m, n)
%RQ economy RQ or QL factorization of a full matrix A.
% No special handling is done for rank-deficient matrices.
%
% [R,Q] = rq (A)
% [R,Q] = rq (A,m,n)
%
% If A is m-by-n with m <= n, then R*Q=A is factorized with R upper triangular
% and m-by-m. Q is m-by-n with orthonormal rows, where Q*Q' = eye (m), but
% Q'*Q is not identity. RQ works quickly when A is upper trapezoidal, but
% also works in the general case. With n=3 and m=5, an upper trapezoidal A:
%
% x x x x x
% . x x x x
% . . x x x
%
% The factorization is R*Q = A where R is upper triangular and m-by-m,
% and Q is m-by-n:
%
% R * Q = A
% x x x x x x x x x x x x x
% . x x x x x x x . x x x x
% . . x x x x x x . . x x x
%
% Q also happens to be upper trapezoidal if A is upper trapezoidal.
% With two optional input arguments (m,n), only A (1:m,1:n) is factorized.
%
% If m > n, then Q*R=A is computed where "R" is lower triangular and Q
% has orthonormal columns (Q'*Q is identity).
%
% Example
%
% A = rand (3,4), [R, Q] = rq (A), norm (R*Q-A), norm (Q*Q'-eye(3))
% C = rand (4,3), [L, Q] = rq (C), norm (Q*L-C), norm (Q'*Q-eye(3))
%
% See also qr.
% Copyright 2011-2012, Timothy A. Davis, http://www.suitesparse.com
if (issparse (A))
% RQ would actually work, but it would be very inefficient since no fill
% reducing ordering is used. That would require a row permutation of R.
error ('FACTORIZE:rq:sparse', 'RQ is not designed for sparse matrices.') ;
end
if (nargin == 1)
[m, n] = size (A) ;
end
if (m <= n)
%---------------------------------------------------------------------------
% RQ factorization of a short-and-fat matrix A
%---------------------------------------------------------------------------
[Q, R] = qr (A (m:-1:1, n:-1:1)', 0) ;
R = R (end:-1:1, end:-1:1)' ;
Q = Q (end:-1:1, end:-1:1)' ;
% Below is a step-by-step working description of the algorithm. Each of
% the error norms will be small. This code will operate correctly if
% uncommented, it will just be slower than the 3 lines of code above.
% (1) The A matrix is transposed and its rows and columns are reversed.
% The row/column reversal can be viewed as multiplication of A by row and
% column permutations, so this operation makes sense in terms of linear
% algebra: H = (Pm*A*Pn)' where Pm and Pn are permutation matrices of
% size m and n, respectively.
% H = A (m:-1:1, n:-1:1)' ;
% H now has the following form. This is good, because qr(H) can exploit
% the 3 zeros in the lower triangular part, to reduce the computation time.
%
% x x x
% x x x
% x x x
% . x x
% . . x
%
% We could instead factorize A', which has the following shape:
%
% x . .
% x x .
% x x x
% x x x
% x x x
%
% but the QR method in MATLAB cannot exploit the zeros in upper triangular
% part A'.
% (2) The QR factorzation of H is computed. QR in MATLAB takes advantage
% of the zeros in H. The resulting Q is n-by-m, R is m-by-m.
% [Q, R] = qr (H, 0) ;
% err = norm (Q*R-H)
% Q * R = H
%
% x x x 1 x x x x x (a "1" denotes the R(1,1) entry,
% x x x . x x x x x so it can be followed in the
% x x x . . x x x x operations below)
% . x x . x x
% . . x . . x
% (3) The columns of R and H are reversed. This is the same as multiplying
% both sides of the equation by the Pm m-by-m permutation matrix on the
% right.
% R = R (:, end:-1:1) ;
% H = H (:, end:-1:1) ;
% err = norm (Q*R-H)
% Q * R = H
%
% x x x x x 1 x x x
% x x x x x . x x x
% x x x x . . x x x
% . x x x x .
% . . x x . .
% (4) Both sides of the equation are transposed.
% H = H' ;
% R = R' ;
% Q = Q' ;
% err = norm (R*Q-H)
% R * Q = H
%
% x x x x x x . . x x x x x
% x x . x x x x . x x x x .
% 1 . . x x x x x x x x . .
% (5) The columns of Q and H are reversed. This is the same as multiplying
% both sides of the equation by the Pn n-by-n permutation matrix on the
% right. H is now equal to A again.
% H = H (:, end:-1:1) ;
% Q = Q (:, end:-1:1) ;
% err = norm (A-H)
% err = norm (R*Q-H)
% R * Q = A
%
% x x x . . x x x x x x x x
% x x . . x x x x . x x x x
% 1 . . x x x x x . . x x x
% (6) The columns of R and rows of Q are reversed. This the same as
% inserting the product of Pm*Pm' = I between R and Q, where Pm is the
% m-by-m permutation matrix.
% R = R (:, end:-1:1) ;
% Q = Q (end:-1:1, :) ;
% err = norm (R*Q-H)
% err = norm (Q*Q' - eye (m))
% R * Q = A
%
% x x x x x x x x x x x x x
% . x x . x x x x . x x x x
% . . 1 . . x x x . . x x x
else
%---------------------------------------------------------------------------
% QL factorization of a tall-and-thin matrix A
%---------------------------------------------------------------------------
[R, Q] = rq (A', n, m) ;
R = R' ;
Q = Q' ;
end
|
