View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.linear;
24  
25  import java.util.Random;
26  
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.junit.Assert;
29  import org.junit.Test;
30  
31  
32  public class QRDecompositionTest {
33      private double[][] testData3x3NonSingular = {
34              { 12, -51, 4 },
35              { 6, 167, -68 },
36              { -4, 24, -41 }, };
37  
38      private double[][] testData3x3Singular = {
39              { 1, 4, 7, },
40              { 2, 5, 8, },
41              { 3, 6, 9, }, };
42  
43      private double[][] testData3x4 = {
44              { 12, -51, 4, 1 },
45              { 6, 167, -68, 2 },
46              { -4, 24, -41, 3 }, };
47  
48      private double[][] testData4x3 = {
49              { 12, -51, 4, },
50              { 6, 167, -68, },
51              { -4, 24, -41, },
52              { -5, 34, 7, }, };
53  
54      private static final double entryTolerance = 10e-16;
55  
56      private static final double normTolerance = 10e-14;
57  
58      /** test dimensions */
59      @Test
60      public void testDimensions() {
61          checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
62  
63          checkDimension(MatrixUtils.createRealMatrix(testData4x3));
64  
65          checkDimension(MatrixUtils.createRealMatrix(testData3x4));
66  
67          Random r = new Random(643895747384642l);
68          int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
69          int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
70          checkDimension(createTestMatrix(r, p, q));
71          checkDimension(createTestMatrix(r, q, p));
72  
73      }
74  
75      private void checkDimension(RealMatrix m) {
76          int rows = m.getRowDimension();
77          int columns = m.getColumnDimension();
78          QRDecomposition qr = new QRDecomposition(m);
79          Assert.assertEquals(rows,    qr.getQ().getRowDimension());
80          Assert.assertEquals(rows,    qr.getQ().getColumnDimension());
81          Assert.assertEquals(rows,    qr.getR().getRowDimension());
82          Assert.assertEquals(columns, qr.getR().getColumnDimension());
83      }
84  
85      /** test A = QR */
86      @Test
87      public void testAEqualQR() {
88          checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
89  
90          checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
91  
92          checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4));
93  
94          checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3));
95  
96          Random r = new Random(643895747384642l);
97          int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
98          int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
99          checkAEqualQR(createTestMatrix(r, p, q));
100 
101         checkAEqualQR(createTestMatrix(r, q, p));
102 
103     }
104 
105     private void checkAEqualQR(RealMatrix m) {
106         QRDecomposition qr = new QRDecomposition(m);
107         double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm1();
108         Assert.assertEquals(0, norm, normTolerance);
109     }
110 
111     /** test the orthogonality of Q */
112     @Test
113     public void testQOrthogonal() {
114         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular));
115 
116         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular));
117 
118         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4));
119 
120         checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3));
121 
122         Random r = new Random(643895747384642l);
123         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
124         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
125         checkQOrthogonal(createTestMatrix(r, p, q));
126 
127         checkQOrthogonal(createTestMatrix(r, q, p));
128 
129     }
130 
131     private void checkQOrthogonal(RealMatrix m) {
132         QRDecomposition qr = new QRDecomposition(m);
133         RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
134         double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm1();
135         Assert.assertEquals(0, norm, normTolerance);
136     }
137 
138     /** test that R is upper triangular */
139     @Test
140     public void testRUpperTriangular() {
141         RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
142         checkUpperTriangular(new QRDecomposition(matrix).getR());
143 
144         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
145         checkUpperTriangular(new QRDecomposition(matrix).getR());
146 
147         matrix = MatrixUtils.createRealMatrix(testData3x4);
148         checkUpperTriangular(new QRDecomposition(matrix).getR());
149 
150         matrix = MatrixUtils.createRealMatrix(testData4x3);
151         checkUpperTriangular(new QRDecomposition(matrix).getR());
152 
153         Random r = new Random(643895747384642l);
154         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
155         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
156         matrix = createTestMatrix(r, p, q);
157         checkUpperTriangular(new QRDecomposition(matrix).getR());
158 
159         matrix = createTestMatrix(r, p, q);
160         checkUpperTriangular(new QRDecomposition(matrix).getR());
161 
162     }
163 
164     private void checkUpperTriangular(RealMatrix m) {
165         m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
166             @Override
167             public void visit(int row, int column, double value) {
168                 if (column < row) {
169                     Assert.assertEquals(0.0, value, entryTolerance);
170                 }
171             }
172         });
173     }
174 
175     /** test that H is trapezoidal */
176     @Test
177     public void testHTrapezoidal() {
178         RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
179         checkTrapezoidal(new QRDecomposition(matrix).getH());
180 
181         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
182         checkTrapezoidal(new QRDecomposition(matrix).getH());
183 
184         matrix = MatrixUtils.createRealMatrix(testData3x4);
185         checkTrapezoidal(new QRDecomposition(matrix).getH());
186 
187         matrix = MatrixUtils.createRealMatrix(testData4x3);
188         checkTrapezoidal(new QRDecomposition(matrix).getH());
189 
190         Random r = new Random(643895747384642l);
191         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
192         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
193         matrix = createTestMatrix(r, p, q);
194         checkTrapezoidal(new QRDecomposition(matrix).getH());
195 
196         matrix = createTestMatrix(r, p, q);
197         checkTrapezoidal(new QRDecomposition(matrix).getH());
198 
199     }
200 
201     private void checkTrapezoidal(RealMatrix m) {
202         m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
203             @Override
204             public void visit(int row, int column, double value) {
205                 if (column > row) {
206                     Assert.assertEquals(0.0, value, entryTolerance);
207                 }
208             }
209         });
210     }
211     /** test matrices values */
212     @Test
213     public void testMatricesValues() {
214         QRDecomposition qr =
215             new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
216         RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] {
217                 { -12.0 / 14.0,   69.0 / 175.0,  -58.0 / 175.0 },
218                 {  -6.0 / 14.0, -158.0 / 175.0,    6.0 / 175.0 },
219                 {   4.0 / 14.0,  -30.0 / 175.0, -165.0 / 175.0 }
220         });
221         RealMatrix rRef = MatrixUtils.createRealMatrix(new double[][] {
222                 { -14.0,  -21.0, 14.0 },
223                 {   0.0, -175.0, 70.0 },
224                 {   0.0,    0.0, 35.0 }
225         });
226         RealMatrix hRef = MatrixUtils.createRealMatrix(new double[][] {
227                 { 26.0 / 14.0, 0.0, 0.0 },
228                 {  6.0 / 14.0, 648.0 / 325.0, 0.0 },
229                 { -4.0 / 14.0,  36.0 / 325.0, 2.0 }
230         });
231 
232         // check values against known references
233         RealMatrix q = qr.getQ();
234         Assert.assertEquals(0, q.subtract(qRef).getNorm1(), 1.0e-13);
235         RealMatrix qT = qr.getQT();
236         Assert.assertEquals(0, qT.subtract(qRef.transpose()).getNorm1(), 1.0e-13);
237         RealMatrix r = qr.getR();
238         Assert.assertEquals(0, r.subtract(rRef).getNorm1(), 1.0e-13);
239         RealMatrix h = qr.getH();
240         Assert.assertEquals(0, h.subtract(hRef).getNorm1(), 1.0e-13);
241 
242         // check the same cached instance is returned the second time
243         Assert.assertTrue(q == qr.getQ());
244         Assert.assertTrue(r == qr.getR());
245         Assert.assertTrue(h == qr.getH());
246 
247     }
248 
249     @Test(expected=MathIllegalArgumentException.class)
250     public void testNonInvertible() {
251         QRDecomposition qr =
252             new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular));
253         qr.getSolver().getInverse();
254     }
255 
256     @Test
257     public void testInvertTallSkinny() {
258         RealMatrix a     = MatrixUtils.createRealMatrix(testData4x3);
259         DecompositionSolver solver = new QRDecomposition(a).getSolver();
260         RealMatrix pinv  = solver.getInverse();
261         Assert.assertEquals(0, pinv.multiply(a).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm1(), 1.0e-6);
262         Assert.assertEquals(testData4x3.length,    solver.getRowDimension());
263         Assert.assertEquals(testData4x3[0].length, solver.getColumnDimension());
264     }
265 
266     @Test
267     public void testInvertShortWide() {
268         RealMatrix a = MatrixUtils.createRealMatrix(testData3x4);
269         DecompositionSolver solver = new QRDecomposition(a).getSolver();
270         RealMatrix pinv  = solver.getInverse();
271         Assert.assertEquals(0, a.multiply(pinv).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm1(), 1.0e-6);
272         Assert.assertEquals(0, pinv.multiply(a).getSubMatrix(0, 2, 0, 2).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm1(), 1.0e-6);
273         Assert.assertEquals(testData3x4.length,    solver.getRowDimension());
274         Assert.assertEquals(testData3x4[0].length, solver.getColumnDimension());
275     }
276 
277     private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
278         RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
279         m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
280             @Override
281             public double visit(int row, int column, double value) {
282                 return 2.0 * r.nextDouble() - 1.0;
283             }
284         });
285         return m;
286     }
287 
288     @Test(expected=MathIllegalArgumentException.class)
289     public void testQRSingular() {
290         final RealMatrix a = MatrixUtils.createRealMatrix(new double[][] {
291             { 1, 6, 4 }, { 2, 4, -1 }, { -1, 2, 5 }
292         });
293         final RealVector b = new ArrayRealVector(new double[]{ 5, 6, 1 });
294         new QRDecomposer(1.0e-15).decompose(a).solve(b);
295     }
296 
297 }