commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From gr...@apache.org
Subject svn commit: r1179935 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/linear/ test/java/org/apache/commons/math/linear/
Date Fri, 07 Oct 2011 05:21:17 GMT
Author: gregs
Date: Fri Oct  7 05:21:17 2011
New Revision: 1179935

URL: http://svn.apache.org/viewvc?rev=1179935&view=rev
Log:
JIRA Math-630 First push of PivotingQRDecomposition

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java

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;
+    }
+}



Mime
View raw message