/*
* 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 junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
import org.apache.commons.math.linear.BlockRealMatrix;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.MatrixVisitorException;
import org.apache.commons.math.linear.QRDecomposition;
import org.apache.commons.math.linear.QRDecompositionImpl;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.ArrayRealVector;
public class QRSolverTest extends TestCase {
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 }
};
public QRSolverTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(QRSolverTest.class);
suite.setName("QRSolver Tests");
return suite;
}
/** test rank */
public void testRank() {
DecompositionSolver solver =
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
assertTrue(solver.isNonSingular());
solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
assertFalse(solver.isNonSingular());
solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
assertTrue(solver.isNonSingular());
solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
assertTrue(solver.isNonSingular());
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
DecompositionSolver solver =
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve rank errors */
public void testSolveRankErrors() {
DecompositionSolver solver =
new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
QRDecomposition decomposition =
new QRDecompositionImpl(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
assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm());
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
final double[] x = solver.solve(b.getColumn(i));
final double error = new ArrayRealVector(x).subtract(xRef.getColumnVector(i)).getNorm();
assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).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();
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();
assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm());
}
}
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 QRDecompositionImpl(a).getSolver().solve(b);
assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p * q);
}
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);
RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b);
// too many equations, the system cannot be solved at all
assertTrue(x.subtract(xRef).getNorm() / (p * q) > 0.01);
// the last unknown should have been set to 0
assertEquals(0.0, x.getSubMatrix(p, q - 1, 0, x.getColumnDimension() - 1).getNorm());
}
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)
throws MatrixVisitorException {
return 2.0 * r.nextDouble() - 1.0;
}
});
return m;
}
}