| /* |
| * 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 |
| * |
| * http://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. |
| */ |
| |
| package org.apache.commons.math4.legacy.linear; |
| |
| import java.util.Arrays; |
| |
| import org.apache.commons.math4.legacy.exception.DimensionMismatchException; |
| import org.apache.commons.math4.core.jdkmath.JdkMath; |
| import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; |
| |
| |
| /** |
| * Calculates the QR-decomposition of a matrix. |
| * <p>The QR-decomposition of a matrix A consists of two matrices Q and R |
| * that satisfy: A = 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 m×n.</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> |
| * |
| * @since 1.2 (changed to concrete class in 3.0) |
| */ |
| public class QRDecomposition { |
| /** |
| * A packed TRANSPOSED representation of the QR decomposition. |
| * <p>The elements BELOW the diagonal are the elements of the UPPER triangular |
| * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors |
| * from which an explicit form of Q can be recomputed if desired.</p> |
| */ |
| private double[][] qrt; |
| /** The diagonal elements of R. */ |
| private double[] rDiag; |
| /** Cached value of Q. */ |
| private RealMatrix cachedQ; |
| /** Cached value of QT. */ |
| private RealMatrix cachedQT; |
| /** Cached value of R. */ |
| private RealMatrix cachedR; |
| /** Cached value of H. */ |
| private RealMatrix cachedH; |
| /** Singularity threshold. */ |
| private final double threshold; |
| |
| /** |
| * Calculates the QR-decomposition of the given matrix. |
| * The singularity threshold defaults to zero. |
| * |
| * @param matrix The matrix to decompose. |
| * |
| * @see #QRDecomposition(RealMatrix,double) |
| */ |
| public QRDecomposition(RealMatrix matrix) { |
| this(matrix, 0d); |
| } |
| |
| /** |
| * Calculates the QR-decomposition of the given matrix. |
| * |
| * @param matrix The matrix to decompose. |
| * @param threshold Singularity threshold. |
| * The matrix will be considered singular if the absolute value of |
| * any of the diagonal elements of the "R" matrix is smaller than |
| * the threshold. |
| */ |
| public QRDecomposition(RealMatrix matrix, |
| double threshold) { |
| this.threshold = threshold; |
| |
| final int m = matrix.getRowDimension(); |
| final int n = matrix.getColumnDimension(); |
| qrt = matrix.transpose().getData(); |
| rDiag = new double[JdkMath.min(m, n)]; |
| cachedQ = null; |
| cachedQT = null; |
| cachedR = null; |
| cachedH = null; |
| |
| decompose(qrt); |
| |
| } |
| |
| /** Decompose matrix. |
| * @param matrix transposed matrix |
| * @since 3.2 |
| */ |
| protected void decompose(double[][] matrix) { |
| for (int minor = 0; minor < JdkMath.min(matrix.length, matrix[0].length); minor++) { |
| performHouseholderReflection(minor, matrix); |
| } |
| } |
| |
| /** Perform Householder reflection for a minor A(minor, minor) of A. |
| * @param minor minor index |
| * @param matrix transposed matrix |
| * @since 3.2 |
| */ |
| protected void performHouseholderReflection(int minor, double[][] matrix) { |
| |
| final double[] qrtMinor = matrix[minor]; |
| |
| /* |
| * Let x be the first column of the minor, and a^2 = |x|^2. |
| * x will be in the positions qr[minor][minor] through qr[m][minor]. |
| * The first column of the transformed minor will be (a,0,0,..)' |
| * The sign of a is chosen to be opposite to the sign of the first |
| * component of x. Let's find a: |
| */ |
| double xNormSqr = 0; |
| for (int row = minor; row < qrtMinor.length; row++) { |
| final double c = qrtMinor[row]; |
| xNormSqr += c * c; |
| } |
| final double a = (qrtMinor[minor] > 0) ? -JdkMath.sqrt(xNormSqr) : JdkMath.sqrt(xNormSqr); |
| rDiag[minor] = a; |
| |
| if (a != 0.0) { |
| |
| /* |
| * Calculate the normalized reflection vector v and transform |
| * the first column. We know the norm of v beforehand: v = x-ae |
| * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> = |
| * a^2+a^2-2a<x,e> = 2a*(a - <x,e>). |
| * Here <x, e> is now qr[minor][minor]. |
| * v = x-ae is stored in the column at qr: |
| */ |
| qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor]) |
| |
| /* |
| * Transform the rest of the columns of the minor: |
| * They will be transformed by the matrix H = I-2vv'/|v|^2. |
| * If x is a column vector of the minor, then |
| * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v. |
| * Therefore the transformation is easily calculated by |
| * subtracting the column vector (2<x,v>/|v|^2)v from x. |
| * |
| * Let 2<x,v>/|v|^2 = alpha. From above we have |
| * |v|^2 = -2a*(qr[minor][minor]), so |
| * alpha = -<x,v>/(a*qr[minor][minor]) |
| */ |
| for (int col = minor+1; col < matrix.length; col++) { |
| final double[] qrtCol = matrix[col]; |
| double alpha = 0; |
| for (int row = minor; row < qrtCol.length; row++) { |
| alpha -= qrtCol[row] * qrtMinor[row]; |
| } |
| alpha /= a * qrtMinor[minor]; |
| |
| // Subtract the column vector alpha*v from x. |
| for (int row = minor; row < qrtCol.length; row++) { |
| qrtCol[row] -= alpha * qrtMinor[row]; |
| } |
| } |
| } |
| } |
| |
| |
| /** |
| * Returns the matrix R of the decomposition. |
| * <p>R is an upper-triangular matrix</p> |
| * @return the R matrix |
| */ |
| public RealMatrix getR() { |
| |
| if (cachedR == null) { |
| |
| // R is supposed to be m x n |
| final int n = qrt.length; |
| final int m = qrt[0].length; |
| double[][] ra = new double[m][n]; |
| // copy the diagonal from rDiag and the upper triangle of qr |
| for (int row = JdkMath.min(m, n) - 1; row >= 0; row--) { |
| ra[row][row] = rDiag[row]; |
| for (int col = row + 1; col < n; col++) { |
| ra[row][col] = qrt[col][row]; |
| } |
| } |
| cachedR = MatrixUtils.createRealMatrix(ra); |
| } |
| |
| // return the cached matrix |
| return cachedR; |
| } |
| |
| /** |
| * Returns the matrix Q of the decomposition. |
| * <p>Q is an orthogonal matrix</p> |
| * @return the Q matrix |
| */ |
| public RealMatrix getQ() { |
| if (cachedQ == null) { |
| cachedQ = getQT().transpose(); |
| } |
| return cachedQ; |
| } |
| |
| /** |
| * Returns the transpose of the matrix Q of the decomposition. |
| * <p>Q is an orthogonal matrix</p> |
| * @return the transpose of the Q matrix, Q<sup>T</sup> |
| */ |
| public RealMatrix getQT() { |
| if (cachedQT == null) { |
| |
| // QT is supposed to be m x m |
| final int n = qrt.length; |
| final int m = qrt[0].length; |
| double[][] qta = new double[m][m]; |
| |
| /* |
| * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then |
| * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in |
| * succession to the result |
| */ |
| for (int minor = m - 1; minor >= JdkMath.min(m, n); minor--) { |
| qta[minor][minor] = 1.0d; |
| } |
| |
| for (int minor = JdkMath.min(m, n)-1; minor >= 0; minor--){ |
| final double[] qrtMinor = qrt[minor]; |
| qta[minor][minor] = 1.0d; |
| if (qrtMinor[minor] != 0.0) { |
| for (int col = minor; col < m; col++) { |
| double alpha = 0; |
| for (int row = minor; row < m; row++) { |
| alpha -= qta[col][row] * qrtMinor[row]; |
| } |
| alpha /= rDiag[minor] * qrtMinor[minor]; |
| |
| for (int row = minor; row < m; row++) { |
| qta[col][row] += -alpha * qrtMinor[row]; |
| } |
| } |
| } |
| } |
| cachedQT = MatrixUtils.createRealMatrix(qta); |
| } |
| |
| // return the cached matrix |
| return cachedQT; |
| } |
| |
| /** |
| * Returns the Householder reflector vectors. |
| * <p>H is a lower trapezoidal matrix whose columns represent |
| * each successive Householder reflector vector. This matrix is used |
| * to compute Q.</p> |
| * @return a matrix containing the Householder reflector vectors |
| */ |
| public RealMatrix getH() { |
| if (cachedH == null) { |
| |
| final int n = qrt.length; |
| final int m = qrt[0].length; |
| double[][] ha = new double[m][n]; |
| for (int i = 0; i < m; ++i) { |
| for (int j = 0; j < JdkMath.min(i + 1, n); ++j) { |
| ha[i][j] = qrt[j][i] / -rDiag[j]; |
| } |
| } |
| cachedH = MatrixUtils.createRealMatrix(ha); |
| } |
| |
| // return the cached matrix |
| return cachedH; |
| } |
| |
| /** |
| * 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 QRDecomposition#QRDecomposition(RealMatrix, |
| * double) construction}, an error will be triggered when |
| * the {@link DecompositionSolver#solve(RealVector) solve} method will be called. |
| * </p> |
| * @return a solver |
| */ |
| public DecompositionSolver getSolver() { |
| return new Solver(qrt, rDiag, threshold); |
| } |
| |
| /** Specialized solver. */ |
| private static final class Solver implements DecompositionSolver { |
| /** |
| * A packed TRANSPOSED representation of the QR decomposition. |
| * <p>The elements BELOW the diagonal are the elements of the UPPER triangular |
| * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors |
| * from which an explicit form of Q can be recomputed if desired.</p> |
| */ |
| private final double[][] qrt; |
| /** The diagonal elements of R. */ |
| private final double[] rDiag; |
| /** Singularity threshold. */ |
| private final double threshold; |
| |
| /** |
| * Build a solver from decomposed matrix. |
| * |
| * @param qrt Packed TRANSPOSED representation of the QR decomposition. |
| * @param rDiag Diagonal elements of R. |
| * @param threshold Singularity threshold. |
| */ |
| private Solver(final double[][] qrt, |
| final double[] rDiag, |
| final double threshold) { |
| this.qrt = qrt; |
| this.rDiag = rDiag; |
| this.threshold = threshold; |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public boolean isNonSingular() { |
| return !checkSingular(rDiag, threshold, false); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public RealVector solve(RealVector b) { |
| final int n = qrt.length; |
| final int m = qrt[0].length; |
| if (b.getDimension() != m) { |
| throw new DimensionMismatchException(b.getDimension(), m); |
| } |
| checkSingular(rDiag, threshold, true); |
| |
| final double[] x = new double[n]; |
| final double[] y = b.toArray(); |
| |
| // apply Householder transforms to solve Q.y = b |
| for (int minor = 0; minor < JdkMath.min(m, n); minor++) { |
| |
| final double[] qrtMinor = qrt[minor]; |
| double dotProduct = 0; |
| for (int row = minor; row < m; row++) { |
| dotProduct += y[row] * qrtMinor[row]; |
| } |
| dotProduct /= rDiag[minor] * qrtMinor[minor]; |
| |
| for (int row = minor; row < m; row++) { |
| y[row] += dotProduct * qrtMinor[row]; |
| } |
| } |
| |
| // solve triangular system R.x = y |
| for (int row = rDiag.length - 1; row >= 0; --row) { |
| y[row] /= rDiag[row]; |
| final double yRow = y[row]; |
| final double[] qrtRow = qrt[row]; |
| x[row] = yRow; |
| for (int i = 0; i < row; i++) { |
| y[i] -= yRow * qrtRow[i]; |
| } |
| } |
| |
| return new ArrayRealVector(x, false); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public RealMatrix solve(RealMatrix b) { |
| final int n = qrt.length; |
| final int m = qrt[0].length; |
| if (b.getRowDimension() != m) { |
| throw new DimensionMismatchException(b.getRowDimension(), m); |
| } |
| checkSingular(rDiag, threshold, true); |
| |
| final int columns = b.getColumnDimension(); |
| final int blockSize = BlockRealMatrix.BLOCK_SIZE; |
| final int cBlocks = (columns + blockSize - 1) / blockSize; |
| final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns); |
| final double[][] y = new double[b.getRowDimension()][blockSize]; |
| final double[] alpha = new double[blockSize]; |
| |
| for (int kBlock = 0; kBlock < cBlocks; ++kBlock) { |
| final int kStart = kBlock * blockSize; |
| final int kEnd = JdkMath.min(kStart + blockSize, columns); |
| final int kWidth = kEnd - kStart; |
| |
| // get the right hand side vector |
| b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y); |
| |
| // apply Householder transforms to solve Q.y = b |
| for (int minor = 0; minor < JdkMath.min(m, n); minor++) { |
| final double[] qrtMinor = qrt[minor]; |
| final double factor = 1.0 / (rDiag[minor] * qrtMinor[minor]); |
| |
| Arrays.fill(alpha, 0, kWidth, 0.0); |
| for (int row = minor; row < m; ++row) { |
| final double d = qrtMinor[row]; |
| final double[] yRow = y[row]; |
| for (int k = 0; k < kWidth; ++k) { |
| alpha[k] += d * yRow[k]; |
| } |
| } |
| for (int k = 0; k < kWidth; ++k) { |
| alpha[k] *= factor; |
| } |
| |
| for (int row = minor; row < m; ++row) { |
| final double d = qrtMinor[row]; |
| final double[] yRow = y[row]; |
| for (int k = 0; k < kWidth; ++k) { |
| yRow[k] += alpha[k] * d; |
| } |
| } |
| } |
| |
| // solve triangular system R.x = y |
| for (int j = rDiag.length - 1; j >= 0; --j) { |
| final int jBlock = j / blockSize; |
| final int jStart = jBlock * blockSize; |
| final double factor = 1.0 / rDiag[j]; |
| final double[] yJ = y[j]; |
| final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock]; |
| int index = (j - jStart) * kWidth; |
| for (int k = 0; k < kWidth; ++k) { |
| yJ[k] *= factor; |
| xBlock[index++] = yJ[k]; |
| } |
| |
| final double[] qrtJ = qrt[j]; |
| for (int i = 0; i < j; ++i) { |
| final double rIJ = qrtJ[i]; |
| final double[] yI = y[i]; |
| for (int k = 0; k < kWidth; ++k) { |
| yI[k] -= yJ[k] * rIJ; |
| } |
| } |
| } |
| } |
| |
| return new BlockRealMatrix(n, columns, xBlocks, false); |
| } |
| |
| /** |
| * {@inheritDoc} |
| * @throws SingularMatrixException if the decomposed matrix is singular. |
| */ |
| @Override |
| public RealMatrix getInverse() { |
| return solve(MatrixUtils.createRealIdentityMatrix(qrt[0].length)); |
| } |
| |
| /** |
| * Check singularity. |
| * |
| * @param diag Diagonal elements of the R matrix. |
| * @param min Singularity threshold. |
| * @param raise Whether to raise a {@link SingularMatrixException} |
| * if any element of the diagonal fails the check. |
| * @return {@code true} if any element of the diagonal is smaller |
| * or equal to {@code min}. |
| * @throws SingularMatrixException if the matrix is singular and |
| * {@code raise} is {@code true}. |
| */ |
| private static boolean checkSingular(double[] diag, |
| double min, |
| boolean raise) { |
| final int len = diag.length; |
| for (int i = 0; i < len; i++) { |
| final double d = diag[i]; |
| if (JdkMath.abs(d) <= min) { |
| if (raise) { |
| final SingularMatrixException e = new SingularMatrixException(); |
| e.getContext().addMessage(LocalizedFormats.NUMBER_TOO_SMALL, d, min); |
| e.getContext().addMessage(LocalizedFormats.INDEX, i); |
| throw e; |
| } else { |
| return true; |
| } |
| } |
| } |
| return false; |
| } |
| } |
| } |