1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.linear;
19
20 import java.util.Random;
21
22 import junit.framework.Test;
23 import junit.framework.TestCase;
24 import junit.framework.TestSuite;
25
26 import org.apache.commons.math.linear.DecompositionSolver;
27 import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
28 import org.apache.commons.math.linear.BlockRealMatrix;
29 import org.apache.commons.math.linear.InvalidMatrixException;
30 import org.apache.commons.math.linear.MatrixUtils;
31 import org.apache.commons.math.linear.MatrixVisitorException;
32 import org.apache.commons.math.linear.QRDecomposition;
33 import org.apache.commons.math.linear.QRDecompositionImpl;
34 import org.apache.commons.math.linear.RealMatrix;
35 import org.apache.commons.math.linear.RealVector;
36 import org.apache.commons.math.linear.ArrayRealVector;
37
38 public class QRSolverTest extends TestCase {
39 double[][] testData3x3NonSingular = {
40 { 12, -51, 4 },
41 { 6, 167, -68 },
42 { -4, 24, -41 }
43 };
44
45 double[][] testData3x3Singular = {
46 { 1, 2, 2 },
47 { 2, 4, 6 },
48 { 4, 8, 12 }
49 };
50
51 double[][] testData3x4 = {
52 { 12, -51, 4, 1 },
53 { 6, 167, -68, 2 },
54 { -4, 24, -41, 3 }
55 };
56
57 double[][] testData4x3 = {
58 { 12, -51, 4 },
59 { 6, 167, -68 },
60 { -4, 24, -41 },
61 { -5, 34, 7 }
62 };
63
64 public QRSolverTest(String name) {
65 super(name);
66 }
67
68 public static Test suite() {
69 TestSuite suite = new TestSuite(QRSolverTest.class);
70 suite.setName("QRSolver Tests");
71 return suite;
72 }
73
74
75 public void testRank() {
76 DecompositionSolver solver =
77 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
78 assertTrue(solver.isNonSingular());
79
80 solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
81 assertFalse(solver.isNonSingular());
82
83 solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
84 assertTrue(solver.isNonSingular());
85
86 solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
87 assertTrue(solver.isNonSingular());
88
89 }
90
91
92 public void testSolveDimensionErrors() {
93 DecompositionSolver solver =
94 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
95 RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
96 try {
97 solver.solve(b);
98 fail("an exception should have been thrown");
99 } catch (IllegalArgumentException iae) {
100
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
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
117 } catch (Exception e) {
118 fail("wrong exception caught");
119 }
120 }
121
122
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
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
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
148 } catch (Exception e) {
149 fail("wrong exception caught");
150 }
151 }
152
153
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
166 assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm());
167
168
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
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
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
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
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
226 assertTrue(x.subtract(xRef).getNorm() / (p * q) > 0.01);
227
228
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 }