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.CalculusFieldElement;
26  import org.hipparchus.Field;
27  import org.hipparchus.analysis.differentiation.Gradient;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.util.Binary64Field;
30  import org.hipparchus.util.MathArrays;
31  import org.junit.jupiter.api.Test;
32  
33  import java.util.Random;
34  
35  import static org.junit.jupiter.api.Assertions.assertEquals;
36  import static org.junit.jupiter.api.Assertions.assertThrows;
37  import static org.junit.jupiter.api.Assertions.assertTrue;
38  
39  
40  class FieldQRDecompositionTest {
41      private double[][] testData3x3NonSingular = {
42              { 12, -51, 4 },
43              { 6, 167, -68 },
44              { -4, 24, -41 }, };
45  
46      private double[][] testData3x3Singular = {
47              { 1, 4, 7, },
48              { 2, 5, 8, },
49              { 3, 6, 9, }, };
50  
51      private double[][] testData3x4 = {
52              { 12, -51, 4, 1 },
53              { 6, 167, -68, 2 },
54              { -4, 24, -41, 3 }, };
55  
56      private double[][] testData4x3 = {
57              { 12, -51, 4, },
58              { 6, 167, -68, },
59              { -4, 24, -41, },
60              { -5, 34, 7, }, };
61  
62      Gradient zero = Gradient.variable(1, 0, 0.);
63      Field<Gradient> GradientField = zero.getField();
64      private static final double entryTolerance = 10e-16;
65  
66      private static final double normTolerance = 10e-14;
67  
68      /**Testing if the dimensions are correct.*/
69      @Test
70      void testDimensions() {
71          doTestDimensions(GradientField);
72          doTestDimensions(Binary64Field.getInstance());   
73      }
74  
75      /**Testing if is impossible to solve QR.*/
76      @Test
77      void testQRSingular(){
78          assertThrows(MathIllegalArgumentException.class, () -> {
79              QRSingular(GradientField);
80              QRSingular(Binary64Field.getInstance());
81          });
82      }
83  
84      /**Testing if Q is orthogonal*/
85      @Test
86      void testQOrthogonal(){
87          QOrthogonal(GradientField);
88          QOrthogonal(Binary64Field.getInstance());
89      }
90  
91      /**Testing if A = Q * R*/
92      @Test
93      void testAEqualQR(){
94          AEqualQR(GradientField);
95          AEqualQR(Binary64Field.getInstance());
96      }
97  
98      /**Testing if R is upper triangular.*/
99      @Test
100     void testRUpperTriangular(){
101         RUpperTriangular(GradientField);
102         RUpperTriangular(Binary64Field.getInstance());
103     }
104 
105     /**Testing if H is trapezoidal.*/
106     @Test
107     void testHTrapezoidal(){
108         HTrapezoidal(GradientField);
109         HTrapezoidal(Binary64Field.getInstance());
110     }
111 
112     /**Testing the values of the matrices.*/
113     @Test
114     void testMatricesValues(){
115         MatricesValues(GradientField);
116         MatricesValues(Binary64Field.getInstance());
117     }
118 
119     /**Testing if there is an error inverting a non invertible matrix.*/
120     @Test
121     void testNonInvertible(){
122         assertThrows(MathIllegalArgumentException.class, () -> {
123             NonInvertible(GradientField);
124             NonInvertible(Binary64Field.getInstance());
125         });
126     }
127 
128     /**Testing to invert a tall and skinny matrix.*/
129     @Test
130     void testInvertTallSkinny(){
131         InvertTallSkinny(GradientField);
132         InvertTallSkinny(Binary64Field.getInstance());
133     }
134 
135     /**Testing to invert a short and wide matrix.*/
136     @Test
137     void testInvertShortWide(){
138         InvertShortWide(GradientField);
139         InvertShortWide(Binary64Field.getInstance());
140     }
141 
142     private <T extends CalculusFieldElement<T>> void doTestDimensions(Field<T> field){
143         T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
144         T[][] data3x4= convert(field, testData3x4            );
145         T[][] data4x3= convert(field, testData4x3            );
146         checkDimension(MatrixUtils.createFieldMatrix( data3x3NS));
147 
148         checkDimension(MatrixUtils.createFieldMatrix( data4x3));
149 
150         checkDimension(MatrixUtils.createFieldMatrix( data3x4));
151 
152         Random r = new Random(643895747384642l);
153         int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
154         int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
155         checkDimension(createTestMatrix(field,r, p, q));
156         checkDimension(createTestMatrix(field, r, q, p));
157 
158     }
159 
160     private <T extends CalculusFieldElement<T>> void checkDimension(FieldMatrix<T> m) {
161         int rows = m.getRowDimension();
162         int columns = m.getColumnDimension();
163         FieldQRDecomposition<T> qr = new FieldQRDecomposition<T>(m);
164         assertEquals(rows,    qr.getQ().getRowDimension());
165         assertEquals(rows,    qr.getQ().getColumnDimension());
166         assertEquals(rows,    qr.getR().getRowDimension());
167         assertEquals(columns, qr.getR().getColumnDimension());
168     }
169 
170     private  <T extends CalculusFieldElement<T>> void AEqualQR(Field<T> field) {
171         T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
172         T[][] data3x3S= convert(field, testData3x3Singular    );
173         T[][] data3x4= convert(field, testData3x4            );
174         T[][] data4x3= convert(field, testData4x3            );
175         checkAEqualQR(MatrixUtils.createFieldMatrix( data3x3NS));
176 
177         checkAEqualQR(MatrixUtils.createFieldMatrix( data3x3S));
178 
179         checkAEqualQR(MatrixUtils.createFieldMatrix( data3x4));
180 
181         checkAEqualQR(MatrixUtils.createFieldMatrix( data4x3));
182 
183         Random r = new Random(643895747384642l);
184         int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
185         int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
186         checkAEqualQR(createTestMatrix(field, r, p, q));
187 
188         checkAEqualQR(createTestMatrix(field, r, q, p));
189 
190     }
191 
192     private  <T extends CalculusFieldElement<T>> void checkAEqualQR(FieldMatrix<T> m) {
193         FieldQRDecomposition<T> qr = new FieldQRDecomposition<>(m);
194         T norm = norm(qr.getQ().multiply(qr.getR()).subtract(m));
195         assertEquals(0, norm.getReal(), normTolerance);
196     }
197 
198     private  <T extends CalculusFieldElement<T>> void QOrthogonal(Field<T> field) {
199         T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
200         T[][] data3x3S= convert(field, testData3x3Singular    );
201         T[][] data3x4= convert(field, testData3x4            );
202         T[][] data4x3= convert(field, testData4x3            );
203         checkQOrthogonal(MatrixUtils.createFieldMatrix( data3x3NS));
204 
205         checkQOrthogonal(MatrixUtils.createFieldMatrix( data3x3S));
206 
207         checkQOrthogonal(MatrixUtils.createFieldMatrix( data3x4));
208 
209         checkQOrthogonal(MatrixUtils.createFieldMatrix( data4x3));
210 
211         Random r = new Random(643895747384642l);
212         int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
213         int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
214         checkQOrthogonal(createTestMatrix(field, r, p, q));
215 
216         checkQOrthogonal(createTestMatrix(field, r, q, p));
217 
218     }
219 
220     private  <T extends CalculusFieldElement<T>> void checkQOrthogonal(FieldMatrix<T> m) {
221         FieldQRDecomposition<T> qr = new FieldQRDecomposition<T>(m);
222         FieldMatrix<T> eye = MatrixUtils.createFieldIdentityMatrix(m.getField(),m.getRowDimension());
223         T norm = norm(qr.getQT().multiply(qr.getQ()).subtract(eye));
224         assertEquals(0, norm.getReal(), normTolerance);
225     }
226 
227     private  <T extends CalculusFieldElement<T>> void RUpperTriangular(Field<T> field) {
228         T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
229         T[][] data3x3S= convert(field, testData3x3Singular    );
230         T[][] data3x4= convert(field, testData3x4            );
231         T[][] data4x3= convert(field, testData4x3            );
232         FieldMatrix<T> matrix = MatrixUtils.createFieldMatrix( data3x3NS);
233         checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
234 
235         matrix = MatrixUtils.createFieldMatrix( data3x3S);
236         checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
237 
238         matrix = MatrixUtils.createFieldMatrix( data3x4);
239         checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
240 
241         matrix = MatrixUtils.createFieldMatrix( data4x3);
242         checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
243 
244         Random r = new Random(643895747384642l);
245         int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
246         int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
247         matrix = createTestMatrix(field, r, p, q);
248         checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
249 
250         matrix = createTestMatrix(field, r, p, q);
251         checkUpperTriangular(new FieldQRDecomposition<T>(matrix).getR());
252 
253     }
254 
255     private  <T extends CalculusFieldElement<T>> void checkUpperTriangular(FieldMatrix<T> m) {
256         m.walkInOptimizedOrder(new DefaultFieldMatrixPreservingVisitor<T>(m.getField().getZero()) {
257             @Override
258             public void visit(int row, int column, T value) {
259                 if (column < row) {
260                     assertEquals(0.0, value.getReal(), entryTolerance);
261                 }
262             }
263         });
264     }
265 
266     private <T extends CalculusFieldElement<T>> void HTrapezoidal(Field<T> field) {
267         T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
268         T[][] data3x3S= convert(field, testData3x3Singular    );
269         T[][] data3x4= convert(field, testData3x4            );
270         T[][] data4x3= convert(field, testData4x3            );
271         FieldMatrix<T> matrix = MatrixUtils.createFieldMatrix( data3x3NS);
272         checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
273 
274         matrix = MatrixUtils.createFieldMatrix( data3x3S);
275         checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
276 
277         matrix = MatrixUtils.createFieldMatrix( data3x4);
278         checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
279 
280         matrix = MatrixUtils.createFieldMatrix( data4x3);
281         checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
282 
283         Random r = new Random(643895747384642l);
284         int    p = (5 * BlockFieldMatrix.BLOCK_SIZE) / 4;
285         int    q = (7 * BlockFieldMatrix.BLOCK_SIZE) / 4;
286         matrix = createTestMatrix(field, r, p, q);
287         checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
288 
289         matrix = createTestMatrix(field, r, p, q);
290         checkTrapezoidal(new FieldQRDecomposition<T>(matrix).getH());
291 
292     }
293 
294     private  <T extends CalculusFieldElement<T>> void checkTrapezoidal(FieldMatrix<T> m) {
295         m.walkInOptimizedOrder(new DefaultFieldMatrixPreservingVisitor<T>(m.getField().getZero()) {
296             @Override
297             public void visit(int row, int column, T value) {
298                 if (column > row) {
299                     assertEquals(0.0, value.getReal(), entryTolerance);
300                 }
301             }
302         });
303     }
304     
305     private  <T extends CalculusFieldElement<T>> void MatricesValues(Field<T> field) {
306         T[][] data3x3NS= convert(field, testData3x3NonSingular ); 
307         FieldQRDecomposition<T> qr =
308             new FieldQRDecomposition<T>(MatrixUtils.createFieldMatrix( data3x3NS));
309         FieldMatrix<T> qRef = MatrixUtils.createFieldMatrix( convert(field,new double[][] {
310                 { -12.0 / 14.0,   69.0 / 175.0,  -58.0 / 175.0 },
311                 {  -6.0 / 14.0, -158.0 / 175.0,    6.0 / 175.0 },
312                 {   4.0 / 14.0,  -30.0 / 175.0, -165.0 / 175.0 }
313         }));
314         FieldMatrix<T> rRef = MatrixUtils.createFieldMatrix( convert(field, new double[][] {
315                 { -14.0,  -21.0, 14.0 },
316                 {   0.0, -175.0, 70.0 },
317                 {   0.0,    0.0, 35.0 }
318         }));
319         FieldMatrix<T> hRef = MatrixUtils.createFieldMatrix( convert(field,new double[][] {
320                 { 26.0 / 14.0, 0.0, 0.0 },
321                 {  6.0 / 14.0, 648.0 / 325.0, 0.0 },
322                 { -4.0 / 14.0,  36.0 / 325.0, 2.0 }
323         }));
324 
325         // check values against known references
326         FieldMatrix<T> q = qr.getQ();
327         assertEquals(0, norm(q.subtract(qRef)).getReal(), 1.0e-13);
328         FieldMatrix<T> qT = qr.getQT();
329         assertEquals(0, norm(qT.subtract(qRef.transpose())).getReal(), 1.0e-13);
330         FieldMatrix<T> r = qr.getR();
331         assertEquals(0, norm(r.subtract(rRef)).getReal(), 1.0e-13);
332         FieldMatrix<T> h = qr.getH();
333         assertEquals(0, norm(h.subtract(hRef)).getReal(), 1.0e-13);
334 
335         // check the same cached instance is returned the second time
336         assertTrue(q == qr.getQ());
337         assertTrue(r == qr.getR());
338         assertTrue(h == qr.getH());
339 
340     }
341 
342     private  <T extends CalculusFieldElement<T>> void NonInvertible(Field<T> field) {
343         T[][] data3x3S= convert(field, testData3x3Singular    );
344         FieldQRDecomposition<T> qr =
345             new FieldQRDecomposition<T>(MatrixUtils.createFieldMatrix( data3x3S));
346         qr.getSolver().getInverse();
347     }
348 
349     private  <T extends CalculusFieldElement<T>> void InvertTallSkinny(Field<T> field) {
350         T[][] data4x3= convert(field, testData4x3            );
351         FieldMatrix<T> a     = MatrixUtils.createFieldMatrix(data4x3);
352         FieldDecompositionSolver<T> solver = new FieldQRDecomposer<>(field.getZero()).decompose(a);
353         FieldMatrix<T> pinv  = solver.getInverse();
354         assertEquals(0, norm(pinv.multiply(a).subtract(MatrixUtils.createFieldIdentityMatrix(field, 3))).getReal(), 1.0e-6);
355         assertEquals(testData4x3.length,    solver.getRowDimension());
356         assertEquals(testData4x3[0].length, solver.getColumnDimension());
357     }
358 
359     private  <T extends CalculusFieldElement<T>> void InvertShortWide(Field<T> field) {
360         T[][] data3x4= convert(field, testData3x4            );
361         FieldMatrix<T> a = MatrixUtils.createFieldMatrix( data3x4);
362         FieldDecompositionSolver<T> solver = new FieldQRDecomposition<T>(a).getSolver();
363         FieldMatrix<T> pinv  = solver.getInverse();
364         assertEquals(0,norm( a.multiply(pinv).subtract(MatrixUtils.createFieldIdentityMatrix(field, 3))).getReal(), 1.0e-6);
365         assertEquals(0,norm( pinv.multiply(a).getSubMatrix(0, 2, 0, 2).subtract(MatrixUtils.createFieldIdentityMatrix(field, 3))).getReal(), 1.0e-6);
366         assertEquals(testData3x4.length,    solver.getRowDimension());
367         assertEquals(testData3x4[0].length, solver.getColumnDimension());
368     }
369 
370     private  <T extends CalculusFieldElement<T>> FieldMatrix<T> createTestMatrix(Field<T> field, final Random r, final int rows, final int columns) {
371         FieldMatrix<T> m = MatrixUtils.createFieldMatrix(field, rows, columns);
372         m.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<T>(field.getOne()){
373             @Override
374             public T visit(int row, int column,T value) {
375                 return field.getZero().add(2.0 * r.nextDouble() - 1.0);
376             }
377         });
378         return m;
379     }
380 
381     private  <T extends CalculusFieldElement<T>> void QRSingular(Field<T> field) {
382         final FieldMatrix<T> a = MatrixUtils.createFieldMatrix(convert( field, new double[][] {
383             { 1, 6, 4 }, { 2, 4, -1 }, { -1, 2, 5 }
384         }));
385         T[] vv = MathArrays.buildArray(field,3);
386         vv[0] = field.getZero().add(5);
387         vv[1] = field.getZero().add(6);
388         vv[2] = field.getZero().add(1);
389         
390         final FieldVector<T> b = new ArrayFieldVector<T>(field, vv);
391         new FieldQRDecomposition<T>(a, field.getZero().add(1.0e-15)).getSolver().solve(b);
392     }
393     
394     private <T extends CalculusFieldElement<T>> T norm(FieldMatrix<T> FM ){
395         return walkInColumnOrder(FM, new FieldMatrixPreservingVisitor<T>() {
396 
397             /** Last row index. */
398             private double endRow;
399 
400             /** Sum of absolute values on one column. */
401             private T columnSum;
402 
403             /** Maximal sum across all columns. */
404             private T maxColSum;
405 
406             /** {@inheritDoc} */
407             @Override
408             public void start(final int rows, final int columns,
409                               final int startRow, final int endRow,
410                               final int startColumn, final int endColumn) {
411                 this.endRow = endRow;
412                 columnSum   = FM.getField().getZero();
413                 maxColSum   = FM.getField().getZero();
414             }
415 
416             /** {@inheritDoc} */
417             @Override
418             public void visit(final int row, final int column, final T value) {
419                 columnSum = columnSum.add(value).abs();
420                 if (row == endRow) {
421                     maxColSum = (maxColSum.getReal() > columnSum.getReal()) ? maxColSum : columnSum ;
422                     columnSum = FM.getField().getZero();
423                 }
424             }
425 
426             /** {@inheritDoc} */
427             @Override
428             public T end() {
429                 return maxColSum;
430             }
431 
432         });
433     }
434     
435     private <T extends CalculusFieldElement<T>> T walkInColumnOrder(FieldMatrix<T> FM, FieldMatrixPreservingVisitor<T> visitor) {
436         final int rows    = FM.getRowDimension();
437         final int columns = FM.getColumnDimension();
438         visitor.start(rows, columns, 0, rows - 1, 0, columns - 1);
439         for (int column = 0; column < columns; ++column) {
440             for (int row = 0; row < rows; ++row) {
441                 visitor.visit(row, column, FM.getEntry(row, column));
442             }
443         }
444         return visitor.end();
445     }
446     
447     private <T extends CalculusFieldElement<T>> T[][] convert(Field<T> field, double[][] value){
448         T[][] res = MathArrays.buildArray(field, value.length, value[0].length);
449         for (int ii = 0; ii < (value.length); ii++){
450             for (int jj = 0; jj < (value[0].length); jj++){
451                 res[ii][jj] = field.getZero().add(value[ii][jj]);
452             }
453         }
454         return res;
455         
456     }
457 
458 
459 }