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