/*
* 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.linear.DefaultRealMatrixChangingVisitor;
import org.apache.commons.math.linear.DefaultRealMatrixPreservingVisitor;
import org.apache.commons.math.linear.BlockRealMatrix;
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 junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class QRDecompositionImplTest extends TestCase {
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;
public QRDecompositionImplTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(QRDecompositionImplTest.class);
suite.setName("QRDecompositionImpl Tests");
return suite;
}
/** test dimensions */
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();
QRDecomposition qr = new QRDecompositionImpl(m);
assertEquals(rows, qr.getQ().getRowDimension());
assertEquals(rows, qr.getQ().getColumnDimension());
assertEquals(rows, qr.getR().getRowDimension());
assertEquals(columns, qr.getR().getColumnDimension());
}
/** test A = QR */
public void testAEqualQR() {
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) {
QRDecomposition qr = new QRDecompositionImpl(m);
double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm();
assertEquals(0, norm, normTolerance);
}
/** test the orthogonality of Q */
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) {
QRDecomposition qr = new QRDecompositionImpl(m);
RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm();
assertEquals(0, norm, normTolerance);
}
/** test that R is upper triangular */
public void testRUpperTriangular() {
RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
matrix = MatrixUtils.createRealMatrix(testData3x4);
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
matrix = MatrixUtils.createRealMatrix(testData4x3);
checkUpperTriangular(new QRDecompositionImpl(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 QRDecompositionImpl(matrix).getR());
matrix = createTestMatrix(r, p, q);
checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
}
private void checkUpperTriangular(RealMatrix m) {
m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
@Override
public void visit(int row, int column, double value) {
if (column < row) {
assertEquals(0.0, value, entryTolerance);
}
}
});
}
/** test that H is trapezoidal */
public void testHTrapezoidal() {
RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
matrix = MatrixUtils.createRealMatrix(testData3x4);
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
matrix = MatrixUtils.createRealMatrix(testData4x3);
checkTrapezoidal(new QRDecompositionImpl(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 QRDecompositionImpl(matrix).getH());
matrix = createTestMatrix(r, p, q);
checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
}
private void checkTrapezoidal(RealMatrix m) {
m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
@Override
public void visit(int row, int column, double value) {
if (column > row) {
assertEquals(0.0, value, entryTolerance);
}
}
});
}
/** test matrices values */
public void testMatricesValues() {
QRDecomposition qr =
new QRDecompositionImpl(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();
assertEquals(0, q.subtract(qRef).getNorm(), 1.0e-13);
RealMatrix qT = qr.getQT();
assertEquals(0, qT.subtract(qRef.transpose()).getNorm(), 1.0e-13);
RealMatrix r = qr.getR();
assertEquals(0, r.subtract(rRef).getNorm(), 1.0e-13);
RealMatrix h = qr.getH();
assertEquals(0, h.subtract(hRef).getNorm(), 1.0e-13);
// check the same cached instance is returned the second time
assertTrue(q == qr.getQ());
assertTrue(r == qr.getR());
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)
throws MatrixVisitorException {
return 2.0 * r.nextDouble() - 1.0;
}
});
return m;
}
}