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 }