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