You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by gr...@apache.org on 2011/10/07 07:21:17 UTC
svn commit: r1179935 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/linear/
test/java/org/apache/commons/math/linear/
Author: gregs
Date: Fri Oct 7 05:21:17 2011
New Revision: 1179935
URL: http://svn.apache.org/viewvc?rev=1179935&view=rev
Log:
JIRA Math-630 First push of PivotingQRDecomposition
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java?rev=1179935&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/PivotingQRDecomposition.java Fri Oct 7 05:21:17 2011
@@ -0,0 +1,421 @@
+/*
+ * Copyright 2011 The Apache Software Foundation.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.linear;
+
+import java.util.Arrays;
+import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.exception.DimensionMismatchException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ *
+ * @author gregsterijevski
+ */
+public class PivotingQRDecomposition {
+
+ private double[][] qr;
+ /** The diagonal elements of R. */
+ private double[] rDiag;
+ /** Cached value of Q. */
+ private RealMatrix cachedQ;
+ /** Cached value of QT. */
+ private RealMatrix cachedQT;
+ /** Cached value of R. */
+ private RealMatrix cachedR;
+ /** Cached value of H. */
+ private RealMatrix cachedH;
+ /** permutation info */
+ private int[] permutation;
+ /** the rank **/
+ private int rank;
+ /** vector of column multipliers */
+ private double[] beta;
+
+ public boolean isSingular() {
+ return rank != qr[0].length;
+ }
+
+ public int getRank() {
+ return rank;
+ }
+
+ public int[] getOrder() {
+ return MathUtils.copyOf(permutation);
+ }
+
+ public PivotingQRDecomposition(RealMatrix matrix) throws ConvergenceException {
+ this(matrix, 1.0e-16, true);
+ }
+
+ public PivotingQRDecomposition(RealMatrix matrix, boolean allowPivot) throws ConvergenceException {
+ this(matrix, 1.0e-16, allowPivot);
+ }
+
+ public PivotingQRDecomposition(RealMatrix matrix, double qrRankingThreshold,
+ boolean allowPivot) throws ConvergenceException {
+ final int rows = matrix.getRowDimension();
+ final int cols = matrix.getColumnDimension();
+ qr = matrix.getData();
+ rDiag = new double[cols];
+ //final double[] norms = new double[cols];
+ this.beta = new double[cols];
+ this.permutation = new int[cols];
+ cachedQ = null;
+ cachedQT = null;
+ cachedR = null;
+ cachedH = null;
+
+ /*- initialize the permutation vector and calculate the norms */
+ for (int k = 0; k < cols; ++k) {
+ permutation[k] = k;
+ }
+ // transform the matrix column after column
+ for (int k = 0; k < cols; ++k) {
+ // select the column with the greatest norm on active components
+ int nextColumn = -1;
+ double ak2 = Double.NEGATIVE_INFINITY;
+ if (allowPivot) {
+ for (int i = k; i < cols; ++i) {
+ double norm2 = 0;
+ for (int j = k; j < rows; ++j) {
+ final double aki = qr[j][permutation[i]];
+ norm2 += aki * aki;
+ }
+ if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
+ throw new ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN,
+ rows, cols);
+ }
+ if (norm2 > ak2) {
+ nextColumn = i;
+ ak2 = norm2;
+ }
+ }
+ } else {
+ nextColumn = k;
+ ak2 = 0.0;
+ for (int j = k; j < rows; ++j) {
+ final double aki = qr[j][k];
+ ak2 += aki * aki;
+ }
+ }
+ if (ak2 <= qrRankingThreshold) {
+ rank = k;
+ for (int i = rank; i < rows; i++) {
+ for (int j = i + 1; j < cols; j++) {
+ qr[i][permutation[j]] = 0.0;
+ }
+ }
+ return;
+ }
+ final int pk = permutation[nextColumn];
+ permutation[nextColumn] = permutation[k];
+ permutation[k] = pk;
+
+ // choose alpha such that Hk.u = alpha ek
+ final double akk = qr[k][pk];
+ final double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2);
+ final double betak = 1.0 / (ak2 - akk * alpha);
+ beta[pk] = betak;
+
+ // transform the current column
+ rDiag[pk] = alpha;
+ qr[k][pk] -= alpha;
+
+ // transform the remaining columns
+ for (int dk = cols - 1 - k; dk > 0; --dk) {
+ double gamma = 0;
+ for (int j = k; j < rows; ++j) {
+ gamma += qr[j][pk] * qr[j][permutation[k + dk]];
+ }
+ gamma *= betak;
+ for (int j = k; j < rows; ++j) {
+ qr[j][permutation[k + dk]] -= gamma * qr[j][pk];
+ }
+ }
+ }
+ rank = cols;
+ return;
+ }
+
+ /**
+ * Returns the matrix Q of the decomposition.
+ * <p>Q is an orthogonal matrix</p>
+ * @return the Q matrix
+ */
+ public RealMatrix getQ() {
+ if (cachedQ == null) {
+ cachedQ = getQT().transpose();
+ }
+ return cachedQ;
+ }
+
+ /**
+ * Returns the transpose of the matrix Q of the decomposition.
+ * <p>Q is an orthogonal matrix</p>
+ * @return the Q matrix
+ */
+ public RealMatrix getQT() {
+ if (cachedQT == null) {
+
+ // QT is supposed to be m x m
+ final int n = qr[0].length;
+ final int m = qr.length;
+ cachedQT = MatrixUtils.createRealMatrix(m, m);
+
+ /*
+ * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then
+ * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in
+ * succession to the result
+ */
+ for (int minor = m - 1; minor >= rank; minor--) {
+ cachedQT.setEntry(minor, minor, 1.0);
+ }
+
+ for (int minor = rank - 1; minor >= 0; minor--) {
+ //final double[] qrtMinor = qrt[minor];
+ final int p_minor = permutation[minor];
+ cachedQT.setEntry(minor, minor, 1.0);
+ //if (qrtMinor[minor] != 0.0) {
+ for (int col = minor; col < m; col++) {
+ double alpha = 0.0;
+ for (int row = minor; row < m; row++) {
+ alpha -= cachedQT.getEntry(col, row) * qr[row][p_minor];
+ }
+ alpha /= rDiag[p_minor] * qr[minor][p_minor];
+ for (int row = minor; row < m; row++) {
+ cachedQT.addToEntry(col, row, -alpha * qr[row][p_minor]);
+ }
+ }
+ //}
+ }
+ }
+ // return the cached matrix
+ return cachedQT;
+ }
+
+ /**
+ * Returns the matrix R of the decomposition.
+ * <p>R is an upper-triangular matrix</p>
+ * @return the R matrix
+ */
+ public RealMatrix getR() {
+ if (cachedR == null) {
+ // R is supposed to be m x n
+ final int n = qr[0].length;
+ final int m = qr.length;
+ cachedR = MatrixUtils.createRealMatrix(m, n);
+ // copy the diagonal from rDiag and the upper triangle of qr
+ for (int row = rank - 1; row >= 0; row--) {
+ cachedR.setEntry(row, row, rDiag[permutation[row]]);
+ for (int col = row + 1; col < n; col++) {
+ cachedR.setEntry(row, col, qr[row][permutation[col]]);
+ }
+ }
+ }
+ // return the cached matrix
+ return cachedR;
+ }
+
+ public RealMatrix getH() {
+ if (cachedH == null) {
+ final int n = qr[0].length;
+ final int m = qr.length;
+ cachedH = MatrixUtils.createRealMatrix(m, n);
+ for (int i = 0; i < m; ++i) {
+ for (int j = 0; j < FastMath.min(i + 1, n); ++j) {
+ final int p_j = permutation[j];
+ cachedH.setEntry(i, j, qr[i][p_j] / -rDiag[p_j]);
+ }
+ }
+ }
+ // return the cached matrix
+ return cachedH;
+ }
+
+ public RealMatrix getPermutationMatrix() {
+ RealMatrix rm = MatrixUtils.createRealMatrix(qr[0].length, qr[0].length);
+ for (int i = 0; i < this.qr[0].length; i++) {
+ rm.setEntry(permutation[i], i, 1.0);
+ }
+ return rm;
+ }
+
+ public DecompositionSolver getSolver() {
+ return new Solver(qr, rDiag, permutation, rank);
+ }
+
+ /** Specialized solver. */
+ private static class Solver implements DecompositionSolver {
+
+ /**
+ * A packed TRANSPOSED representation of the QR decomposition.
+ * <p>The elements BELOW the diagonal are the elements of the UPPER triangular
+ * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors
+ * from which an explicit form of Q can be recomputed if desired.</p>
+ */
+ private final double[][] qr;
+ /** The diagonal elements of R. */
+ private final double[] rDiag;
+ /** The rank of the matrix */
+ private final int rank;
+ /** The permutation matrix */
+ private final int[] perm;
+
+ /**
+ * Build a solver from decomposed matrix.
+ * @param qrt packed TRANSPOSED representation of the QR decomposition
+ * @param rDiag diagonal elements of R
+ */
+ private Solver(final double[][] qr, final double[] rDiag, int[] perm, int rank) {
+ this.qr = qr;
+ this.rDiag = rDiag;
+ this.perm = perm;
+ this.rank = rank;
+ }
+
+ /** {@inheritDoc} */
+ public boolean isNonSingular() {
+ if (qr.length >= qr[0].length) {
+ return rank == qr[0].length;
+ } else { //qr.length < qr[0].length
+ return rank == qr.length;
+ }
+ }
+
+ /** {@inheritDoc} */
+ public RealVector solve(RealVector b) {
+ final int n = qr[0].length;
+ final int m = qr.length;
+ if (b.getDimension() != m) {
+ throw new DimensionMismatchException(b.getDimension(), m);
+ }
+ if (!isNonSingular()) {
+ throw new SingularMatrixException();
+ }
+
+ final double[] x = new double[n];
+ final double[] y = b.toArray();
+
+ // apply Householder transforms to solve Q.y = b
+ for (int minor = 0; minor < rank; minor++) {
+ final int m_idx = perm[minor];
+ double dotProduct = 0;
+ for (int row = minor; row < m; row++) {
+ dotProduct += y[row] * qr[row][m_idx];
+ }
+ dotProduct /= rDiag[m_idx] * qr[minor][m_idx];
+ for (int row = minor; row < m; row++) {
+ y[row] += dotProduct * qr[row][m_idx];
+ }
+ }
+ // solve triangular system R.x = y
+ for (int row = rank - 1; row >= 0; --row) {
+ final int m_row = perm[row];
+ y[row] /= rDiag[m_row];
+ final double yRow = y[row];
+ //final double[] qrtRow = qrt[row];
+ x[perm[row]] = yRow;
+ for (int i = 0; i < row; i++) {
+ y[i] -= yRow * qr[i][m_row];
+ }
+ }
+ return new ArrayRealVector(x, false);
+ }
+
+ /** {@inheritDoc} */
+ public RealMatrix solve(RealMatrix b) {
+ final int cols = qr[0].length;
+ final int rows = qr.length;
+ if (b.getRowDimension() != rows) {
+ throw new DimensionMismatchException(b.getRowDimension(), rows);
+ }
+ if (!isNonSingular()) {
+ throw new SingularMatrixException();
+ }
+
+ final int columns = b.getColumnDimension();
+ final int blockSize = BlockRealMatrix.BLOCK_SIZE;
+ final int cBlocks = (columns + blockSize - 1) / blockSize;
+ final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(cols, columns);
+ final double[][] y = new double[b.getRowDimension()][blockSize];
+ final double[] alpha = new double[blockSize];
+ //final BlockRealMatrix result = new BlockRealMatrix(cols, columns, xBlocks, false);
+ for (int kBlock = 0; kBlock < cBlocks; ++kBlock) {
+ final int kStart = kBlock * blockSize;
+ final int kEnd = FastMath.min(kStart + blockSize, columns);
+ final int kWidth = kEnd - kStart;
+ // get the right hand side vector
+ b.copySubMatrix(0, rows - 1, kStart, kEnd - 1, y);
+
+ // apply Householder transforms to solve Q.y = b
+ for (int minor = 0; minor < rank; minor++) {
+ final int m_idx = perm[minor];
+ final double factor = 1.0 / (rDiag[m_idx] * qr[minor][m_idx]);
+
+ Arrays.fill(alpha, 0, kWidth, 0.0);
+ for (int row = minor; row < rows; ++row) {
+ final double d = qr[row][m_idx];
+ final double[] yRow = y[row];
+ for (int k = 0; k < kWidth; ++k) {
+ alpha[k] += d * yRow[k];
+ }
+ }
+ for (int k = 0; k < kWidth; ++k) {
+ alpha[k] *= factor;
+ }
+
+ for (int row = minor; row < rows; ++row) {
+ final double d = qr[row][m_idx];
+ final double[] yRow = y[row];
+ for (int k = 0; k < kWidth; ++k) {
+ yRow[k] += alpha[k] * d;
+ }
+ }
+ }
+
+ // solve triangular system R.x = y
+ for (int j = rank - 1; j >= 0; --j) {
+ final int jBlock = perm[j] / blockSize; //which block
+ final int jStart = jBlock * blockSize; // idx of top corner of block in my coord
+ final double factor = 1.0 / rDiag[perm[j]];
+ final double[] yJ = y[j];
+ final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock];
+ int index = (perm[j] - jStart) * kWidth; //to local (block) coordinates
+ for (int k = 0; k < kWidth; ++k) {
+ yJ[k] *= factor;
+ xBlock[index++] = yJ[k];
+ }
+ for (int i = 0; i < j; ++i) {
+ final double rIJ = qr[i][perm[j]];
+ final double[] yI = y[i];
+ for (int k = 0; k < kWidth; ++k) {
+ yI[k] -= yJ[k] * rIJ;
+ }
+ }
+ }
+ }
+ //return result;
+ return new BlockRealMatrix(cols, columns, xBlocks, false);
+ }
+
+ /** {@inheritDoc} */
+ public RealMatrix getInverse() {
+ return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
+ }
+ }
+}
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java?rev=1179935&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRDecompositionTest.java Fri Oct 7 05:21:17 2011
@@ -0,0 +1,257 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.linear;
+
+import java.util.Random;
+
+
+import org.apache.commons.math.ConvergenceException;
+import org.junit.Assert;
+import org.junit.Test;
+
+
+public class PivotingQRDecompositionTest {
+ double[][] testData3x3NonSingular = {
+ { 12, -51, 4 },
+ { 6, 167, -68 },
+ { -4, 24, -41 }, };
+
+ double[][] testData3x3Singular = {
+ { 1, 4, 7, },
+ { 2, 5, 8, },
+ { 3, 6, 9, }, };
+
+ double[][] testData3x4 = {
+ { 12, -51, 4, 1 },
+ { 6, 167, -68, 2 },
+ { -4, 24, -41, 3 }, };
+
+ double[][] testData4x3 = {
+ { 12, -51, 4, },
+ { 6, 167, -68, },
+ { -4, 24, -41, },
+ { -5, 34, 7, }, };
+
+ private static final double entryTolerance = 10e-16;
+
+ private static final double normTolerance = 10e-14;
+
+ /** test dimensions */
+ @Test
+ public void testDimensions() throws ConvergenceException {
+ checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+
+ checkDimension(MatrixUtils.createRealMatrix(testData4x3));
+
+ checkDimension(MatrixUtils.createRealMatrix(testData3x4));
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ checkDimension(createTestMatrix(r, p, q));
+ checkDimension(createTestMatrix(r, q, p));
+
+ }
+
+ private void checkDimension(RealMatrix m) throws ConvergenceException {
+ int rows = m.getRowDimension();
+ int columns = m.getColumnDimension();
+ PivotingQRDecomposition qr = new PivotingQRDecomposition(m);
+ Assert.assertEquals(rows, qr.getQ().getRowDimension());
+ Assert.assertEquals(rows, qr.getQ().getColumnDimension());
+ Assert.assertEquals(rows, qr.getR().getRowDimension());
+ Assert.assertEquals(columns, qr.getR().getColumnDimension());
+ }
+
+ /** test A = QR */
+ @Test
+ public void testAEqualQR() throws ConvergenceException {
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
+
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4));
+
+ checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3));
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ checkAEqualQR(createTestMatrix(r, p, q));
+
+ checkAEqualQR(createTestMatrix(r, q, p));
+
+ }
+
+ private void checkAEqualQR(RealMatrix m) throws ConvergenceException {
+ PivotingQRDecomposition qr = new PivotingQRDecomposition(m);
+ RealMatrix prod = qr.getQ().multiply(qr.getR()).multiply(qr.getPermutationMatrix().transpose());
+ double norm = prod.subtract(m).getNorm();
+ Assert.assertEquals(0, norm, normTolerance);
+ }
+
+ /** test the orthogonality of Q */
+ @Test
+ public void testQOrthogonal() throws ConvergenceException{
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular));
+
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4));
+
+ checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3));
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ checkQOrthogonal(createTestMatrix(r, p, q));
+
+ checkQOrthogonal(createTestMatrix(r, q, p));
+
+ }
+
+ private void checkQOrthogonal(RealMatrix m) throws ConvergenceException{
+ PivotingQRDecomposition qr = new PivotingQRDecomposition(m);
+ RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
+ double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm();
+ Assert.assertEquals(0, norm, normTolerance);
+ }
+//
+ /** test that R is upper triangular */
+ @Test
+ public void testRUpperTriangular() throws ConvergenceException{
+ RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x4);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = MatrixUtils.createRealMatrix(testData4x3);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ matrix = createTestMatrix(r, p, q);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ matrix = createTestMatrix(r, p, q);
+ checkUpperTriangular(new PivotingQRDecomposition(matrix).getR());
+
+ }
+
+ private void checkUpperTriangular(RealMatrix m) {
+ m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
+ @Override
+ public void visit(int row, int column, double value) {
+ if (column < row) {
+ Assert.assertEquals(0.0, value, entryTolerance);
+ }
+ }
+ });
+ }
+
+ /** test that H is trapezoidal */
+ @Test
+ public void testHTrapezoidal() throws ConvergenceException{
+ RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = MatrixUtils.createRealMatrix(testData3x4);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = MatrixUtils.createRealMatrix(testData4x3);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ Random r = new Random(643895747384642l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ matrix = createTestMatrix(r, p, q);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ matrix = createTestMatrix(r, p, q);
+ checkTrapezoidal(new PivotingQRDecomposition(matrix).getH());
+
+ }
+
+ private void checkTrapezoidal(RealMatrix m) {
+ m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
+ @Override
+ public void visit(int row, int column, double value) {
+ if (column > row) {
+ Assert.assertEquals(0.0, value, entryTolerance);
+ }
+ }
+ });
+ }
+ /** test matrices values */
+ @Test
+ public void testMatricesValues() throws ConvergenceException{
+ PivotingQRDecomposition qr =
+ new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular),false);
+ RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] {
+ { -12.0 / 14.0, 69.0 / 175.0, -58.0 / 175.0 },
+ { -6.0 / 14.0, -158.0 / 175.0, 6.0 / 175.0 },
+ { 4.0 / 14.0, -30.0 / 175.0, -165.0 / 175.0 }
+ });
+ RealMatrix rRef = MatrixUtils.createRealMatrix(new double[][] {
+ { -14.0, -21.0, 14.0 },
+ { 0.0, -175.0, 70.0 },
+ { 0.0, 0.0, 35.0 }
+ });
+ RealMatrix hRef = MatrixUtils.createRealMatrix(new double[][] {
+ { 26.0 / 14.0, 0.0, 0.0 },
+ { 6.0 / 14.0, 648.0 / 325.0, 0.0 },
+ { -4.0 / 14.0, 36.0 / 325.0, 2.0 }
+ });
+
+ // check values against known references
+ RealMatrix q = qr.getQ();
+ Assert.assertEquals(0, q.subtract(qRef).getNorm(), 1.0e-13);
+ RealMatrix qT = qr.getQT();
+ Assert.assertEquals(0, qT.subtract(qRef.transpose()).getNorm(), 1.0e-13);
+ RealMatrix r = qr.getR();
+ Assert.assertEquals(0, r.subtract(rRef).getNorm(), 1.0e-13);
+ RealMatrix h = qr.getH();
+ Assert.assertEquals(0, h.subtract(hRef).getNorm(), 1.0e-13);
+
+ // check the same cached instance is returned the second time
+ Assert.assertTrue(q == qr.getQ());
+ Assert.assertTrue(r == qr.getR());
+ Assert.assertTrue(h == qr.getH());
+
+ }
+
+ private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
+ RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
+ m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
+ @Override
+ public double visit(int row, int column, double value) {
+ return 2.0 * r.nextDouble() - 1.0;
+ }
+ });
+ return m;
+ }
+
+}
Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java?rev=1179935&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/PivotingQRSolverTest.java Fri Oct 7 05:21:17 2011
@@ -0,0 +1,201 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.linear;
+
+import java.util.Random;
+
+import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.exception.MathIllegalArgumentException;
+
+import org.junit.Test;
+import org.junit.Assert;
+
+public class PivotingQRSolverTest {
+ double[][] testData3x3NonSingular = {
+ { 12, -51, 4 },
+ { 6, 167, -68 },
+ { -4, 24, -41 }
+ };
+
+ double[][] testData3x3Singular = {
+ { 1, 2, 2 },
+ { 2, 4, 6 },
+ { 4, 8, 12 }
+ };
+
+ double[][] testData3x4 = {
+ { 12, -51, 4, 1 },
+ { 6, 167, -68, 2 },
+ { -4, 24, -41, 3 }
+ };
+
+ double[][] testData4x3 = {
+ { 12, -51, 4 },
+ { 6, 167, -68 },
+ { -4, 24, -41 },
+ { -5, 34, 7 }
+ };
+
+ /** test rank */
+ @Test
+ public void testRank() throws ConvergenceException {
+ DecompositionSolver solver =
+ new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ solver = new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
+ Assert.assertFalse(solver.isNonSingular());
+
+ solver = new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ solver = new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ }
+
+ /** test solve dimension errors */
+ @Test
+ public void testSolveDimensionErrors() throws ConvergenceException {
+ DecompositionSolver solver =
+ new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
+ RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
+ try {
+ solver.solve(b);
+ Assert.fail("an exception should have been thrown");
+ } catch (MathIllegalArgumentException iae) {
+ // expected behavior
+ }
+ try {
+ solver.solve(b.getColumnVector(0));
+ Assert.fail("an exception should have been thrown");
+ } catch (MathIllegalArgumentException iae) {
+ // expected behavior
+ }
+ }
+
+ /** test solve rank errors */
+ @Test
+ public void testSolveRankErrors() throws ConvergenceException {
+ DecompositionSolver solver =
+ new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
+ RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]);
+ try {
+ solver.solve(b);
+ Assert.fail("an exception should have been thrown");
+ } catch (SingularMatrixException iae) {
+ // expected behavior
+ }
+ try {
+ solver.solve(b.getColumnVector(0));
+ Assert.fail("an exception should have been thrown");
+ } catch (SingularMatrixException iae) {
+ // expected behavior
+ }
+ }
+
+ /** test solve */
+ @Test
+ public void testSolve() throws ConvergenceException {
+ PivotingQRDecomposition decomposition =
+ new PivotingQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+ DecompositionSolver solver = decomposition.getSolver();
+ RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
+ { -102, 12250 }, { 544, 24500 }, { 167, -36750 }
+ });
+
+ RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
+ { 1, 2515 }, { 2, 422 }, { -3, 898 }
+ });
+
+ // using RealMatrix
+ Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 2.0e-14 * xRef.getNorm());
+
+ // using ArrayRealVector
+ for (int i = 0; i < b.getColumnDimension(); ++i) {
+ final RealVector x = solver.solve(b.getColumnVector(i));
+ final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
+ Assert.assertEquals(0, error, 3.0e-14 * xRef.getColumnVector(i).getNorm());
+ }
+
+ // using RealVector with an alternate implementation
+ for (int i = 0; i < b.getColumnDimension(); ++i) {
+ ArrayRealVectorTest.RealVectorTestImpl v =
+ new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
+ final RealVector x = solver.solve(v);
+ final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
+ Assert.assertEquals(0, error, 3.0e-14 * xRef.getColumnVector(i).getNorm());
+ }
+
+ }
+
+ @Test
+ public void testOverdetermined() throws ConvergenceException {
+ final Random r = new Random(5559252868205245l);
+ int p = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ RealMatrix a = createTestMatrix(r, p, q);
+ RealMatrix xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);
+
+ // build a perturbed system: A.X + noise = B
+ RealMatrix b = a.multiply(xRef);
+ final double noise = 0.001;
+ b.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
+ @Override
+ public double visit(int row, int column, double value) {
+ return value * (1.0 + noise * (2 * r.nextDouble() - 1));
+ }
+ });
+
+ // despite perturbation, the least square solution should be pretty good
+ RealMatrix x = new PivotingQRDecomposition(a).getSolver().solve(b);
+ Assert.assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p * q);
+
+ }
+
+ @Test
+ public void testUnderdetermined() throws ConvergenceException {
+ final Random r = new Random(42185006424567123l);
+ int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
+ RealMatrix a = createTestMatrix(r, p, q);
+ RealMatrix xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);
+ RealMatrix b = a.multiply(xRef);
+ PivotingQRDecomposition pqr = new PivotingQRDecomposition(a);
+ RealMatrix x = pqr.getSolver().solve(b);
+ Assert.assertTrue(x.subtract(xRef).getNorm() / (p * q) > 0.01);
+ int count=0;
+ for( int i = 0 ; i < q; i++){
+ if( x.getRowVector(i).getNorm() == 0.0 ){
+ ++count;
+ }
+ }
+ Assert.assertEquals("Zeroed rows", q-p, count);
+ }
+
+ private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
+ RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
+ m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
+ @Override
+ public double visit(int row, int column, double value) {
+ return 2.0 * r.nextDouble() - 1.0;
+ }
+ });
+ return m;
+ }
+}
Re: svn commit: r1179935 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/linear/ test/java/org/apache/commons/math/linear/
Posted by Greg Sterijevski <gs...@gmail.com>.
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 <gs...@gmail.com>wrote:
> Will do. My aplogies! -Greg
>
>
> On Fri, Oct 7, 2011 at 3:55 AM, Luc Maisonobe <Lu...@free.fr>wrote:
>
>> 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<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<de...@commons.apache.org>
>> For additional commands, e-mail: dev-help@commons.apache.org
>>
>>
>
Re: svn commit: r1179935 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/linear/ test/java/org/apache/commons/math/linear/
Posted by Greg Sterijevski <gs...@gmail.com>.
Will do. My aplogies! -Greg
On Fri, Oct 7, 2011 at 3:55 AM, Luc Maisonobe <Lu...@free.fr> wrote:
> 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<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<de...@commons.apache.org>
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
Re: svn commit: r1179935 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/linear/
test/java/org/apache/commons/math/linear/
Posted by Luc Maisonobe <Lu...@free.fr>.
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