/* * 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.math.linear; import java.util.Arrays; import org.apache.commons.math.MathRuntimeException; import org.apache.commons.math.exception.util.LocalizedFormats; import org.apache.commons.math.util.FastMath; /** * Calculates the QR-decomposition of a matrix. *
The QR-decomposition of a matrix A consists of two matrices Q and R * that satisfy: A = QR, Q is orthogonal (QTQ = I), and R is * upper triangular. If A is m×n, Q is m×m and R m×n.
*This class compute the decomposition using Householder reflectors.
*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.
* * @see MathWorld * @see Wikipedia * * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ * @since 1.2 */ public class QRDecompositionImpl implements QRDecomposition { /** * A packed TRANSPOSED representation of the QR decomposition. *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.
*/ 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; /** * Calculates the QR-decomposition of the given matrix. * @param matrix The matrix to decompose. */ public QRDecompositionImpl(RealMatrix matrix) { final int m = matrix.getRowDimension(); final int n = matrix.getColumnDimension(); qrt = matrix.transpose().getData(); rDiag = new double[FastMath.min(m, n)]; cachedQ = null; cachedQT = null; cachedR = null; cachedH = null; /* * The QR decomposition of a matrix A is calculated using Householder * reflectors by repeating the following operations to each minor * A(minor,minor) of A: */ for (int minor = 0; minor < FastMath.min(m, n); minor++) { final double[] qrtMinor = qrt[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 < m; row++) { final double c = qrtMinor[row]; xNormSqr += c * c; } final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.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 =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.
*/ private final double[][] qrt; /** The diagonal elements of R. */ private final double[] rDiag; /** * Build a solver from decomposed matrix. * @param qrt packed TRANSPOSED representation of the QR decomposition * @param rDiag diagonal elements of R */ private Solver(final double[][] qrt, final double[] rDiag) { this.qrt = qrt; this.rDiag = rDiag; } /** {@inheritDoc} */ public boolean isNonSingular() { for (double diag : rDiag) { if (diag == 0) { return false; } } return true; } /** {@inheritDoc} */ public double[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException { final int n = qrt.length; final int m = qrt[0].length; if (b.length != m) { throw MathRuntimeException.createIllegalArgumentException( LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.length, m); } if (!isNonSingular()) { throw new SingularMatrixException(); } final double[] x = new double[n]; final double[] y = b.clone(); // apply Householder transforms to solve Q.y = b for (int minor = 0; minor < FastMath.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 x; } /** {@inheritDoc} */ public RealVector solve(RealVector b) throws IllegalArgumentException, InvalidMatrixException { try { return solve((ArrayRealVector) b); } catch (ClassCastException cce) { return new ArrayRealVector(solve(b.getData()), false); } } /** Solve the linear equation A × X = B. *The A matrix is implicit here. It is
* @param b right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B * @throws IllegalArgumentException if matrices dimensions don't match * @throws InvalidMatrixException if decomposed matrix is singular */ public ArrayRealVector solve(ArrayRealVector b) throws IllegalArgumentException, InvalidMatrixException { return new ArrayRealVector(solve(b.getDataRef()), false); } /** {@inheritDoc} */ public RealMatrix solve(RealMatrix b) throws IllegalArgumentException, InvalidMatrixException { final int n = qrt.length; final int m = qrt[0].length; if (b.getRowDimension() != m) { throw MathRuntimeException.createIllegalArgumentException( LocalizedFormats.DIMENSIONS_MISMATCH_2x2, b.getRowDimension(), b.getColumnDimension(), m, "n"); } if (!isNonSingular()) { throw new SingularMatrixException(); } 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 = FastMath.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 < FastMath.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} */ public RealMatrix getInverse() throws InvalidMatrixException { return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length)); } } }