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 junit.framework.Test;
023    import junit.framework.TestCase;
024    import junit.framework.TestSuite;
025    
026    import org.apache.commons.math.linear.DecompositionSolver;
027    import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
028    import org.apache.commons.math.linear.BlockRealMatrix;
029    import org.apache.commons.math.linear.InvalidMatrixException;
030    import org.apache.commons.math.linear.MatrixUtils;
031    import org.apache.commons.math.linear.MatrixVisitorException;
032    import org.apache.commons.math.linear.QRDecomposition;
033    import org.apache.commons.math.linear.QRDecompositionImpl;
034    import org.apache.commons.math.linear.RealMatrix;
035    import org.apache.commons.math.linear.RealVector;
036    import org.apache.commons.math.linear.ArrayRealVector;
037    
038    public class QRSolverTest extends TestCase {
039        double[][] testData3x3NonSingular = { 
040                { 12, -51,   4 }, 
041                {  6, 167, -68 },
042                { -4,  24, -41 }
043        };
044    
045        double[][] testData3x3Singular = { 
046                { 1, 2,  2 }, 
047                { 2, 4,  6 },
048                { 4, 8, 12 }
049        };
050    
051        double[][] testData3x4 = { 
052                { 12, -51,   4, 1 }, 
053                {  6, 167, -68, 2 },
054                { -4,  24, -41, 3 }
055        };
056    
057        double[][] testData4x3 = { 
058                { 12, -51,   4 }, 
059                {  6, 167, -68 },
060                { -4,  24, -41 }, 
061                { -5,  34,   7 }
062        };
063    
064        public QRSolverTest(String name) {
065            super(name);
066        }
067    
068        public static Test suite() {
069            TestSuite suite = new TestSuite(QRSolverTest.class);
070            suite.setName("QRSolver Tests");
071            return suite;
072        }
073    
074        /** test rank */
075        public void testRank() {
076            DecompositionSolver solver =
077                new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
078            assertTrue(solver.isNonSingular());
079    
080            solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
081            assertFalse(solver.isNonSingular());
082    
083            solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
084            assertTrue(solver.isNonSingular());
085    
086            solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
087            assertTrue(solver.isNonSingular());
088    
089        }
090    
091        /** test solve dimension errors */
092        public void testSolveDimensionErrors() {
093            DecompositionSolver solver =
094                new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
095            RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
096            try {
097                solver.solve(b);
098                fail("an exception should have been thrown");
099            } catch (IllegalArgumentException iae) {
100                // expected behavior
101            } catch (Exception e) {
102                fail("wrong exception caught");
103            }
104            try {
105                solver.solve(b.getColumn(0));
106                fail("an exception should have been thrown");
107            } catch (IllegalArgumentException iae) {
108                // expected behavior
109            } catch (Exception e) {
110                fail("wrong exception caught");
111            }
112            try {
113                solver.solve(b.getColumnVector(0));
114                fail("an exception should have been thrown");
115            } catch (IllegalArgumentException iae) {
116                // expected behavior
117            } catch (Exception e) {
118                fail("wrong exception caught");
119            }
120        }
121    
122        /** test solve rank errors */
123        public void testSolveRankErrors() {
124            DecompositionSolver solver =
125                new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
126            RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]);
127            try {
128                solver.solve(b);
129                fail("an exception should have been thrown");
130            } catch (InvalidMatrixException iae) {
131                // expected behavior
132            } catch (Exception e) {
133                fail("wrong exception caught");
134            }
135            try {
136                solver.solve(b.getColumn(0));
137                fail("an exception should have been thrown");
138            } catch (InvalidMatrixException iae) {
139                // expected behavior
140            } catch (Exception e) {
141                fail("wrong exception caught");
142            }
143            try {
144                solver.solve(b.getColumnVector(0));
145                fail("an exception should have been thrown");
146            } catch (InvalidMatrixException iae) {
147                // expected behavior
148            } catch (Exception e) {
149                fail("wrong exception caught");
150            }
151        }
152    
153        /** test solve */
154        public void testSolve() {
155            QRDecomposition decomposition =
156                new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular));
157            DecompositionSolver solver = decomposition.getSolver();
158            RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
159                    { -102, 12250 }, { 544, 24500 }, { 167, -36750 }
160            });
161            RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
162                    { 1, 2515 }, { 2, 422 }, { -3, 898 }
163            });
164    
165            // using RealMatrix
166            assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 2.0e-16 * xRef.getNorm());
167    
168            // using double[]
169            for (int i = 0; i < b.getColumnDimension(); ++i) {
170                final double[] x = solver.solve(b.getColumn(i));
171                final double error = new ArrayRealVector(x).subtract(xRef.getColumnVector(i)).getNorm();
172                assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm());
173            }
174    
175            // using ArrayRealVector
176            for (int i = 0; i < b.getColumnDimension(); ++i) {
177                final RealVector x = solver.solve(b.getColumnVector(i));
178                final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
179                assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm());
180            }
181    
182            // using RealVector with an alternate implementation
183            for (int i = 0; i < b.getColumnDimension(); ++i) {
184                ArrayRealVectorTest.RealVectorTestImpl v =
185                    new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
186                final RealVector x = solver.solve(v);
187                final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
188                assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm());
189            }
190    
191        }
192    
193        public void testOverdetermined() {
194            final Random r    = new Random(5559252868205245l);
195            int          p    = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
196            int          q    = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
197            RealMatrix   a    = createTestMatrix(r, p, q);
198            RealMatrix   xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);
199    
200            // build a perturbed system: A.X + noise = B
201            RealMatrix b = a.multiply(xRef);
202            final double noise = 0.001;
203            b.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
204                @Override
205                public double visit(int row, int column, double value) {
206                    return value * (1.0 + noise * (2 * r.nextDouble() - 1));
207                }
208            });
209    
210            // despite perturbation, the least square solution should be pretty good
211            RealMatrix x = new QRDecompositionImpl(a).getSolver().solve(b);
212            assertEquals(0, x.subtract(xRef).getNorm(), 0.01 * noise * p * q);
213    
214        }
215    
216        public void testUnderdetermined() {
217            final Random r    = new Random(42185006424567123l);
218            int          p    = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
219            int          q    = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
220            RealMatrix   a    = createTestMatrix(r, p, q);
221            RealMatrix   xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);
222            RealMatrix   b    = a.multiply(xRef);
223            RealMatrix   x = new QRDecompositionImpl(a).getSolver().solve(b);
224    
225            // too many equations, the system cannot be solved at all
226            assertTrue(x.subtract(xRef).getNorm() / (p * q) > 0.01);
227    
228            // the last unknown should have been set to 0
229            assertEquals(0.0, x.getSubMatrix(p, q - 1, 0, x.getColumnDimension() - 1).getNorm());
230    
231        }
232    
233        private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
234            RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
235            m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
236                @Override
237                public double visit(int row, int column, double value)
238                    throws MatrixVisitorException {
239                    return 2.0 * r.nextDouble() - 1.0;
240                }
241            });
242            return m;
243        }
244    }