    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    }