|
/Users/lyon/j4p/src/math/linearAlgebra/QRDecomposition.java
|
1 package math.linearAlgebra;
2
3 import math.linearAlgebra.util.Maths;
4
5 /** QR Decomposition.
6 <P>
7 For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
8 orthogonal matrix Q and an n-by-n upper triangular matrix R so that
9 A = Q*R.
10 <P>
11 The QR decompostion always exists, even if the matrix does not have
12 full rank, so the constructor will never fail. The primary use of the
13 QR decomposition is in the least squares solution of nonsquare systems
14 of simultaneous linear equations. This will fail if isFullRank()
15 returns false.
16 */
17
18 public class QRDecomposition implements java.io.Serializable {
19
20 /* ------------------------
21 Class variables
22 * ------------------------ */
23
24 /** Array for internal storage of decomposition.
25 @serial internal array storage.
26 */
27 private double[][] QR;
28
29 /** Row and column dimensions.
30 @serial column dimension.
31 @serial row dimension.
32 */
33 private int m, n;
34
35 /** Array for internal storage of diagonal of R.
36 @serial diagonal of R.
37 */
38 private double[] Rdiag;
39
40 /* ------------------------
41 Constructor
42 * ------------------------ */
43
44 /** QR Decomposition, computed by Householder reflections.
45 @param A Rectangular matrix
46 Structure to access R and the Householder vectors and compute Q.
47 */
48
49 public QRDecomposition(Matrix A) {
50 // Initialize.
51 QR = A.getArrayCopy();
52 m = A.getRowDimension();
53 n = A.getColumnDimension();
54 Rdiag = new double[n];
55
56 // Main loop.
57 for (int k = 0; k < n; k++) {
58 // Compute 2-norm of k-th column without under/overflow.
59 double nrm = 0;
60 for (int i = k; i < m; i++) {
61 nrm = Maths.hypot(nrm, QR[i][k]);
62 }
63
64 if (nrm != 0.0) {
65 // Form k-th Householder vector.
66 if (QR[k][k] < 0) {
67 nrm = -nrm;
68 }
69 for (int i = k; i < m; i++) {
70 QR[i][k] /= nrm;
71 }
72 QR[k][k] += 1.0;
73
74 // Apply transformation to remaining columns.
75 for (int j = k + 1; j < n; j++) {
76 double s = 0.0;
77 for (int i = k; i < m; i++) {
78 s += QR[i][k] * QR[i][j];
79 }
80 s = -s / QR[k][k];
81 for (int i = k; i < m; i++) {
82 QR[i][j] += s * QR[i][k];
83 }
84 }
85 }
86 Rdiag[k] = -nrm;
87 }
88 }
89
90 /* ------------------------
91 Public Methods
92 * ------------------------ */
93
94 /** Is the matrix full rank?
95 @return true if R, and hence A, has full rank.
96 */
97
98 public boolean isFullRank() {
99 for (int j = 0; j < n; j++) {
100 if (Rdiag[j] == 0)
101 return false;
102 }
103 return true;
104 }
105
106 /** Return the Householder vectors
107 @return Lower trapezoidal matrix whose columns define the reflections
108 */
109
110 public Matrix getH() {
111 Matrix X = new Matrix(m, n);
112 double[][] H = X.getArray();
113 for (int i = 0; i < m; i++) {
114 for (int j = 0; j < n; j++) {
115 if (i >= j) {
116 H[i][j] = QR[i][j];
117 } else {
118 H[i][j] = 0.0;
119 }
120 }
121 }
122 return X;
123 }
124
125 /** Return the upper triangular factor
126 @return R
127 */
128
129 public Matrix getR() {
130 Matrix X = new Matrix(n, n);
131 double[][] R = X.getArray();
132 for (int i = 0; i < n; i++) {
133 for (int j = 0; j < n; j++) {
134 if (i < j) {
135 R[i][j] = QR[i][j];
136 } else if (i == j) {
137 R[i][j] = Rdiag[i];
138 } else {
139 R[i][j] = 0.0;
140 }
141 }
142 }
143 return X;
144 }
145
146 /** Generate and return the (economy-sized) orthogonal factor
147 @return Q
148 */
149
150 public Matrix getQ() {
151 Matrix X = new Matrix(m, n);
152 double[][] Q = X.getArray();
153 for (int k = n - 1; k >= 0; k--) {
154 for (int i = 0; i < m; i++) {
155 Q[i][k] = 0.0;
156 }
157 Q[k][k] = 1.0;
158 for (int j = k; j < n; j++) {
159 if (QR[k][k] != 0) {
160 double s = 0.0;
161 for (int i = k; i < m; i++) {
162 s += QR[i][k] * Q[i][j];
163 }
164 s = -s / QR[k][k];
165 for (int i = k; i < m; i++) {
166 Q[i][j] += s * QR[i][k];
167 }
168 }
169 }
170 }
171 return X;
172 }
173
174 /** Least squares solution of A*X = B
175 @param B A Matrix with as many rows as A and any number of columns.
176 @return X that minimizes the two norm of Q*R*X-B.
177 @exception IllegalArgumentException Matrix row dimensions must agree.
178 @exception RuntimeException Matrix is rank deficient.
179 */
180
181 public Matrix solve(Matrix B) {
182 if (B.getRowDimension() != m) {
183 throw new IllegalArgumentException("Matrix row dimensions must agree.");
184 }
185 if (!this.isFullRank()) {
186 throw new RuntimeException("Matrix is rank deficient.");
187 }
188
189 // Copy right hand side
190 int nx = B.getColumnDimension();
191 double[][] X = B.getArrayCopy();
192
193 // Compute Y = transpose(Q)*B
194 for (int k = 0; k < n; k++) {
195 for (int j = 0; j < nx; j++) {
196 double s = 0.0;
197 for (int i = k; i < m; i++) {
198 s += QR[i][k] * X[i][j];
199 }
200 s = -s / QR[k][k];
201 for (int i = k; i < m; i++) {
202 X[i][j] += s * QR[i][k];
203 }
204 }
205 }
206 // Solve R*X = Y;
207 for (int k = n - 1; k >= 0; k--) {
208 for (int j = 0; j < nx; j++) {
209 X[k][j] /= Rdiag[k];
210 }
211 for (int i = 0; i < k; i++) {
212 for (int j = 0; j < nx; j++) {
213 X[i][j] -= X[k][j] * QR[i][k];
214 }
215 }
216 }
217 return (new Matrix(X, n, nx).getMatrix(0, n - 1, 0, nx - 1));
218 }
219 }
220