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    *      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.math4.legacy.linear;
19  
20  import java.util.Random;
21  
22  import org.junit.Assert;
23  import org.junit.Test;
24  
25  
26  public class QRDecompositionTest {
27      private double[][] testData3x3NonSingular = {
28              { 12, -51, 4 },
29              { 6, 167, -68 },
30              { -4, 24, -41 }, };
31  
32      private double[][] testData3x3Singular = {
33              { 1, 4, 7, },
34              { 2, 5, 8, },
35              { 3, 6, 9, }, };
36  
37      private double[][] testData3x4 = {
38              { 12, -51, 4, 1 },
39              { 6, 167, -68, 2 },
40              { -4, 24, -41, 3 }, };
41  
42      private double[][] testData4x3 = {
43              { 12, -51, 4, },
44              { 6, 167, -68, },
45              { -4, 24, -41, },
46              { -5, 34, 7, }, };
47  
48      private static final double entryTolerance = 10e-16;
49  
50      private static final double normTolerance = 10e-14;
51  
52      /** test dimensions */
53      @Test
54      public void testDimensions() {
55          checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
56  
57          checkDimension(MatrixUtils.createRealMatrix(testData4x3));
58  
59          checkDimension(MatrixUtils.createRealMatrix(testData3x4));
60  
61          Random r = new Random(643895747384642L);
62          int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
63          int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
64          checkDimension(createTestMatrix(r, p, q));
65          checkDimension(createTestMatrix(r, q, p));
66      }
67  
68      private void checkDimension(RealMatrix m) {
69          int rows = m.getRowDimension();
70          int columns = m.getColumnDimension();
71          QRDecomposition qr = new QRDecomposition(m);
72          Assert.assertEquals(rows,    qr.getQ().getRowDimension());
73          Assert.assertEquals(rows,    qr.getQ().getColumnDimension());
74          Assert.assertEquals(rows,    qr.getR().getRowDimension());
75          Assert.assertEquals(columns, qr.getR().getColumnDimension());
76      }
77  
78      /** test A = QR */
79      @Test
80      public void testAEqualQR() {
81          checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
82  
83          checkAEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
84  
85          checkAEqualQR(MatrixUtils.createRealMatrix(testData3x4));
86  
87          checkAEqualQR(MatrixUtils.createRealMatrix(testData4x3));
88  
89          Random r = new Random(643895747384642L);
90          int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
91          int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
92          checkAEqualQR(createTestMatrix(r, p, q));
93  
94          checkAEqualQR(createTestMatrix(r, q, p));
95      }
96  
97      private void checkAEqualQR(RealMatrix m) {
98          QRDecomposition qr = new QRDecomposition(m);
99          double norm = qr.getQ().multiply(qr.getR()).subtract(m).getNorm();
100         Assert.assertEquals(0, norm, normTolerance);
101     }
102 
103     /** test the orthogonality of Q */
104     @Test
105     public void testQOrthogonal() {
106         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3NonSingular));
107 
108         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x3Singular));
109 
110         checkQOrthogonal(MatrixUtils.createRealMatrix(testData3x4));
111 
112         checkQOrthogonal(MatrixUtils.createRealMatrix(testData4x3));
113 
114         Random r = new Random(643895747384642L);
115         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
116         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
117         checkQOrthogonal(createTestMatrix(r, p, q));
118 
119         checkQOrthogonal(createTestMatrix(r, q, p));
120     }
121 
122     private void checkQOrthogonal(RealMatrix m) {
123         QRDecomposition qr = new QRDecomposition(m);
124         RealMatrix eye = MatrixUtils.createRealIdentityMatrix(m.getRowDimension());
125         double norm = qr.getQT().multiply(qr.getQ()).subtract(eye).getNorm();
126         Assert.assertEquals(0, norm, normTolerance);
127     }
128 
129     /** test that R is upper triangular */
130     @Test
131     public void testRUpperTriangular() {
132         RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
133         checkUpperTriangular(new QRDecomposition(matrix).getR());
134 
135         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
136         checkUpperTriangular(new QRDecomposition(matrix).getR());
137 
138         matrix = MatrixUtils.createRealMatrix(testData3x4);
139         checkUpperTriangular(new QRDecomposition(matrix).getR());
140 
141         matrix = MatrixUtils.createRealMatrix(testData4x3);
142         checkUpperTriangular(new QRDecomposition(matrix).getR());
143 
144         Random r = new Random(643895747384642L);
145         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
146         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
147         matrix = createTestMatrix(r, p, q);
148         checkUpperTriangular(new QRDecomposition(matrix).getR());
149 
150         matrix = createTestMatrix(r, p, q);
151         checkUpperTriangular(new QRDecomposition(matrix).getR());
152     }
153 
154     private void checkUpperTriangular(RealMatrix m) {
155         m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
156             @Override
157             public void visit(int row, int column, double value) {
158                 if (column < row) {
159                     Assert.assertEquals(0.0, value, entryTolerance);
160                 }
161             }
162         });
163     }
164 
165     /** test that H is trapezoidal */
166     @Test
167     public void testHTrapezoidal() {
168         RealMatrix matrix = MatrixUtils.createRealMatrix(testData3x3NonSingular);
169         checkTrapezoidal(new QRDecomposition(matrix).getH());
170 
171         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
172         checkTrapezoidal(new QRDecomposition(matrix).getH());
173 
174         matrix = MatrixUtils.createRealMatrix(testData3x4);
175         checkTrapezoidal(new QRDecomposition(matrix).getH());
176 
177         matrix = MatrixUtils.createRealMatrix(testData4x3);
178         checkTrapezoidal(new QRDecomposition(matrix).getH());
179 
180         Random r = new Random(643895747384642L);
181         int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
182         int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
183         matrix = createTestMatrix(r, p, q);
184         checkTrapezoidal(new QRDecomposition(matrix).getH());
185 
186         matrix = createTestMatrix(r, p, q);
187         checkTrapezoidal(new QRDecomposition(matrix).getH());
188     }
189 
190     private void checkTrapezoidal(RealMatrix m) {
191         m.walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
192             @Override
193             public void visit(int row, int column, double value) {
194                 if (column > row) {
195                     Assert.assertEquals(0.0, value, entryTolerance);
196                 }
197             }
198         });
199     }
200     /** test matrices values */
201     @Test
202     public void testMatricesValues() {
203         QRDecomposition qr =
204             new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3NonSingular));
205         RealMatrix qRef = MatrixUtils.createRealMatrix(new double[][] {
206                 { -12.0 / 14.0,   69.0 / 175.0,  -58.0 / 175.0 },
207                 {  -6.0 / 14.0, -158.0 / 175.0,    6.0 / 175.0 },
208                 {   4.0 / 14.0,  -30.0 / 175.0, -165.0 / 175.0 }
209         });
210         RealMatrix rRef = MatrixUtils.createRealMatrix(new double[][] {
211                 { -14.0,  -21.0, 14.0 },
212                 {   0.0, -175.0, 70.0 },
213                 {   0.0,    0.0, 35.0 }
214         });
215         RealMatrix hRef = MatrixUtils.createRealMatrix(new double[][] {
216                 { 26.0 / 14.0, 0.0, 0.0 },
217                 {  6.0 / 14.0, 648.0 / 325.0, 0.0 },
218                 { -4.0 / 14.0,  36.0 / 325.0, 2.0 }
219         });
220 
221         // check values against known references
222         RealMatrix q = qr.getQ();
223         Assert.assertEquals(0, q.subtract(qRef).getNorm(), 1.0e-13);
224         RealMatrix qT = qr.getQT();
225         Assert.assertEquals(0, qT.subtract(qRef.transpose()).getNorm(), 1.0e-13);
226         RealMatrix r = qr.getR();
227         Assert.assertEquals(0, r.subtract(rRef).getNorm(), 1.0e-13);
228         RealMatrix h = qr.getH();
229         Assert.assertEquals(0, h.subtract(hRef).getNorm(), 1.0e-13);
230 
231         // check the same cached instance is returned the second time
232         Assert.assertSame(q, qr.getQ());
233         Assert.assertSame(r, qr.getR());
234         Assert.assertSame(h, qr.getH());
235     }
236 
237     @Test(expected=SingularMatrixException.class)
238     public void testNonInvertible() {
239         QRDecomposition qr =
240             new QRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular));
241         qr.getSolver().getInverse();
242     }
243 
244     @Test
245     public void testInvertTallSkinny() {
246         RealMatrix a     = MatrixUtils.createRealMatrix(testData4x3);
247         RealMatrix pinv  = new QRDecomposition(a).getSolver().getInverse();
248         Assert.assertEquals(0, pinv.multiply(a).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm(), 1.0e-6);
249     }
250 
251     @Test
252     public void testInvertShortWide() {
253         RealMatrix a = MatrixUtils.createRealMatrix(testData3x4);
254         RealMatrix pinv  = new QRDecomposition(a).getSolver().getInverse();
255         Assert.assertEquals(0, a.multiply(pinv).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm(), 1.0e-6);
256         Assert.assertEquals(0, pinv.multiply(a).getSubMatrix(0, 2, 0, 2).subtract(MatrixUtils.createRealIdentityMatrix(3)).getNorm(), 1.0e-6);
257     }
258 
259     private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
260         RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
261         m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
262             @Override
263             public double visit(int row, int column, double value) {
264                 return 2.0 * r.nextDouble() - 1.0;
265             }
266         });
267         return m;
268     }
269 
270     @Test(expected=SingularMatrixException.class)
271     public void testQRSingular() {
272         final RealMatrix a = MatrixUtils.createRealMatrix(new double[][] {
273             { 1, 6, 4 }, { 2, 4, -1 }, { -1, 2, 5 }
274         });
275         final RealVector b = new ArrayRealVector(new double[]{ 5, 6, 1 });
276         new QRDecomposition(a, 1.0e-15).getSolver().solve(b);
277     }
278 }