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