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 org.hipparchus.util.FastMath;
26  
27  
28  /**
29   * Calculates the rank-revealing QR-decomposition of a matrix, with column pivoting.
30   * <p>The rank-revealing QR-decomposition of a matrix A consists of three matrices Q,
31   * R and P such that AP=QR.  Q is orthogonal (Q<sup>T</sup>Q = I), and R is upper triangular.
32   * If A is m&times;n, Q is m&times;m and R is m&times;n and P is n&times;n.</p>
33   * <p>QR decomposition with column pivoting produces a rank-revealing QR
34   * decomposition and the {@link #getRank(double)} method may be used to return the rank of the
35   * input matrix A.</p>
36   * <p>This class compute the decomposition using Householder reflectors.</p>
37   * <p>For efficiency purposes, the decomposition in packed form is transposed.
38   * This allows inner loop to iterate inside rows, which is much more cache-efficient
39   * in Java.</p>
40   * <p>This class is based on the class with similar name from the
41   * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
42   * following changes:</p>
43   * <ul>
44   *   <li>a {@link #getQT() getQT} method has been added,</li>
45   *   <li>the {@code solve} and {@code isFullRank} methods have been replaced
46   *   by a {@link #getSolver() getSolver} method and the equivalent methods
47   *   provided by the returned {@link DecompositionSolver}.</li>
48   * </ul>
49   *
50   * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
51   * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
52   *
53   */
54  public class RRQRDecomposition extends QRDecomposition {
55  
56      /** An array to record the column pivoting for later creation of P. */
57      private int[] p;
58  
59      /** Cached value of P. */
60      private RealMatrix cachedP;
61  
62  
63      /**
64       * Calculates the QR-decomposition of the given matrix.
65       * The singularity threshold defaults to zero.
66       *
67       * @param matrix The matrix to decompose.
68       *
69       * @see #RRQRDecomposition(RealMatrix, double)
70       */
71      public RRQRDecomposition(RealMatrix matrix) {
72          this(matrix, 0d);
73      }
74  
75     /**
76       * Calculates the QR-decomposition of the given matrix.
77       *
78       * @param matrix The matrix to decompose.
79       * @param threshold Singularity threshold.
80       * @see #RRQRDecomposition(RealMatrix)
81       */
82      public RRQRDecomposition(RealMatrix matrix,  double threshold) {
83          super(matrix, threshold);
84      }
85  
86      /** Decompose matrix.
87       * @param qrt transposed matrix
88       */
89      @Override
90      protected void decompose(double[][] qrt) {
91          p = new int[qrt.length];
92          for (int i = 0; i < p.length; i++) {
93              p[i] = i;
94          }
95          super.decompose(qrt);
96      }
97  
98      /** Perform Householder reflection for a minor A(minor, minor) of A.
99       * @param minor minor index
100      * @param qrt transposed matrix
101      */
102     @Override
103     protected void performHouseholderReflection(int minor, double[][] qrt) {
104 
105         double l2NormSquaredMax = 0;
106         // Find the unreduced column with the greatest L2-Norm
107         int l2NormSquaredMaxIndex = minor;
108         for (int i = minor; i < qrt.length; i++) {
109             double l2NormSquared = 0;
110             for (int j = 0; j < qrt[i].length; j++) {
111                 l2NormSquared += qrt[i][j] * qrt[i][j];
112             }
113             if (l2NormSquared > l2NormSquaredMax) {
114                 l2NormSquaredMax = l2NormSquared;
115                 l2NormSquaredMaxIndex = i;
116             }
117         }
118         // swap the current column with that with the greated L2-Norm and record in p
119         if (l2NormSquaredMaxIndex != minor) {
120             double[] tmp1 = qrt[minor];
121             qrt[minor] = qrt[l2NormSquaredMaxIndex];
122             qrt[l2NormSquaredMaxIndex] = tmp1;
123             int tmp2 = p[minor];
124             p[minor] = p[l2NormSquaredMaxIndex];
125             p[l2NormSquaredMaxIndex] = tmp2;
126         }
127 
128         super.performHouseholderReflection(minor, qrt);
129 
130     }
131 
132 
133     /**
134      * Returns the pivot matrix, P, used in the QR Decomposition of matrix A such that AP = QR.
135      *
136      * If no pivoting is used in this decomposition then P is equal to the identity matrix.
137      *
138      * @return a permutation matrix.
139      */
140     public RealMatrix getP() {
141         if (cachedP == null) {
142             int n = p.length;
143             cachedP = MatrixUtils.createRealMatrix(n,n);
144             for (int i = 0; i < n; i++) {
145                 cachedP.setEntry(p[i], i, 1);
146             }
147         }
148         return cachedP ;
149     }
150 
151     /**
152      * Return the effective numerical matrix rank.
153      * <p>The effective numerical rank is the number of non-negligible
154      * singular values.</p>
155      * <p>This implementation looks at Frobenius norms of the sequence of
156      * bottom right submatrices.  When a large fall in norm is seen,
157      * the rank is returned. The drop is computed as:</p>
158      * <pre>
159      *   (thisNorm/lastNorm) * rNorm &lt; dropThreshold
160      * </pre>
161      * <p>
162      * where thisNorm is the Frobenius norm of the current submatrix,
163      * lastNorm is the Frobenius norm of the previous submatrix,
164      * rNorm is is the Frobenius norm of the complete matrix
165      * </p>
166      *
167      * @param dropThreshold threshold triggering rank computation
168      * @return effective numerical matrix rank
169      */
170     public int getRank(final double dropThreshold) {
171         RealMatrix r    = getR();
172         int rows        = r.getRowDimension();
173         int columns     = r.getColumnDimension();
174         int rank        = 1;
175         double lastNorm = r.getFrobeniusNorm();
176         double rNorm    = lastNorm;
177         while (rank < FastMath.min(rows, columns)) {
178             double thisNorm = r.getSubMatrix(rank, rows - 1, rank, columns - 1).getFrobeniusNorm();
179             if (thisNorm == 0 || (thisNorm / lastNorm) * rNorm < dropThreshold) {
180                 break;
181             }
182             lastNorm = thisNorm;
183             rank++;
184         }
185         return rank;
186     }
187 
188     /**
189      * Get a solver for finding the A &times; X = B solution in least square sense.
190      * <p>
191      * Least Square sense means a solver can be computed for an overdetermined system,
192      * (i.e. a system with more equations than unknowns, which corresponds to a tall A
193      * matrix with more rows than columns). In any case, if the matrix is singular
194      * within the tolerance set at {@link RRQRDecomposition#RRQRDecomposition(RealMatrix,
195      * double) construction}, an error will be triggered when
196      * the {@link DecompositionSolver#solve(RealVector) solve} method will be called.
197      * </p>
198      * @return a solver
199      */
200     @Override
201     public DecompositionSolver getSolver() {
202         return new Solver(super.getSolver(), this.getP());
203     }
204 
205     /** Specialized solver. */
206     private static class Solver implements DecompositionSolver {
207 
208         /** Upper level solver. */
209         private final DecompositionSolver upper;
210 
211         /** A permutation matrix for the pivots used in the QR decomposition */
212         private RealMatrix p;
213 
214         /**
215          * Build a solver from decomposed matrix.
216          *
217          * @param upper upper level solver.
218          * @param p permutation matrix
219          */
220         private Solver(final DecompositionSolver upper, final RealMatrix p) {
221             this.upper = upper;
222             this.p     = p;
223         }
224 
225         /** {@inheritDoc} */
226         @Override
227         public boolean isNonSingular() {
228             return upper.isNonSingular();
229         }
230 
231         /** {@inheritDoc} */
232         @Override
233         public RealVector solve(RealVector b) {
234             return p.operate(upper.solve(b));
235         }
236 
237         /** {@inheritDoc} */
238         @Override
239         public RealMatrix solve(RealMatrix b) {
240             return p.multiply(upper.solve(b));
241         }
242 
243         /**
244          * {@inheritDoc}
245          * @throws org.hipparchus.exception.MathIllegalArgumentException
246          * if the decomposed matrix is singular.
247          */
248         @Override
249         public RealMatrix getInverse() {
250             return solve(MatrixUtils.createRealIdentityMatrix(p.getRowDimension()));
251         }
252         /** {@inheritDoc} */
253         @Override
254         public int getRowDimension() {
255             return upper.getRowDimension();
256         }
257 
258         /** {@inheritDoc} */
259         @Override
260         public int getColumnDimension() {
261             return upper.getColumnDimension();
262         }
263 
264     }
265 }