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 org.hipparchus.exception.MathIllegalArgumentException;
26  import org.junit.jupiter.api.Test;
27  
28  import java.util.Random;
29  
30  import static org.junit.jupiter.api.Assertions.assertEquals;
31  import static org.junit.jupiter.api.Assertions.assertFalse;
32  import static org.junit.jupiter.api.Assertions.assertTrue;
33  import static org.junit.jupiter.api.Assertions.fail;
34  
35  class RRQRSolverTest {
36      double[][] testData3x3NonSingular = {
37              { 12, -51,   4 },
38              {  6, 167, -68 },
39              { -4,  24, -41 }
40      };
41  
42      double[][] testData3x3Singular = {
43              { 1, 2,  2 },
44              { 2, 4,  6 },
45              { 4, 8, 12 }
46      };
47  
48      double[][] testData3x4 = {
49              { 12, -51,   4, 1 },
50              {  6, 167, -68, 2 },
51              { -4,  24, -41, 3 }
52      };
53  
54      double[][] testData4x3 = {
55              { 12, -51,   4 },
56              {  6, 167, -68 },
57              { -4,  24, -41 },
58              { -5,  34,   7 }
59      };
60  
61      /** test rank */
62      @Test
63      void testRank() {
64          DecompositionSolver solver =
65              new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular), 1.0e-16).getSolver();
66          assertTrue(solver.isNonSingular());
67  
68          solver = new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular), 1.0e-16).getSolver();
69          assertFalse(solver.isNonSingular());
70  
71          solver = new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x4), 1.0e-16).getSolver();
72          assertTrue(solver.isNonSingular());
73  
74          solver = new RRQRDecomposition(MatrixUtils.createRealMatrix(testData4x3), 1.0e-16).getSolver();
75          assertTrue(solver.isNonSingular());
76  
77      }
78  
79      /** test solve dimension errors */
80      @Test
81      void testSolveDimensionErrors() {
82          DecompositionSolver solver =
83              new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular)).getSolver();
84          RealMatrix b = MatrixUtils.createRealMatrix(new double[2][2]);
85          try {
86              solver.solve(b);
87              fail("an exception should have been thrown");
88          } catch (MathIllegalArgumentException iae) {
89              // expected behavior
90          }
91          try {
92              solver.solve(b.getColumnVector(0));
93              fail("an exception should have been thrown");
94          } catch (MathIllegalArgumentException iae) {
95              // expected behavior
96          }
97  
98      }
99  
100     /** test solve rank errors */
101     @Test
102     void testSolveRankErrors() {
103         DecompositionSolver solver =
104             new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular), 1.0e-16).getSolver();
105         RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]);
106         try {
107             solver.solve(b);
108             fail("an exception should have been thrown");
109         } catch (MathIllegalArgumentException iae) {
110             // expected behavior
111         }
112         try {
113             solver.solve(b.getColumnVector(0));
114             fail("an exception should have been thrown");
115         } catch (MathIllegalArgumentException iae) {
116             // expected behavior
117         }
118 
119     }
120 
121     /** test solve */
122     @Test
123     void testSolve() {
124         RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
125                 { -102, 12250 }, { 544, 24500 }, { 167, -36750 }
126         });
127         RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
128                 { 1, 2515 }, { 2, 422 }, { -3, 898 }
129         });
130 
131 
132         RRQRDecomposition decomposition = new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
133         DecompositionSolver solver = decomposition.getSolver();
134         assertEquals(testData3x3NonSingular.length, solver.getRowDimension());
135         assertEquals(testData3x3NonSingular[0].length, solver.getColumnDimension());
136 
137         // using RealMatrix
138         assertEquals(0, solver.solve(b).subtract(xRef).getNorm1(), 4.0e-16 * xRef.getNorm1());
139 
140         // using ArrayRealVector
141         for (int i = 0; i < b.getColumnDimension(); ++i) {
142             final RealVector x = solver.solve(b.getColumnVector(i));
143             final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
144             assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm());
145         }
146 
147         // using RealVector with an alternate implementation
148         for (int i = 0; i < b.getColumnDimension(); ++i) {
149             ArrayRealVectorTest.RealVectorTestImpl v =
150                 new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
151             final RealVector x = solver.solve(v);
152             final double error = x.subtract(xRef.getColumnVector(i)).getNorm();
153             assertEquals(0, error, 3.0e-16 * xRef.getColumnVector(i).getNorm());
154         }
155 
156     }
157 
158     @Test
159     void testOverdetermined() {
160         final Random r    = new Random(5559252868205245l);
161         int          p    = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
162         int          q    = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
163         RealMatrix   a    = createTestMatrix(r, p, q);
164         RealMatrix   xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);
165 
166         // build a perturbed system: A.X + noise = B
167         RealMatrix b = a.multiply(xRef);
168         final double noise = 0.001;
169         b.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
170             @Override
171             public double visit(int row, int column, double value) {
172                 return value * (1.0 + noise * (2 * r.nextDouble() - 1));
173             }
174         });
175 
176         // despite perturbation, the least square solution should be pretty good
177         RealMatrix x = new RRQRDecomposition(a).getSolver().solve(b);
178         assertEquals(0, x.subtract(xRef).getNorm1(), 0.01 * noise * p * q);
179 
180     }
181 
182     @Test
183     void testUnderdetermined() {
184         final Random r    = new Random(42185006424567123l);
185         int          p    = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
186         int          q    = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
187         RealMatrix   a    = createTestMatrix(r, p, q);
188         RealMatrix   xRef = createTestMatrix(r, q, BlockRealMatrix.BLOCK_SIZE + 3);
189         RealMatrix   b    = a.multiply(xRef);
190         RRQRDecomposition rrqrd = new RRQRDecomposition(a);
191         RealMatrix   x = rrqrd.getSolver().solve(b);
192 
193         // too many equations, the system cannot be solved at all
194         assertTrue(x.subtract(xRef).getNorm1() / (p * q) > 0.01);
195 
196         // the last permuted unknown should have been set to 0
197         RealMatrix permuted = rrqrd.getP().transposeMultiply(x);
198         assertEquals(0.0, permuted.getSubMatrix(p, q - 1, 0, permuted.getColumnDimension() - 1).getNorm1(), 0);
199 
200     }
201 
202     private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
203         RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
204         m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
205                 @Override
206                     public double visit(int row, int column, double value) {
207                     return 2.0 * r.nextDouble() - 1.0;
208                 }
209             });
210         return m;
211     }
212 }