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.UnitTestUtils;
28  import org.hipparchus.complex.Complex;
29  import org.hipparchus.exception.LocalizedCoreFormats;
30  import org.hipparchus.exception.MathIllegalArgumentException;
31  import org.hipparchus.fraction.Fraction;
32  import org.hipparchus.fraction.FractionField;
33  import org.hipparchus.util.Binary64Field;
34  import org.hipparchus.util.MathArrays;
35  import org.junit.Assert;
36  import org.junit.Test;
37  
38  public class FieldLUDecompositionTest {
39      private Fraction[][] testData = {
40              { new Fraction(1), new Fraction(2), new Fraction(3)},
41              { new Fraction(2), new Fraction(5), new Fraction(3)},
42              { new Fraction(1), new Fraction(0), new Fraction(8)}
43      };
44      private Fraction[][] testDataMinus = {
45              { new Fraction(-1), new Fraction(-2), new Fraction(-3)},
46              { new Fraction(-2), new Fraction(-5), new Fraction(-3)},
47              { new Fraction(-1),  new Fraction(0), new Fraction(-8)}
48      };
49      private Fraction[][] luData = {
50              { new Fraction(2), new Fraction(3), new Fraction(3) },
51              { new Fraction(0), new Fraction(5), new Fraction(7) },
52              { new Fraction(6), new Fraction(9), new Fraction(8) }
53      };
54  
55      // singular matrices
56      private Fraction[][] singular = {
57              { new Fraction(2), new Fraction(3) },
58              { new Fraction(2), new Fraction(3) }
59      };
60      private Fraction[][] bigSingular = {
61              { new Fraction(1), new Fraction(2),   new Fraction(3),    new Fraction(4) },
62              { new Fraction(2), new Fraction(5),   new Fraction(3),    new Fraction(4) },
63              { new Fraction(7), new Fraction(3), new Fraction(256), new Fraction(1930) },
64              { new Fraction(3), new Fraction(7),   new Fraction(6),    new Fraction(8) }
65      }; // 4th row = 1st + 2nd
66  
67      /** test dimensions */
68      @Test
69      public void testDimensions() {
70          FieldMatrix<Fraction> matrix =
71              new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
72          FieldLUDecomposition<Fraction> LU = new FieldLUDecomposition<Fraction>(matrix);
73          Assert.assertEquals(testData.length, LU.getL().getRowDimension());
74          Assert.assertEquals(testData.length, LU.getL().getColumnDimension());
75          Assert.assertEquals(testData.length, LU.getU().getRowDimension());
76          Assert.assertEquals(testData.length, LU.getU().getColumnDimension());
77          Assert.assertEquals(testData.length, LU.getP().getRowDimension());
78          Assert.assertEquals(testData.length, LU.getP().getColumnDimension());
79  
80      }
81  
82      /** test non-square matrix */
83      @Test
84      public void testNonSquare() {
85          try {
86              // we don't use FractionField.getInstance() for testing purposes
87              new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(new Fraction[][] {
88                      { Fraction.ZERO, Fraction.ZERO },
89                      { Fraction.ZERO, Fraction.ZERO },
90                      { Fraction.ZERO, Fraction.ZERO }
91              }));
92              Assert.fail("Expected MathIllegalArgumentException");
93          } catch (MathIllegalArgumentException ime) {
94              Assert.assertEquals(LocalizedCoreFormats.NON_SQUARE_MATRIX, ime.getSpecifier());
95          }
96      }
97  
98      /** test PA = LU */
99      @Test
100     public void testPAEqualLU() {
101         FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
102         FieldLUDecomposition<Fraction> lu = new FieldLUDecomposition<Fraction>(matrix);
103         FieldMatrix<Fraction> l = lu.getL();
104         FieldMatrix<Fraction> u = lu.getU();
105         FieldMatrix<Fraction> p = lu.getP();
106         UnitTestUtils.assertEquals(p.multiply(matrix), l.multiply(u));
107 
108         matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testDataMinus);
109         lu = new FieldLUDecomposition<Fraction>(matrix);
110         l = lu.getL();
111         u = lu.getU();
112         p = lu.getP();
113         UnitTestUtils.assertEquals(p.multiply(matrix), l.multiply(u));
114 
115         matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), 17, 17);
116         for (int i = 0; i < matrix.getRowDimension(); ++i) {
117             matrix.setEntry(i, i, Fraction.ONE);
118         }
119         lu = new FieldLUDecomposition<Fraction>(matrix);
120         l = lu.getL();
121         u = lu.getU();
122         p = lu.getP();
123         UnitTestUtils.assertEquals(p.multiply(matrix), l.multiply(u));
124 
125         matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), singular);
126         lu = new FieldLUDecomposition<Fraction>(matrix);
127         Assert.assertFalse(lu.getSolver().isNonSingular());
128         Assert.assertNull(lu.getL());
129         Assert.assertNull(lu.getU());
130         Assert.assertNull(lu.getP());
131 
132         matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), bigSingular);
133         lu = new FieldLUDecomposition<Fraction>(matrix);
134         Assert.assertFalse(lu.getSolver().isNonSingular());
135         Assert.assertNull(lu.getL());
136         Assert.assertNull(lu.getU());
137         Assert.assertNull(lu.getP());
138 
139     }
140 
141     /** test that L is lower triangular with unit diagonal */
142     @Test
143     public void testLLowerTriangular() {
144         FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
145         FieldMatrix<Fraction> l = new FieldLUDecomposition<Fraction>(matrix).getL();
146         for (int i = 0; i < l.getRowDimension(); i++) {
147             Assert.assertEquals(Fraction.ONE, l.getEntry(i, i));
148             for (int j = i + 1; j < l.getColumnDimension(); j++) {
149                 Assert.assertEquals(Fraction.ZERO, l.getEntry(i, j));
150             }
151         }
152     }
153 
154     /** test that U is upper triangular */
155     @Test
156     public void testUUpperTriangular() {
157         FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
158         FieldMatrix<Fraction> u = new FieldLUDecomposition<Fraction>(matrix).getU();
159         for (int i = 0; i < u.getRowDimension(); i++) {
160             for (int j = 0; j < i; j++) {
161                 Assert.assertEquals(Fraction.ZERO, u.getEntry(i, j));
162             }
163         }
164     }
165 
166     /** test that P is a permutation matrix */
167     @Test
168     public void testPPermutation() {
169         FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
170         FieldMatrix<Fraction> p   = new FieldLUDecomposition<Fraction>(matrix).getP();
171 
172         FieldMatrix<Fraction> ppT = p.multiply(p.transpose());
173         FieldMatrix<Fraction> id  =
174             new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(),
175                                           p.getRowDimension(), p.getRowDimension());
176         for (int i = 0; i < id.getRowDimension(); ++i) {
177             id.setEntry(i, i, Fraction.ONE);
178         }
179         UnitTestUtils.assertEquals(id, ppT);
180 
181         for (int i = 0; i < p.getRowDimension(); i++) {
182             int zeroCount  = 0;
183             int oneCount   = 0;
184             int otherCount = 0;
185             for (int j = 0; j < p.getColumnDimension(); j++) {
186                 final Fraction e = p.getEntry(i, j);
187                 if (e.equals(Fraction.ZERO)) {
188                     ++zeroCount;
189                 } else if (e.equals(Fraction.ONE)) {
190                     ++oneCount;
191                 } else {
192                     ++otherCount;
193                 }
194             }
195             Assert.assertEquals(p.getColumnDimension() - 1, zeroCount);
196             Assert.assertEquals(1, oneCount);
197             Assert.assertEquals(0, otherCount);
198         }
199 
200         for (int j = 0; j < p.getColumnDimension(); j++) {
201             int zeroCount  = 0;
202             int oneCount   = 0;
203             int otherCount = 0;
204             for (int i = 0; i < p.getRowDimension(); i++) {
205                 final Fraction e = p.getEntry(i, j);
206                 if (e.equals(Fraction.ZERO)) {
207                     ++zeroCount;
208                 } else if (e.equals(Fraction.ONE)) {
209                     ++oneCount;
210                 } else {
211                     ++otherCount;
212                 }
213             }
214             Assert.assertEquals(p.getRowDimension() - 1, zeroCount);
215             Assert.assertEquals(1, oneCount);
216             Assert.assertEquals(0, otherCount);
217         }
218 
219     }
220 
221 
222     /** test singular */
223     @Test
224     public void testSingular() {
225         final FieldMatrix<Fraction> m = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
226         FieldLUDecomposition<Fraction> lu = new FieldLUDecomposition<Fraction>(m);
227         Assert.assertTrue(lu.getSolver().isNonSingular());
228         Assert.assertEquals(new Fraction(-1, 1), lu.getDeterminant());
229         lu = new FieldLUDecomposition<>(m.getSubMatrix(0, 1, 0, 1));
230         Assert.assertTrue(lu.getSolver().isNonSingular());
231         Assert.assertEquals(new Fraction(+1, 1), lu.getDeterminant());
232         lu = new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), singular));
233         Assert.assertFalse(lu.getSolver().isNonSingular());
234         Assert.assertEquals(new Fraction(0, 1), lu.getDeterminant());
235         lu = new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), bigSingular));
236         Assert.assertFalse(lu.getSolver().isNonSingular());
237         Assert.assertEquals(new Fraction(0, 1), lu.getDeterminant());
238         try {
239             lu.getSolver().solve(new ArrayFieldVector<>(new Fraction[] { Fraction.ONE, Fraction.ONE, Fraction.ONE, Fraction.ONE }));
240             Assert.fail("an exception should have been thrown");
241         } catch (MathIllegalArgumentException miae) {
242             Assert.assertEquals(LocalizedCoreFormats.SINGULAR_MATRIX, miae.getSpecifier());
243         }
244     }
245 
246     /** test matrices values */
247     @Test
248     public void testMatricesValues1() {
249        FieldLUDecomposition<Fraction> lu =
250             new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData));
251         FieldMatrix<Fraction> lRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
252                 { new Fraction(1),   new Fraction(0),   new Fraction(0) },
253                 { new Fraction(0.5), new Fraction(1),   new Fraction(0) },
254                 { new Fraction(0.5), new Fraction(0.2), new Fraction(1) }
255         });
256         FieldMatrix<Fraction> uRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
257                 { new Fraction(2), new Fraction(5),    new Fraction(3)   },
258                 { new Fraction(0), new Fraction(-2.5), new Fraction(6.5) },
259                 { new Fraction(0), new Fraction(0),    new Fraction(0.2) }
260         });
261         FieldMatrix<Fraction> pRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
262                 { new Fraction(0), new Fraction(1), new Fraction(0) },
263                 { new Fraction(0), new Fraction(0), new Fraction(1) },
264                 { new Fraction(1), new Fraction(0), new Fraction(0) }
265         });
266         int[] pivotRef = { 1, 2, 0 };
267 
268         // check values against known references
269         FieldMatrix<Fraction> l = lu.getL();
270         UnitTestUtils.assertEquals(lRef, l);
271         FieldMatrix<Fraction> u = lu.getU();
272         UnitTestUtils.assertEquals(uRef, u);
273         FieldMatrix<Fraction> p = lu.getP();
274         UnitTestUtils.assertEquals(pRef, p);
275         int[] pivot = lu.getPivot();
276         for (int i = 0; i < pivotRef.length; ++i) {
277             Assert.assertEquals(pivotRef[i], pivot[i]);
278         }
279 
280         // check the same cached instance is returned the second time
281         Assert.assertTrue(l == lu.getL());
282         Assert.assertTrue(u == lu.getU());
283         Assert.assertTrue(p == lu.getP());
284 
285     }
286 
287     /** test matrices values */
288     @Test
289     public void testMatricesValues2() {
290        FieldLUDecomposition<Fraction> lu =
291             new FieldLUDecomposition<Fraction>(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), luData));
292         FieldMatrix<Fraction> lRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
293                 { new Fraction(1),         new Fraction(0), new Fraction(0) },
294                 { new Fraction(0),         new Fraction(1), new Fraction(0) },
295                 { new Fraction(1.0 / 3.0), new Fraction(0), new Fraction(1) }
296         });
297         FieldMatrix<Fraction> uRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
298                 { new Fraction(6), new Fraction(9), new Fraction(8)         },
299                 { new Fraction(0), new Fraction(5), new Fraction(7)         },
300                 { new Fraction(0), new Fraction(0), new Fraction(1.0 / 3.0) }
301         });
302         FieldMatrix<Fraction> pRef = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), new Fraction[][] {
303                 { new Fraction(0), new Fraction(0), new Fraction(1) },
304                 { new Fraction(0), new Fraction(1), new Fraction(0) },
305                 { new Fraction(1), new Fraction(0), new Fraction(0) }
306         });
307         int[] pivotRef = { 2, 1, 0 };
308 
309         // check values against known references
310         FieldMatrix<Fraction> l = lu.getL();
311         UnitTestUtils.assertEquals(lRef, l);
312         FieldMatrix<Fraction> u = lu.getU();
313         UnitTestUtils.assertEquals(uRef, u);
314         FieldMatrix<Fraction> p = lu.getP();
315         UnitTestUtils.assertEquals(pRef, p);
316         int[] pivot = lu.getPivot();
317         for (int i = 0; i < pivotRef.length; ++i) {
318             Assert.assertEquals(pivotRef[i], pivot[i]);
319         }
320 
321         // check the same cached instance is returned the second time
322         Assert.assertTrue(l == lu.getL());
323         Assert.assertTrue(u == lu.getU());
324         Assert.assertTrue(p == lu.getP());
325     }
326 
327     @Test
328     public void testSignedZeroPivot() {
329         FieldMatrix<Complex> m = new Array2DRowFieldMatrix<>(new Complex[][] {
330             { new Complex(-0.0, 0.0), Complex.ONE },
331             { Complex.ONE, Complex.ZERO }
332         });
333         FieldVector<Complex> v = new ArrayFieldVector<>(new Complex[] {
334             new Complex(2, 0),
335             new Complex(0, 2)
336         });
337         FieldDecompositionSolver<Complex> solver = new FieldLUDecomposition<>(m).getSolver();
338         FieldVector<Complex> u = solver.solve(v);
339         Assert.assertEquals(u.getEntry(0), new Complex(0, 2));
340         Assert.assertEquals(u.getEntry(1), new Complex(2, 0));
341     }
342 
343     @Test
344     public void testSolve() {
345         FieldDecompositionSolver<Fraction> solver =
346                         new FieldLUDecomposer<Fraction>(e -> e.isZero()).
347                         decompose(new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(),
348                                                                       testData));
349         FieldVector<Fraction> solution = solver.solve(new ArrayFieldVector<>(new Fraction[] {
350             new Fraction(1, 2), new Fraction(2, 3), new Fraction(3,4)
351         }));
352         Assert.assertEquals(testData.length, solution.getDimension());
353         Assert.assertEquals(new Fraction(-31, 12), solution.getEntry(0));
354         Assert.assertEquals(new Fraction( 11, 12), solution.getEntry(1));
355         Assert.assertEquals(new Fraction(  5, 12), solution.getEntry(2));
356         Assert.assertEquals(testData.length,    solver.getRowDimension());
357         Assert.assertEquals(testData[0].length, solver.getColumnDimension());
358         try {
359             solver.solve(new ArrayFieldVector<>(new Fraction[] { Fraction.ONE, Fraction.ONE }));
360             Assert.fail("an exception should have been thrown");
361         } catch (MathIllegalArgumentException miae) {
362             Assert.assertEquals(LocalizedCoreFormats.DIMENSIONS_MISMATCH, miae.getSpecifier());
363         }
364     }
365 
366     @Test
367     public void testComparisonWithReal() {
368         doTestComparisonWithReal(Binary64Field.getInstance());
369     }
370 
371     private <T extends CalculusFieldElement<T>> void doTestComparisonWithReal(final Field<T> field) {
372 
373         ////////////
374         // Test with a real version
375         ////////////
376 
377         final double[][] jacobianReal = new double[][] {
378             {-1.8079069467383695, -0.5276841137999425, -0.06927544502469293, 575.7094908176842, -1864.684268657213, -820.8524955582242},
379             {1.121475385353888E-7, 4.3674817490819154E-8, 9.740062323061996E-9, -6.304893098996501E-5, 2.48921714502984E-4, 1.083365579991483E-4},
380             {4.395254068576291E-8, -1.0258110498819202E-7, -5.5724863389796155E-8, -1.4063182668276462E-4, -1.609675082956865E-5, 6.006390276284299E-6},
381             {-8.949490536614748E-10, 4.4866323700855295E-9, -1.0819706965376411E-8, -1.09340948914072E-5, 5.481570585126429E-5, -1.32190432709699E-4},
382             {2.423811020572752E-8, -1.2151249212880152E-7, 2.9303260196492917E-7, -2.6617404304907148E-6, 1.334405654416438E-5, -3.2179766387795136E-5},
383             {-5.564319851994915E-8, 2.1983585343061848E-7, -2.2238994423695564E-7, 2.768985626446657E-4, 6.781392777371218E-5, 4.0155285354156046E-5}
384         };
385 
386         final RealMatrix matrixReal = MatrixUtils.createRealMatrix(jacobianReal);
387 
388         final DecompositionSolver solverReal = new LUDecomposition(matrixReal).getSolver();
389         final RealMatrix inverseReal = solverReal.getInverse();
390 
391         ////////////
392         // Test with a field version
393         ////////////
394 
395         final T[][] jacobianField = MathArrays.buildArray(field, 6, 6);
396         for (int row = 0; row < matrixReal.getRowDimension(); row++) {
397             for (int column = 0; column < matrixReal.getColumnDimension(); column++) {
398                 jacobianField[row][column] = field.getZero().add(jacobianReal[row][column]);
399             }
400         }
401 
402         final FieldMatrix<T> matrixField = MatrixUtils.createFieldMatrix(jacobianField);
403 
404         final FieldDecompositionSolver<T> solverField = new FieldLUDecomposition<>(matrixField).getSolver();
405         final FieldMatrix<T> inverseField = solverField.getInverse();
406 
407         // Verify
408         for (int row = 0; row < inverseReal.getRowDimension(); row++) {
409             for (int column = 0; column < inverseReal.getColumnDimension(); column++) {
410                Assert.assertEquals(inverseReal.getEntry(row, column), inverseField.getEntry(row, column).getReal(), 1.0e-15); 
411             }
412         }
413 
414     }
415 
416     @Test
417     public void testIssue134() {
418 
419         FieldMatrix<Fraction> matrix = new Array2DRowFieldMatrix<Fraction>(FractionField.getInstance(), testData);
420         FieldLUDecomposition<Fraction> lu = new FieldLUDecomposition<Fraction>(matrix, e -> e.isZero(), false);
421 
422         // L
423         final FieldMatrix<Fraction> l = lu.getL();
424         for (int i = 0; i < l.getRowDimension(); i++) {
425             Assert.assertEquals(Fraction.ONE, l.getEntry(i, i));
426             for (int j = i + 1; j < l.getColumnDimension(); j++) {
427                 Assert.assertEquals(Fraction.ZERO, l.getEntry(i, j));
428             }
429         }
430 
431         // U
432         final FieldMatrix<Fraction> u = lu.getU();
433         for (int i = 0; i < u.getRowDimension(); i++) {
434             for (int j = 0; j < i; j++) {
435                 Assert.assertEquals(Fraction.ZERO, u.getEntry(i, j));
436             }
437         }
438 
439     }
440 
441 }