RRQRDecomposition.java

  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.  * This is not the original file distributed by the Apache Software Foundation
  19.  * It has been modified by the Hipparchus project
  20.  */

  21. package org.hipparchus.linear;

  22. import org.hipparchus.util.FastMath;


  23. /**
  24.  * Calculates the rank-revealing QR-decomposition of a matrix, with column pivoting.
  25.  * <p>The rank-revealing QR-decomposition of a matrix A consists of three matrices Q,
  26.  * R and P such that AP=QR.  Q is orthogonal (Q<sup>T</sup>Q = I), and R is upper triangular.
  27.  * If A is m&times;n, Q is m&times;m and R is m&times;n and P is n&times;n.</p>
  28.  * <p>QR decomposition with column pivoting produces a rank-revealing QR
  29.  * decomposition and the {@link #getRank(double)} method may be used to return the rank of the
  30.  * input matrix A.</p>
  31.  * <p>This class compute the decomposition using Householder reflectors.</p>
  32.  * <p>For efficiency purposes, the decomposition in packed form is transposed.
  33.  * This allows inner loop to iterate inside rows, which is much more cache-efficient
  34.  * in Java.</p>
  35.  * <p>This class is based on the class with similar name from the
  36.  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
  37.  * following changes:</p>
  38.  * <ul>
  39.  *   <li>a {@link #getQT() getQT} method has been added,</li>
  40.  *   <li>the {@code solve} and {@code isFullRank} methods have been replaced
  41.  *   by a {@link #getSolver() getSolver} method and the equivalent methods
  42.  *   provided by the returned {@link DecompositionSolver}.</li>
  43.  * </ul>
  44.  *
  45.  * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
  46.  * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
  47.  *
  48.  */
  49. public class RRQRDecomposition extends QRDecomposition {

  50.     /** An array to record the column pivoting for later creation of P. */
  51.     private int[] p;

  52.     /** Cached value of P. */
  53.     private RealMatrix cachedP;


  54.     /**
  55.      * Calculates the QR-decomposition of the given matrix.
  56.      * The singularity threshold defaults to zero.
  57.      *
  58.      * @param matrix The matrix to decompose.
  59.      *
  60.      * @see #RRQRDecomposition(RealMatrix, double)
  61.      */
  62.     public RRQRDecomposition(RealMatrix matrix) {
  63.         this(matrix, 0d);
  64.     }

  65.    /**
  66.      * Calculates the QR-decomposition of the given matrix.
  67.      *
  68.      * @param matrix The matrix to decompose.
  69.      * @param threshold Singularity threshold.
  70.      * @see #RRQRDecomposition(RealMatrix)
  71.      */
  72.     public RRQRDecomposition(RealMatrix matrix,  double threshold) {
  73.         super(matrix, threshold);
  74.     }

  75.     /** Decompose matrix.
  76.      * @param qrt transposed matrix
  77.      */
  78.     @Override
  79.     protected void decompose(double[][] qrt) {
  80.         p = new int[qrt.length];
  81.         for (int i = 0; i < p.length; i++) {
  82.             p[i] = i;
  83.         }
  84.         super.decompose(qrt);
  85.     }

  86.     /** Perform Householder reflection for a minor A(minor, minor) of A.
  87.      * @param minor minor index
  88.      * @param qrt transposed matrix
  89.      */
  90.     @Override
  91.     protected void performHouseholderReflection(int minor, double[][] qrt) {

  92.         double l2NormSquaredMax = 0;
  93.         // Find the unreduced column with the greatest L2-Norm
  94.         int l2NormSquaredMaxIndex = minor;
  95.         for (int i = minor; i < qrt.length; i++) {
  96.             double l2NormSquared = 0;
  97.             for (int j = 0; j < qrt[i].length; j++) {
  98.                 l2NormSquared += qrt[i][j] * qrt[i][j];
  99.             }
  100.             if (l2NormSquared > l2NormSquaredMax) {
  101.                 l2NormSquaredMax = l2NormSquared;
  102.                 l2NormSquaredMaxIndex = i;
  103.             }
  104.         }
  105.         // swap the current column with that with the greated L2-Norm and record in p
  106.         if (l2NormSquaredMaxIndex != minor) {
  107.             double[] tmp1 = qrt[minor];
  108.             qrt[minor] = qrt[l2NormSquaredMaxIndex];
  109.             qrt[l2NormSquaredMaxIndex] = tmp1;
  110.             int tmp2 = p[minor];
  111.             p[minor] = p[l2NormSquaredMaxIndex];
  112.             p[l2NormSquaredMaxIndex] = tmp2;
  113.         }

  114.         super.performHouseholderReflection(minor, qrt);

  115.     }


  116.     /**
  117.      * Returns the pivot matrix, P, used in the QR Decomposition of matrix A such that AP = QR.
  118.      *
  119.      * If no pivoting is used in this decomposition then P is equal to the identity matrix.
  120.      *
  121.      * @return a permutation matrix.
  122.      */
  123.     public RealMatrix getP() {
  124.         if (cachedP == null) {
  125.             int n = p.length;
  126.             cachedP = MatrixUtils.createRealMatrix(n,n);
  127.             for (int i = 0; i < n; i++) {
  128.                 cachedP.setEntry(p[i], i, 1);
  129.             }
  130.         }
  131.         return cachedP ;
  132.     }

  133.     /**
  134.      * Return the effective numerical matrix rank.
  135.      * <p>The effective numerical rank is the number of non-negligible
  136.      * singular values.</p>
  137.      * <p>This implementation looks at Frobenius norms of the sequence of
  138.      * bottom right submatrices.  When a large fall in norm is seen,
  139.      * the rank is returned. The drop is computed as:</p>
  140.      * <pre>
  141.      *   (thisNorm/lastNorm) * rNorm &lt; dropThreshold
  142.      * </pre>
  143.      * <p>
  144.      * where thisNorm is the Frobenius norm of the current submatrix,
  145.      * lastNorm is the Frobenius norm of the previous submatrix,
  146.      * rNorm is is the Frobenius norm of the complete matrix
  147.      * </p>
  148.      *
  149.      * @param dropThreshold threshold triggering rank computation
  150.      * @return effective numerical matrix rank
  151.      */
  152.     public int getRank(final double dropThreshold) {
  153.         RealMatrix r    = getR();
  154.         int rows        = r.getRowDimension();
  155.         int columns     = r.getColumnDimension();
  156.         int rank        = 1;
  157.         double lastNorm = r.getFrobeniusNorm();
  158.         double rNorm    = lastNorm;
  159.         while (rank < FastMath.min(rows, columns)) {
  160.             double thisNorm = r.getSubMatrix(rank, rows - 1, rank, columns - 1).getFrobeniusNorm();
  161.             if (thisNorm == 0 || (thisNorm / lastNorm) * rNorm < dropThreshold) {
  162.                 break;
  163.             }
  164.             lastNorm = thisNorm;
  165.             rank++;
  166.         }
  167.         return rank;
  168.     }

  169.     /**
  170.      * Get a solver for finding the A &times; X = B solution in least square sense.
  171.      * <p>
  172.      * Least Square sense means a solver can be computed for an overdetermined system,
  173.      * (i.e. a system with more equations than unknowns, which corresponds to a tall A
  174.      * matrix with more rows than columns). In any case, if the matrix is singular
  175.      * within the tolerance set at {@link RRQRDecomposition#RRQRDecomposition(RealMatrix,
  176.      * double) construction}, an error will be triggered when
  177.      * the {@link DecompositionSolver#solve(RealVector) solve} method will be called.
  178.      * </p>
  179.      * @return a solver
  180.      */
  181.     @Override
  182.     public DecompositionSolver getSolver() {
  183.         return new Solver(super.getSolver(), this.getP());
  184.     }

  185.     /** Specialized solver. */
  186.     private static class Solver implements DecompositionSolver {

  187.         /** Upper level solver. */
  188.         private final DecompositionSolver upper;

  189.         /** A permutation matrix for the pivots used in the QR decomposition */
  190.         private RealMatrix p;

  191.         /**
  192.          * Build a solver from decomposed matrix.
  193.          *
  194.          * @param upper upper level solver.
  195.          * @param p permutation matrix
  196.          */
  197.         private Solver(final DecompositionSolver upper, final RealMatrix p) {
  198.             this.upper = upper;
  199.             this.p     = p;
  200.         }

  201.         /** {@inheritDoc} */
  202.         @Override
  203.         public boolean isNonSingular() {
  204.             return upper.isNonSingular();
  205.         }

  206.         /** {@inheritDoc} */
  207.         @Override
  208.         public RealVector solve(RealVector b) {
  209.             return p.operate(upper.solve(b));
  210.         }

  211.         /** {@inheritDoc} */
  212.         @Override
  213.         public RealMatrix solve(RealMatrix b) {
  214.             return p.multiply(upper.solve(b));
  215.         }

  216.         /**
  217.          * {@inheritDoc}
  218.          * @throws org.hipparchus.exception.MathIllegalArgumentException
  219.          * if the decomposed matrix is singular.
  220.          */
  221.         @Override
  222.         public RealMatrix getInverse() {
  223.             return solve(MatrixUtils.createRealIdentityMatrix(p.getRowDimension()));
  224.         }
  225.         /** {@inheritDoc} */
  226.         @Override
  227.         public int getRowDimension() {
  228.             return upper.getRowDimension();
  229.         }

  230.         /** {@inheritDoc} */
  231.         @Override
  232.         public int getColumnDimension() {
  233.             return upper.getColumnDimension();
  234.         }

  235.     }
  236. }