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.math3.linear;
019
020 import java.util.Arrays;
021
022 import org.apache.commons.math3.exception.DimensionMismatchException;
023 import org.apache.commons.math3.util.FastMath;
024
025
026 /**
027 * Calculates the QR-decomposition of a matrix.
028 * <p>The QR-decomposition of a matrix A consists of two matrices Q and R
029 * that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is
030 * upper triangular. If A is m×n, Q is m×m and R m×n.</p>
031 * <p>This class compute the decomposition using Householder reflectors.</p>
032 * <p>For efficiency purposes, the decomposition in packed form is transposed.
033 * This allows inner loop to iterate inside rows, which is much more cache-efficient
034 * in Java.</p>
035 * <p>This class is based on the class with similar name from the
036 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
037 * following changes:</p>
038 * <ul>
039 * <li>a {@link #getQT() getQT} method has been added,</li>
040 * <li>the {@code solve} and {@code isFullRank} methods have been replaced
041 * by a {@link #getSolver() getSolver} method and the equivalent methods
042 * provided by the returned {@link DecompositionSolver}.</li>
043 * </ul>
044 *
045 * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
046 * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
047 *
048 * @version $Id: QRDecomposition.java 1416643 2012-12-03 19:37:14Z tn $
049 * @since 1.2 (changed to concrete class in 3.0)
050 */
051 public class QRDecomposition {
052 /**
053 * A packed TRANSPOSED representation of the QR decomposition.
054 * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
055 * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
056 * from which an explicit form of Q can be recomputed if desired.</p>
057 */
058 private double[][] qrt;
059 /** The diagonal elements of R. */
060 private double[] rDiag;
061 /** Cached value of Q. */
062 private RealMatrix cachedQ;
063 /** Cached value of QT. */
064 private RealMatrix cachedQT;
065 /** Cached value of R. */
066 private RealMatrix cachedR;
067 /** Cached value of H. */
068 private RealMatrix cachedH;
069 /** Singularity threshold. */
070 private final double threshold;
071
072 /**
073 * Calculates the QR-decomposition of the given matrix.
074 * The singularity threshold defaults to zero.
075 *
076 * @param matrix The matrix to decompose.
077 *
078 * @see #QRDecomposition(RealMatrix,double)
079 */
080 public QRDecomposition(RealMatrix matrix) {
081 this(matrix, 0d);
082 }
083
084 /**
085 * Calculates the QR-decomposition of the given matrix.
086 *
087 * @param matrix The matrix to decompose.
088 * @param threshold Singularity threshold.
089 */
090 public QRDecomposition(RealMatrix matrix,
091 double threshold) {
092 this.threshold = threshold;
093
094 final int m = matrix.getRowDimension();
095 final int n = matrix.getColumnDimension();
096 qrt = matrix.transpose().getData();
097 rDiag = new double[FastMath.min(m, n)];
098 cachedQ = null;
099 cachedQT = null;
100 cachedR = null;
101 cachedH = null;
102
103 /*
104 * The QR decomposition of a matrix A is calculated using Householder
105 * reflectors by repeating the following operations to each minor
106 * A(minor,minor) of A:
107 */
108 for (int minor = 0; minor < FastMath.min(m, n); minor++) {
109
110 final double[] qrtMinor = qrt[minor];
111
112 /*
113 * Let x be the first column of the minor, and a^2 = |x|^2.
114 * x will be in the positions qr[minor][minor] through qr[m][minor].
115 * The first column of the transformed minor will be (a,0,0,..)'
116 * The sign of a is chosen to be opposite to the sign of the first
117 * component of x. Let's find a:
118 */
119 double xNormSqr = 0;
120 for (int row = minor; row < m; row++) {
121 final double c = qrtMinor[row];
122 xNormSqr += c * c;
123 }
124 final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
125 rDiag[minor] = a;
126
127 if (a != 0.0) {
128
129 /*
130 * Calculate the normalized reflection vector v and transform
131 * the first column. We know the norm of v beforehand: v = x-ae
132 * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
133 * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
134 * Here <x, e> is now qr[minor][minor].
135 * v = x-ae is stored in the column at qr:
136 */
137 qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
138
139 /*
140 * Transform the rest of the columns of the minor:
141 * They will be transformed by the matrix H = I-2vv'/|v|^2.
142 * If x is a column vector of the minor, then
143 * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
144 * Therefore the transformation is easily calculated by
145 * subtracting the column vector (2<x,v>/|v|^2)v from x.
146 *
147 * Let 2<x,v>/|v|^2 = alpha. From above we have
148 * |v|^2 = -2a*(qr[minor][minor]), so
149 * alpha = -<x,v>/(a*qr[minor][minor])
150 */
151 for (int col = minor+1; col < n; col++) {
152 final double[] qrtCol = qrt[col];
153 double alpha = 0;
154 for (int row = minor; row < m; row++) {
155 alpha -= qrtCol[row] * qrtMinor[row];
156 }
157 alpha /= a * qrtMinor[minor];
158
159 // Subtract the column vector alpha*v from x.
160 for (int row = minor; row < m; row++) {
161 qrtCol[row] -= alpha * qrtMinor[row];
162 }
163 }
164 }
165 }
166 }
167
168 /**
169 * Returns the matrix R of the decomposition.
170 * <p>R is an upper-triangular matrix</p>
171 * @return the R matrix
172 */
173 public RealMatrix getR() {
174
175 if (cachedR == null) {
176
177 // R is supposed to be m x n
178 final int n = qrt.length;
179 final int m = qrt[0].length;
180 double[][] ra = new double[m][n];
181 // copy the diagonal from rDiag and the upper triangle of qr
182 for (int row = FastMath.min(m, n) - 1; row >= 0; row--) {
183 ra[row][row] = rDiag[row];
184 for (int col = row + 1; col < n; col++) {
185 ra[row][col] = qrt[col][row];
186 }
187 }
188 cachedR = MatrixUtils.createRealMatrix(ra);
189 }
190
191 // return the cached matrix
192 return cachedR;
193 }
194
195 /**
196 * Returns the matrix Q of the decomposition.
197 * <p>Q is an orthogonal matrix</p>
198 * @return the Q matrix
199 */
200 public RealMatrix getQ() {
201 if (cachedQ == null) {
202 cachedQ = getQT().transpose();
203 }
204 return cachedQ;
205 }
206
207 /**
208 * Returns the transpose of the matrix Q of the decomposition.
209 * <p>Q is an orthogonal matrix</p>
210 * @return the transpose of the Q matrix, Q<sup>T</sup>
211 */
212 public RealMatrix getQT() {
213 if (cachedQT == null) {
214
215 // QT is supposed to be m x m
216 final int n = qrt.length;
217 final int m = qrt[0].length;
218 double[][] qta = new double[m][m];
219
220 /*
221 * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
222 * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
223 * succession to the result
224 */
225 for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) {
226 qta[minor][minor] = 1.0d;
227 }
228
229 for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){
230 final double[] qrtMinor = qrt[minor];
231 qta[minor][minor] = 1.0d;
232 if (qrtMinor[minor] != 0.0) {
233 for (int col = minor; col < m; col++) {
234 double alpha = 0;
235 for (int row = minor; row < m; row++) {
236 alpha -= qta[col][row] * qrtMinor[row];
237 }
238 alpha /= rDiag[minor] * qrtMinor[minor];
239
240 for (int row = minor; row < m; row++) {
241 qta[col][row] += -alpha * qrtMinor[row];
242 }
243 }
244 }
245 }
246 cachedQT = MatrixUtils.createRealMatrix(qta);
247 }
248
249 // return the cached matrix
250 return cachedQT;
251 }
252
253 /**
254 * Returns the Householder reflector vectors.
255 * <p>H is a lower trapezoidal matrix whose columns represent
256 * each successive Householder reflector vector. This matrix is used
257 * to compute Q.</p>
258 * @return a matrix containing the Householder reflector vectors
259 */
260 public RealMatrix getH() {
261 if (cachedH == null) {
262
263 final int n = qrt.length;
264 final int m = qrt[0].length;
265 double[][] ha = new double[m][n];
266 for (int i = 0; i < m; ++i) {
267 for (int j = 0; j < FastMath.min(i + 1, n); ++j) {
268 ha[i][j] = qrt[j][i] / -rDiag[j];
269 }
270 }
271 cachedH = MatrixUtils.createRealMatrix(ha);
272 }
273
274 // return the cached matrix
275 return cachedH;
276 }
277
278 /**
279 * Get a solver for finding the A × X = B solution in least square sense.
280 * @return a solver
281 */
282 public DecompositionSolver getSolver() {
283 return new Solver(qrt, rDiag, threshold);
284 }
285
286 /** Specialized solver. */
287 private static class Solver implements DecompositionSolver {
288 /**
289 * A packed TRANSPOSED representation of the QR decomposition.
290 * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
291 * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
292 * from which an explicit form of Q can be recomputed if desired.</p>
293 */
294 private final double[][] qrt;
295 /** The diagonal elements of R. */
296 private final double[] rDiag;
297 /** Singularity threshold. */
298 private final double threshold;
299
300 /**
301 * Build a solver from decomposed matrix.
302 *
303 * @param qrt Packed TRANSPOSED representation of the QR decomposition.
304 * @param rDiag Diagonal elements of R.
305 * @param threshold Singularity threshold.
306 */
307 private Solver(final double[][] qrt,
308 final double[] rDiag,
309 final double threshold) {
310 this.qrt = qrt;
311 this.rDiag = rDiag;
312 this.threshold = threshold;
313 }
314
315 /** {@inheritDoc} */
316 public boolean isNonSingular() {
317 for (double diag : rDiag) {
318 if (FastMath.abs(diag) <= threshold) {
319 return false;
320 }
321 }
322 return true;
323 }
324
325 /** {@inheritDoc} */
326 public RealVector solve(RealVector b) {
327 final int n = qrt.length;
328 final int m = qrt[0].length;
329 if (b.getDimension() != m) {
330 throw new DimensionMismatchException(b.getDimension(), m);
331 }
332 if (!isNonSingular()) {
333 throw new SingularMatrixException();
334 }
335
336 final double[] x = new double[n];
337 final double[] y = b.toArray();
338
339 // apply Householder transforms to solve Q.y = b
340 for (int minor = 0; minor < FastMath.min(m, n); minor++) {
341
342 final double[] qrtMinor = qrt[minor];
343 double dotProduct = 0;
344 for (int row = minor; row < m; row++) {
345 dotProduct += y[row] * qrtMinor[row];
346 }
347 dotProduct /= rDiag[minor] * qrtMinor[minor];
348
349 for (int row = minor; row < m; row++) {
350 y[row] += dotProduct * qrtMinor[row];
351 }
352 }
353
354 // solve triangular system R.x = y
355 for (int row = rDiag.length - 1; row >= 0; --row) {
356 y[row] /= rDiag[row];
357 final double yRow = y[row];
358 final double[] qrtRow = qrt[row];
359 x[row] = yRow;
360 for (int i = 0; i < row; i++) {
361 y[i] -= yRow * qrtRow[i];
362 }
363 }
364
365 return new ArrayRealVector(x, false);
366 }
367
368 /** {@inheritDoc} */
369 public RealMatrix solve(RealMatrix b) {
370 final int n = qrt.length;
371 final int m = qrt[0].length;
372 if (b.getRowDimension() != m) {
373 throw new DimensionMismatchException(b.getRowDimension(), m);
374 }
375 if (!isNonSingular()) {
376 throw new SingularMatrixException();
377 }
378
379 final int columns = b.getColumnDimension();
380 final int blockSize = BlockRealMatrix.BLOCK_SIZE;
381 final int cBlocks = (columns + blockSize - 1) / blockSize;
382 final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns);
383 final double[][] y = new double[b.getRowDimension()][blockSize];
384 final double[] alpha = new double[blockSize];
385
386 for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
387 final int kStart = kBlock * blockSize;
388 final int kEnd = FastMath.min(kStart + blockSize, columns);
389 final int kWidth = kEnd - kStart;
390
391 // get the right hand side vector
392 b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y);
393
394 // apply Householder transforms to solve Q.y = b
395 for (int minor = 0; minor < FastMath.min(m, n); minor++) {
396 final double[] qrtMinor = qrt[minor];
397 final double factor = 1.0 / (rDiag[minor] * qrtMinor[minor]);
398
399 Arrays.fill(alpha, 0, kWidth, 0.0);
400 for (int row = minor; row < m; ++row) {
401 final double d = qrtMinor[row];
402 final double[] yRow = y[row];
403 for (int k = 0; k < kWidth; ++k) {
404 alpha[k] += d * yRow[k];
405 }
406 }
407 for (int k = 0; k < kWidth; ++k) {
408 alpha[k] *= factor;
409 }
410
411 for (int row = minor; row < m; ++row) {
412 final double d = qrtMinor[row];
413 final double[] yRow = y[row];
414 for (int k = 0; k < kWidth; ++k) {
415 yRow[k] += alpha[k] * d;
416 }
417 }
418 }
419
420 // solve triangular system R.x = y
421 for (int j = rDiag.length - 1; j >= 0; --j) {
422 final int jBlock = j / blockSize;
423 final int jStart = jBlock * blockSize;
424 final double factor = 1.0 / rDiag[j];
425 final double[] yJ = y[j];
426 final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
427 int index = (j - jStart) * kWidth;
428 for (int k = 0; k < kWidth; ++k) {
429 yJ[k] *= factor;
430 xBlock[index++] = yJ[k];
431 }
432
433 final double[] qrtJ = qrt[j];
434 for (int i = 0; i < j; ++i) {
435 final double rIJ = qrtJ[i];
436 final double[] yI = y[i];
437 for (int k = 0; k < kWidth; ++k) {
438 yI[k] -= yJ[k] * rIJ;
439 }
440 }
441 }
442 }
443
444 return new BlockRealMatrix(n, columns, xBlocks, false);
445 }
446
447 /** {@inheritDoc} */
448 public RealMatrix getInverse() {
449 return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
450 }
451 }
452 }