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 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  
33  
34  class RRQRDecompositionTest {
35      private double[][] testData3x3NonSingular = {
36              { 12, -51, 4 },
37              { 6, 167, -68 },
38              { -4, 24, -41 }, };
39  
40      private double[][] testData3x3Singular = {
41              { 1, 4, 7, },
42              { 2, 5, 8, },
43              { 3, 6, 9, }, };
44  
45      private double[][] testData3x4 = {
46              { 12, -51, 4, 1 },
47              { 6, 167, -68, 2 },
48              { -4, 24, -41, 3 }, };
49  
50      private double[][] testData4x3 = {
51              { 12, -51, 4, },
52              { 6, 167, -68, },
53              { -4, 24, -41, },
54              { -5, 34, 7, }, };
55  
56      private static final double entryTolerance = 10e-16;
57  
58      private static final double normTolerance = 10e-14;
59  
60      /** test dimensions */
61      @Test
62      void testDimensions() {
63          checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
64  
65          checkDimension(MatrixUtils.createRealMatrix(testData4x3));
66  
67          checkDimension(MatrixUtils.createRealMatrix(testData3x4));
68  
69          Random r = new Random(643895747384642l);
70          int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
71          int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
72          checkDimension(createTestMatrix(r, p, q));
73          checkDimension(createTestMatrix(r, q, p));
74  
75      }
76  
77      private void checkDimension(RealMatrix m) {
78          int rows = m.getRowDimension();
79          int columns = m.getColumnDimension();
80          RRQRDecomposition qr = new RRQRDecomposition(m);
81          assertEquals(rows,    qr.getQ().getRowDimension());
82          assertEquals(rows,    qr.getQ().getColumnDimension());
83          assertEquals(rows,    qr.getR().getRowDimension());
84          assertEquals(columns, qr.getR().getColumnDimension());
85      }
86  
87      /** test AP = QR */
88      @Test
89      void testAPEqualQR() {
90          checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
91  
92          checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
93  
94          checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x4));
95  
96          checkAPEqualQR(MatrixUtils.createRealMatrix(testData4x3));
97  
98          Random r = new Random(643895747384642l);
99          int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
100         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
101         checkAPEqualQR(createTestMatrix(r, p, q));
102 
103         checkAPEqualQR(createTestMatrix(r, q, p));
104 
105     }
106 
107     private void checkAPEqualQR(RealMatrix m) {
108         RRQRDecomposition rrqr = new RRQRDecomposition(m);
109         double norm = rrqr.getQ().multiply(rrqr.getR()).subtract(m.multiply(rrqr.getP())).getNorm1();
110         assertEquals(0, norm, normTolerance);
111     }
112 
113     /** test the orthogonality of Q */
114     @Test
115     void testQOrthogonal() {
116         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular));
117 
118         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular));
119 
120         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4));
121 
122         checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3));
123 
124         Random r = new Random(643895747384642l);
125         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
126         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
127         checkQOrthogonal(createTestMatrix(r, p, q));
128 
129         checkQOrthogonal(createTestMatrix(r, q, p));
130 
131     }
132 
133     private void checkQOrthogonal(RealMatrix m) {
134         RRQRDecomposition qr = new RRQRDecomposition(m);
135         RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
136         double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm1();
137         assertEquals(0, norm, normTolerance);
138     }
139 
140     /** test that R is upper triangular */
141     @Test
142     void testRUpperTriangular() {
143         RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
144         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
145 
146         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
147         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
148 
149         matrix = MatrixUtils.createRealMatrix(testData3x4);
150         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
151 
152         matrix = MatrixUtils.createRealMatrix(testData4x3);
153         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
154 
155         Random r = new Random(643895747384642l);
156         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
157         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
158         matrix = createTestMatrix(r, p, q);
159         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
160 
161         matrix = createTestMatrix(r, p, q);
162         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
163 
164     }
165 
166     private void checkUpperTriangular(RealMatrix m) {
167         m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
168             @Override
169             public void visit(int row, int column, double value) {
170                 if (column < row) {
171                     assertEquals(0.0, value, entryTolerance);
172                 }
173             }
174         });
175     }
176 
177     /** test that H is trapezoidal */
178     @Test
179     void testHTrapezoidal() {
180         RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
181         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
182 
183         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
184         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
185 
186         matrix = MatrixUtils.createRealMatrix(testData3x4);
187         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
188 
189         matrix = MatrixUtils.createRealMatrix(testData4x3);
190         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
191 
192         Random r = new Random(643895747384642l);
193         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
194         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
195         matrix = createTestMatrix(r, p, q);
196         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
197 
198         matrix = createTestMatrix(r, p, q);
199         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
200 
201     }
202 
203     private void checkTrapezoidal(RealMatrix m) {
204         m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
205             @Override
206             public void visit(int row, int column, double value) {
207                 if (column > row) {
208                     assertEquals(0.0, value, entryTolerance);
209                 }
210             }
211         });
212     }
213 
214     @Test
215     void testNonInvertible() {
216         assertThrows(MathIllegalArgumentException.class, () -> {
217             RRQRDecomposition qr =
218                 new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular), 3.0e-16);
219             qr.getSolver().getInverse();
220         });
221     }
222 
223     private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
224         RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
225         m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
226             @Override
227             public double visit(int row, int column, double value) {
228                 return 2.0 * r.nextDouble() - 1.0;
229             }
230         });
231         return m;
232     }
233 
234     /** test the rank is returned correctly */
235     @Test
236     void testRank() {
237         double[][] d = { { 1, 1, 1 }, { 0, 0, 0 }, { 1, 2, 3 } };
238         RealMatrix m = new Array2DRowRealMatrix(d);
239         RRQRDecomposition qr = new RRQRDecomposition(m);
240         assertEquals(2, qr.getRank(1.0e-16));
241     }
242 
243 }