1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.linear;
24
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.junit.jupiter.api.Test;
27
28 import java.util.Random;
29
30 import static org.junit.jupiter.api.Assertions.assertEquals;
31 import static org.junit.jupiter.api.Assertions.assertThrows;
32 import static org.junit.jupiter.api.Assertions.assertTrue;
33
34
35 class QRDecompositionTest {
36 private double[][] testData3x3NonSingular = {
37 { 12, -51, 4 },
38 { 6, 167, -68 },
39 { -4, 24, -41 }, };
40
41 private double[][] testData3x3Singular = {
42 { 1, 4, 7, },
43 { 2, 5, 8, },
44 { 3, 6, 9, }, };
45
46 private double[][] testData3x4 = {
47 { 12, -51, 4, 1 },
48 { 6, 167, -68, 2 },
49 { -4, 24, -41, 3 }, };
50
51 private 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
62 @Test
63 void testDimensions() {
64 checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
65
66 checkDimension(MatrixUtils.createRealMatrix(testData4x3));
67
68 checkDimension(MatrixUtils.createRealMatrix(testData3x4));
69
70 Random r = new Random(643895747384642l);
71 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
72 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
73 checkDimension(createTestMatrix(r, p, q));
74 checkDimension(createTestMatrix(r, q, p));
75
76 }
77
78 private void checkDimension(RealMatrix m) {
79 int rows = m.getRowDimension();
80 int columns = m.getColumnDimension();
81 QRDecomposition qr = new QRDecomposition(m);
82 assertEquals(rows, qr.getQ().getRowDimension());
83 assertEquals(rows, qr.getQ().getColumnDimension());
84 assertEquals(rows, qr.getR().getRowDimension());
85 assertEquals(columns, qr.getR().getColumnDimension());
86 }
87
88
89 @Test
90 void testAEqualQR() {
91 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
92
93 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
94
95 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4));
96
97 checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3));
98
99 Random r = new Random(643895747384642l);
100 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
101 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
102 checkAEqualQR(createTestMatrix(r, p, q));
103
104 checkAEqualQR(createTestMatrix(r, q, p));
105
106 }
107
108 private void checkAEqualQR(RealMatrix m) {
109 QRDecomposition qr = new QRDecomposition(m);
110 double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm1();
111 assertEquals(0, norm, normTolerance);
112 }
113
114
115 @Test
116 void testQOrthogonal() {
117 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular));
118
119 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular));
120
121 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4));
122
123 checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3));
124
125 Random r = new Random(643895747384642l);
126 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
127 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
128 checkQOrthogonal(createTestMatrix(r, p, q));
129
130 checkQOrthogonal(createTestMatrix(r, q, p));
131
132 }
133
134 private void checkQOrthogonal(RealMatrix m) {
135 QRDecomposition qr = new QRDecomposition(m);
136 RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
137 double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm1();
138 assertEquals(0, norm, normTolerance);
139 }
140
141
142 @Test
143 void testRUpperTriangular() {
144 RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
145 checkUpperTriangular(new QRDecomposition(matrix).getR());
146
147 matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
148 checkUpperTriangular(new QRDecomposition(matrix).getR());
149
150 matrix = MatrixUtils.createRealMatrix(testData3x4);
151 checkUpperTriangular(new QRDecomposition(matrix).getR());
152
153 matrix = MatrixUtils.createRealMatrix(testData4x3);
154 checkUpperTriangular(new QRDecomposition(matrix).getR());
155
156 Random r = new Random(643895747384642l);
157 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
158 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
159 matrix = createTestMatrix(r, p, q);
160 checkUpperTriangular(new QRDecomposition(matrix).getR());
161
162 matrix = createTestMatrix(r, p, q);
163 checkUpperTriangular(new QRDecomposition(matrix).getR());
164
165 }
166
167 private void checkUpperTriangular(RealMatrix m) {
168 m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
169 @Override
170 public void visit(int row, int column, double value) {
171 if (column < row) {
172 assertEquals(0.0, value, entryTolerance);
173 }
174 }
175 });
176 }
177
178
179 @Test
180 void testHTrapezoidal() {
181 RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
182 checkTrapezoidal(new QRDecomposition(matrix).getH());
183
184 matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
185 checkTrapezoidal(new QRDecomposition(matrix).getH());
186
187 matrix = MatrixUtils.createRealMatrix(testData3x4);
188 checkTrapezoidal(new QRDecomposition(matrix).getH());
189
190 matrix = MatrixUtils.createRealMatrix(testData4x3);
191 checkTrapezoidal(new QRDecomposition(matrix).getH());
192
193 Random r = new Random(643895747384642l);
194 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
195 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
196 matrix = createTestMatrix(r, p, q);
197 checkTrapezoidal(new QRDecomposition(matrix).getH());
198
199 matrix = createTestMatrix(r, p, q);
200 checkTrapezoidal(new QRDecomposition(matrix).getH());
201
202 }
203
204 private void checkTrapezoidal(RealMatrix m) {
205 m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
206 @Override
207 public void visit(int row, int column, double value) {
208 if (column > row) {
209 assertEquals(0.0, value, entryTolerance);
210 }
211 }
212 });
213 }
214
215
216 @Test
217 void testMatricesValues() {
218 QRDecomposition qr =
219 new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
220 RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] {
221 { -12.0 / 14.0, 69.0 / 175.0, -58.0 / 175.0 },
222 { -6.0 / 14.0, -158.0 / 175.0, 6.0 / 175.0 },
223 { 4.0 / 14.0, -30.0 / 175.0, -165.0 / 175.0 }
224 });
225 RealMatrix rRef = MatrixUtils.createRealMatrix(new double[][] {
226 { -14.0, -21.0, 14.0 },
227 { 0.0, -175.0, 70.0 },
228 { 0.0, 0.0, 35.0 }
229 });
230 RealMatrix hRef = MatrixUtils.createRealMatrix(new double[][] {
231 { 26.0 / 14.0, 0.0, 0.0 },
232 { 6.0 / 14.0, 648.0 / 325.0, 0.0 },
233 { -4.0 / 14.0, 36.0 / 325.0, 2.0 }
234 });
235
236
237 RealMatrix q = qr.getQ();
238 assertEquals(0, q.subtract(qRef).getNorm1(), 1.0e-13);
239 RealMatrix qT = qr.getQT();
240 assertEquals(0, qT.subtract(qRef.transpose()).getNorm1(), 1.0e-13);
241 RealMatrix r = qr.getR();
242 assertEquals(0, r.subtract(rRef).getNorm1(), 1.0e-13);
243 RealMatrix h = qr.getH();
244 assertEquals(0, h.subtract(hRef).getNorm1(), 1.0e-13);
245
246
247 assertTrue(q == qr.getQ());
248 assertTrue(r == qr.getR());
249 assertTrue(h == qr.getH());
250
251 }
252
253 @Test
254 void testNonInvertible() {
255 assertThrows(MathIllegalArgumentException.class, () -> {
256 QRDecomposition qr =
257 new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular));
258 qr.getSolver().getInverse();
259 });
260 }
261
262 @Test
263 void testInvertTallSkinny() {
264 RealMatrix a = MatrixUtils.createRealMatrix(testData4x3);
265 DecompositionSolver solver = new QRDecomposition(a).getSolver();
266 RealMatrix pinv = solver.getInverse();
267 assertEquals(0, pinv.multiply(a).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm1(), 1.0e-6);
268 assertEquals(testData4x3.length, solver.getRowDimension());
269 assertEquals(testData4x3[0].length, solver.getColumnDimension());
270 }
271
272 @Test
273 void testInvertShortWide() {
274 RealMatrix a = MatrixUtils.createRealMatrix(testData3x4);
275 DecompositionSolver solver = new QRDecomposition(a).getSolver();
276 RealMatrix pinv = solver.getInverse();
277 assertEquals(0, a.multiply(pinv).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm1(), 1.0e-6);
278 assertEquals(0, pinv.multiply(a).getSubMatrix(0, 2, 0, 2).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm1(), 1.0e-6);
279 assertEquals(testData3x4.length, solver.getRowDimension());
280 assertEquals(testData3x4[0].length, solver.getColumnDimension());
281 }
282
283 private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
284 RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
285 m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
286 @Override
287 public double visit(int row, int column, double value) {
288 return 2.0 * r.nextDouble() - 1.0;
289 }
290 });
291 return m;
292 }
293
294 @Test
295 void testQRSingular() {
296 assertThrows(MathIllegalArgumentException.class, () -> {
297 final RealMatrix a = MatrixUtils.createRealMatrix(new double[][]{
298 {1, 6, 4}, {2, 4, -1}, {-1, 2, 5}
299 });
300 final RealVector b = new ArrayRealVector(new double[]{5, 6, 1});
301 new QRDecomposer(1.0e-15).decompose(a).solve(b);
302 });
303 }
304
305 }