Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java?rev=1179935&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
Fri Oct 7 05:21:17 2011
@@ -0,0 +1,421 @@
+/*
+ * Copyright 2011 The Apache Software Foundation.
+ *
+ * Licensed 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.util.MathUtils;
+import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.exception.DimensionMismatchException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ *
+ * @author gregsterijevski
+ */
+public class PivotingQRDecomposition {
+
+ private double[][] qr;
+ /** 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;
+ /** permutation info */
+ private int[] permutation;
+ /** the rank **/
+ private int rank;
+ /** vector of column multipliers */
+ private double[] beta;
+
+ public boolean isSingular() {
+ return rank != qr[0].length;
+ }
+
+ public int getRank() {
+ return rank;
+ }
+
+ public int[] getOrder() {
+ return MathUtils.copyOf(permutation);
+ }
+
+ public PivotingQRDecomposition(RealMatrix matrix) throws
ConvergenceException {
+ this(matrix, 1.0e-16, true);
+ }
+
+ public PivotingQRDecomposition(RealMatrix matrix, boolean allowPivot)
throws ConvergenceException {
+ this(matrix, 1.0e-16, allowPivot);
+ }
+
+ public PivotingQRDecomposition(RealMatrix matrix, double
qrRankingThreshold,
+ boolean allowPivot) throws ConvergenceException {
+ final int rows = matrix.getRowDimension();
+ final int cols = matrix.getColumnDimension();
+ qr = matrix.getData();
+ rDiag = new double[cols];
+ //final double[] norms = new double[cols];
+ this.beta = new double[cols];
+ this.permutation = new int[cols];
+ cachedQ = null;
+ cachedQT = null;
+ cachedR = null;
+ cachedH = null;
+
+ /*- initialize the permutation vector and calculate the norms */
+ for (int k = 0; k< cols; ++k) {
+ permutation[k] = k;
+ }
+ // transform the matrix column after column
+ for (int k = 0; k< cols; ++k) {
+ // select the column with the greatest norm on active components
+ int nextColumn = -1;
+ double ak2 = Double.NEGATIVE_INFINITY;
+ if (allowPivot) {
+ for (int i = k; i< cols; ++i) {
+ double norm2 = 0;
+ for (int j = k; j< rows; ++j) {
+ final double aki = qr[j][permutation[i]];
+ norm2 += aki * aki;
+ }
+ if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
+ throw new
ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN,
+ rows, cols);
+ }
+ if (norm2> ak2) {
+ nextColumn = i;
+ ak2 = norm2;
+ }
+ }
+ } else {
+ nextColumn = k;
+ ak2 = 0.0;
+ for (int j = k; j< rows; ++j) {
+ final double aki = qr[j][k];
+ ak2 += aki * aki;
+ }
+ }
+ if (ak2<= qrRankingThreshold) {
+ rank = k;
+ for (int i = rank; i< rows; i++) {
+ for (int j = i + 1; j< cols; j++) {
+ qr[i][permutation[j]] = 0.0;
+ }
+ }
+ return;
+ }
+ final int pk = permutation[nextColumn];
+ permutation[nextColumn] = permutation[k];
+ permutation[k] = pk;
+
+ // choose alpha such that Hk.u = alpha ek
+ final double akk = qr[k][pk];
+ final double alpha = (akk> 0) ? -FastMath.sqrt(ak2) :
FastMath.sqrt(ak2);
+ final double betak = 1.0 / (ak2 - akk * alpha);
+ beta[pk] = betak;
+
+ // transform the current column
+ rDiag[pk] = alpha;
+ qr[k][pk] -= alpha;
+
+ // transform the remaining columns
+ for (int dk = cols - 1 - k; dk> 0; --dk) {
+ double gamma = 0;
+ for (int j = k; j< rows; ++j) {
+ gamma += qr[j][pk] * qr[j][permutation[k + dk]];
+ }
+ gamma *= betak;
+ for (int j = k; j< rows; ++j) {
+ qr[j][permutation[k + dk]] -= gamma * qr[j][pk];
+ }
+ }
+ }
+ rank = cols;
+ return;
+ }
+
+ /**
+ * 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 Q matrix
+ */
+ public RealMatrix getQT() {
+ if (cachedQT == null) {
+
+ // QT is supposed to be m x m
+ final int n = qr[0].length;
+ final int m = qr.length;
+ cachedQT = MatrixUtils.createRealMatrix(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>= rank; minor--) {
+ cachedQT.setEntry(minor, minor, 1.0);
+ }
+
+ for (int minor = rank - 1; minor>= 0; minor--) {
+ //final double[] qrtMinor = qrt[minor];
+ final int p_minor = permutation[minor];
+ cachedQT.setEntry(minor, minor, 1.0);
+ //if (qrtMinor[minor] != 0.0) {
+ for (int col = minor; col< m; col++) {
+ double alpha = 0.0;
+ for (int row = minor; row< m; row++) {
+ alpha -= cachedQT.getEntry(col, row) *
qr[row][p_minor];
+ }
+ alpha /= rDiag[p_minor] * qr[minor][p_minor];
+ for (int row = minor; row< m; row++) {
+ cachedQT.addToEntry(col, row, -alpha *
qr[row][p_minor]);
+ }
+ }
+ //}
+ }
+ }
+ // return the cached matrix
+ return cachedQT;
+ }
+
+ /**
+ * 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 = qr[0].length;
+ final int m = qr.length;
+ cachedR = MatrixUtils.createRealMatrix(m, n);
+ // copy the diagonal from rDiag and the upper triangle of qr
+ for (int row = rank - 1; row>= 0; row--) {
+ cachedR.setEntry(row, row, rDiag[permutation[row]]);
+ for (int col = row + 1; col< n; col++) {
+ cachedR.setEntry(row, col, qr[row][permutation[col]]);
+ }
+ }
+ }
+ // return the cached matrix
+ return cachedR;
+ }
+
+ public RealMatrix getH() {
+ if (cachedH == null) {
+ final int n = qr[0].length;
+ final int m = qr.length;
+ cachedH = MatrixUtils.createRealMatrix(m, n);
+ for (int i = 0; i< m; ++i) {
+ for (int j = 0; j< FastMath.min(i + 1, n); ++j) {
+ final int p_j = permutation[j];
+ cachedH.setEntry(i, j, qr[i][p_j] / -rDiag[p_j]);
+ }
+ }
+ }
+ // return the cached matrix
+ return cachedH;
+ }
+
+ public RealMatrix getPermutationMatrix() {
+ RealMatrix rm = MatrixUtils.createRealMatrix(qr[0].length,
qr[0].length);
+ for (int i = 0; i< this.qr[0].length; i++) {
+ rm.setEntry(permutation[i], i, 1.0);
+ }
+ return rm;
+ }
+
+ public DecompositionSolver getSolver() {
+ return new Solver(qr, rDiag, permutation, rank);
+ }
+
+ /** Specialized solver. */
+ private static 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[][] qr;
+ /** The diagonal elements of R. */
+ private final double[] rDiag;
+ /** The rank of the matrix */
+ private final int rank;
+ /** The permutation matrix */
+ private final int[] perm;
+
+ /**
+ * 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[][] qr, final double[] rDiag, int[] perm,
int rank) {
+ this.qr = qr;
+ this.rDiag = rDiag;
+ this.perm = perm;
+ this.rank = rank;
+ }
+
+ /** {@inheritDoc} */
+ public boolean isNonSingular() {
+ if (qr.length>= qr[0].length) {
+ return rank == qr[0].length;
+ } else { //qr.length< qr[0].length
+ return rank == qr.length;
+ }
+ }
+
+ /** {@inheritDoc} */
+ public RealVector solve(RealVector b) {
+ final int n = qr[0].length;
+ final int m = qr.length;
+ if (b.getDimension() != m) {
+ throw new DimensionMismatchException(b.getDimension(), m);
+ }
+ if (!isNonSingular()) {
+ throw new SingularMatrixException();
+ }
+
+ final double[] x = new double[n];
+ final double[] y = b.toArray();
+
+ // apply Householder transforms to solve Q.y = b
+ for (int minor = 0; minor< rank; minor++) {
+ final int m_idx = perm[minor];
+ double dotProduct = 0;
+ for (int row = minor; row< m; row++) {
+ dotProduct += y[row] * qr[row][m_idx];
+ }
+ dotProduct /= rDiag[m_idx] * qr[minor][m_idx];
+ for (int row = minor; row< m; row++) {
+ y[row] += dotProduct * qr[row][m_idx];
+ }
+ }
+ // solve triangular system R.x = y
+ for (int row = rank - 1; row>= 0; --row) {
+ final int m_row = perm[row];
+ y[row] /= rDiag[m_row];
+ final double yRow = y[row];
+ //final double[] qrtRow = qrt[row];
+ x[perm[row]] = yRow;
+ for (int i = 0; i< row; i++) {
+ y[i] -= yRow * qr[i][m_row];
+ }
+ }
+ return new ArrayRealVector(x, false);
+ }
+
+ /** {@inheritDoc} */
+ public RealMatrix solve(RealMatrix b) {
+ final int cols = qr[0].length;
+ final int rows = qr.length;
+ if (b.getRowDimension() != rows) {
+ throw new DimensionMismatchException(b.getRowDimension(),
rows);
+ }
+ 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(cols, columns);
+ final double[][] y = new double[b.getRowDimension()][blockSize];
+ final double[] alpha = new double[blockSize];
+ //final BlockRealMatrix result = new BlockRealMatrix(cols,
columns, xBlocks, false);
+ 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, rows - 1, kStart, kEnd - 1, y);
+
+ // apply Householder transforms to solve Q.y = b
+ for (int minor = 0; minor< rank; minor++) {
+ final int m_idx = perm[minor];
+ final double factor = 1.0 / (rDiag[m_idx] *
qr[minor][m_idx]);
+
+ Arrays.fill(alpha, 0, kWidth, 0.0);
+ for (int row = minor; row< rows; ++row) {
+ final double d = qr[row][m_idx];
+ 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< rows; ++row) {
+ final double d = qr[row][m_idx];
+ 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 = rank - 1; j>= 0; --j) {
+ final int jBlock = perm[j] / blockSize; //which block
+ final int jStart = jBlock * blockSize; // idx of top
corner of block in my coord
+ final double factor = 1.0 / rDiag[perm[j]];
+ final double[] yJ = y[j];
+ final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
+ int index = (perm[j] - jStart) * kWidth; //to local
(block) coordinates
+ for (int k = 0; k< kWidth; ++k) {
+ yJ[k] *= factor;
+ xBlock[index++] = yJ[k];
+ }
+ for (int i = 0; i< j; ++i) {
+ final double rIJ = qr[i][perm[j]];
+ final double[] yI = y[i];
+ for (int k = 0; k< kWidth; ++k) {
+ yI[k] -= yJ[k] * rIJ;
+ }
+ }
+ }
+ }
+ //return result;
+ return new BlockRealMatrix(cols, columns, xBlocks, false);
+ }
+
+ /** {@inheritDoc} */
+ public RealMatrix getInverse() {
+ return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
+ }
+ }
+}
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java?rev=1179935&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java
Fri Oct 7 05:21:17 2011
@@ -0,0 +1,257 @@
+/*
+ * 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.Random;
+
+
+import org.apache.commons.math.ConvergenceException;
+import org.junit.Assert;
+import org.junit.Test;
+
+
+public class PivotingQRDecompositionTest {
+ double[][] testData3x3NonSingular = {
+ { 12, -51, 4 },
+ { 6, 167, -68 },
+ { -4, 24, -41 }, };
+
+ double[][] testData3x3Singular = {
+ { 1, 4, 7, },
+ { 2, 5, 8, },
+ { 3, 6, 9, }, };
+
+ double[][] testData3x4 = {
+ { 12, -51, 4, 1 },
+ { 6, 167, -68, 2 },
+ { -4, 24, -41, 3 }, };
+
+ double[][] testData4x3 = {
+ { 12, -51, 4, },
+ { 6, 167, -68, },
+ { -4, 24, -41, },
+ { -5, 34, 7, }, };
+
+ private static final double entryTolerance = 10e-16;
+
+ private static final double normTolerance = 10e-14;
+
+ /** test dimensions */
+ @Test
+ public void testDimensions() throws ConvergenceException {
+ checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+
+ checkDimension(MatrixUtils.createRealMatrix(testData4x3));
+
+ checkDimension(MatrixUtils.createRealMatrix(testData3x4));
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ checkDimension(createTestMatrix(r, p, q));
+ checkDimension(createTestMatrix(r, q, p));
+
+ }
+
+ private void checkDimension(RealMatrix m) throws ConvergenceException {
+ int rows = m.getRowDimension();
+ int columns = m.getColumnDimension();
+ PivotingQRDecomposition qr = new PivotingQRDecomposition(m);
+ Assert.assertEquals(rows, qr.getQ().getRowDimension());
+ Assert.assertEquals(rows, qr.getQ().getColumnDimension());
+ Assert.assertEquals(rows, qr.getR().getRowDimension());
+ Assert.assertEquals(columns, qr.getR().getColumnDimension());
+ }
+
+ /** test A = QR */
+ @Test
+ public void testAEqualQR() throws ConvergenceException {
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
+
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4));
+
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3));
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ checkAEqualQR(createTestMatrix(r, p, q));
+
+ checkAEqualQR(createTestMatrix(r, q, p));
+
+ }
+
+ private void checkAEqualQR(RealMatrix m) throws ConvergenceException {
+ PivotingQRDecomposition qr = new PivotingQRDecomposition(m);
+ RealMatrix prod =
qr.getQ().multiply(qr.getR()).multiply(qr.getPermutationMatrix().transpose());
+ double norm = prod.subtract(m).getNorm();
+ Assert.assertEquals(0, norm, normTolerance);
+ }
+
+ /** test the orthogonality of Q */
+ @Test
+ public void testQOrthogonal() throws ConvergenceException{
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular));
+
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4));
+
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3));
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ checkQOrthogonal(createTestMatrix(r, p, q));
+
+ checkQOrthogonal(createTestMatrix(r, q, p));
+
+ }
+
+ private void checkQOrthogonal(RealMatrix m) throws ConvergenceException{
+ PivotingQRDecomposition qr = new PivotingQRDecomposition(m);
+ RealMatrix eye =
MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
+ double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm();
+ Assert.assertEquals(0, norm, normTolerance);
+ }
+//
+ /** test that R is upper triangular */
+ @Test
+ public void testRUpperTriangular() throws ConvergenceException{
+ RealMatrix matrix =
MatrixUtils.createRealMatrix(testData3x3NonSingular);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x4);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = MatrixUtils.createRealMatrix(testData4x3);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ matrix = createTestMatrix(r, p, q);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = createTestMatrix(r, p, q);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ }
+
+ private void checkUpperTriangular(RealMatrix m) {
+ m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
+ @Override
+ public void visit(int row, int column, double value) {
+ if (column< row) {
+ Assert.assertEquals(0.0, value, entryTolerance);
+ }
+ }
+ });
+ }
+
+ /** test that H is trapezoidal */
+ @Test
+ public void testHTrapezoidal() throws ConvergenceException{
+ RealMatrix matrix =
MatrixUtils.createRealMatrix(testData3x3NonSingular);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x4);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = MatrixUtils.createRealMatrix(testData4x3);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ matrix = createTestMatrix(r, p, q);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = createTestMatrix(r, p, q);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ }
+
+ private void checkTrapezoidal(RealMatrix m) {
+ m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
+ @Override
+ public void visit(int row, int column, double value) {
+ if (column> row) {
+ Assert.assertEquals(0.0, value, entryTolerance);
+ }
+ }
+ });
+ }
+ /** test matrices values */
+ @Test
+ public void testMatricesValues() throws ConvergenceException{
+ PivotingQRDecomposition qr =
+ new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular),false);
+ RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] {
+ { -12.0 / 14.0, 69.0 / 175.0, -58.0 / 175.0 },
+ { -6.0 / 14.0, -158.0 / 175.0, 6.0 / 175.0 },
+ { 4.0 / 14.0, -30.0 / 175.0, -165.0 / 175.0 }
+ });
+ RealMatrix rRef = MatrixUtils.createRealMatrix(new double[][] {
+ { -14.0, -21.0, 14.0 },
+ { 0.0, -175.0, 70.0 },
+ { 0.0, 0.0, 35.0 }
+ });
+ RealMatrix hRef = MatrixUtils.createRealMatrix(new double[][] {
+ { 26.0 / 14.0, 0.0, 0.0 },
+ { 6.0 / 14.0, 648.0 / 325.0, 0.0 },
+ { -4.0 / 14.0, 36.0 / 325.0, 2.0 }
+ });
+
+ // check values against known references
+ RealMatrix q = qr.getQ();
+ Assert.assertEquals(0, q.subtract(qRef).getNorm(), 1.0e-13);
+ RealMatrix qT = qr.getQT();
+ Assert.assertEquals(0, qT.subtract(qRef.transpose()).getNorm(),
1.0e-13);
+ RealMatrix r = qr.getR();
+ Assert.assertEquals(0, r.subtract(rRef).getNorm(), 1.0e-13);
+ RealMatrix h = qr.getH();
+ Assert.assertEquals(0, h.subtract(hRef).getNorm(), 1.0e-13);
+
+ // check the same cached instance is returned the second time
+ Assert.assertTrue(q == qr.getQ());
+ Assert.assertTrue(r == qr.getR());
+ Assert.assertTrue(h == qr.getH());
+
+ }
+
+ private RealMatrix createTestMatrix(final Random r, final int rows, final
int columns) {
+ RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
+ m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
+ @Override
+ public double visit(int row, int column, double value) {
+ return 2.0 * r.nextDouble() - 1.0;
+ }
+ });
+ return m;
+ }
+
+}
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java?rev=1179935&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java
Fri Oct 7 05:21:17 2011
@@ -0,0 +1,201 @@
+/*
+ * 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.Random;
+
+import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.exception.MathIllegalArgumentException;
+
+import org.junit.Test;
+import org.junit.Assert;
+
+public class PivotingQRSolverTest {
+ double[][] testData3x3NonSingular = {
+ { 12, -51, 4 },
+ { 6, 167, -68 },
+ { -4, 24, -41 }
+ };
+
+ double[][] testData3x3Singular = {
+ { 1, 2, 2 },
+ { 2, 4, 6 },
+ { 4, 8, 12 }
+ };
+
+ double[][] testData3x4 = {
+ { 12, -51, 4, 1 },
+ { 6, 167, -68, 2 },
+ { -4, 24, -41, 3 }
+ };
+
+ double[][] testData4x3 = {
+ { 12, -51, 4 },
+ { 6, 167, -68 },
+ { -4, 24, -41 },
+ { -5, 34, 7 }
+ };
+
+ /** test rank */
+ @Test
+ public void testRank() throws ConvergenceException {
+ DecompositionSolver solver =
+ new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ solver = new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
+ Assert.assertFalse(solver.isNonSingular());
+
+ solver = new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ solver = new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ }
+
+ /** test solve dimension errors */
+ @Test
+ public void testSolveDimensionErrors() throws ConvergenceException {
+ DecompositionSolver solver =
+ new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
+ RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
+ try {
+ solver.solve(b);
+ Assert.fail("an exception should have been thrown");
+ } catch (MathIllegalArgumentException iae) {
+ // expected behavior
+ }
+ try {
+ solver.solve(b.getColumnVector(0));
+ Assert.fail("an exception should have been thrown");
+ } catch (MathIllegalArgumentException iae) {
+ // expected behavior
+ }
+ }
+
+ /** test solve rank errors */
+ @Test
+ public void testSolveRankErrors() throws ConvergenceException {
+ DecompositionSolver solver =
+ new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
+ RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]);
+ try {
+ solver.solve(b);
+ Assert.fail("an exception should have been thrown");
+ } catch (SingularMatrixException iae) {
+ // expected behavior
+ }
+ try {
+ solver.solve(b.getColumnVector(0));
+ Assert.fail("an exception should have been thrown");
+ } catch (SingularMatrixException iae) {
+ // expected behavior
+ }
+ }
+
+ /** test solve */
+ @Test
+ public void testSolve() throws ConvergenceException {
+ PivotingQRDecomposition decomposition =
+ new
PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+ DecompositionSolver solver = decomposition.getSolver();
+ RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
+ { -102, 12250 }, { 544, 24500 }, { 167, -36750 }
+ });
+
+ RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
+ { 1, 2515 }, { 2, 422 }, { -3, 898 }
+ });
+
+ // using RealMatrix
+ Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(),
2.0e-14 * xRef.getNorm());
+
+ // using ArrayRealVector
+ for (int i = 0; i< b.getColumnDimension(); ++i) {
+ final RealVector x = solver.solve(b.getColumnVector(i));
+ final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
+ Assert.assertEquals(0, error, 3.0e-14 *
xRef.getColumnVector(i).getNorm());
+ }
+
+ // using RealVector with an alternate implementation
+ for (int i = 0; i< b.getColumnDimension(); ++i) {
+ ArrayRealVectorTest.RealVectorTestImpl v =
+ new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
+ final RealVector x = solver.solve(v);
+ final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
+ Assert.assertEquals(0, error, 3.0e-14 *
xRef.getColumnVector(i).getNorm());
+ }
+
+ }
+
+ @Test
+ public void testOverdetermined() throws ConvergenceException {
+ final Random r = new Random(5559252868205245l);
+ int p = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ RealMatrix a = createTestMatrix(r, p, q);
+ RealMatrix xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE
+ 3);
+
+ // build a perturbed system: A.X + noise = B
+ RealMatrix b = a.multiply(xRef);
+ final double noise = 0.001;
+ b.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
+ @Override
+ public double visit(int row, int column, double value) {
+ return value * (1.0 + noise * (2 * r.nextDouble() - 1));
+ }
+ });
+
+ // despite perturbation, the least square solution should be pretty
good
+ RealMatrix x = new PivotingQRDecomposition(a).getSolver().solve(b);
+ Assert.assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p *
q);
+
+ }
+
+ @Test
+ public void testUnderdetermined() throws ConvergenceException {
+ final Random r = new Random(42185006424567123l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ RealMatrix a = createTestMatrix(r, p, q);
+ RealMatrix xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE
+ 3);
+ RealMatrix b = a.multiply(xRef);
+ PivotingQRDecomposition pqr = new PivotingQRDecomposition(a);
+ RealMatrix x = pqr.getSolver().solve(b);
+ Assert.assertTrue(x.subtract(xRef).getNorm() / (p * q)> 0.01);
+ int count=0;
+ for( int i = 0 ; i< q; i++){
+ if( x.getRowVector(i).getNorm() == 0.0 ){
+ ++count;
+ }
+ }
+ Assert.assertEquals("Zeroed rows", q-p, count);
+ }
+
+ private RealMatrix createTestMatrix(final Random r, final int rows, final
int columns) {
+ RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
+ m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
+ @Override
+ public double visit(int row, int column, double value) {
+ return 2.0 * r.nextDouble() - 1.0;
+ }
+ });
+ return m;
+ }
+}
