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 org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
23 import org.apache.commons.math.linear.DefaultRealMatrixPreservingVisitor;
24 import org.apache.commons.math.linear.BlockRealMatrix;
25 import org.apache.commons.math.linear.MatrixUtils;
26 import org.apache.commons.math.linear.MatrixVisitorException;
27 import org.apache.commons.math.linear.QRDecomposition;
28 import org.apache.commons.math.linear.QRDecompositionImpl;
29 import org.apache.commons.math.linear.RealMatrix;
30
31 import junit.framework.Test;
32 import junit.framework.TestCase;
33 import junit.framework.TestSuite;
34
35 public class QRDecompositionImplTest extends TestCase {
36 double[][] testData3x3NonSingular = {
37 { 12, -51, 4 },
38 { 6, 167, -68 },
39 { -4, 24, -41 }, };
40
41 double[][] testData3x3Singular = {
42 { 1, 4, 7, },
43 { 2, 5, 8, },
44 { 3, 6, 9, }, };
45
46 double[][] testData3x4 = {
47 { 12, -51, 4, 1 },
48 { 6, 167, -68, 2 },
49 { -4, 24, -41, 3 }, };
50
51 double[][] testData4x3 = {
52 { 12, -51, 4, },
53 { 6, 167, -68, },
54 { -4, 24, -41, },
55 { -5, 34, 7, }, };
56
57 private static final double entryTolerance = 10e-16;
58
59 private static final double normTolerance = 10e-14;
60
61 public QRDecompositionImplTest(String name) {
62 super(name);
63 }
64
65 public static Test suite() {
66 TestSuite suite = new TestSuite(QRDecompositionImplTest.class);
67 suite.setName("QRDecompositionImpl Tests");
68 return suite;
69 }
70
71
72 public void testDimensions() {
73 checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
74
75 checkDimension(MatrixUtils.createRealMatrix(testData4x3));
76
77 checkDimension(MatrixUtils.createRealMatrix(testData3x4));
78
79 Random r = new Random(643895747384642l);
80 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
81 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
82 checkDimension(createTestMatrix(r, p, q));
83 checkDimension(createTestMatrix(r, q, p));
84
85 }
86
87 private void checkDimension(RealMatrix m) {
88 int rows = m.getRowDimension();
89 int columns = m.getColumnDimension();
90 QRDecomposition qr = new QRDecompositionImpl(m);
91 assertEquals(rows, qr.getQ().getRowDimension());
92 assertEquals(rows, qr.getQ().getColumnDimension());
93 assertEquals(rows, qr.getR().getRowDimension());
94 assertEquals(columns, qr.getR().getColumnDimension());
95 }
96
97
98 public void testAEqualQR() {
99 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
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
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
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
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
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
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 }