commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Luc Maisonobe <Luc.Maison...@free.fr>
Subject Re: 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 08:55:08 GMT
Le 07/10/2011 07:21, gregs@apache.org a écrit :
> 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

Hello Greg,

It seems the files do not have the right subversion properties.
Could you check your global subversion settings and make sure 
[auto-props] is set correctly ?

Thanks
Luc

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


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message