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&times;n, Q is m&times;m and R m&times;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 &times; 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    }