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 /** 21 * Calculates the QR-decomposition of a matrix. In the QR-decomposition of 22 * a matrix A consists of two matrices Q and R that satisfy: A = QR, Q is 23 * orthogonal (Q<sup>T</sup>Q = I), and R is upper triangular. If A is 24 * m×n, Q is m×m and R m×n. 25 * <p> 26 * Implemented using Householder reflectors.</p> 27 * 28 * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a> 29 * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a> 30 * 31 * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $ 32 * @since 1.2 33 */ 34 public class QRDecompositionImpl implements QRDecomposition { 35 36 /** 37 * A packed representation of the QR decomposition. The elements above the 38 * diagonal are the elements of R, and the columns of the lower triangle 39 * are the Householder reflector vectors of which an explicit form of Q can 40 * be calculated. 41 */ 42 private double[][] qr; 43 44 /** 45 * The diagonal elements of R. 46 */ 47 private double[] rDiag; 48 49 /** 50 * The row dimension of the given matrix. The size of Q will be m x m, the 51 * size of R will be m x n. 52 */ 53 private int m; 54 55 /** 56 * The column dimension of the given matrix. The size of R will be m x n. 57 */ 58 private int n; 59 60 /** 61 * Calculates the QR decomposition of the given matrix. 62 * 63 * @param matrix The matrix to decompose. 64 */ 65 public QRDecompositionImpl(RealMatrix matrix) { 66 m = matrix.getRowDimension(); 67 n = matrix.getColumnDimension(); 68 qr = matrix.getData(); 69 rDiag = new double[n]; 70 71 /* 72 * The QR decomposition of a matrix A is calculated using Householder 73 * reflectors by repeating the following operations to each minor 74 * A(minor,minor) of A: 75 */ 76 for (int minor = 0; minor < Math.min(m, n); minor++) { 77 /* 78 * Let x be the first column of the minor, and a^2 = |x|^2. 79 * x will be in the positions qr[minor][minor] through qr[m][minor]. 80 * The first column of the transformed minor will be (a,0,0,..)' 81 * The sign of a is chosen to be opposite to the sign of the first 82 * component of x. Let's find a: 83 */ 84 double xNormSqr = 0; 85 for (int row = minor; row < m; row++) { 86 xNormSqr += qr[row][minor]*qr[row][minor]; 87 } 88 double a = Math.sqrt(xNormSqr); 89 if (qr[minor][minor] > 0) a = -a; 90 rDiag[minor] = a; 91 92 if (a != 0.0) { 93 94 /* 95 * Calculate the normalized reflection vector v and transform 96 * the first column. We know the norm of v beforehand: v = x-ae 97 * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> = 98 * a^2+a^2-2a<x,e> = 2a*(a - <x,e>). 99 * Here <x, e> is now qr[minor][minor]. 100 * v = x-ae is stored in the column at qr: 101 */ 102 qr[minor][minor] -= a; // now |v|^2 = -2a*(qr[minor][minor]) 103 104 /* 105 * Transform the rest of the columns of the minor: 106 * They will be transformed by the matrix H = I-2vv'/|v|^2. 107 * If x is a column vector of the minor, then 108 * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v. 109 * Therefore the transformation is easily calculated by 110 * subtracting the column vector (2<x,v>/|v|^2)v from x. 111 * 112 * Let 2<x,v>/|v|^2 = alpha. From above we have 113 * |v|^2 = -2a*(qr[minor][minor]), so 114 * alpha = -<x,v>/(a*qr[minor][minor]) 115 */ 116 for (int col = minor+1; col < n; col++) { 117 double alpha = 0; 118 for (int row = minor; row < m; row++) { 119 alpha -= qr[row][col]*qr[row][minor]; 120 } 121 alpha /= a*qr[minor][minor]; 122 123 // Subtract the column vector alpha*v from x. 124 for (int row = minor; row < m; row++) { 125 qr[row][col] -= alpha*qr[row][minor]; 126 } 127 } 128 } 129 } 130 } 131 132 /** 133 * Returns the matrix R of the QR-decomposition. 134 * 135 * @return the R matrix 136 */ 137 public RealMatrix getR() 138 { 139 // R is supposed to be m x n 140 RealMatrixImpl ret = new RealMatrixImpl(m,n); 141 double[][] r = ret.getDataRef(); 142 143 // copy the diagonal from rDiag and the upper triangle of qr 144 for (int row = Math.min(m,n)-1; row >= 0; row--) { 145 r[row][row] = rDiag[row]; 146 for (int col = row+1; col < n; col++) { 147 r[row][col] = qr[row][col]; 148 } 149 } 150 return ret; 151 } 152 153 /** 154 * Returns the matrix Q of the QR-decomposition. 155 * 156 * @return the Q matrix 157 */ 158 public RealMatrix getQ() 159 { 160 // Q is supposed to be m x m 161 RealMatrixImpl ret = new RealMatrixImpl(m,m); 162 double[][] Q = ret.getDataRef(); 163 164 /* 165 * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then 166 * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in 167 * succession to the result 168 */ 169 for (int minor = m-1; minor >= Math.min(m,n); minor--) { 170 Q[minor][minor]=1; 171 } 172 173 for (int minor = Math.min(m,n)-1; minor >= 0; minor--){ 174 Q[minor][minor] = 1; 175 if (qr[minor][minor] != 0.0) { 176 for (int col = minor; col < m; col++) { 177 double alpha = 0; 178 for (int row = minor; row < m; row++) { 179 alpha -= Q[row][col] * qr[row][minor]; 180 } 181 alpha /= rDiag[minor]*qr[minor][minor]; 182 183 for (int row = minor; row < m; row++) { 184 Q[row][col] -= alpha*qr[row][minor]; 185 } 186 } 187 } 188 } 189 190 return ret; 191 } 192 }