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.math.linear;
19  
20  import java.util.Random;
21  
22  import org.apache.commons.math.linear.MatrixUtils;
23  import org.apache.commons.math.linear.RealMatrix;
24  import org.apache.commons.math.linear.SingularValueDecomposition;
25  import org.apache.commons.math.linear.SingularValueDecompositionImpl;
26  
27  import junit.framework.Test;
28  import junit.framework.TestCase;
29  import junit.framework.TestSuite;
30  
31  public class SingularValueDecompositionImplTest extends TestCase {
32  
33      private double[][] testSquare = {
34              { 24.0 / 25.0, 43.0 / 25.0 },
35              { 57.0 / 25.0, 24.0 / 25.0 }
36      };
37  
38      private double[][] testNonSquare = {
39          {  -540.0 / 625.0,  963.0 / 625.0, -216.0 / 625.0 },
40          { -1730.0 / 625.0, -744.0 / 625.0, 1008.0 / 625.0 },
41          {  -720.0 / 625.0, 1284.0 / 625.0, -288.0 / 625.0 },
42          {  -360.0 / 625.0,  192.0 / 625.0, 1756.0 / 625.0 },
43      };
44  
45      private static final double normTolerance = 10e-14;
46  
47      public SingularValueDecompositionImplTest(String name) {
48          super(name);
49      }
50  
51      public static Test suite() {
52          TestSuite suite = new TestSuite(SingularValueDecompositionImplTest.class);
53          suite.setName("SingularValueDecompositionImpl Tests");
54          return suite;
55      }
56  
57      public void testMoreRows() {
58          final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
59          final int rows    = singularValues.length + 2;
60          final int columns = singularValues.length;
61          Random r = new Random(15338437322523l);
62          SingularValueDecomposition svd =
63              new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
64          double[] computedSV = svd.getSingularValues();
65          assertEquals(singularValues.length, computedSV.length);
66          for (int i = 0; i < singularValues.length; ++i) {
67              assertEquals(singularValues[i], computedSV[i], 1.0e-10);
68          }
69      }
70  
71      public void testMoreColumns() {
72          final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
73          final int rows    = singularValues.length;
74          final int columns = singularValues.length + 2;
75          Random r = new Random(732763225836210l);
76          SingularValueDecomposition svd =
77              new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
78          double[] computedSV = svd.getSingularValues();
79          assertEquals(singularValues.length, computedSV.length);
80          for (int i = 0; i < singularValues.length; ++i) {
81              assertEquals(singularValues[i], computedSV[i], 1.0e-10);
82          }
83      }
84  
85      /** test dimensions */
86      public void testDimensions() {
87          RealMatrix matrix = MatrixUtils.createRealMatrix(testSquare);
88          final int m = matrix.getRowDimension();
89          final int n = matrix.getColumnDimension();
90          SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
91          assertEquals(m, svd.getU().getRowDimension());
92          assertEquals(m, svd.getU().getColumnDimension());
93          assertEquals(m, svd.getS().getColumnDimension());
94          assertEquals(n, svd.getS().getColumnDimension());
95          assertEquals(n, svd.getV().getRowDimension());
96          assertEquals(n, svd.getV().getColumnDimension());
97  
98      }
99  
100     /** Test based on a dimension 4 Hadamard matrix. */
101     public void testHadamard() {
102         RealMatrix matrix = new Array2DRowRealMatrix(new double[][] {
103                 {15.0 / 2.0,  5.0 / 2.0,  9.0 / 2.0,  3.0 / 2.0 },
104                 { 5.0 / 2.0, 15.0 / 2.0,  3.0 / 2.0,  9.0 / 2.0 },
105                 { 9.0 / 2.0,  3.0 / 2.0, 15.0 / 2.0,  5.0 / 2.0 },
106                 { 3.0 / 2.0,  9.0 / 2.0,  5.0 / 2.0, 15.0 / 2.0 }
107         }, false);
108         SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
109         assertEquals(16.0, svd.getSingularValues()[0], 1.0e-14);
110         assertEquals( 8.0, svd.getSingularValues()[1], 1.0e-14);
111         assertEquals( 4.0, svd.getSingularValues()[2], 1.0e-14);
112         assertEquals( 2.0, svd.getSingularValues()[3], 1.0e-14);
113 
114         RealMatrix fullCovariance = new Array2DRowRealMatrix(new double[][] {
115                 {  85.0 / 1024, -51.0 / 1024, -75.0 / 1024,  45.0 / 1024 },
116                 { -51.0 / 1024,  85.0 / 1024,  45.0 / 1024, -75.0 / 1024 },
117                 { -75.0 / 1024,  45.0 / 1024,  85.0 / 1024, -51.0 / 1024 },
118                 {  45.0 / 1024, -75.0 / 1024, -51.0 / 1024,  85.0 / 1024 }
119         }, false);
120         assertEquals(0.0,
121                      fullCovariance.subtract(svd.getCovariance(0.0)).getNorm(),
122                      1.0e-14);
123 
124         RealMatrix halfCovariance = new Array2DRowRealMatrix(new double[][] {
125                 {   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024 },
126                 {  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024 },
127                 {   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024 },
128                 {  -3.0 / 1024,   5.0 / 1024,  -3.0 / 1024,   5.0 / 1024 }
129         }, false);
130         assertEquals(0.0,
131                      halfCovariance.subtract(svd.getCovariance(6.0)).getNorm(),
132                      1.0e-14);
133 
134     }
135 
136     /** test A = USVt */
137     public void testAEqualUSVt() {
138         checkAEqualUSVt(MatrixUtils.createRealMatrix(testSquare));
139         checkAEqualUSVt(MatrixUtils.createRealMatrix(testNonSquare));
140         checkAEqualUSVt(MatrixUtils.createRealMatrix(testNonSquare).transpose());
141     }
142 
143     public void checkAEqualUSVt(final RealMatrix matrix) {
144         SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
145         RealMatrix u = svd.getU();
146         RealMatrix s = svd.getS();
147         RealMatrix v = svd.getV();
148         double norm = u.multiply(s).multiply(v.transpose()).subtract(matrix).getNorm();
149         assertEquals(0, norm, normTolerance);
150 
151     }
152 
153     /** test that U is orthogonal */
154     public void testUOrthogonal() {
155         checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getU());
156         checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare)).getU());
157         checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getU());
158     }
159 
160     /** test that V is orthogonal */
161     public void testVOrthogonal() {
162         checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getV());
163         checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare)).getV());
164         checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getV());
165     }
166 
167     public void checkOrthogonal(final RealMatrix m) {
168         RealMatrix mTm = m.transpose().multiply(m);
169         RealMatrix id  = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
170         assertEquals(0, mTm.subtract(id).getNorm(), normTolerance);
171     }
172 
173     /** test matrices values */
174     public void testMatricesValues1() {
175        SingularValueDecomposition svd =
176             new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
177         RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
178                 { 3.0 / 5.0, -4.0 / 5.0 },
179                 { 4.0 / 5.0,  3.0 / 5.0 }
180         });
181         RealMatrix sRef = MatrixUtils.createRealMatrix(new double[][] {
182                 { 3.0, 0.0 },
183                 { 0.0, 1.0 }
184         });
185         RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
186                 { 4.0 / 5.0,  3.0 / 5.0 },
187                 { 3.0 / 5.0, -4.0 / 5.0 }
188         });
189 
190         // check values against known references
191         RealMatrix u = svd.getU();
192         assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
193         RealMatrix s = svd.getS();
194         assertEquals(0, s.subtract(sRef).getNorm(), normTolerance);
195         RealMatrix v = svd.getV();
196         assertEquals(0, v.subtract(vRef).getNorm(), normTolerance);
197 
198         // check the same cached instance is returned the second time
199         assertTrue(u == svd.getU());
200         assertTrue(s == svd.getS());
201         assertTrue(v == svd.getV());
202         
203     }
204 
205     /** test matrices values */
206     public void testMatricesValues2() {
207 
208         RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
209             {  0.0 / 5.0,  3.0 / 5.0,  0.0 / 5.0 },
210             { -4.0 / 5.0,  0.0 / 5.0, -3.0 / 5.0 },
211             {  0.0 / 5.0,  4.0 / 5.0,  0.0 / 5.0 },
212             { -3.0 / 5.0,  0.0 / 5.0,  4.0 / 5.0 }
213         });
214         RealMatrix sRef = MatrixUtils.createRealMatrix(new double[][] {
215             { 4.0, 0.0, 0.0 },
216             { 0.0, 3.0, 0.0 },
217             { 0.0, 0.0, 2.0 }
218         });
219         RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
220             {  80.0 / 125.0,  -60.0 / 125.0, 75.0 / 125.0 },
221             {  24.0 / 125.0,  107.0 / 125.0, 60.0 / 125.0 },
222             { -93.0 / 125.0,  -24.0 / 125.0, 80.0 / 125.0 }
223         });
224 
225         // check values against known references
226         SingularValueDecomposition svd =
227             new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare));
228         RealMatrix u = svd.getU();
229         assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
230         RealMatrix s = svd.getS();
231         assertEquals(0, s.subtract(sRef).getNorm(), normTolerance);
232         RealMatrix v = svd.getV();
233         assertEquals(0, v.subtract(vRef).getNorm(), normTolerance);
234 
235         // check the same cached instance is returned the second time
236         assertTrue(u == svd.getU());
237         assertTrue(s == svd.getS());
238         assertTrue(v == svd.getV());
239 
240     }
241 
242     /** test condition number */
243     public void testConditionNumber() {
244         SingularValueDecompositionImpl svd =
245             new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
246         assertEquals(3.0, svd.getConditionNumber(), 1.0e-15);
247     }
248 
249     private RealMatrix createTestMatrix(final Random r, final int rows, final int columns,
250                                         final double[] singularValues) {
251         final RealMatrix u =
252             EigenDecompositionImplTest.createOrthogonalMatrix(r, rows);
253         final RealMatrix d =
254             EigenDecompositionImplTest.createDiagonalMatrix(singularValues, rows, columns);
255         final RealMatrix v =
256             EigenDecompositionImplTest.createOrthogonalMatrix(r, columns);
257         return u.multiply(d).multiply(v);
258     }
259 
260 }