View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math3.linear;
19  
20  import java.util.Arrays;
21  
22  import org.apache.commons.math3.exception.DimensionMismatchException;
23  import org.apache.commons.math3.util.FastMath;
24  
25  
26  /**
27   * Calculates the QR-decomposition of a matrix.
28   * <p>The QR-decomposition of a matrix A consists of two matrices Q and R
29   * that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is
30   * upper triangular. If A is m&times;n, Q is m&times;m and R m&times;n.</p>
31   * <p>This class compute the decomposition using Householder reflectors.</p>
32   * <p>For efficiency purposes, the decomposition in packed form is transposed.
33   * This allows inner loop to iterate inside rows, which is much more cache-efficient
34   * in Java.</p>
35   * <p>This class is based on the class with similar name from the
36   * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
37   * following changes:</p>
38   * <ul>
39   *   <li>a {@link #getQT() getQT} method has been added,</li>
40   *   <li>the {@code solve} and {@code isFullRank} methods have been replaced
41   *   by a {@link #getSolver() getSolver} method and the equivalent methods
42   *   provided by the returned {@link DecompositionSolver}.</li>
43   * </ul>
44   *
45   * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
46   * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
47   *
48   * @version $Id: QRDecomposition.java 1462423 2013-03-29 07:25:18Z luc $
49   * @since 1.2 (changed to concrete class in 3.0)
50   */
51  public class QRDecomposition {
52      /**
53       * A packed TRANSPOSED representation of the QR decomposition.
54       * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
55       * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
56       * from which an explicit form of Q can be recomputed if desired.</p>
57       */
58      private double[][] qrt;
59      /** The diagonal elements of R. */
60      private double[] rDiag;
61      /** Cached value of Q. */
62      private RealMatrix cachedQ;
63      /** Cached value of QT. */
64      private RealMatrix cachedQT;
65      /** Cached value of R. */
66      private RealMatrix cachedR;
67      /** Cached value of H. */
68      private RealMatrix cachedH;
69      /** Singularity threshold. */
70      private final double threshold;
71  
72      /**
73       * Calculates the QR-decomposition of the given matrix.
74       * The singularity threshold defaults to zero.
75       *
76       * @param matrix The matrix to decompose.
77       *
78       * @see #QRDecomposition(RealMatrix,double)
79       */
80      public QRDecomposition(RealMatrix matrix) {
81          this(matrix, 0d);
82      }
83  
84      /**
85       * Calculates the QR-decomposition of the given matrix.
86       *
87       * @param matrix The matrix to decompose.
88       * @param threshold Singularity threshold.
89       */
90      public QRDecomposition(RealMatrix matrix,
91                             double threshold) {
92          this.threshold = threshold;
93  
94          final int m = matrix.getRowDimension();
95          final int n = matrix.getColumnDimension();
96          qrt = matrix.transpose().getData();
97          rDiag = new double[FastMath.min(m, n)];
98          cachedQ  = null;
99          cachedQT = null;
100         cachedR  = null;
101         cachedH  = null;
102 
103         decompose(qrt);
104 
105     }
106 
107     /** Decompose matrix.
108      * @param matrix transposed matrix
109      * @since 3.2
110      */
111     protected void decompose(double[][] matrix) {
112         for (int minor = 0; minor < FastMath.min(qrt.length, qrt[0].length); minor++) {
113             performHouseholderReflection(minor, qrt);
114         }
115     }
116 
117     /** Perform Householder reflection for a minor A(minor, minor) of A.
118      * @param minor minor index
119      * @param matrix transposed matrix
120      * @since 3.2
121      */
122     protected void performHouseholderReflection(int minor, double[][] matrix) {
123 
124         final double[] qrtMinor = qrt[minor];
125 
126         /*
127          * Let x be the first column of the minor, and a^2 = |x|^2.
128          * x will be in the positions qr[minor][minor] through qr[m][minor].
129          * The first column of the transformed minor will be (a,0,0,..)'
130          * The sign of a is chosen to be opposite to the sign of the first
131          * component of x. Let's find a:
132          */
133         double xNormSqr = 0;
134         for (int row = minor; row < qrtMinor.length; row++) {
135             final double c = qrtMinor[row];
136             xNormSqr += c * c;
137         }
138         final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr);
139         rDiag[minor] = a;
140 
141         if (a != 0.0) {
142 
143             /*
144              * Calculate the normalized reflection vector v and transform
145              * the first column. We know the norm of v beforehand: v = x-ae
146              * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
147              * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
148              * Here <x, e> is now qr[minor][minor].
149              * v = x-ae is stored in the column at qr:
150              */
151             qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
152 
153             /*
154              * Transform the rest of the columns of the minor:
155              * They will be transformed by the matrix H = I-2vv'/|v|^2.
156              * If x is a column vector of the minor, then
157              * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
158              * Therefore the transformation is easily calculated by
159              * subtracting the column vector (2<x,v>/|v|^2)v from x.
160              *
161              * Let 2<x,v>/|v|^2 = alpha. From above we have
162              * |v|^2 = -2a*(qr[minor][minor]), so
163              * alpha = -<x,v>/(a*qr[minor][minor])
164              */
165             for (int col = minor+1; col < qrt.length; col++) {
166                 final double[] qrtCol = qrt[col];
167                 double alpha = 0;
168                 for (int row = minor; row < qrtCol.length; row++) {
169                     alpha -= qrtCol[row] * qrtMinor[row];
170                 }
171                 alpha /= a * qrtMinor[minor];
172 
173                 // Subtract the column vector alpha*v from x.
174                 for (int row = minor; row < qrtCol.length; row++) {
175                     qrtCol[row] -= alpha * qrtMinor[row];
176                 }
177             }
178         }
179     }
180 
181 
182     /**
183      * Returns the matrix R of the decomposition.
184      * <p>R is an upper-triangular matrix</p>
185      * @return the R matrix
186      */
187     public RealMatrix getR() {
188 
189         if (cachedR == null) {
190 
191             // R is supposed to be m x n
192             final int n = qrt.length;
193             final int m = qrt[0].length;
194             double[][] ra = new double[m][n];
195             // copy the diagonal from rDiag and the upper triangle of qr
196             for (int row = FastMath.min(m, n) - 1; row >= 0; row--) {
197                 ra[row][row] = rDiag[row];
198                 for (int col = row + 1; col < n; col++) {
199                     ra[row][col] = qrt[col][row];
200                 }
201             }
202             cachedR = MatrixUtils.createRealMatrix(ra);
203         }
204 
205         // return the cached matrix
206         return cachedR;
207     }
208 
209     /**
210      * Returns the matrix Q of the decomposition.
211      * <p>Q is an orthogonal matrix</p>
212      * @return the Q matrix
213      */
214     public RealMatrix getQ() {
215         if (cachedQ == null) {
216             cachedQ = getQT().transpose();
217         }
218         return cachedQ;
219     }
220 
221     /**
222      * Returns the transpose of the matrix Q of the decomposition.
223      * <p>Q is an orthogonal matrix</p>
224      * @return the transpose of the Q matrix, Q<sup>T</sup>
225      */
226     public RealMatrix getQT() {
227         if (cachedQT == null) {
228 
229             // QT is supposed to be m x m
230             final int n = qrt.length;
231             final int m = qrt[0].length;
232             double[][] qta = new double[m][m];
233 
234             /*
235              * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
236              * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
237              * succession to the result
238              */
239             for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) {
240                 qta[minor][minor] = 1.0d;
241             }
242 
243             for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){
244                 final double[] qrtMinor = qrt[minor];
245                 qta[minor][minor] = 1.0d;
246                 if (qrtMinor[minor] != 0.0) {
247                     for (int col = minor; col < m; col++) {
248                         double alpha = 0;
249                         for (int row = minor; row < m; row++) {
250                             alpha -= qta[col][row] * qrtMinor[row];
251                         }
252                         alpha /= rDiag[minor] * qrtMinor[minor];
253 
254                         for (int row = minor; row < m; row++) {
255                             qta[col][row] += -alpha * qrtMinor[row];
256                         }
257                     }
258                 }
259             }
260             cachedQT = MatrixUtils.createRealMatrix(qta);
261         }
262 
263         // return the cached matrix
264         return cachedQT;
265     }
266 
267     /**
268      * Returns the Householder reflector vectors.
269      * <p>H is a lower trapezoidal matrix whose columns represent
270      * each successive Householder reflector vector. This matrix is used
271      * to compute Q.</p>
272      * @return a matrix containing the Householder reflector vectors
273      */
274     public RealMatrix getH() {
275         if (cachedH == null) {
276 
277             final int n = qrt.length;
278             final int m = qrt[0].length;
279             double[][] ha = new double[m][n];
280             for (int i = 0; i < m; ++i) {
281                 for (int j = 0; j < FastMath.min(i + 1, n); ++j) {
282                     ha[i][j] = qrt[j][i] / -rDiag[j];
283                 }
284             }
285             cachedH = MatrixUtils.createRealMatrix(ha);
286         }
287 
288         // return the cached matrix
289         return cachedH;
290     }
291 
292     /**
293      * Get a solver for finding the A &times; X = B solution in least square sense.
294      * @return a solver
295      */
296     public DecompositionSolver getSolver() {
297         return new Solver(qrt, rDiag, threshold);
298     }
299 
300     /** Specialized solver. */
301     private static class Solver implements DecompositionSolver {
302         /**
303          * A packed TRANSPOSED representation of the QR decomposition.
304          * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
305          * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
306          * from which an explicit form of Q can be recomputed if desired.</p>
307          */
308         private final double[][] qrt;
309         /** The diagonal elements of R. */
310         private final double[] rDiag;
311         /** Singularity threshold. */
312         private final double threshold;
313 
314         /**
315          * Build a solver from decomposed matrix.
316          *
317          * @param qrt Packed TRANSPOSED representation of the QR decomposition.
318          * @param rDiag Diagonal elements of R.
319          * @param threshold Singularity threshold.
320          */
321         private Solver(final double[][] qrt,
322                        final double[] rDiag,
323                        final double threshold) {
324             this.qrt   = qrt;
325             this.rDiag = rDiag;
326             this.threshold = threshold;
327         }
328 
329         /** {@inheritDoc} */
330         public boolean isNonSingular() {
331             for (double diag : rDiag) {
332                 if (FastMath.abs(diag) <= threshold) {
333                     return false;
334                 }
335             }
336             return true;
337         }
338 
339         /** {@inheritDoc} */
340         public RealVector solve(RealVector b) {
341             final int n = qrt.length;
342             final int m = qrt[0].length;
343             if (b.getDimension() != m) {
344                 throw new DimensionMismatchException(b.getDimension(), m);
345             }
346             if (!isNonSingular()) {
347                 throw new SingularMatrixException();
348             }
349 
350             final double[] x = new double[n];
351             final double[] y = b.toArray();
352 
353             // apply Householder transforms to solve Q.y = b
354             for (int minor = 0; minor < FastMath.min(m, n); minor++) {
355 
356                 final double[] qrtMinor = qrt[minor];
357                 double dotProduct = 0;
358                 for (int row = minor; row < m; row++) {
359                     dotProduct += y[row] * qrtMinor[row];
360                 }
361                 dotProduct /= rDiag[minor] * qrtMinor[minor];
362 
363                 for (int row = minor; row < m; row++) {
364                     y[row] += dotProduct * qrtMinor[row];
365                 }
366             }
367 
368             // solve triangular system R.x = y
369             for (int row = rDiag.length - 1; row >= 0; --row) {
370                 y[row] /= rDiag[row];
371                 final double yRow = y[row];
372                 final double[] qrtRow = qrt[row];
373                 x[row] = yRow;
374                 for (int i = 0; i < row; i++) {
375                     y[i] -= yRow * qrtRow[i];
376                 }
377             }
378 
379             return new ArrayRealVector(x, false);
380         }
381 
382         /** {@inheritDoc} */
383         public RealMatrix solve(RealMatrix b) {
384             final int n = qrt.length;
385             final int m = qrt[0].length;
386             if (b.getRowDimension() != m) {
387                 throw new DimensionMismatchException(b.getRowDimension(), m);
388             }
389             if (!isNonSingular()) {
390                 throw new SingularMatrixException();
391             }
392 
393             final int columns        = b.getColumnDimension();
394             final int blockSize      = BlockRealMatrix.BLOCK_SIZE;
395             final int cBlocks        = (columns + blockSize - 1) / blockSize;
396             final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns);
397             final double[][] y       = new double[b.getRowDimension()][blockSize];
398             final double[]   alpha   = new double[blockSize];
399 
400             for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
401                 final int kStart = kBlock * blockSize;
402                 final int kEnd   = FastMath.min(kStart + blockSize, columns);
403                 final int kWidth = kEnd - kStart;
404 
405                 // get the right hand side vector
406                 b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y);
407 
408                 // apply Householder transforms to solve Q.y = b
409                 for (int minor = 0; minor < FastMath.min(m, n); minor++) {
410                     final double[] qrtMinor = qrt[minor];
411                     final double factor     = 1.0 / (rDiag[minor] * qrtMinor[minor]);
412 
413                     Arrays.fill(alpha, 0, kWidth, 0.0);
414                     for (int row = minor; row < m; ++row) {
415                         final double   d    = qrtMinor[row];
416                         final double[] yRow = y[row];
417                         for (int k = 0; k < kWidth; ++k) {
418                             alpha[k] += d * yRow[k];
419                         }
420                     }
421                     for (int k = 0; k < kWidth; ++k) {
422                         alpha[k] *= factor;
423                     }
424 
425                     for (int row = minor; row < m; ++row) {
426                         final double   d    = qrtMinor[row];
427                         final double[] yRow = y[row];
428                         for (int k = 0; k < kWidth; ++k) {
429                             yRow[k] += alpha[k] * d;
430                         }
431                     }
432                 }
433 
434                 // solve triangular system R.x = y
435                 for (int j = rDiag.length - 1; j >= 0; --j) {
436                     final int      jBlock = j / blockSize;
437                     final int      jStart = jBlock * blockSize;
438                     final double   factor = 1.0 / rDiag[j];
439                     final double[] yJ     = y[j];
440                     final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
441                     int index = (j - jStart) * kWidth;
442                     for (int k = 0; k < kWidth; ++k) {
443                         yJ[k]          *= factor;
444                         xBlock[index++] = yJ[k];
445                     }
446 
447                     final double[] qrtJ = qrt[j];
448                     for (int i = 0; i < j; ++i) {
449                         final double rIJ  = qrtJ[i];
450                         final double[] yI = y[i];
451                         for (int k = 0; k < kWidth; ++k) {
452                             yI[k] -= yJ[k] * rIJ;
453                         }
454                     }
455                 }
456             }
457 
458             return new BlockRealMatrix(n, columns, xBlocks, false);
459         }
460 
461         /** {@inheritDoc} */
462         public RealMatrix getInverse() {
463             return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
464         }
465     }
466 }