RRQRDecomposition.java
- /*
- * Licensed to the Apache Software Foundation (ASF) under one or more
- * contributor license agreements. See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The ASF licenses this file to You under the Apache License, Version 2.0
- * (the "License"); you may not use this file except in compliance with
- * the License. You may obtain a copy of the License at
- *
- * https://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
- /*
- * This is not the original file distributed by the Apache Software Foundation
- * It has been modified by the Hipparchus project
- */
- package org.hipparchus.linear;
- import org.hipparchus.util.FastMath;
- /**
- * Calculates the rank-revealing QR-decomposition of a matrix, with column pivoting.
- * <p>The rank-revealing QR-decomposition of a matrix A consists of three matrices Q,
- * R and P such that AP=QR. Q is orthogonal (Q<sup>T</sup>Q = I), and R is upper triangular.
- * If A is m×n, Q is m×m and R is m×n and P is n×n.</p>
- * <p>QR decomposition with column pivoting produces a rank-revealing QR
- * decomposition and the {@link #getRank(double)} method may be used to return the rank of the
- * input matrix A.</p>
- * <p>This class compute the decomposition using Householder reflectors.</p>
- * <p>For efficiency purposes, the decomposition in packed form is transposed.
- * This allows inner loop to iterate inside rows, which is much more cache-efficient
- * in Java.</p>
- * <p>This class is based on the class with similar name from the
- * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
- * following changes:</p>
- * <ul>
- * <li>a {@link #getQT() getQT} method has been added,</li>
- * <li>the {@code solve} and {@code isFullRank} methods have been replaced
- * by a {@link #getSolver() getSolver} method and the equivalent methods
- * provided by the returned {@link DecompositionSolver}.</li>
- * </ul>
- *
- * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a>
- * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a>
- *
- */
- public class RRQRDecomposition extends QRDecomposition {
- /** An array to record the column pivoting for later creation of P. */
- private int[] p;
- /** Cached value of P. */
- private RealMatrix cachedP;
- /**
- * Calculates the QR-decomposition of the given matrix.
- * The singularity threshold defaults to zero.
- *
- * @param matrix The matrix to decompose.
- *
- * @see #RRQRDecomposition(RealMatrix, double)
- */
- public RRQRDecomposition(RealMatrix matrix) {
- this(matrix, 0d);
- }
- /**
- * Calculates the QR-decomposition of the given matrix.
- *
- * @param matrix The matrix to decompose.
- * @param threshold Singularity threshold.
- * @see #RRQRDecomposition(RealMatrix)
- */
- public RRQRDecomposition(RealMatrix matrix, double threshold) {
- super(matrix, threshold);
- }
- /** Decompose matrix.
- * @param qrt transposed matrix
- */
- @Override
- protected void decompose(double[][] qrt) {
- p = new int[qrt.length];
- for (int i = 0; i < p.length; i++) {
- p[i] = i;
- }
- super.decompose(qrt);
- }
- /** Perform Householder reflection for a minor A(minor, minor) of A.
- * @param minor minor index
- * @param qrt transposed matrix
- */
- @Override
- protected void performHouseholderReflection(int minor, double[][] qrt) {
- double l2NormSquaredMax = 0;
- // Find the unreduced column with the greatest L2-Norm
- int l2NormSquaredMaxIndex = minor;
- for (int i = minor; i < qrt.length; i++) {
- double l2NormSquared = 0;
- for (int j = 0; j < qrt[i].length; j++) {
- l2NormSquared += qrt[i][j] * qrt[i][j];
- }
- if (l2NormSquared > l2NormSquaredMax) {
- l2NormSquaredMax = l2NormSquared;
- l2NormSquaredMaxIndex = i;
- }
- }
- // swap the current column with that with the greated L2-Norm and record in p
- if (l2NormSquaredMaxIndex != minor) {
- double[] tmp1 = qrt[minor];
- qrt[minor] = qrt[l2NormSquaredMaxIndex];
- qrt[l2NormSquaredMaxIndex] = tmp1;
- int tmp2 = p[minor];
- p[minor] = p[l2NormSquaredMaxIndex];
- p[l2NormSquaredMaxIndex] = tmp2;
- }
- super.performHouseholderReflection(minor, qrt);
- }
- /**
- * Returns the pivot matrix, P, used in the QR Decomposition of matrix A such that AP = QR.
- *
- * If no pivoting is used in this decomposition then P is equal to the identity matrix.
- *
- * @return a permutation matrix.
- */
- public RealMatrix getP() {
- if (cachedP == null) {
- int n = p.length;
- cachedP = MatrixUtils.createRealMatrix(n,n);
- for (int i = 0; i < n; i++) {
- cachedP.setEntry(p[i], i, 1);
- }
- }
- return cachedP ;
- }
- /**
- * Return the effective numerical matrix rank.
- * <p>The effective numerical rank is the number of non-negligible
- * singular values.</p>
- * <p>This implementation looks at Frobenius norms of the sequence of
- * bottom right submatrices. When a large fall in norm is seen,
- * the rank is returned. The drop is computed as:</p>
- * <pre>
- * (thisNorm/lastNorm) * rNorm < dropThreshold
- * </pre>
- * <p>
- * where thisNorm is the Frobenius norm of the current submatrix,
- * lastNorm is the Frobenius norm of the previous submatrix,
- * rNorm is is the Frobenius norm of the complete matrix
- * </p>
- *
- * @param dropThreshold threshold triggering rank computation
- * @return effective numerical matrix rank
- */
- public int getRank(final double dropThreshold) {
- RealMatrix r = getR();
- int rows = r.getRowDimension();
- int columns = r.getColumnDimension();
- int rank = 1;
- double lastNorm = r.getFrobeniusNorm();
- double rNorm = lastNorm;
- while (rank < FastMath.min(rows, columns)) {
- double thisNorm = r.getSubMatrix(rank, rows - 1, rank, columns - 1).getFrobeniusNorm();
- if (thisNorm == 0 || (thisNorm / lastNorm) * rNorm < dropThreshold) {
- break;
- }
- lastNorm = thisNorm;
- rank++;
- }
- return rank;
- }
- /**
- * Get a solver for finding the A × X = B solution in least square sense.
- * <p>
- * Least Square sense means a solver can be computed for an overdetermined system,
- * (i.e. a system with more equations than unknowns, which corresponds to a tall A
- * matrix with more rows than columns). In any case, if the matrix is singular
- * within the tolerance set at {@link RRQRDecomposition#RRQRDecomposition(RealMatrix,
- * double) construction}, an error will be triggered when
- * the {@link DecompositionSolver#solve(RealVector) solve} method will be called.
- * </p>
- * @return a solver
- */
- @Override
- public DecompositionSolver getSolver() {
- return new Solver(super.getSolver(), this.getP());
- }
- /** Specialized solver. */
- private static class Solver implements DecompositionSolver {
- /** Upper level solver. */
- private final DecompositionSolver upper;
- /** A permutation matrix for the pivots used in the QR decomposition */
- private RealMatrix p;
- /**
- * Build a solver from decomposed matrix.
- *
- * @param upper upper level solver.
- * @param p permutation matrix
- */
- private Solver(final DecompositionSolver upper, final RealMatrix p) {
- this.upper = upper;
- this.p = p;
- }
- /** {@inheritDoc} */
- @Override
- public boolean isNonSingular() {
- return upper.isNonSingular();
- }
- /** {@inheritDoc} */
- @Override
- public RealVector solve(RealVector b) {
- return p.operate(upper.solve(b));
- }
- /** {@inheritDoc} */
- @Override
- public RealMatrix solve(RealMatrix b) {
- return p.multiply(upper.solve(b));
- }
- /**
- * {@inheritDoc}
- * @throws org.hipparchus.exception.MathIllegalArgumentException
- * if the decomposed matrix is singular.
- */
- @Override
- public RealMatrix getInverse() {
- return solve(MatrixUtils.createRealIdentityMatrix(p.getRowDimension()));
- }
- /** {@inheritDoc} */
- @Override
- public int getRowDimension() {
- return upper.getRowDimension();
- }
- /** {@inheritDoc} */
- @Override
- public int getColumnDimension() {
- return upper.getColumnDimension();
- }
- }
- }