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