001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.linear;
019    
020    import java.util.Random;
021    
022    import org.apache.commons.math.linear.MatrixUtils;
023    import org.apache.commons.math.linear.RealMatrix;
024    import org.apache.commons.math.linear.SingularValueDecomposition;
025    import org.apache.commons.math.linear.SingularValueDecompositionImpl;
026    
027    import junit.framework.Test;
028    import junit.framework.TestCase;
029    import junit.framework.TestSuite;
030    
031    public class SingularValueDecompositionImplTest extends TestCase {
032    
033        private double[][] testSquare = {
034                { 24.0 / 25.0, 43.0 / 25.0 },
035                { 57.0 / 25.0, 24.0 / 25.0 }
036        };
037    
038        private double[][] testNonSquare = {
039            {  -540.0 / 625.0,  963.0 / 625.0, -216.0 / 625.0 },
040            { -1730.0 / 625.0, -744.0 / 625.0, 1008.0 / 625.0 },
041            {  -720.0 / 625.0, 1284.0 / 625.0, -288.0 / 625.0 },
042            {  -360.0 / 625.0,  192.0 / 625.0, 1756.0 / 625.0 },
043        };
044    
045        private static final double normTolerance = 10e-14;
046    
047        public SingularValueDecompositionImplTest(String name) {
048            super(name);
049        }
050    
051        public static Test suite() {
052            TestSuite suite = new TestSuite(SingularValueDecompositionImplTest.class);
053            suite.setName("SingularValueDecompositionImpl Tests");
054            return suite;
055        }
056    
057        public void testMoreRows() {
058            final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
059            final int rows    = singularValues.length + 2;
060            final int columns = singularValues.length;
061            Random r = new Random(15338437322523l);
062            SingularValueDecomposition svd =
063                new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
064            double[] computedSV = svd.getSingularValues();
065            assertEquals(singularValues.length, computedSV.length);
066            for (int i = 0; i < singularValues.length; ++i) {
067                assertEquals(singularValues[i], computedSV[i], 1.0e-10);
068            }
069        }
070    
071        public void testMoreColumns() {
072            final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
073            final int rows    = singularValues.length;
074            final int columns = singularValues.length + 2;
075            Random r = new Random(732763225836210l);
076            SingularValueDecomposition svd =
077                new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
078            double[] computedSV = svd.getSingularValues();
079            assertEquals(singularValues.length, computedSV.length);
080            for (int i = 0; i < singularValues.length; ++i) {
081                assertEquals(singularValues[i], computedSV[i], 1.0e-10);
082            }
083        }
084    
085        /** test dimensions */
086        public void testDimensions() {
087            RealMatrix matrix = MatrixUtils.createRealMatrix(testSquare);
088            final int m = matrix.getRowDimension();
089            final int n = matrix.getColumnDimension();
090            SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
091            assertEquals(m, svd.getU().getRowDimension());
092            assertEquals(m, svd.getU().getColumnDimension());
093            assertEquals(m, svd.getS().getColumnDimension());
094            assertEquals(n, svd.getS().getColumnDimension());
095            assertEquals(n, svd.getV().getRowDimension());
096            assertEquals(n, svd.getV().getColumnDimension());
097    
098        }
099    
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    }