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 junit.framework.Test; 023 import junit.framework.TestCase; 024 import junit.framework.TestSuite; 025 026 import org.apache.commons.math.linear.DecompositionSolver; 027 import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor; 028 import org.apache.commons.math.linear.BlockRealMatrix; 029 import org.apache.commons.math.linear.InvalidMatrixException; 030 import org.apache.commons.math.linear.MatrixUtils; 031 import org.apache.commons.math.linear.MatrixVisitorException; 032 import org.apache.commons.math.linear.QRDecomposition; 033 import org.apache.commons.math.linear.QRDecompositionImpl; 034 import org.apache.commons.math.linear.RealMatrix; 035 import org.apache.commons.math.linear.RealVector; 036 import org.apache.commons.math.linear.ArrayRealVector; 037 038 public class QRSolverTest extends TestCase { 039 double[][] testData3x3NonSingular = { 040 { 12, -51, 4 }, 041 { 6, 167, -68 }, 042 { -4, 24, -41 } 043 }; 044 045 double[][] testData3x3Singular = { 046 { 1, 2, 2 }, 047 { 2, 4, 6 }, 048 { 4, 8, 12 } 049 }; 050 051 double[][] testData3x4 = { 052 { 12, -51, 4, 1 }, 053 { 6, 167, -68, 2 }, 054 { -4, 24, -41, 3 } 055 }; 056 057 double[][] testData4x3 = { 058 { 12, -51, 4 }, 059 { 6, 167, -68 }, 060 { -4, 24, -41 }, 061 { -5, 34, 7 } 062 }; 063 064 public QRSolverTest(String name) { 065 super(name); 066 } 067 068 public static Test suite() { 069 TestSuite suite = new TestSuite(QRSolverTest.class); 070 suite.setName("QRSolver Tests"); 071 return suite; 072 } 073 074 /** test rank */ 075 public void testRank() { 076 DecompositionSolver solver = 077 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver(); 078 assertTrue(solver.isNonSingular()); 079 080 solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver(); 081 assertFalse(solver.isNonSingular()); 082 083 solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x4)).getSolver(); 084 assertTrue(solver.isNonSingular()); 085 086 solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData4x3)).getSolver(); 087 assertTrue(solver.isNonSingular()); 088 089 } 090 091 /** test solve dimension errors */ 092 public void testSolveDimensionErrors() { 093 DecompositionSolver solver = 094 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver(); 095 RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]); 096 try { 097 solver.solve(b); 098 fail("an exception should have been thrown"); 099 } catch (IllegalArgumentException iae) { 100 // expected behavior 101 } catch (Exception e) { 102 fail("wrong exception caught"); 103 } 104 try { 105 solver.solve(b.getColumn(0)); 106 fail("an exception should have been thrown"); 107 } catch (IllegalArgumentException iae) { 108 // expected behavior 109 } catch (Exception e) { 110 fail("wrong exception caught"); 111 } 112 try { 113 solver.solve(b.getColumnVector(0)); 114 fail("an exception should have been thrown"); 115 } catch (IllegalArgumentException iae) { 116 // expected behavior 117 } catch (Exception e) { 118 fail("wrong exception caught"); 119 } 120 } 121 122 /** test solve rank errors */ 123 public void testSolveRankErrors() { 124 DecompositionSolver solver = 125 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver(); 126 RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]); 127 try { 128 solver.solve(b); 129 fail("an exception should have been thrown"); 130 } catch (InvalidMatrixException iae) { 131 // expected behavior 132 } catch (Exception e) { 133 fail("wrong exception caught"); 134 } 135 try { 136 solver.solve(b.getColumn(0)); 137 fail("an exception should have been thrown"); 138 } catch (InvalidMatrixException iae) { 139 // expected behavior 140 } catch (Exception e) { 141 fail("wrong exception caught"); 142 } 143 try { 144 solver.solve(b.getColumnVector(0)); 145 fail("an exception should have been thrown"); 146 } catch (InvalidMatrixException iae) { 147 // expected behavior 148 } catch (Exception e) { 149 fail("wrong exception caught"); 150 } 151 } 152 153 /** test solve */ 154 public void testSolve() { 155 QRDecomposition decomposition = 156 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)); 157 DecompositionSolver solver = decomposition.getSolver(); 158 RealMatrix b = MatrixUtils.createRealMatrix(new double[][] { 159 { -102, 12250 }, { 544, 24500 }, { 167, -36750 } 160 }); 161 RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] { 162 { 1, 2515 }, { 2, 422 }, { -3, 898 } 163 }); 164 165 // using RealMatrix 166 assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm()); 167 168 // using double[] 169 for (int i = 0; i < b.getColumnDimension(); ++i) { 170 final double[] x = solver.solve(b.getColumn(i)); 171 final double error = new ArrayRealVector(x).subtract(xRef.getColumnVector(i)).getNorm(); 172 assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm()); 173 } 174 175 // using ArrayRealVector 176 for (int i = 0; i < b.getColumnDimension(); ++i) { 177 final RealVector x = solver.solve(b.getColumnVector(i)); 178 final double error = x.subtract(xRef.getColumnVector(i)).getNorm(); 179 assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm()); 180 } 181 182 // using RealVector with an alternate implementation 183 for (int i = 0; i < b.getColumnDimension(); ++i) { 184 ArrayRealVectorTest.RealVectorTestImpl v = 185 new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i)); 186 final RealVector x = solver.solve(v); 187 final double error = x.subtract(xRef.getColumnVector(i)).getNorm(); 188 assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm()); 189 } 190 191 } 192 193 public void testOverdetermined() { 194 final Random r = new Random(5559252868205245l); 195 int p = (7 * BlockRealMatrix.BLOCK_SIZE) / 4; 196 int q = (5 * BlockRealMatrix.BLOCK_SIZE) / 4; 197 RealMatrix a = createTestMatrix(r, p, q); 198 RealMatrix xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3); 199 200 // build a perturbed system: A.X + noise = B 201 RealMatrix b = a.multiply(xRef); 202 final double noise = 0.001; 203 b.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() { 204 @Override 205 public double visit(int row, int column, double value) { 206 return value * (1.0 + noise * (2 * r.nextDouble() - 1)); 207 } 208 }); 209 210 // despite perturbation, the least square solution should be pretty good 211 RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b); 212 assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p * q); 213 214 } 215 216 public void testUnderdetermined() { 217 final Random r = new Random(42185006424567123l); 218 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4; 219 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4; 220 RealMatrix a = createTestMatrix(r, p, q); 221 RealMatrix xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3); 222 RealMatrix b = a.multiply(xRef); 223 RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b); 224 225 // too many equations, the system cannot be solved at all 226 assertTrue(x.subtract(xRef).getNorm() / (p * q) > 0.01); 227 228 // the last unknown should have been set to 0 229 assertEquals(0.0, x.getSubMatrix(p, q - 1, 0, x.getColumnDimension() - 1).getNorm()); 230 231 } 232 233 private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) { 234 RealMatrix m = MatrixUtils.createRealMatrix(rows, columns); 235 m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){ 236 @Override 237 public double visit(int row, int column, double value) 238 throws MatrixVisitorException { 239 return 2.0 * r.nextDouble() - 1.0; 240 } 241 }); 242 return m; 243 } 244 }