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.DefaultRealMatrixChangingVisitor;
023 import org.apache.commons.math.linear.DefaultRealMatrixPreservingVisitor;
024 import org.apache.commons.math.linear.BlockRealMatrix;
025 import org.apache.commons.math.linear.MatrixUtils;
026 import org.apache.commons.math.linear.MatrixVisitorException;
027 import org.apache.commons.math.linear.QRDecomposition;
028 import org.apache.commons.math.linear.QRDecompositionImpl;
029 import org.apache.commons.math.linear.RealMatrix;
030
031 import junit.framework.Test;
032 import junit.framework.TestCase;
033 import junit.framework.TestSuite;
034
035 public class QRDecompositionImplTest extends TestCase {
036 double[][] testData3x3NonSingular = {
037 { 12, -51, 4 },
038 { 6, 167, -68 },
039 { -4, 24, -41 }, };
040
041 double[][] testData3x3Singular = {
042 { 1, 4, 7, },
043 { 2, 5, 8, },
044 { 3, 6, 9, }, };
045
046 double[][] testData3x4 = {
047 { 12, -51, 4, 1 },
048 { 6, 167, -68, 2 },
049 { -4, 24, -41, 3 }, };
050
051 double[][] testData4x3 = {
052 { 12, -51, 4, },
053 { 6, 167, -68, },
054 { -4, 24, -41, },
055 { -5, 34, 7, }, };
056
057 private static final double entryTolerance = 10e-16;
058
059 private static final double normTolerance = 10e-14;
060
061 public QRDecompositionImplTest(String name) {
062 super(name);
063 }
064
065 public static Test suite() {
066 TestSuite suite = new TestSuite(QRDecompositionImplTest.class);
067 suite.setName("QRDecompositionImpl Tests");
068 return suite;
069 }
070
071 /** test dimensions */
072 public void testDimensions() {
073 checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
074
075 checkDimension(MatrixUtils.createRealMatrix(testData4x3));
076
077 checkDimension(MatrixUtils.createRealMatrix(testData3x4));
078
079 Random r = new Random(643895747384642l);
080 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
081 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
082 checkDimension(createTestMatrix(r, p, q));
083 checkDimension(createTestMatrix(r, q, p));
084
085 }
086
087 private void checkDimension(RealMatrix m) {
088 int rows = m.getRowDimension();
089 int columns = m.getColumnDimension();
090 QRDecomposition qr = new QRDecompositionImpl(m);
091 assertEquals(rows, qr.getQ().getRowDimension());
092 assertEquals(rows, qr.getQ().getColumnDimension());
093 assertEquals(rows, qr.getR().getRowDimension());
094 assertEquals(columns, qr.getR().getColumnDimension());
095 }
096
097 /** test A = QR */
098 public void testAEqualQR() {
099 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
100
101 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
102
103 checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4));
104
105 checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3));
106
107 Random r = new Random(643895747384642l);
108 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
109 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
110 checkAEqualQR(createTestMatrix(r, p, q));
111
112 checkAEqualQR(createTestMatrix(r, q, p));
113
114 }
115
116 private void checkAEqualQR(RealMatrix m) {
117 QRDecomposition qr = new QRDecompositionImpl(m);
118 double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm();
119 assertEquals(0, norm, normTolerance);
120 }
121
122 /** test the orthogonality of Q */
123 public void testQOrthogonal() {
124 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular));
125
126 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular));
127
128 checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4));
129
130 checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3));
131
132 Random r = new Random(643895747384642l);
133 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
134 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
135 checkQOrthogonal(createTestMatrix(r, p, q));
136
137 checkQOrthogonal(createTestMatrix(r, q, p));
138
139 }
140
141 private void checkQOrthogonal(RealMatrix m) {
142 QRDecomposition qr = new QRDecompositionImpl(m);
143 RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
144 double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm();
145 assertEquals(0, norm, normTolerance);
146 }
147
148 /** test that R is upper triangular */
149 public void testRUpperTriangular() {
150 RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
151 checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
152
153 matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
154 checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
155
156 matrix = MatrixUtils.createRealMatrix(testData3x4);
157 checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
158
159 matrix = MatrixUtils.createRealMatrix(testData4x3);
160 checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
161
162 Random r = new Random(643895747384642l);
163 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
164 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
165 matrix = createTestMatrix(r, p, q);
166 checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
167
168 matrix = createTestMatrix(r, p, q);
169 checkUpperTriangular(new QRDecompositionImpl(matrix).getR());
170
171 }
172
173 private void checkUpperTriangular(RealMatrix m) {
174 m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
175 @Override
176 public void visit(int row, int column, double value) {
177 if (column < row) {
178 assertEquals(0.0, value, entryTolerance);
179 }
180 }
181 });
182 }
183
184 /** test that H is trapezoidal */
185 public void testHTrapezoidal() {
186 RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
187 checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
188
189 matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
190 checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
191
192 matrix = MatrixUtils.createRealMatrix(testData3x4);
193 checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
194
195 matrix = MatrixUtils.createRealMatrix(testData4x3);
196 checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
197
198 Random r = new Random(643895747384642l);
199 int p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
200 int q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
201 matrix = createTestMatrix(r, p, q);
202 checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
203
204 matrix = createTestMatrix(r, p, q);
205 checkTrapezoidal(new QRDecompositionImpl(matrix).getH());
206
207 }
208
209 private void checkTrapezoidal(RealMatrix m) {
210 m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
211 @Override
212 public void visit(int row, int column, double value) {
213 if (column > row) {
214 assertEquals(0.0, value, entryTolerance);
215 }
216 }
217 });
218 }
219 /** test matrices values */
220 public void testMatricesValues() {
221 QRDecomposition qr =
222 new QRDecompositionImpl(MatrixUtils.createRealMatrix(testData3x3NonSingular));
223 RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] {
224 { -12.0 / 14.0, 69.0 / 175.0, -58.0 / 175.0 },
225 { -6.0 / 14.0, -158.0 / 175.0, 6.0 / 175.0 },
226 { 4.0 / 14.0, -30.0 / 175.0, -165.0 / 175.0 }
227 });
228 RealMatrix rRef = MatrixUtils.createRealMatrix(new double[][] {
229 { -14.0, -21.0, 14.0 },
230 { 0.0, -175.0, 70.0 },
231 { 0.0, 0.0, 35.0 }
232 });
233 RealMatrix hRef = MatrixUtils.createRealMatrix(new double[][] {
234 { 26.0 / 14.0, 0.0, 0.0 },
235 { 6.0 / 14.0, 648.0 / 325.0, 0.0 },
236 { -4.0 / 14.0, 36.0 / 325.0, 2.0 }
237 });
238
239 // check values against known references
240 RealMatrix q = qr.getQ();
241 assertEquals(0, q.subtract(qRef).getNorm(), 1.0e-13);
242 RealMatrix qT = qr.getQT();
243 assertEquals(0, qT.subtract(qRef.transpose()).getNorm(), 1.0e-13);
244 RealMatrix r = qr.getR();
245 assertEquals(0, r.subtract(rRef).getNorm(), 1.0e-13);
246 RealMatrix h = qr.getH();
247 assertEquals(0, h.subtract(hRef).getNorm(), 1.0e-13);
248
249 // check the same cached instance is returned the second time
250 assertTrue(q == qr.getQ());
251 assertTrue(r == qr.getR());
252 assertTrue(h == qr.getH());
253
254 }
255
256 private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
257 RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
258 m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
259 @Override
260 public double visit(int row, int column, double value)
261 throws MatrixVisitorException {
262 return 2.0 * r.nextDouble() - 1.0;
263 }
264 });
265 return m;
266 }
267
268 }