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
018package org.apache.commons.math4.linear;
019
020import org.apache.commons.math4.util.FastMath;
021
022
023/**
024 * Calculates the rank-revealing QR-decomposition of a matrix, with column pivoting.
025 * <p>The rank-revealing QR-decomposition of a matrix A consists of three matrices Q,
026 * R and P such that AP=QR.  Q is orthogonal (Q<sup>T</sup>Q = I), and R is upper triangular.
027 * If A is m&times;n, Q is m&times;m and R is m&times;n and P is n&times;n.</p>
028 * <p>QR decomposition with column pivoting produces a rank-revealing QR
029 * decomposition and the {@link #getRank(double)} method may be used to return the rank of the
030 * input matrix A.</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 * @since 3.2
049 */
050public class RRQRDecomposition extends QRDecomposition {
051
052    /** An array to record the column pivoting for later creation of P. */
053    private int[] p;
054
055    /** Cached value of P. */
056    private RealMatrix cachedP;
057
058
059    /**
060     * Calculates the QR-decomposition of the given matrix.
061     * The singularity threshold defaults to zero.
062     *
063     * @param matrix The matrix to decompose.
064     *
065     * @see #RRQRDecomposition(RealMatrix, double)
066     */
067    public RRQRDecomposition(RealMatrix matrix) {
068        this(matrix, 0d);
069    }
070
071    /**
072     * Calculates the QR-decomposition of the given matrix.
073     *
074     * @param matrix The matrix to decompose.
075     * @param threshold Singularity threshold.
076     * @see #RRQRDecomposition(RealMatrix)
077     */
078    public RRQRDecomposition(RealMatrix matrix,  double threshold) {
079        super(matrix, threshold);
080    }
081
082    /** Decompose matrix.
083     * @param qrt transposed matrix
084     */
085    @Override
086    protected void decompose(double[][] qrt) {
087        p = new int[qrt.length];
088        for (int i = 0; i < p.length; i++) {
089            p[i] = i;
090        }
091        super.decompose(qrt);
092    }
093
094    /** Perform Householder reflection for a minor A(minor, minor) of A.
095     *
096     * @param minor minor index
097     * @param qrt transposed matrix
098     */
099    @Override
100    protected void performHouseholderReflection(int minor, double[][] qrt) {
101        double l2NormSquaredMax = 0;
102        // Find the unreduced column with the greatest L2-Norm
103        int l2NormSquaredMaxIndex = minor;
104        for (int i = minor; i < qrt.length; i++) {
105            double l2NormSquared = 0;
106            for (int j = minor; j < qrt[i].length; j++) {
107                l2NormSquared += qrt[i][j] * qrt[i][j];
108            }
109            if (l2NormSquared > l2NormSquaredMax) {
110                l2NormSquaredMax = l2NormSquared;
111                l2NormSquaredMaxIndex = i;
112            }
113        }
114        // swap the current column with that with the greated L2-Norm and record in p
115        if (l2NormSquaredMaxIndex != minor) {
116            double[] tmp1 = qrt[minor];
117            qrt[minor] = qrt[l2NormSquaredMaxIndex];
118            qrt[l2NormSquaredMaxIndex] = tmp1;
119            int tmp2 = p[minor];
120            p[minor] = p[l2NormSquaredMaxIndex];
121            p[l2NormSquaredMaxIndex] = tmp2;
122        }
123
124        super.performHouseholderReflection(minor, qrt);
125    }
126
127
128    /**
129     * Returns the pivot matrix, P, used in the QR Decomposition of matrix A such that AP = QR.
130     *
131     * If no pivoting is used in this decomposition then P is equal to the identity matrix.
132     *
133     * @return a permutation matrix.
134     */
135    public RealMatrix getP() {
136        if (cachedP == null) {
137            int n = p.length;
138            cachedP = MatrixUtils.createRealMatrix(n,n);
139            for (int i = 0; i < n; i++) {
140                cachedP.setEntry(p[i], i, 1);
141            }
142        }
143        return cachedP ;
144    }
145
146    /**
147     * Return the effective numerical matrix rank.
148     * <p>The effective numerical rank is the number of non-negligible
149     * singular values.</p>
150     * <p>This implementation looks at Frobenius norms of the sequence of
151     * bottom right submatrices.  When a large fall in norm is seen,
152     * the rank is returned. The drop is computed as:</p>
153     * <pre>
154     *   (thisNorm/lastNorm) * rNorm &lt; dropThreshold
155     * </pre>
156     * <p>
157     * where thisNorm is the Frobenius norm of the current submatrix,
158     * lastNorm is the Frobenius norm of the previous submatrix,
159     * rNorm is is the Frobenius norm of the complete matrix
160     * </p>
161     *
162     * @param dropThreshold threshold triggering rank computation
163     * @return effective numerical matrix rank
164     */
165    public int getRank(final double dropThreshold) {
166        RealMatrix r    = getR();
167        int rows        = r.getRowDimension();
168        int columns     = r.getColumnDimension();
169        int rank        = 1;
170        double lastNorm = r.getFrobeniusNorm();
171        double rNorm    = lastNorm;
172        while (rank < FastMath.min(rows, columns)) {
173            double thisNorm = r.getSubMatrix(rank, rows - 1, rank, columns - 1).getFrobeniusNorm();
174            if (thisNorm == 0 || (thisNorm / lastNorm) * rNorm < dropThreshold) {
175                break;
176            }
177            lastNorm = thisNorm;
178            rank++;
179        }
180        return rank;
181    }
182
183    /**
184     * Get a solver for finding the A &times; X = B solution in least square sense.
185     * <p>
186     * Least Square sense means a solver can be computed for an overdetermined system,
187     * (i.e. a system with more equations than unknowns, which corresponds to a tall A
188     * matrix with more rows than columns). In any case, if the matrix is singular
189     * within the tolerance set at {@link RRQRDecomposition#RRQRDecomposition(RealMatrix,
190     * double) construction}, an error will be triggered when
191     * the {@link DecompositionSolver#solve(RealVector) solve} method will be called.
192     * </p>
193     * @return a solver
194     */
195    @Override
196    public DecompositionSolver getSolver() {
197        return new Solver(super.getSolver(), this.getP());
198    }
199
200    /** Specialized solver. */
201    private static class Solver implements DecompositionSolver {
202
203        /** Upper level solver. */
204        private final DecompositionSolver upper;
205
206        /** A permutation matrix for the pivots used in the QR decomposition */
207        private RealMatrix p;
208
209        /**
210         * Build a solver from decomposed matrix.
211         *
212         * @param upper upper level solver.
213         * @param p permutation matrix
214         */
215        private Solver(final DecompositionSolver upper, final RealMatrix p) {
216            this.upper = upper;
217            this.p     = p;
218        }
219
220        /** {@inheritDoc} */
221        @Override
222        public boolean isNonSingular() {
223            return upper.isNonSingular();
224        }
225
226        /** {@inheritDoc} */
227        @Override
228        public RealVector solve(RealVector b) {
229            return p.operate(upper.solve(b));
230        }
231
232        /** {@inheritDoc} */
233        @Override
234        public RealMatrix solve(RealMatrix b) {
235            return p.multiply(upper.solve(b));
236        }
237
238        /**
239         * {@inheritDoc}
240         * @throws SingularMatrixException if the decomposed matrix is singular.
241         */
242        @Override
243        public RealMatrix getInverse() {
244            return solve(MatrixUtils.createRealIdentityMatrix(p.getRowDimension()));
245        }
246    }
247}