001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.linear;
019
020 import java.util.Arrays;
021
022 import org.apache.commons.math.MathRuntimeException;
023
024
025 /**
026 * Calculates the QR-decomposition of a matrix.
027 * <p>The QR-decomposition of a matrix A consists of two matrices Q and R
028 * that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is
029 * upper triangular. If A is m×n, Q is m×m and R m×n.</p>
030 * <p>This class compute the decomposition using Householder reflectors.</p>
031 * <p>For efficiency purposes, the decomposition in packed form is transposed.
032 * This allows inner loop to iterate inside rows, which is much more cache-efficient
033 * in Java.</p>
034 *
035 * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
036 * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
037 *
038 * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $
039 * @since 1.2
040 */
041 public class QRDecompositionImpl implements QRDecomposition {
042
043 /**
044 * A packed TRANSPOSED representation of the QR decomposition.
045 * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
046 * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
047 * from which an explicit form of Q can be recomputed if desired.</p>
048 */
049 private double[][] qrt;
050
051 /** The diagonal elements of R. */
052 private double[] rDiag;
053
054 /** Cached value of Q. */
055 private RealMatrix cachedQ;
056
057 /** Cached value of QT. */
058 private RealMatrix cachedQT;
059
060 /** Cached value of R. */
061 private RealMatrix cachedR;
062
063 /** Cached value of H. */
064 private RealMatrix cachedH;
065
066 /**
067 * Calculates the QR-decomposition of the given matrix.
068 * @param matrix The matrix to decompose.
069 */
070 public QRDecompositionImpl(RealMatrix matrix) {
071
072 final int m = matrix.getRowDimension();
073 final int n = matrix.getColumnDimension();
074 qrt = matrix.transpose().getData();
075 rDiag = new double[Math.min(m, n)];
076 cachedQ = null;
077 cachedQT = null;
078 cachedR = null;
079 cachedH = null;
080
081 /*
082 * The QR decomposition of a matrix A is calculated using Householder
083 * reflectors by repeating the following operations to each minor
084 * A(minor,minor) of A:
085 */
086 for (int minor = 0; minor < Math.min(m, n); minor++) {
087
088 final double[] qrtMinor = qrt[minor];
089
090 /*
091 * Let x be the first column of the minor, and a^2 = |x|^2.
092 * x will be in the positions qr[minor][minor] through qr[m][minor].
093 * The first column of the transformed minor will be (a,0,0,..)'
094 * The sign of a is chosen to be opposite to the sign of the first
095 * component of x. Let's find a:
096 */
097 double xNormSqr = 0;
098 for (int row = minor; row < m; row++) {
099 final double c = qrtMinor[row];
100 xNormSqr += c * c;
101 }
102 final double a = (qrtMinor[minor] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
103 rDiag[minor] = a;
104
105 if (a != 0.0) {
106
107 /*
108 * Calculate the normalized reflection vector v and transform
109 * the first column. We know the norm of v beforehand: v = x-ae
110 * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
111 * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
112 * Here <x, e> is now qr[minor][minor].
113 * v = x-ae is stored in the column at qr:
114 */
115 qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
116
117 /*
118 * Transform the rest of the columns of the minor:
119 * They will be transformed by the matrix H = I-2vv'/|v|^2.
120 * If x is a column vector of the minor, then
121 * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
122 * Therefore the transformation is easily calculated by
123 * subtracting the column vector (2<x,v>/|v|^2)v from x.
124 *
125 * Let 2<x,v>/|v|^2 = alpha. From above we have
126 * |v|^2 = -2a*(qr[minor][minor]), so
127 * alpha = -<x,v>/(a*qr[minor][minor])
128 */
129 for (int col = minor+1; col < n; col++) {
130 final double[] qrtCol = qrt[col];
131 double alpha = 0;
132 for (int row = minor; row < m; row++) {
133 alpha -= qrtCol[row] * qrtMinor[row];
134 }
135 alpha /= a * qrtMinor[minor];
136
137 // Subtract the column vector alpha*v from x.
138 for (int row = minor; row < m; row++) {
139 qrtCol[row] -= alpha * qrtMinor[row];
140 }
141 }
142 }
143 }
144 }
145
146 /** {@inheritDoc} */
147 public RealMatrix getR() {
148
149 if (cachedR == null) {
150
151 // R is supposed to be m x n
152 final int n = qrt.length;
153 final int m = qrt[0].length;
154 cachedR = MatrixUtils.createRealMatrix(m, n);
155
156 // copy the diagonal from rDiag and the upper triangle of qr
157 for (int row = Math.min(m, n) - 1; row >= 0; row--) {
158 cachedR.setEntry(row, row, rDiag[row]);
159 for (int col = row + 1; col < n; col++) {
160 cachedR.setEntry(row, col, qrt[col][row]);
161 }
162 }
163
164 }
165
166 // return the cached matrix
167 return cachedR;
168
169 }
170
171 /** {@inheritDoc} */
172 public RealMatrix getQ() {
173 if (cachedQ == null) {
174 cachedQ = getQT().transpose();
175 }
176 return cachedQ;
177 }
178
179 /** {@inheritDoc} */
180 public RealMatrix getQT() {
181
182 if (cachedQT == null) {
183
184 // QT is supposed to be m x m
185 final int n = qrt.length;
186 final int m = qrt[0].length;
187 cachedQT = MatrixUtils.createRealMatrix(m, m);
188
189 /*
190 * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
191 * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
192 * succession to the result
193 */
194 for (int minor = m - 1; minor >= Math.min(m, n); minor--) {
195 cachedQT.setEntry(minor, minor, 1.0);
196 }
197
198 for (int minor = Math.min(m, n)-1; minor >= 0; minor--){
199 final double[] qrtMinor = qrt[minor];
200 cachedQT.setEntry(minor, minor, 1.0);
201 if (qrtMinor[minor] != 0.0) {
202 for (int col = minor; col < m; col++) {
203 double alpha = 0;
204 for (int row = minor; row < m; row++) {
205 alpha -= cachedQT.getEntry(col, row) * qrtMinor[row];
206 }
207 alpha /= rDiag[minor] * qrtMinor[minor];
208
209 for (int row = minor; row < m; row++) {
210 cachedQT.addToEntry(col, row, -alpha * qrtMinor[row]);
211 }
212 }
213 }
214 }
215
216 }
217
218 // return the cached matrix
219 return cachedQT;
220
221 }
222
223 /** {@inheritDoc} */
224 public RealMatrix getH() {
225
226 if (cachedH == null) {
227
228 final int n = qrt.length;
229 final int m = qrt[0].length;
230 cachedH = MatrixUtils.createRealMatrix(m, n);
231 for (int i = 0; i < m; ++i) {
232 for (int j = 0; j < Math.min(i + 1, n); ++j) {
233 cachedH.setEntry(i, j, qrt[j][i] / -rDiag[j]);
234 }
235 }
236
237 }
238
239 // return the cached matrix
240 return cachedH;
241
242 }
243
244 /** {@inheritDoc} */
245 public DecompositionSolver getSolver() {
246 return new Solver(qrt, rDiag);
247 }
248
249 /** Specialized solver. */
250 private static class Solver implements DecompositionSolver {
251
252 /**
253 * A packed TRANSPOSED representation of the QR decomposition.
254 * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
255 * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
256 * from which an explicit form of Q can be recomputed if desired.</p>
257 */
258 private final double[][] qrt;
259
260 /** The diagonal elements of R. */
261 private final double[] rDiag;
262
263 /**
264 * Build a solver from decomposed matrix.
265 * @param qrt packed TRANSPOSED representation of the QR decomposition
266 * @param rDiag diagonal elements of R
267 */
268 private Solver(final double[][] qrt, final double[] rDiag) {
269 this.qrt = qrt;
270 this.rDiag = rDiag;
271 }
272
273 /** {@inheritDoc} */
274 public boolean isNonSingular() {
275
276 for (double diag : rDiag) {
277 if (diag == 0) {
278 return false;
279 }
280 }
281 return true;
282
283 }
284
285 /** {@inheritDoc} */
286 public double[] solve(double[] b)
287 throws IllegalArgumentException, InvalidMatrixException {
288
289 final int n = qrt.length;
290 final int m = qrt[0].length;
291 if (b.length != m) {
292 throw MathRuntimeException.createIllegalArgumentException(
293 "vector length mismatch: got {0} but expected {1}",
294 b.length, m);
295 }
296 if (!isNonSingular()) {
297 throw new SingularMatrixException();
298 }
299
300 final double[] x = new double[n];
301 final double[] y = b.clone();
302
303 // apply Householder transforms to solve Q.y = b
304 for (int minor = 0; minor < Math.min(m, n); minor++) {
305
306 final double[] qrtMinor = qrt[minor];
307 double dotProduct = 0;
308 for (int row = minor; row < m; row++) {
309 dotProduct += y[row] * qrtMinor[row];
310 }
311 dotProduct /= rDiag[minor] * qrtMinor[minor];
312
313 for (int row = minor; row < m; row++) {
314 y[row] += dotProduct * qrtMinor[row];
315 }
316
317 }
318
319 // solve triangular system R.x = y
320 for (int row = rDiag.length - 1; row >= 0; --row) {
321 y[row] /= rDiag[row];
322 final double yRow = y[row];
323 final double[] qrtRow = qrt[row];
324 x[row] = yRow;
325 for (int i = 0; i < row; i++) {
326 y[i] -= yRow * qrtRow[i];
327 }
328 }
329
330 return x;
331
332 }
333
334 /** {@inheritDoc} */
335 public RealVector solve(RealVector b)
336 throws IllegalArgumentException, InvalidMatrixException {
337 try {
338 return solve((ArrayRealVector) b);
339 } catch (ClassCastException cce) {
340 return new ArrayRealVector(solve(b.getData()), false);
341 }
342 }
343
344 /** Solve the linear equation A × X = B.
345 * <p>The A matrix is implicit here. It is </p>
346 * @param b right-hand side of the equation A × X = B
347 * @return a vector X that minimizes the two norm of A × X - B
348 * @throws IllegalArgumentException if matrices dimensions don't match
349 * @throws InvalidMatrixException if decomposed matrix is singular
350 */
351 public ArrayRealVector solve(ArrayRealVector b)
352 throws IllegalArgumentException, InvalidMatrixException {
353 return new ArrayRealVector(solve(b.getDataRef()), false);
354 }
355
356 /** {@inheritDoc} */
357 public RealMatrix solve(RealMatrix b)
358 throws IllegalArgumentException, InvalidMatrixException {
359
360 final int n = qrt.length;
361 final int m = qrt[0].length;
362 if (b.getRowDimension() != m) {
363 throw MathRuntimeException.createIllegalArgumentException(
364 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
365 b.getRowDimension(), b.getColumnDimension(), m, "n");
366 }
367 if (!isNonSingular()) {
368 throw new SingularMatrixException();
369 }
370
371 final int columns = b.getColumnDimension();
372 final int blockSize = BlockRealMatrix.BLOCK_SIZE;
373 final int cBlocks = (columns + blockSize - 1) / blockSize;
374 final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns);
375 final double[][] y = new double[b.getRowDimension()][blockSize];
376 final double[] alpha = new double[blockSize];
377
378 for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
379 final int kStart = kBlock * blockSize;
380 final int kEnd = Math.min(kStart + blockSize, columns);
381 final int kWidth = kEnd - kStart;
382
383 // get the right hand side vector
384 b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y);
385
386 // apply Householder transforms to solve Q.y = b
387 for (int minor = 0; minor < Math.min(m, n); minor++) {
388 final double[] qrtMinor = qrt[minor];
389 final double factor = 1.0 / (rDiag[minor] * qrtMinor[minor]);
390
391 Arrays.fill(alpha, 0, kWidth, 0.0);
392 for (int row = minor; row < m; ++row) {
393 final double d = qrtMinor[row];
394 final double[] yRow = y[row];
395 for (int k = 0; k < kWidth; ++k) {
396 alpha[k] += d * yRow[k];
397 }
398 }
399 for (int k = 0; k < kWidth; ++k) {
400 alpha[k] *= factor;
401 }
402
403 for (int row = minor; row < m; ++row) {
404 final double d = qrtMinor[row];
405 final double[] yRow = y[row];
406 for (int k = 0; k < kWidth; ++k) {
407 yRow[k] += alpha[k] * d;
408 }
409 }
410
411 }
412
413 // solve triangular system R.x = y
414 for (int j = rDiag.length - 1; j >= 0; --j) {
415 final int jBlock = j / blockSize;
416 final int jStart = jBlock * blockSize;
417 final double factor = 1.0 / rDiag[j];
418 final double[] yJ = y[j];
419 final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
420 int index = (j - jStart) * kWidth;
421 for (int k = 0; k < kWidth; ++k) {
422 yJ[k] *= factor;
423 xBlock[index++] = yJ[k];
424 }
425
426 final double[] qrtJ = qrt[j];
427 for (int i = 0; i < j; ++i) {
428 final double rIJ = qrtJ[i];
429 final double[] yI = y[i];
430 for (int k = 0; k < kWidth; ++k) {
431 yI[k] -= yJ[k] * rIJ;
432 }
433 }
434
435 }
436
437 }
438
439 return new BlockRealMatrix(n, columns, xBlocks, false);
440
441 }
442
443 /** {@inheritDoc} */
444 public RealMatrix getInverse()
445 throws InvalidMatrixException {
446 return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
447 }
448
449 }
450
451 }