Will do. My aplogies! -Greg On Fri, Oct 7, 2011 at 3:55 AM, Luc Maisonobe <luc.maison...@free.fr> wrote:
> Le 07/10/2011 07:21, gr...@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<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<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<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<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<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 DefaultRealMatrixPreservingVis**itor() >> { >> + @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 DefaultRealMatrixPreservingVis**itor() >> { >> + @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 DefaultRealMatrixChangingVisit**or(){ >> + @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<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<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 DefaultRealMatrixChangingVisit**or() >> { >> + @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 DefaultRealMatrixChangingVisit**or() >> { >> + @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<dev-unsubscr...@commons.apache.org> > For additional commands, e-mail: dev-h...@commons.apache.org > >
