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    }