Author: luc
Date: Tue Mar 12 17:08:52 2013
New Revision: 1455627
URL: http://svn.apache.org/r1455627
Log:
Added rank revealing QR decomposition.
Patch applied after conversion to current status and slight adaptations.
JIRA: MATH-630
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RRQRDecomposition.java
(with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRDecompositionTest.java
- copied, changed from r1455260,
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRSolverTest.java
(with props)
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/QRDecomposition.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1455627&r1=1455626&r2=1455627&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Tue Mar 12 17:08:52 2013
@@ -55,6 +55,9 @@ This is a minor release: It combines bug
Changes to existing features were made in a backwards-compatible
way such as to allow drop-in replacement of the v3.1[.1] JAR file.
">
+ <action dev="luc" type="fix" issue="MATH-630" due-to="Christopher Nix" >
+ Added rank revealing QR decomposition.
+ </action>
<action dev="luc" type="fix" issue="MATH-570" due-to="Arne Plöse" >
ArrayFieldVector can now be constructed from any FieldVector.
</action>
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/QRDecomposition.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/QRDecomposition.java?rev=1455627&r1=1455626&r2=1455627&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/QRDecomposition.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/QRDecomposition.java
Tue Mar 12 17:08:52 2013
@@ -100,71 +100,83 @@ public class QRDecomposition {
cachedR = null;
cachedH = null;
+ decompose(qrt);
+
+ }
+
+ /** Decompose matrix.
+ * @param qrt transposed matrix
+ */
+ protected void decompose(double[][] qrt) {
+ for (int minor = 0; minor < FastMath.min(qrt.length, qrt[0].length);
minor++) {
+ performHouseholderReflection(minor, qrt);
+ }
+ }
+
+ /** Perform Householder reflection for a minor A(minor, minor) of A.
+ * @param minor minor index
+ * @param qrt transposed matrix
+ */
+ protected void performHouseholderReflection(int minor, double[][] qrt) {
+
+ final double[] qrtMinor = qrt[minor];
+
/*
- * The QR decomposition of a matrix A is calculated using Householder
- * reflectors by repeating the following operations to each minor
- * A(minor,minor) of A:
+ * Let x be the first column of the minor, and a^2 = |x|^2.
+ * x will be in the positions qr[minor][minor] through qr[m][minor].
+ * The first column of the transformed minor will be (a,0,0,..)'
+ * The sign of a is chosen to be opposite to the sign of the first
+ * component of x. Let's find a:
*/
- for (int minor = 0; minor < FastMath.min(m, n); minor++) {
+ double xNormSqr = 0;
+ for (int row = minor; row < qrtMinor.length; row++) {
+ final double c = qrtMinor[row];
+ xNormSqr += c * c;
+ }
+ final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) :
FastMath.sqrt(xNormSqr);
+ rDiag[minor] = a;
- final double[] qrtMinor = qrt[minor];
+ if (a != 0.0) {
/*
- * Let x be the first column of the minor, and a^2 = |x|^2.
- * x will be in the positions qr[minor][minor] through
qr[m][minor].
- * The first column of the transformed minor will be (a,0,0,..)'
- * The sign of a is chosen to be opposite to the sign of the first
- * component of x. Let's find a:
+ * Calculate the normalized reflection vector v and transform
+ * the first column. We know the norm of v beforehand: v = x-ae
+ * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
+ * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
+ * Here <x, e> is now qr[minor][minor].
+ * v = x-ae is stored in the column at qr:
*/
- double xNormSqr = 0;
- for (int row = minor; row < m; row++) {
- final double c = qrtMinor[row];
- xNormSqr += c * c;
- }
- final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr)
: FastMath.sqrt(xNormSqr);
- rDiag[minor] = a;
-
- if (a != 0.0) {
-
- /*
- * Calculate the normalized reflection vector v and transform
- * the first column. We know the norm of v beforehand: v = x-ae
- * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
- * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
- * Here <x, e> is now qr[minor][minor].
- * v = x-ae is stored in the column at qr:
- */
- qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
-
- /*
- * Transform the rest of the columns of the minor:
- * They will be transformed by the matrix H = I-2vv'/|v|^2.
- * If x is a column vector of the minor, then
- * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
- * Therefore the transformation is easily calculated by
- * subtracting the column vector (2<x,v>/|v|^2)v from x.
- *
- * Let 2<x,v>/|v|^2 = alpha. From above we have
- * |v|^2 = -2a*(qr[minor][minor]), so
- * alpha = -<x,v>/(a*qr[minor][minor])
- */
- for (int col = minor+1; col < n; col++) {
- final double[] qrtCol = qrt[col];
- double alpha = 0;
- for (int row = minor; row < m; row++) {
- alpha -= qrtCol[row] * qrtMinor[row];
- }
- alpha /= a * qrtMinor[minor];
+ qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor])
- // Subtract the column vector alpha*v from x.
- for (int row = minor; row < m; row++) {
- qrtCol[row] -= alpha * qrtMinor[row];
- }
+ /*
+ * Transform the rest of the columns of the minor:
+ * They will be transformed by the matrix H = I-2vv'/|v|^2.
+ * If x is a column vector of the minor, then
+ * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v.
+ * Therefore the transformation is easily calculated by
+ * subtracting the column vector (2<x,v>/|v|^2)v from x.
+ *
+ * Let 2<x,v>/|v|^2 = alpha. From above we have
+ * |v|^2 = -2a*(qr[minor][minor]), so
+ * alpha = -<x,v>/(a*qr[minor][minor])
+ */
+ for (int col = minor+1; col < qrt.length; col++) {
+ final double[] qrtCol = qrt[col];
+ double alpha = 0;
+ for (int row = minor; row < qrtCol.length; row++) {
+ alpha -= qrtCol[row] * qrtMinor[row];
+ }
+ alpha /= a * qrtMinor[minor];
+
+ // Subtract the column vector alpha*v from x.
+ for (int row = minor; row < qrtCol.length; row++) {
+ qrtCol[row] -= alpha * qrtMinor[row];
}
}
}
}
+
/**
* Returns the matrix R of the decomposition.
* <p>R is an upper-triangular matrix</p>
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RRQRDecomposition.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RRQRDecomposition.java?rev=1455627&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RRQRDecomposition.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RRQRDecomposition.java
Tue Mar 12 17:08:52 2013
@@ -0,0 +1,230 @@
+/*
+ * 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.math3.linear;
+
+import org.apache.commons.math3.util.FastMath;
+
+
+/**
+ * Calculates the rank-revealing QR-decomposition of a matrix, with column
pivoting.
+ * <p>The rank-revealing QR-decomposition of a matrix A consists of three
matrices Q,
+ * R and P such that AP=QR. Q is orthogonal (Q<sup>T</sup>Q = I), and R is
upper triangular.
+ * If A is m×n, Q is m×m and R is m×n and P is n×n.</p>
+ * <p>QR decomposition with column pivoting produces a rank-revealing QR
+ * decomposition and the {@link #getRank(double)} method may be used to return
the rank of the
+ * input matrix A.</p>
+ * <p>This class compute the decomposition using Householder reflectors.</p>
+ * <p>For efficiency purposes, the decomposition in packed form is transposed.
+ * This allows inner loop to iterate inside rows, which is much more
cache-efficient
+ * in Java.</p>
+ * <p>This class is based on the class with similar name from the
+ * <a href="http://math.nist.gov/javanumerics/jama/";>JAMA</a> library, with the
+ * following changes:</p>
+ * <ul>
+ * <li>a {@link #getQT() getQT} method has been added,</li>
+ * <li>the {@code solve} and {@code isFullRank} methods have been replaced
+ * by a {@link #getSolver() getSolver} method and the equivalent methods
+ * provided by the returned {@link DecompositionSolver}.</li>
+ * </ul>
+ *
+ * @see <a
href="http://mathworld.wolfram.com/QRDecomposition.html";>MathWorld</a>
+ * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition";>Wikipedia</a>
+ *
+ * @version $Id$
+ * @since 3.2
+ */
+public class RRQRDecomposition extends QRDecomposition {
+
+ /** An array to record the column pivoting for later creation of P. */
+ private int[] p;
+
+ /** Cached value of P. */
+ private RealMatrix cachedP;
+
+
+ /**
+ * Calculates the QR-decomposition of the given matrix.
+ * The singularity threshold defaults to zero.
+ *
+ * @param matrix The matrix to decompose.
+ *
+ * @see #QRDecomposition(RealMatrix,double)
+ */
+ public RRQRDecomposition(RealMatrix matrix) {
+ this(matrix, 0d);
+ }
+
+ /**
+ * Calculates the QR-decomposition of the given matrix.
+ *
+ * @param matrix The matrix to decompose.
+ * @param threshold Singularity threshold.
+ */
+ public RRQRDecomposition(RealMatrix matrix, double threshold) {
+ super(matrix, threshold);
+ }
+
+ /** Decompose matrix.
+ * @param qrt transposed matrix
+ */
+ protected void decompose(double[][] qrt) {
+ p = new int[qrt.length];
+ for (int i = 0; i < p.length; i++) {
+ p[i] = i;
+ }
+ super.decompose(qrt);
+ }
+
+ /** Perform Householder reflection for a minor A(minor, minor) of A.
+ * @param minor minor index
+ * @param qrt transposed matrix
+ */
+ protected void performHouseholderReflection(int minor, double[][] qrt) {
+
+ double l2NormSquaredMax = 0;
+ // Find the unreduced column with the greatest L2-Norm
+ int l2NormSquaredMaxIndex = minor;
+ for (int i = minor; i < qrt.length; i++) {
+ double l2NormSquared = 0;
+ for (int j = 0; j < qrt[i].length; j++) {
+ l2NormSquared += qrt[i][j] * qrt[i][j];
+ }
+ if (l2NormSquared > l2NormSquaredMax) {
+ l2NormSquaredMax = l2NormSquared;
+ l2NormSquaredMaxIndex = i;
+ }
+ }
+ // swap the current column with that with the greated L2-Norm and
record in p
+ if (l2NormSquaredMaxIndex != minor) {
+ double[] tmp1 = qrt[minor];
+ qrt[minor] = qrt[l2NormSquaredMaxIndex];
+ qrt[l2NormSquaredMaxIndex] = tmp1;
+ int tmp2 = p[minor];
+ p[minor] = p[l2NormSquaredMaxIndex];
+ p[l2NormSquaredMaxIndex] = tmp2;
+ }
+
+ super.performHouseholderReflection(minor, qrt);
+
+ }
+
+
+ /**
+ * Returns the pivot matrix, P, used in the QR Decomposition of matrix A
such that AP = QR.
+ *
+ * If no pivoting is used in this decomposition then P is equal to the
identity matrix.
+ *
+ * @return a permutation matrix.
+ */
+ public RealMatrix getP() {
+ if (cachedP == null) {
+ int n = p.length;
+ cachedP = MatrixUtils.createRealMatrix(n,n);
+ for (int i = 0; i < n; i++) {
+ cachedP.setEntry(p[i], i, 1);
+ }
+ }
+ return cachedP ;
+ }
+
+ /**
+ * Return the effective numerical matrix rank.
+ * <p>The effective numerical rank is the number of non-negligible
+ * singular values.</p>
+ * <p>This implementation looks at Frobenius norms of the sequence of
+ * bottom right submatrices. When a large fall in norm is seen,
+ * the rank is returned. The drop is computed as:</p>
+ * <pre>
+ * (thisNorm/lastNorm) * rNorm < dropThreshold
+ * </pre>
+ * <p>
+ * where thisNorm is the Frobenius norm of the current submatrix,
+ * lastNorm is the Frobenius norm of the previous submatrix,
+ * rNorm is is the Frobenius norm of the complete matrix
+ * </p>
+ *
+ * @param dropThreshold threshold triggering rank computation
+ * @return effective numerical matrix rank
+ */
+ public int getRank(final double dropThreshold) {
+ RealMatrix r = getR();
+ int rows = r.getRowDimension();
+ int columns = r.getColumnDimension();
+ int rank = 1;
+ double lastNorm = r.getFrobeniusNorm();
+ double rNorm = lastNorm;
+ while (rank < FastMath.min(rows, columns)) {
+ double thisNorm = r.getSubMatrix(rank, rows - 1, rank, columns -
1).getFrobeniusNorm();
+ if (thisNorm == 0 || (thisNorm / lastNorm) * rNorm <
dropThreshold) {
+ break;
+ }
+ lastNorm = thisNorm;
+ rank++;
+ }
+ return rank;
+ }
+
+ /**
+ * Get a solver for finding the A × X = B solution in least square
sense.
+ * @return a solver
+ */
+ public DecompositionSolver getSolver() {
+ return new Solver(super.getSolver(), this.getP());
+ }
+
+ /** Specialized solver. */
+ private static class Solver implements DecompositionSolver {
+
+ /** Upper level solver. */
+ private final DecompositionSolver upper;
+
+ /** A permutation matrix for the pivots used in the QR decomposition */
+ private RealMatrix p;
+
+ /**
+ * Build a solver from decomposed matrix.
+ *
+ * @param upper upper level solver.
+ * @param p permutation matrix
+ */
+ private Solver(final DecompositionSolver upper, final RealMatrix p) {
+ this.upper = upper;
+ this.p = p;
+ }
+
+ /** {@inheritDoc} */
+ public boolean isNonSingular() {
+ return upper.isNonSingular();
+ }
+
+ /** {@inheritDoc} */
+ public RealVector solve(RealVector b) {
+ return p.operate(upper.solve(b));
+ }
+
+ /** {@inheritDoc} */
+ public RealMatrix solve(RealMatrix b) {
+ return p.multiply(upper.solve(b));
+ }
+
+ /** {@inheritDoc} */
+ public RealMatrix getInverse() {
+ return
solve(MatrixUtils.createRealIdentityMatrix(p.getRowDimension()));
+ }
+ }
+}
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RRQRDecomposition.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RRQRDecomposition.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java?rev=1455627&r1=1455626&r2=1455627&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java
Tue Mar 12 17:08:52 2013
@@ -245,8 +245,7 @@ public class QRDecompositionTest {
public void testNonInvertible() {
QRDecomposition qr =
new
QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular));
-
- final RealMatrix inv = qr.getSolver().getInverse();
+ qr.getSolver().getInverse();
}
private RealMatrix createTestMatrix(final Random r, final int rows, final
int columns) {
Copied:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRDecompositionTest.java
(from r1455260,
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java)
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRDecompositionTest.java?p2=commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRDecompositionTest.java&p1=commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java&r1=1455260&r2=1455627&rev=1455627&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/QRDecompositionTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRDecompositionTest.java
Tue Mar 12 17:08:52 2013
@@ -18,13 +18,12 @@
package org.apache.commons.math3.linear;
import java.util.Random;
-import org.apache.commons.math3.linear.SingularMatrixException;
import org.junit.Assert;
import org.junit.Test;
-public class QRDecompositionTest {
+public class RRQRDecompositionTest {
private double[][] testData3x3NonSingular = {
{ 12, -51, 4 },
{ 6, 167, -68 },
@@ -70,36 +69,36 @@ public class QRDecompositionTest {
private void checkDimension(RealMatrix m) {
int rows = m.getRowDimension();
int columns = m.getColumnDimension();
- QRDecomposition qr = new QRDecomposition(m);
+ RRQRDecomposition qr = new RRQRDecomposition(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 AP = QR */
@Test
- public void testAEqualQR() {
- checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+ public void testAPEqualQR() {
+ checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
- checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
+ checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
- checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4));
+ checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x4));
- checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3));
+ checkAPEqualQR(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));
+ checkAPEqualQR(createTestMatrix(r, p, q));
- checkAEqualQR(createTestMatrix(r, q, p));
+ checkAPEqualQR(createTestMatrix(r, q, p));
}
- private void checkAEqualQR(RealMatrix m) {
- QRDecomposition qr = new QRDecomposition(m);
- double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm();
+ private void checkAPEqualQR(RealMatrix m) {
+ RRQRDecomposition rrqr = new RRQRDecomposition(m);
+ double norm =
rrqr.getQ().multiply(rrqr.getR()).subtract(m.multiply(rrqr.getP())).getNorm();
Assert.assertEquals(0, norm, normTolerance);
}
@@ -124,7 +123,7 @@ public class QRDecompositionTest {
}
private void checkQOrthogonal(RealMatrix m) {
- QRDecomposition qr = new QRDecomposition(m);
+ RRQRDecomposition qr = new RRQRDecomposition(m);
RealMatrix eye =
MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm();
Assert.assertEquals(0, norm, normTolerance);
@@ -134,25 +133,25 @@ public class QRDecompositionTest {
@Test
public void testRUpperTriangular() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(testData3x3NonSingular);
- checkUpperTriangular(new QRDecomposition(matrix).getR());
+ checkUpperTriangular(new RRQRDecomposition(matrix).getR());
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
- checkUpperTriangular(new QRDecomposition(matrix).getR());
+ checkUpperTriangular(new RRQRDecomposition(matrix).getR());
matrix = MatrixUtils.createRealMatrix(testData3x4);
- checkUpperTriangular(new QRDecomposition(matrix).getR());
+ checkUpperTriangular(new RRQRDecomposition(matrix).getR());
matrix = MatrixUtils.createRealMatrix(testData4x3);
- checkUpperTriangular(new QRDecomposition(matrix).getR());
+ checkUpperTriangular(new RRQRDecomposition(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 QRDecomposition(matrix).getR());
+ checkUpperTriangular(new RRQRDecomposition(matrix).getR());
matrix = createTestMatrix(r, p, q);
- checkUpperTriangular(new QRDecomposition(matrix).getR());
+ checkUpperTriangular(new RRQRDecomposition(matrix).getR());
}
@@ -171,25 +170,25 @@ public class QRDecompositionTest {
@Test
public void testHTrapezoidal() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(testData3x3NonSingular);
- checkTrapezoidal(new QRDecomposition(matrix).getH());
+ checkTrapezoidal(new RRQRDecomposition(matrix).getH());
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
- checkTrapezoidal(new QRDecomposition(matrix).getH());
+ checkTrapezoidal(new RRQRDecomposition(matrix).getH());
matrix = MatrixUtils.createRealMatrix(testData3x4);
- checkTrapezoidal(new QRDecomposition(matrix).getH());
+ checkTrapezoidal(new RRQRDecomposition(matrix).getH());
matrix = MatrixUtils.createRealMatrix(testData4x3);
- checkTrapezoidal(new QRDecomposition(matrix).getH());
+ checkTrapezoidal(new RRQRDecomposition(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 QRDecomposition(matrix).getH());
+ checkTrapezoidal(new RRQRDecomposition(matrix).getH());
matrix = createTestMatrix(r, p, q);
- checkTrapezoidal(new QRDecomposition(matrix).getH());
+ checkTrapezoidal(new RRQRDecomposition(matrix).getH());
}
@@ -203,50 +202,12 @@ public class QRDecompositionTest {
}
});
}
- /** test matrices values */
- @Test
- public void testMatricesValues() {
- QRDecomposition qr =
- new
QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
- 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());
-
- }
@Test(expected=SingularMatrixException.class)
public void testNonInvertible() {
- QRDecomposition qr =
- new
QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular));
-
- final RealMatrix inv = qr.getSolver().getInverse();
+ RRQRDecomposition qr =
+ new
RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular), 3.0e-16);
+ qr.getSolver().getInverse();
}
private RealMatrix createTestMatrix(final Random r, final int rows, final
int columns) {
@@ -260,4 +221,13 @@ public class QRDecompositionTest {
return m;
}
+ /** test the rank is returned correctly */
+ @Test
+ public void testRank() {
+ double[][] d = { { 1, 1, 1 }, { 0, 0, 0 }, { 1, 2, 3 } };
+ RealMatrix m = new Array2DRowRealMatrix(d);
+ RRQRDecomposition qr = new RRQRDecomposition(m);
+ Assert.assertEquals(2, qr.getRank(1.0e-16));
+ }
+
}
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRSolverTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRSolverTest.java?rev=1455627&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRSolverTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRSolverTest.java
Tue Mar 12 17:08:52 2013
@@ -0,0 +1,202 @@
+/*
+ * 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.math3.linear;
+
+import java.util.Random;
+
+import org.apache.commons.math3.exception.MathIllegalArgumentException;
+
+import org.junit.Test;
+import org.junit.Assert;
+
+public class RRQRSolverTest {
+ 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() {
+ DecompositionSolver solver =
+ new
RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular),
1.0e-16).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ solver = new
RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular),
1.0e-16).getSolver();
+ Assert.assertFalse(solver.isNonSingular());
+
+ solver = new
RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x4),
1.0e-16).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ solver = new
RRQRDecomposition(MatrixUtils.createRealMatrix(testData4x3),
1.0e-16).getSolver();
+ Assert.assertTrue(solver.isNonSingular());
+
+ }
+
+ /** test solve dimension errors */
+ @Test
+ public void testSolveDimensionErrors() {
+ DecompositionSolver solver =
+ new
RRQRDecomposition(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() {
+ DecompositionSolver solver =
+ new
RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular),
1.0e-16).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() {
+ RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
+ { -102, 12250 }, { 544, 24500 }, { 167, -36750 }
+ });
+ RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
+ { 1, 2515 }, { 2, 422 }, { -3, 898 }
+ });
+
+
+ RRQRDecomposition decomposition = new
RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
+ DecompositionSolver solver = decomposition.getSolver();
+
+ // using RealMatrix
+ Assert.assertEquals(0, solver.solve(b).subtract(xRef).getNorm(),
3.0e-16 * 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-16 *
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-16 *
xRef.getColumnVector(i).getNorm());
+ }
+
+ }
+
+ @Test
+ public void testOverdetermined() {
+ 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 RRQRDecomposition(a).getSolver().solve(b);
+ Assert.assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p *
q);
+
+ }
+
+ @Test
+ public void testUnderdetermined() {
+ 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);
+ RRQRDecomposition rrqrd = new RRQRDecomposition(a);
+ RealMatrix x = rrqrd.getSolver().solve(b);
+
+ // too many equations, the system cannot be solved at all
+ Assert.assertTrue(x.subtract(xRef).getNorm() / (p * q) > 0.01);
+
+ // the last permuted unknown should have been set to 0
+ RealMatrix permuted = rrqrd.getP().transpose().multiply(x);
+ Assert.assertEquals(0.0, permuted.getSubMatrix(p, q - 1, 0,
permuted.getColumnDimension() - 1).getNorm(), 0);
+
+ }
+
+ 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;
+ }
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRSolverTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RRQRSolverTest.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
