Added rank revealing QR decomposition.
Patch applied after conversion to current status and slight adaptations. JIRA: MATH-630 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1455627 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
62f86d9580
commit
e1f1860894
|
|
@ -55,6 +55,9 @@ This is a minor release: It combines bug fixes and new features.
|
||||||
Changes to existing features were made in a backwards-compatible
|
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.
|
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" >
|
<action dev="luc" type="fix" issue="MATH-570" due-to="Arne Plöse" >
|
||||||
ArrayFieldVector can now be constructed from any FieldVector.
|
ArrayFieldVector can now be constructed from any FieldVector.
|
||||||
</action>
|
</action>
|
||||||
|
|
|
||||||
|
|
@ -100,12 +100,24 @@ public class QRDecomposition {
|
||||||
cachedR = null;
|
cachedR = null;
|
||||||
cachedH = null;
|
cachedH = null;
|
||||||
|
|
||||||
/*
|
decompose(qrt);
|
||||||
* 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:
|
|
||||||
|
/** Decompose matrix.
|
||||||
|
* @param qrt transposed matrix
|
||||||
*/
|
*/
|
||||||
for (int minor = 0; minor < FastMath.min(m, n); minor++) {
|
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];
|
final double[] qrtMinor = qrt[minor];
|
||||||
|
|
||||||
|
|
@ -117,7 +129,7 @@ public class QRDecomposition {
|
||||||
* component of x. Let's find a:
|
* component of x. Let's find a:
|
||||||
*/
|
*/
|
||||||
double xNormSqr = 0;
|
double xNormSqr = 0;
|
||||||
for (int row = minor; row < m; row++) {
|
for (int row = minor; row < qrtMinor.length; row++) {
|
||||||
final double c = qrtMinor[row];
|
final double c = qrtMinor[row];
|
||||||
xNormSqr += c * c;
|
xNormSqr += c * c;
|
||||||
}
|
}
|
||||||
|
|
@ -148,22 +160,22 @@ public class QRDecomposition {
|
||||||
* |v|^2 = -2a*(qr[minor][minor]), so
|
* |v|^2 = -2a*(qr[minor][minor]), so
|
||||||
* alpha = -<x,v>/(a*qr[minor][minor])
|
* alpha = -<x,v>/(a*qr[minor][minor])
|
||||||
*/
|
*/
|
||||||
for (int col = minor+1; col < n; col++) {
|
for (int col = minor+1; col < qrt.length; col++) {
|
||||||
final double[] qrtCol = qrt[col];
|
final double[] qrtCol = qrt[col];
|
||||||
double alpha = 0;
|
double alpha = 0;
|
||||||
for (int row = minor; row < m; row++) {
|
for (int row = minor; row < qrtCol.length; row++) {
|
||||||
alpha -= qrtCol[row] * qrtMinor[row];
|
alpha -= qrtCol[row] * qrtMinor[row];
|
||||||
}
|
}
|
||||||
alpha /= a * qrtMinor[minor];
|
alpha /= a * qrtMinor[minor];
|
||||||
|
|
||||||
// Subtract the column vector alpha*v from x.
|
// Subtract the column vector alpha*v from x.
|
||||||
for (int row = minor; row < m; row++) {
|
for (int row = minor; row < qrtCol.length; row++) {
|
||||||
qrtCol[row] -= alpha * qrtMinor[row];
|
qrtCol[row] -= alpha * qrtMinor[row];
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Returns the matrix R of the decomposition.
|
* Returns the matrix R of the decomposition.
|
||||||
|
|
|
||||||
|
|
@ -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()));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
@ -245,8 +245,7 @@ public class QRDecompositionTest {
|
||||||
public void testNonInvertible() {
|
public void testNonInvertible() {
|
||||||
QRDecomposition qr =
|
QRDecomposition qr =
|
||||||
new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular));
|
new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular));
|
||||||
|
qr.getSolver().getInverse();
|
||||||
final RealMatrix inv = qr.getSolver().getInverse();
|
|
||||||
}
|
}
|
||||||
|
|
||||||
private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
|
private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,233 @@
|
||||||
|
/*
|
||||||
|
* 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.junit.Assert;
|
||||||
|
import org.junit.Test;
|
||||||
|
|
||||||
|
|
||||||
|
public class RRQRDecompositionTest {
|
||||||
|
private double[][] testData3x3NonSingular = {
|
||||||
|
{ 12, -51, 4 },
|
||||||
|
{ 6, 167, -68 },
|
||||||
|
{ -4, 24, -41 }, };
|
||||||
|
|
||||||
|
private double[][] testData3x3Singular = {
|
||||||
|
{ 1, 4, 7, },
|
||||||
|
{ 2, 5, 8, },
|
||||||
|
{ 3, 6, 9, }, };
|
||||||
|
|
||||||
|
private double[][] testData3x4 = {
|
||||||
|
{ 12, -51, 4, 1 },
|
||||||
|
{ 6, 167, -68, 2 },
|
||||||
|
{ -4, 24, -41, 3 }, };
|
||||||
|
|
||||||
|
private 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() {
|
||||||
|
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) {
|
||||||
|
int rows = m.getRowDimension();
|
||||||
|
int columns = m.getColumnDimension();
|
||||||
|
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 AP = QR */
|
||||||
|
@Test
|
||||||
|
public void testAPEqualQR() {
|
||||||
|
checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
|
||||||
|
|
||||||
|
checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
|
||||||
|
|
||||||
|
checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x4));
|
||||||
|
|
||||||
|
checkAPEqualQR(MatrixUtils.createRealMatrix(testData4x3));
|
||||||
|
|
||||||
|
Random r = new Random(643895747384642l);
|
||||||
|
int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
|
||||||
|
int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
|
||||||
|
checkAPEqualQR(createTestMatrix(r, p, q));
|
||||||
|
|
||||||
|
checkAPEqualQR(createTestMatrix(r, q, p));
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
|
||||||
|
/** test the orthogonality of Q */
|
||||||
|
@Test
|
||||||
|
public void testQOrthogonal() {
|
||||||
|
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) {
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
|
||||||
|
/** test that R is upper triangular */
|
||||||
|
@Test
|
||||||
|
public void testRUpperTriangular() {
|
||||||
|
RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
|
||||||
|
checkUpperTriangular(new RRQRDecomposition(matrix).getR());
|
||||||
|
|
||||||
|
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
|
||||||
|
checkUpperTriangular(new RRQRDecomposition(matrix).getR());
|
||||||
|
|
||||||
|
matrix = MatrixUtils.createRealMatrix(testData3x4);
|
||||||
|
checkUpperTriangular(new RRQRDecomposition(matrix).getR());
|
||||||
|
|
||||||
|
matrix = MatrixUtils.createRealMatrix(testData4x3);
|
||||||
|
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 RRQRDecomposition(matrix).getR());
|
||||||
|
|
||||||
|
matrix = createTestMatrix(r, p, q);
|
||||||
|
checkUpperTriangular(new RRQRDecomposition(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() {
|
||||||
|
RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
|
||||||
|
checkTrapezoidal(new RRQRDecomposition(matrix).getH());
|
||||||
|
|
||||||
|
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
|
||||||
|
checkTrapezoidal(new RRQRDecomposition(matrix).getH());
|
||||||
|
|
||||||
|
matrix = MatrixUtils.createRealMatrix(testData3x4);
|
||||||
|
checkTrapezoidal(new RRQRDecomposition(matrix).getH());
|
||||||
|
|
||||||
|
matrix = MatrixUtils.createRealMatrix(testData4x3);
|
||||||
|
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 RRQRDecomposition(matrix).getH());
|
||||||
|
|
||||||
|
matrix = createTestMatrix(r, p, q);
|
||||||
|
checkTrapezoidal(new RRQRDecomposition(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(expected=SingularMatrixException.class)
|
||||||
|
public void testNonInvertible() {
|
||||||
|
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) {
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
|
||||||
|
/** 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));
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
@ -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;
|
||||||
|
}
|
||||||
|
}
|
||||||
Loading…
Reference in New Issue