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 RRQRDecompositionTest {
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          RRQRDecomposition qr = new RRQRDecomposition(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 AP = QR */
79      @Test
80      public void testAPEqualQR() {
81          checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3NonSingular));
82  
83          checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x3Singular));
84  
85          checkAPEqualQR(MatrixUtils.createRealMatrix(testData3x4));
86  
87          checkAPEqualQR(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          checkAPEqualQR(createTestMatrix(r, p, q));
93  
94          checkAPEqualQR(createTestMatrix(r, q, p));
95      }
96  
97      private void checkAPEqualQR(RealMatrix m) {
98          RRQRDecomposition rrqr = new RRQRDecomposition(m);
99          double norm = rrqr.getQ().multiply(rrqr.getR()).subtract(m.multiply(rrqr.getP())).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         RRQRDecomposition qr = new RRQRDecomposition(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 RRQRDecomposition(matrix).getR());
134 
135         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
136         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
137 
138         matrix = MatrixUtils.createRealMatrix(testData3x4);
139         checkUpperTriangular(new RRQRDecomposition(matrix).getR());
140 
141         matrix = MatrixUtils.createRealMatrix(testData4x3);
142         checkUpperTriangular(new RRQRDecomposition(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 RRQRDecomposition(matrix).getR());
149 
150         matrix = createTestMatrix(r, p, q);
151         checkUpperTriangular(new RRQRDecomposition(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 RRQRDecomposition(matrix).getH());
170 
171         matrix = MatrixUtils.createRealMatrix(testData3x3Singular);
172         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
173 
174         matrix = MatrixUtils.createRealMatrix(testData3x4);
175         checkTrapezoidal(new RRQRDecomposition(matrix).getH());
176 
177         matrix = MatrixUtils.createRealMatrix(testData4x3);
178         checkTrapezoidal(new RRQRDecomposition(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 RRQRDecomposition(matrix).getH());
185 
186         matrix = createTestMatrix(r, p, q);
187         checkTrapezoidal(new RRQRDecomposition(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 
201     @Test(expected=SingularMatrixException.class)
202     public void testNonInvertible() {
203         RRQRDecomposition qr =
204             new RRQRDecomposition(MatrixUtils.createRealMatrix(testData3x3Singular), 3.0e-16);
205         qr.getSolver().getInverse();
206     }
207 
208     private RealMatrix createTestMatrix(final Random r, final int rows, final int columns) {
209         RealMatrix m = MatrixUtils.createRealMatrix(rows, columns);
210         m.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor(){
211             @Override
212             public double visit(int row, int column, double value) {
213                 return 2.0 * r.nextDouble() - 1.0;
214             }
215         });
216         return m;
217     }
218 
219     /** test the rank is returned correctly */
220     @Test
221     public void testRank() {
222         double[][] d = { { 1, 1, 1 }, { 0, 0, 0 }, { 1, 2, 3 } };
223         RealMatrix m = new Array2DRowRealMatrix(d);
224         RRQRDecomposition qr = new RRQRDecomposition(m);
225         Assert.assertEquals(2, qr.getRank(0));
226     }
227     @Test
228     public void testRank2() {
229         double[][] d = { { 1, 1, 1 }, { 2, 3, 4 }, { 1, 2, 3 } };
230         RealMatrix m = new Array2DRowRealMatrix(d);
231         RRQRDecomposition qr = new RRQRDecomposition(m);
232         Assert.assertEquals(2, qr.getRank(1e-14));
233     }
234 
235     // MATH-1417
236     @Test
237     public void testRank3() {
238         double[][] d = {
239             {0, 0, 0, 0, 0, 0, 0, 0, 0},
240             {0, 1, 0, 0, 0, 0, 0, 0, 0},
241             {0, 0, 1, 0, 0, 0, 0, 0, 0},
242             {0, 0, 1, 0, 0, 0, 0, 0, 0},
243             {0, 0, 1, 0, 0, 0, 0, 0, 0},
244             {0, 0, 0, 1, 0, 0, 0, 0, 0},
245             {0, 0, 0, 0, 0, 0, 1, 0, 0},
246             {0, 0, 0, 0, 0, 0, 0, 0, 0}
247         };
248         RealMatrix m = new Array2DRowRealMatrix(d);
249         RRQRDecomposition qr = new RRQRDecomposition(m.transpose());
250         Assert.assertEquals(4, qr.getRank(1e-14));
251     }
252 }