Sorry Luc, I had the file set on my linux box. When it died, I moved to the mac and neglected to make the changes.
On Fri, Oct 7, 2011 at 10:21 AM, Greg Sterijevski <gsterijev...@gmail.com>wrote: > 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 >> >> >
