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    *      http://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  package org.apache.commons.math4.legacy.linear;
19  
20  import org.junit.Test;
21  import org.junit.Assert;
22  
23  public class LUDecompositionTest {
24      private double[][] testData = {
25              { 1.0, 2.0, 3.0},
26              { 2.0, 5.0, 3.0},
27              { 1.0, 0.0, 8.0}
28      };
29      private double[][] testDataMinus = {
30              { -1.0, -2.0, -3.0},
31              { -2.0, -5.0, -3.0},
32              { -1.0,  0.0, -8.0}
33      };
34      private double[][] luData = {
35              { 2.0, 3.0, 3.0 },
36              { 0.0, 5.0, 7.0 },
37              { 6.0, 9.0, 8.0 }
38      };
39  
40      // singular matrices
41      private double[][] singular = {
42              { 2.0, 3.0 },
43              { 2.0, 3.0 }
44      };
45      private double[][] bigSingular = {
46              { 1.0, 2.0,   3.0,    4.0 },
47              { 2.0, 5.0,   3.0,    4.0 },
48              { 7.0, 3.0, 256.0, 1930.0 },
49              { 3.0, 7.0,   6.0,    8.0 }
50      }; // 4th row = 1st + 2nd
51  
52      private static final double entryTolerance = 10e-16;
53  
54      private static final double normTolerance = 10e-14;
55  
56      /** test dimensions */
57      @Test
58      public void testDimensions() {
59          RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
60          LUDecomposition LU = new LUDecomposition(matrix);
61          Assert.assertEquals(testData.length, LU.getL().getRowDimension());
62          Assert.assertEquals(testData.length, LU.getL().getColumnDimension());
63          Assert.assertEquals(testData.length, LU.getU().getRowDimension());
64          Assert.assertEquals(testData.length, LU.getU().getColumnDimension());
65          Assert.assertEquals(testData.length, LU.getP().getRowDimension());
66          Assert.assertEquals(testData.length, LU.getP().getColumnDimension());
67      }
68  
69      /** test non-square matrix */
70      @Test
71      public void testNonSquare() {
72          try {
73              new LUDecomposition(MatrixUtils.createRealMatrix(new double[3][2]));
74              Assert.fail("Expecting NonSquareMatrixException");
75          } catch (NonSquareMatrixException ime) {
76              // expected behavior
77          }
78      }
79  
80      /** test PA = LU */
81      @Test
82      public void testPAEqualLU() {
83          RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
84          LUDecomposition lu = new LUDecomposition(matrix);
85          RealMatrix l = lu.getL();
86          RealMatrix u = lu.getU();
87          RealMatrix p = lu.getP();
88          double norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
89          Assert.assertEquals(0, norm, normTolerance);
90  
91          matrix = MatrixUtils.createRealMatrix(testDataMinus);
92          lu = new LUDecomposition(matrix);
93          l = lu.getL();
94          u = lu.getU();
95          p = lu.getP();
96          norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
97          Assert.assertEquals(0, norm, normTolerance);
98  
99          matrix = MatrixUtils.createRealIdentityMatrix(17);
100         lu = new LUDecomposition(matrix);
101         l = lu.getL();
102         u = lu.getU();
103         p = lu.getP();
104         norm = l.multiply(u).subtract(p.multiply(matrix)).getNorm();
105         Assert.assertEquals(0, norm, normTolerance);
106 
107         matrix = MatrixUtils.createRealMatrix(singular);
108         lu = new LUDecomposition(matrix);
109         Assert.assertFalse(lu.getSolver().isNonSingular());
110         Assert.assertNull(lu.getL());
111         Assert.assertNull(lu.getU());
112         Assert.assertNull(lu.getP());
113 
114         matrix = MatrixUtils.createRealMatrix(bigSingular);
115         lu = new LUDecomposition(matrix);
116         Assert.assertFalse(lu.getSolver().isNonSingular());
117         Assert.assertNull(lu.getL());
118         Assert.assertNull(lu.getU());
119         Assert.assertNull(lu.getP());
120     }
121 
122     /** test that L is lower triangular with unit diagonal */
123     @Test
124     public void testLLowerTriangular() {
125         RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
126         RealMatrix l = new LUDecomposition(matrix).getL();
127         for (int i = 0; i < l.getRowDimension(); i++) {
128             Assert.assertEquals(l.getEntry(i, i), 1, entryTolerance);
129             for (int j = i + 1; j < l.getColumnDimension(); j++) {
130                 Assert.assertEquals(l.getEntry(i, j), 0, entryTolerance);
131             }
132         }
133     }
134 
135     /** test that U is upper triangular */
136     @Test
137     public void testUUpperTriangular() {
138         RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
139         RealMatrix u = new LUDecomposition(matrix).getU();
140         for (int i = 0; i < u.getRowDimension(); i++) {
141             for (int j = 0; j < i; j++) {
142                 Assert.assertEquals(u.getEntry(i, j), 0, entryTolerance);
143             }
144         }
145     }
146 
147     /** test that P is a permutation matrix */
148     @Test
149     public void testPPermutation() {
150         RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
151         RealMatrix p   = new LUDecomposition(matrix).getP();
152 
153         RealMatrix ppT = p.multiply(p.transpose());
154         RealMatrix id  = MatrixUtils.createRealIdentityMatrix(p.getRowDimension());
155         Assert.assertEquals(0, ppT.subtract(id).getNorm(), normTolerance);
156 
157         for (int i = 0; i < p.getRowDimension(); i++) {
158             int zeroCount  = 0;
159             int oneCount   = 0;
160             int otherCount = 0;
161             for (int j = 0; j < p.getColumnDimension(); j++) {
162                 final double e = p.getEntry(i, j);
163                 if (e == 0) {
164                     ++zeroCount;
165                 } else if (e == 1) {
166                     ++oneCount;
167                 } else {
168                     ++otherCount;
169                 }
170             }
171             Assert.assertEquals(p.getColumnDimension() - 1, zeroCount);
172             Assert.assertEquals(1, oneCount);
173             Assert.assertEquals(0, otherCount);
174         }
175 
176         for (int j = 0; j < p.getColumnDimension(); j++) {
177             int zeroCount  = 0;
178             int oneCount   = 0;
179             int otherCount = 0;
180             for (int i = 0; i < p.getRowDimension(); i++) {
181                 final double e = p.getEntry(i, j);
182                 if (e == 0) {
183                     ++zeroCount;
184                 } else if (e == 1) {
185                     ++oneCount;
186                 } else {
187                     ++otherCount;
188                 }
189             }
190             Assert.assertEquals(p.getRowDimension() - 1, zeroCount);
191             Assert.assertEquals(1, oneCount);
192             Assert.assertEquals(0, otherCount);
193         }
194     }
195 
196     /** test singular */
197     @Test
198     public void testSingular() {
199         LUDecomposition lu =
200             new LUDecomposition(MatrixUtils.createRealMatrix(testData));
201         Assert.assertTrue(lu.getSolver().isNonSingular());
202         lu = new LUDecomposition(MatrixUtils.createRealMatrix(singular));
203         Assert.assertFalse(lu.getSolver().isNonSingular());
204         lu = new LUDecomposition(MatrixUtils.createRealMatrix(bigSingular));
205         Assert.assertFalse(lu.getSolver().isNonSingular());
206     }
207 
208     /** test matrices values */
209     @Test
210     public void testMatricesValues1() {
211        LUDecomposition lu =
212             new LUDecomposition(MatrixUtils.createRealMatrix(testData));
213         RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
214                 { 1.0, 0.0, 0.0 },
215                 { 0.5, 1.0, 0.0 },
216                 { 0.5, 0.2, 1.0 }
217         });
218         RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
219                 { 2.0,  5.0, 3.0 },
220                 { 0.0, -2.5, 6.5 },
221                 { 0.0,  0.0, 0.2 }
222         });
223         RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
224                 { 0.0, 1.0, 0.0 },
225                 { 0.0, 0.0, 1.0 },
226                 { 1.0, 0.0, 0.0 }
227         });
228         int[] pivotRef = { 1, 2, 0 };
229 
230         // check values against known references
231         RealMatrix l = lu.getL();
232         Assert.assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
233         RealMatrix u = lu.getU();
234         Assert.assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
235         RealMatrix p = lu.getP();
236         Assert.assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
237         int[] pivot = lu.getPivot();
238         for (int i = 0; i < pivotRef.length; ++i) {
239             Assert.assertEquals(pivotRef[i], pivot[i]);
240         }
241 
242         // check the same cached instance is returned the second time
243         Assert.assertSame(l, lu.getL());
244         Assert.assertSame(u, lu.getU());
245         Assert.assertSame(p, lu.getP());
246     }
247 
248     /** test matrices values */
249     @Test
250     public void testMatricesValues2() {
251        LUDecomposition lu =
252             new LUDecomposition(MatrixUtils.createRealMatrix(luData));
253         RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
254                 {    1.0,    0.0, 0.0 },
255                 {    0.0,    1.0, 0.0 },
256                 { 1.0 / 3.0, 0.0, 1.0 }
257         });
258         RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
259                 { 6.0, 9.0,    8.0    },
260                 { 0.0, 5.0,    7.0    },
261                 { 0.0, 0.0, 1.0 / 3.0 }
262         });
263         RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
264                 { 0.0, 0.0, 1.0 },
265                 { 0.0, 1.0, 0.0 },
266                 { 1.0, 0.0, 0.0 }
267         });
268         int[] pivotRef = { 2, 1, 0 };
269 
270         // check values against known references
271         RealMatrix l = lu.getL();
272         Assert.assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
273         RealMatrix u = lu.getU();
274         Assert.assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
275         RealMatrix p = lu.getP();
276         Assert.assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
277         int[] pivot = lu.getPivot();
278         for (int i = 0; i < pivotRef.length; ++i) {
279             Assert.assertEquals(pivotRef[i], pivot[i]);
280         }
281 
282         // check the same cached instance is returned the second time
283         Assert.assertSame(l, lu.getL());
284         Assert.assertSame(u, lu.getU());
285         Assert.assertSame(p, lu.getP());
286     }
287 }