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 org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
23  import org.apache.commons.math.linear.DefaultRealMatrixPreservingVisitor;
24  import org.apache.commons.math.linear.BlockRealMatrix;
25  import org.apache.commons.math.linear.MatrixUtils;
26  import org.apache.commons.math.linear.MatrixVisitorException;
27  import org.apache.commons.math.linear.QRDecomposition;
28  import org.apache.commons.math.linear.QRDecompositionImpl;
29  import org.apache.commons.math.linear.RealMatrix;
30  
31  import junit.framework.Test;
32  import junit.framework.TestCase;
33  import junit.framework.TestSuite;
34  
35  public class QRDecompositionImplTest extends TestCase {
36      double[][] testData3x3NonSingular = { 
37              { 12, -51, 4 }, 
38              { 6, 167, -68 },
39              { -4, 24, -41 }, };
40  
41      double[][] testData3x3Singular = { 
42              { 1, 4, 7, }, 
43              { 2, 5, 8, },
44              { 3, 6, 9, }, };
45  
46      double[][] testData3x4 = { 
47              { 12, -51, 4, 1 }, 
48              { 6, 167, -68, 2 },
49              { -4, 24, -41, 3 }, };
50  
51      double[][] testData4x3 = { 
52              { 12, -51, 4, }, 
53              { 6, 167, -68, },
54              { -4, 24, -41, }, 
55              { -5, 34, 7, }, };
56  
57      private static final double entryTolerance = 10e-16;
58  
59      private static final double normTolerance = 10e-14;
60  
61      public QRDecompositionImplTest(String name) {
62          super(name);
63      }
64  
65      public static Test suite() {
66          TestSuite suite = new TestSuite(QRDecompositionImplTest.class);
67          suite.setName("QRDecompositionImpl Tests");
68          return suite;
69      }
70  
71      /** test dimensions */
72      public void testDimensions() {
73          checkDimension(MatrixUtils.createRealMatrix(testData3x3NonSingular));
74  
75          checkDimension(MatrixUtils.createRealMatrix(testData4x3));
76  
77          checkDimension(MatrixUtils.createRealMatrix(testData3x4));
78  
79          Random r = new Random(643895747384642l);
80          int    p = (5 * BlockRealMatrix.BLOCK_SIZE) / 4;
81          int    q = (7 * BlockRealMatrix.BLOCK_SIZE) / 4;
82          checkDimension(createTestMatrix(r, p, q));
83          checkDimension(createTestMatrix(r, q, p));
84  
85      }
86  
87      private void checkDimension(RealMatrix m) {
88          int rows = m.getRowDimension();
89          int columns = m.getColumnDimension();
90          QRDecomposition qr = new QRDecompositionImpl(m);
91          assertEquals(rows,    qr.getQ().getRowDimension());
92          assertEquals(rows,    qr.getQ().getColumnDimension());
93          assertEquals(rows,    qr.getR().getRowDimension());
94          assertEquals(columns, qr.getR().getColumnDimension());        
95      }
96  
97      /** test A = QR */
98      public void testAEqualQR() {
99          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 }