001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.linear; 019 020 import java.util.Random; 021 022 import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor; 023 import org.apache.commons.math.linear.DefaultRealMatrixPreservingVisitor; 024 import org.apache.commons.math.linear.BlockRealMatrix; 025 import org.apache.commons.math.linear.MatrixUtils; 026 import org.apache.commons.math.linear.MatrixVisitorException; 027 import org.apache.commons.math.linear.QRDecomposition; 028 import org.apache.commons.math.linear.QRDecompositionImpl; 029 import org.apache.commons.math.linear.RealMatrix; 030 031 import junit.framework.Test; 032 import junit.framework.TestCase; 033 import junit.framework.TestSuite; 034 035 public class QRDecompositionImplTest extends TestCase { 036 double[][] testData3x3NonSingular = { 037 { 12, -51, 4 }, 038 { 6, 167, -68 }, 039 { -4, 24, -41 }, }; 040 041 double[][] testData3x3Singular = { 042 { 1, 4, 7, }, 043 { 2, 5, 8, }, 044 { 3, 6, 9, }, }; 045 046 double[][] testData3x4 = { 047 { 12, -51, 4, 1 }, 048 { 6, 167, -68, 2 }, 049 { -4, 24, -41, 3 }, }; 050 051 double[][] testData4x3 = { 052 { 12, -51, 4, }, 053 { 6, 167, -68, }, 054 { -4, 24, -41, }, 055 { -5, 34, 7, }, }; 056 057 private static final double entryTolerance = 10e-16; 058 059 private static final double normTolerance = 10e-14; 060 061 public QRDecompositionImplTest(String name) { 062 super(name); 063 } 064 065 public static Test suite() { 066 TestSuite suite = new TestSuite(QRDecompositionImplTest.class); 067 suite.setName("QRDecompositionImpl Tests"); 068 return suite; 069 } 070 071 /** test dimensions */ 072 public void testDimensions() { 073 checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular)); 074 075 checkDimension(MatrixUtils.createRealMatrix(testData4x3)); 076 077 checkDimension(MatrixUtils.createRealMatrix(testData3x4)); 078 079 Random r = new Random(643895747384642l); 080 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4; 081 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4; 082 checkDimension(createTestMatrix(r, p, q)); 083 checkDimension(createTestMatrix(r, q, p)); 084 085 } 086 087 private void checkDimension(RealMatrix m) { 088 int rows = m.getRowDimension(); 089 int columns = m.getColumnDimension(); 090 QRDecomposition qr = new QRDecompositionImpl(m); 091 assertEquals(rows, qr.getQ().getRowDimension()); 092 assertEquals(rows, qr.getQ().getColumnDimension()); 093 assertEquals(rows, qr.getR().getRowDimension()); 094 assertEquals(columns, qr.getR().getColumnDimension()); 095 } 096 097 /** test A = QR */ 098 public void testAEqualQR() { 099 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular)); 100 101 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular)); 102 103 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4)); 104 105 checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3)); 106 107 Random r = new Random(643895747384642l); 108 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4; 109 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4; 110 checkAEqualQR(createTestMatrix(r, p, q)); 111 112 checkAEqualQR(createTestMatrix(r, q, p)); 113 114 } 115 116 private void checkAEqualQR(RealMatrix m) { 117 QRDecomposition qr = new QRDecompositionImpl(m); 118 double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm(); 119 assertEquals(0, norm, normTolerance); 120 } 121 122 /** test the orthogonality of Q */ 123 public void testQOrthogonal() { 124 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular)); 125 126 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular)); 127 128 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4)); 129 130 checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3)); 131 132 Random r = new Random(643895747384642l); 133 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4; 134 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4; 135 checkQOrthogonal(createTestMatrix(r, p, q)); 136 137 checkQOrthogonal(createTestMatrix(r, q, p)); 138 139 } 140 141 private void checkQOrthogonal(RealMatrix m) { 142 QRDecomposition qr = new QRDecompositionImpl(m); 143 RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension()); 144 double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm(); 145 assertEquals(0, norm, normTolerance); 146 } 147 148 /** test that R is upper triangular */ 149 public void testRUpperTriangular() { 150 RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular); 151 checkUpperTriangular(new QRDecompositionImpl(matrix).getR()); 152 153 matrix = MatrixUtils.createRealMatrix(testData3x3Singular); 154 checkUpperTriangular(new QRDecompositionImpl(matrix).getR()); 155 156 matrix = MatrixUtils.createRealMatrix(testData3x4); 157 checkUpperTriangular(new QRDecompositionImpl(matrix).getR()); 158 159 matrix = MatrixUtils.createRealMatrix(testData4x3); 160 checkUpperTriangular(new QRDecompositionImpl(matrix).getR()); 161 162 Random r = new Random(643895747384642l); 163 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4; 164 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4; 165 matrix = createTestMatrix(r, p, q); 166 checkUpperTriangular(new QRDecompositionImpl(matrix).getR()); 167 168 matrix = createTestMatrix(r, p, q); 169 checkUpperTriangular(new QRDecompositionImpl(matrix).getR()); 170 171 } 172 173 private void checkUpperTriangular(RealMatrix m) { 174 m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() { 175 @Override 176 public void visit(int row, int column, double value) { 177 if (column < row) { 178 assertEquals(0.0, value, entryTolerance); 179 } 180 } 181 }); 182 } 183 184 /** test that H is trapezoidal */ 185 public void testHTrapezoidal() { 186 RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular); 187 checkTrapezoidal(new QRDecompositionImpl(matrix).getH()); 188 189 matrix = MatrixUtils.createRealMatrix(testData3x3Singular); 190 checkTrapezoidal(new QRDecompositionImpl(matrix).getH()); 191 192 matrix = MatrixUtils.createRealMatrix(testData3x4); 193 checkTrapezoidal(new QRDecompositionImpl(matrix).getH()); 194 195 matrix = MatrixUtils.createRealMatrix(testData4x3); 196 checkTrapezoidal(new QRDecompositionImpl(matrix).getH()); 197 198 Random r = new Random(643895747384642l); 199 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4; 200 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4; 201 matrix = createTestMatrix(r, p, q); 202 checkTrapezoidal(new QRDecompositionImpl(matrix).getH()); 203 204 matrix = createTestMatrix(r, p, q); 205 checkTrapezoidal(new QRDecompositionImpl(matrix).getH()); 206 207 } 208 209 private void checkTrapezoidal(RealMatrix m) { 210 m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() { 211 @Override 212 public void visit(int row, int column, double value) { 213 if (column > row) { 214 assertEquals(0.0, value, entryTolerance); 215 } 216 } 217 }); 218 } 219 /** test matrices values */ 220 public void testMatricesValues() { 221 QRDecomposition qr = 222 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)); 223 RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] { 224 { -12.0 / 14.0, 69.0 / 175.0, -58.0 / 175.0 }, 225 { -6.0 / 14.0, -158.0 / 175.0, 6.0 / 175.0 }, 226 { 4.0 / 14.0, -30.0 / 175.0, -165.0 / 175.0 } 227 }); 228 RealMatrix rRef = MatrixUtils.createRealMatrix(new double[][] { 229 { -14.0, -21.0, 14.0 }, 230 { 0.0, -175.0, 70.0 }, 231 { 0.0, 0.0, 35.0 } 232 }); 233 RealMatrix hRef = MatrixUtils.createRealMatrix(new double[][] { 234 { 26.0 / 14.0, 0.0, 0.0 }, 235 { 6.0 / 14.0, 648.0 / 325.0, 0.0 }, 236 { -4.0 / 14.0, 36.0 / 325.0, 2.0 } 237 }); 238 239 // check values against known references 240 RealMatrix q = qr.getQ(); 241 assertEquals(0, q.subtract(qRef).getNorm(), 1.0e-13); 242 RealMatrix qT = qr.getQT(); 243 assertEquals(0, qT.subtract(qRef.transpose()).getNorm(), 1.0e-13); 244 RealMatrix r = qr.getR(); 245 assertEquals(0, r.subtract(rRef).getNorm(), 1.0e-13); 246 RealMatrix h = qr.getH(); 247 assertEquals(0, h.subtract(hRef).getNorm(), 1.0e-13); 248 249 // check the same cached instance is returned the second time 250 assertTrue(q == qr.getQ()); 251 assertTrue(r == qr.getR()); 252 assertTrue(h == qr.getH()); 253 254 } 255 256 private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) { 257 RealMatrix m = MatrixUtils.createRealMatrix(rows, columns); 258 m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){ 259 @Override 260 public double visit(int row, int column, double value) 261 throws MatrixVisitorException { 262 return 2.0 * r.nextDouble() - 1.0; 263 } 264 }); 265 return m; 266 } 267 268 }