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.io.BufferedReader;
26  import java.io.DataInputStream;
27  import java.io.IOException;
28  import java.io.InputStreamReader;
29  import java.util.Random;
30  
31  import org.junit.Assert;
32  import org.junit.Test;
33  
34  
35  public class SingularValueDecompositionTest {
36  
37      private double[][] testSquare = {
38              { 24.0 / 25.0, 43.0 / 25.0 },
39              { 57.0 / 25.0, 24.0 / 25.0 }
40      };
41  
42      private double[][] testNonSquare = {
43          {  -540.0 / 625.0,  963.0 / 625.0, -216.0 / 625.0 },
44          { -1730.0 / 625.0, -744.0 / 625.0, 1008.0 / 625.0 },
45          {  -720.0 / 625.0, 1284.0 / 625.0, -288.0 / 625.0 },
46          {  -360.0 / 625.0,  192.0 / 625.0, 1756.0 / 625.0 },
47      };
48  
49      private static final double normTolerance = 10e-14;
50  
51      @Test
52      public void testMoreRows() {
53          final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
54          final int rows    = singularValues.length + 2;
55          final int columns = singularValues.length;
56          Random r = new Random(15338437322523l);
57          SingularValueDecomposition svd =
58              new SingularValueDecomposition(createTestMatrix(r, rows, columns, singularValues));
59          double[] computedSV = svd.getSingularValues();
60          Assert.assertEquals(singularValues.length, computedSV.length);
61          for (int i = 0; i < singularValues.length; ++i) {
62              Assert.assertEquals(singularValues[i], computedSV[i], 1.0e-10);
63          }
64      }
65  
66      @Test
67      public void testMoreColumns() {
68          final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
69          final int rows    = singularValues.length;
70          final int columns = singularValues.length + 2;
71          Random r = new Random(732763225836210l);
72          SingularValueDecomposition svd =
73              new SingularValueDecomposition(createTestMatrix(r, rows, columns, singularValues));
74          double[] computedSV = svd.getSingularValues();
75          Assert.assertEquals(singularValues.length, computedSV.length);
76          for (int i = 0; i < singularValues.length; ++i) {
77              Assert.assertEquals(singularValues[i], computedSV[i], 1.0e-10);
78          }
79      }
80  
81      /** test dimensions */
82      @Test
83      public void testDimensions() {
84          RealMatrix matrix = MatrixUtils.createRealMatrix(testSquare);
85          final int m = matrix.getRowDimension();
86          final int n = matrix.getColumnDimension();
87          SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
88          Assert.assertEquals(m, svd.getU().getRowDimension());
89          Assert.assertEquals(m, svd.getU().getColumnDimension());
90          Assert.assertEquals(m, svd.getS().getColumnDimension());
91          Assert.assertEquals(n, svd.getS().getColumnDimension());
92          Assert.assertEquals(n, svd.getV().getRowDimension());
93          Assert.assertEquals(n, svd.getV().getColumnDimension());
94  
95      }
96  
97      @Test
98      public void testDecomposer() {
99          MatrixDecomposer decomposer = new SingularValueDecomposer();
100         Assert.assertTrue(decomposer.decompose(MatrixUtils.createRealMatrix(testSquare)).isNonSingular());
101         Assert.assertFalse(decomposer.decompose(MatrixUtils.createRealMatrix(testNonSquare)).isNonSingular());
102     }
103 
104     /** Test based on a dimension 4 Hadamard matrix. */
105     @Test
106     public void testHadamard() {
107         RealMatrix matrix = new Array2DRowRealMatrix(new double[][] {
108                 {15.0 / 2.0,  5.0 / 2.0,  9.0 / 2.0,  3.0 / 2.0 },
109                 { 5.0 / 2.0, 15.0 / 2.0,  3.0 / 2.0,  9.0 / 2.0 },
110                 { 9.0 / 2.0,  3.0 / 2.0, 15.0 / 2.0,  5.0 / 2.0 },
111                 { 3.0 / 2.0,  9.0 / 2.0,  5.0 / 2.0, 15.0 / 2.0 }
112         }, false);
113         SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
114         Assert.assertEquals(16.0, svd.getSingularValues()[0], 1.0e-14);
115         Assert.assertEquals( 8.0, svd.getSingularValues()[1], 1.0e-14);
116         Assert.assertEquals( 4.0, svd.getSingularValues()[2], 1.0e-14);
117         Assert.assertEquals( 2.0, svd.getSingularValues()[3], 1.0e-14);
118 
119         RealMatrix fullCovariance = new Array2DRowRealMatrix(new double[][] {
120                 {  85.0 / 1024, -51.0 / 1024, -75.0 / 1024,  45.0 / 1024 },
121                 { -51.0 / 1024,  85.0 / 1024,  45.0 / 1024, -75.0 / 1024 },
122                 { -75.0 / 1024,  45.0 / 1024,  85.0 / 1024, -51.0 / 1024 },
123                 {  45.0 / 1024, -75.0 / 1024, -51.0 / 1024,  85.0 / 1024 }
124         }, false);
125         Assert.assertEquals(0.0,
126                      fullCovariance.subtract(svd.getCovariance(0.0)).getNorm1(),
127                      1.0e-14);
128 
129         RealMatrix halfCovariance = new Array2DRowRealMatrix(new double[][] {
130                 {   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024 },
131                 {  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024 },
132                 {   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024 },
133                 {  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024 }
134         }, false);
135         Assert.assertEquals(0.0,
136                      halfCovariance.subtract(svd.getCovariance(6.0)).getNorm1(),
137                      1.0e-14);
138 
139     }
140 
141     /** test A = USVt */
142     @Test
143     public void testAEqualUSVt() {
144         checkAEqualUSVt(MatrixUtils.createRealMatrix(testSquare));
145         checkAEqualUSVt(MatrixUtils.createRealMatrix(testNonSquare));
146         checkAEqualUSVt(MatrixUtils.createRealMatrix(testNonSquare).transpose());
147     }
148 
149     public void checkAEqualUSVt(final RealMatrix matrix) {
150         SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
151         RealMatrix u = svd.getU();
152         RealMatrix s = svd.getS();
153         RealMatrix v = svd.getV();
154         double norm = u.multiply(s).multiplyTransposed(v).subtract(matrix).getNorm1();
155         Assert.assertEquals(0, norm, normTolerance);
156 
157     }
158 
159     /** test that U is orthogonal */
160     @Test
161     public void testUOrthogonal() {
162         checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare)).getU());
163         checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare)).getU());
164         checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getU());
165     }
166 
167     /** test that V is orthogonal */
168     @Test
169     public void testVOrthogonal() {
170         checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare)).getV());
171         checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare)).getV());
172         checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getV());
173     }
174 
175     public void checkOrthogonal(final RealMatrix m) {
176         RealMatrix mTm = m.transposeMultiply(m);
177         RealMatrix id  = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
178         Assert.assertEquals(0, mTm.subtract(id).getNorm1(), normTolerance);
179     }
180 
181     /** test matrices values */
182     // This test is useless since whereas the columns of U and V are linked
183     // together, the actual triplet (U,S,V) is not uniquely defined.
184     public void testMatricesValues1() {
185        SingularValueDecomposition svd =
186             new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare));
187         RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
188                 { 3.0 / 5.0, -4.0 / 5.0 },
189                 { 4.0 / 5.0,  3.0 / 5.0 }
190         });
191         RealMatrix sRef = MatrixUtils.createRealMatrix(new double[][] {
192                 { 3.0, 0.0 },
193                 { 0.0, 1.0 }
194         });
195         RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
196                 { 4.0 / 5.0,  3.0 / 5.0 },
197                 { 3.0 / 5.0, -4.0 / 5.0 }
198         });
199 
200         // check values against known references
201         RealMatrix u = svd.getU();
202         Assert.assertEquals(0, u.subtract(uRef).getNorm1(), normTolerance);
203         RealMatrix s = svd.getS();
204         Assert.assertEquals(0, s.subtract(sRef).getNorm1(), normTolerance);
205         RealMatrix v = svd.getV();
206         Assert.assertEquals(0, v.subtract(vRef).getNorm1(), normTolerance);
207 
208         // check the same cached instance is returned the second time
209         Assert.assertTrue(u == svd.getU());
210         Assert.assertTrue(s == svd.getS());
211         Assert.assertTrue(v == svd.getV());
212 
213     }
214 
215     /** test matrices values */
216     // This test is useless since whereas the columns of U and V are linked
217     // together, the actual triplet (U,S,V) is not uniquely defined.
218     public void useless_testMatricesValues2() {
219 
220         RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
221             {  0.0 / 5.0,  3.0 / 5.0,  0.0 / 5.0 },
222             { -4.0 / 5.0,  0.0 / 5.0, -3.0 / 5.0 },
223             {  0.0 / 5.0,  4.0 / 5.0,  0.0 / 5.0 },
224             { -3.0 / 5.0,  0.0 / 5.0,  4.0 / 5.0 }
225         });
226         RealMatrix sRef = MatrixUtils.createRealMatrix(new double[][] {
227             { 4.0, 0.0, 0.0 },
228             { 0.0, 3.0, 0.0 },
229             { 0.0, 0.0, 2.0 }
230         });
231         RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
232             {  80.0 / 125.0,  -60.0 / 125.0, 75.0 / 125.0 },
233             {  24.0 / 125.0,  107.0 / 125.0, 60.0 / 125.0 },
234             { -93.0 / 125.0,  -24.0 / 125.0, 80.0 / 125.0 }
235         });
236 
237         // check values against known references
238         SingularValueDecomposition svd =
239             new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare));
240         RealMatrix u = svd.getU();
241         Assert.assertEquals(0, u.subtract(uRef).getNorm1(), normTolerance);
242         RealMatrix s = svd.getS();
243         Assert.assertEquals(0, s.subtract(sRef).getNorm1(), normTolerance);
244         RealMatrix v = svd.getV();
245         Assert.assertEquals(0, v.subtract(vRef).getNorm1(), normTolerance);
246 
247         // check the same cached instance is returned the second time
248         Assert.assertTrue(u == svd.getU());
249         Assert.assertTrue(s == svd.getS());
250         Assert.assertTrue(v == svd.getV());
251 
252     }
253 
254      /** test MATH-465 */
255     @Test
256     public void testRank() {
257         double[][] d = { { 1, 1, 1 }, { 0, 0, 0 }, { 1, 2, 3 } };
258         RealMatrix m = new Array2DRowRealMatrix(d);
259         SingularValueDecomposition svd = new SingularValueDecomposition(m);
260         Assert.assertEquals(2, svd.getRank());
261     }
262 
263     /** test MATH-583 */
264     @Test
265     public void testStability1() {
266         RealMatrix m = new Array2DRowRealMatrix(201, 201);
267         loadRealMatrix(m,"matrix1.csv");
268         try {
269             new SingularValueDecomposition(m);
270         } catch (Exception e) {
271             Assert.fail("Exception whilst constructing SVD");
272         }
273     }
274 
275     /** test MATH-327 */
276     @Test
277     public void testStability2() {
278         RealMatrix m = new Array2DRowRealMatrix(7, 168);
279         loadRealMatrix(m,"matrix2.csv");
280         try {
281             new SingularValueDecomposition(m);
282         } catch (Throwable e) {
283             Assert.fail("Exception whilst constructing SVD");
284         }
285     }
286 
287     private void loadRealMatrix(RealMatrix m, String resourceName) {
288         try {
289             DataInputStream in = new DataInputStream(getClass().getResourceAsStream(resourceName));
290             BufferedReader br = new BufferedReader(new InputStreamReader(in));
291             String strLine;
292             int row = 0;
293             while ((strLine = br.readLine()) != null) {
294                 if (!strLine.startsWith("#")) {
295                     int col = 0;
296                     for (String entry : strLine.split(",")) {
297                         m.setEntry(row, col++, Double.parseDouble(entry));
298                     }
299                     row++;
300                 }
301             }
302             in.close();
303         } catch (IOException e) {}
304     }
305 
306     /** test condition number */
307     @Test
308     public void testConditionNumber() {
309         SingularValueDecomposition svd =
310             new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare));
311         // replace 1.0e-15 with 1.5e-15
312         Assert.assertEquals(3.0, svd.getConditionNumber(), 1.5e-15);
313     }
314 
315     @Test
316     public void testInverseConditionNumber() {
317         SingularValueDecomposition svd =
318             new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare));
319         Assert.assertEquals(1.0/3.0, svd.getInverseConditionNumber(), 1.5e-15);
320     }
321 
322     private RealMatrix createTestMatrix(final Random r, final int rows, final int columns,
323                                         final double[] singularValues) {
324         final RealMatrix u = EigenDecompositionSymmetricTest.createOrthogonalMatrix(r, rows);
325         final RealMatrix d = new Array2DRowRealMatrix(rows, columns);
326         d.setSubMatrix(MatrixUtils.createRealDiagonalMatrix(singularValues).getData(), 0, 0);
327         final RealMatrix v = EigenDecompositionSymmetricTest.createOrthogonalMatrix(r, columns);
328         return u.multiply(d).multiply(v);
329     }
330 
331     @Test
332     public void testIssue947() {
333         double[][] nans = new double[][] {
334             { Double.NaN, Double.NaN },
335             { Double.NaN, Double.NaN }
336         };
337         RealMatrix m = new Array2DRowRealMatrix(nans, false);
338         SingularValueDecomposition svd = new SingularValueDecomposition(m);
339         Assert.assertTrue(Double.isNaN(svd.getSingularValues()[0]));
340         Assert.assertTrue(Double.isNaN(svd.getSingularValues()[1]));
341     }
342 
343 }