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 junit.framework.Test;
23  import junit.framework.TestCase;
24  import junit.framework.TestSuite;
25  
26  import org.apache.commons.math.linear.DecompositionSolver;
27  import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
28  import org.apache.commons.math.linear.BlockRealMatrix;
29  import org.apache.commons.math.linear.InvalidMatrixException;
30  import org.apache.commons.math.linear.MatrixUtils;
31  import org.apache.commons.math.linear.MatrixVisitorException;
32  import org.apache.commons.math.linear.QRDecomposition;
33  import org.apache.commons.math.linear.QRDecompositionImpl;
34  import org.apache.commons.math.linear.RealMatrix;
35  import org.apache.commons.math.linear.RealVector;
36  import org.apache.commons.math.linear.ArrayRealVector;
37  
38  public class QRSolverTest extends TestCase {
39      double[][] testData3x3NonSingular = { 
40              { 12, -51,   4 }, 
41              {  6, 167, -68 },
42              { -4,  24, -41 }
43      };
44  
45      double[][] testData3x3Singular = { 
46              { 1, 2,  2 }, 
47              { 2, 4,  6 },
48              { 4, 8, 12 }
49      };
50  
51      double[][] testData3x4 = { 
52              { 12, -51,   4, 1 }, 
53              {  6, 167, -68, 2 },
54              { -4,  24, -41, 3 }
55      };
56  
57      double[][] testData4x3 = { 
58              { 12, -51,   4 }, 
59              {  6, 167, -68 },
60              { -4,  24, -41 }, 
61              { -5,  34,   7 }
62      };
63  
64      public QRSolverTest(String name) {
65          super(name);
66      }
67  
68      public static Test suite() {
69          TestSuite suite = new TestSuite(QRSolverTest.class);
70          suite.setName("QRSolver Tests");
71          return suite;
72      }
73  
74      /** test rank */
75      public void testRank() {
76          DecompositionSolver solver =
77              new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
78          assertTrue(solver.isNonSingular());
79  
80          solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3Singular)).getSolver();
81          assertFalse(solver.isNonSingular());
82  
83          solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x4)).getSolver();
84          assertTrue(solver.isNonSingular());
85  
86          solver = new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData4x3)).getSolver();
87          assertTrue(solver.isNonSingular());
88  
89      }
90  
91      /** test solve dimension errors */
92      public void testSolveDimensionErrors() {
93          DecompositionSolver solver =
94              new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
95          RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
96          try {
97              solver.solve(b);
98              fail("an exception should have been thrown");
99          } 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 }