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.math3.linear; 019 020 import java.util.Arrays; 021 022 import org.apache.commons.math3.exception.DimensionMismatchException; 023 import org.apache.commons.math3.util.FastMath; 024 025 026 /** 027 * Calculates the QR-decomposition of a matrix. 028 * <p>The QR-decomposition of a matrix A consists of two matrices Q and R 029 * that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is 030 * upper triangular. If A is m×n, Q is m×m and R m×n.</p> 031 * <p>This class compute the decomposition using Householder reflectors.</p> 032 * <p>For efficiency purposes, the decomposition in packed form is transposed. 033 * This allows inner loop to iterate inside rows, which is much more cache-efficient 034 * in Java.</p> 035 * <p>This class is based on the class with similar name from the 036 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the 037 * following changes:</p> 038 * <ul> 039 * <li>a {@link #getQT() getQT} method has been added,</li> 040 * <li>the {@code solve} and {@code isFullRank} methods have been replaced 041 * by a {@link #getSolver() getSolver} method and the equivalent methods 042 * provided by the returned {@link DecompositionSolver}.</li> 043 * </ul> 044 * 045 * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a> 046 * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a> 047 * 048 * @version $Id: QRDecomposition.java 1416643 2012-12-03 19:37:14Z tn $ 049 * @since 1.2 (changed to concrete class in 3.0) 050 */ 051 public class QRDecomposition { 052 /** 053 * A packed TRANSPOSED representation of the QR decomposition. 054 * <p>The elements BELOW the diagonal are the elements of the UPPER triangular 055 * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors 056 * from which an explicit form of Q can be recomputed if desired.</p> 057 */ 058 private double[][] qrt; 059 /** The diagonal elements of R. */ 060 private double[] rDiag; 061 /** Cached value of Q. */ 062 private RealMatrix cachedQ; 063 /** Cached value of QT. */ 064 private RealMatrix cachedQT; 065 /** Cached value of R. */ 066 private RealMatrix cachedR; 067 /** Cached value of H. */ 068 private RealMatrix cachedH; 069 /** Singularity threshold. */ 070 private final double threshold; 071 072 /** 073 * Calculates the QR-decomposition of the given matrix. 074 * The singularity threshold defaults to zero. 075 * 076 * @param matrix The matrix to decompose. 077 * 078 * @see #QRDecomposition(RealMatrix,double) 079 */ 080 public QRDecomposition(RealMatrix matrix) { 081 this(matrix, 0d); 082 } 083 084 /** 085 * Calculates the QR-decomposition of the given matrix. 086 * 087 * @param matrix The matrix to decompose. 088 * @param threshold Singularity threshold. 089 */ 090 public QRDecomposition(RealMatrix matrix, 091 double threshold) { 092 this.threshold = threshold; 093 094 final int m = matrix.getRowDimension(); 095 final int n = matrix.getColumnDimension(); 096 qrt = matrix.transpose().getData(); 097 rDiag = new double[FastMath.min(m, n)]; 098 cachedQ = null; 099 cachedQT = null; 100 cachedR = null; 101 cachedH = null; 102 103 /* 104 * The QR decomposition of a matrix A is calculated using Householder 105 * reflectors by repeating the following operations to each minor 106 * A(minor,minor) of A: 107 */ 108 for (int minor = 0; minor < FastMath.min(m, n); minor++) { 109 110 final double[] qrtMinor = qrt[minor]; 111 112 /* 113 * Let x be the first column of the minor, and a^2 = |x|^2. 114 * x will be in the positions qr[minor][minor] through qr[m][minor]. 115 * The first column of the transformed minor will be (a,0,0,..)' 116 * The sign of a is chosen to be opposite to the sign of the first 117 * component of x. Let's find a: 118 */ 119 double xNormSqr = 0; 120 for (int row = minor; row < m; row++) { 121 final double c = qrtMinor[row]; 122 xNormSqr += c * c; 123 } 124 final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); 125 rDiag[minor] = a; 126 127 if (a != 0.0) { 128 129 /* 130 * Calculate the normalized reflection vector v and transform 131 * the first column. We know the norm of v beforehand: v = x-ae 132 * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> = 133 * a^2+a^2-2a<x,e> = 2a*(a - <x,e>). 134 * Here <x, e> is now qr[minor][minor]. 135 * v = x-ae is stored in the column at qr: 136 */ 137 qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor]) 138 139 /* 140 * Transform the rest of the columns of the minor: 141 * They will be transformed by the matrix H = I-2vv'/|v|^2. 142 * If x is a column vector of the minor, then 143 * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v. 144 * Therefore the transformation is easily calculated by 145 * subtracting the column vector (2<x,v>/|v|^2)v from x. 146 * 147 * Let 2<x,v>/|v|^2 = alpha. From above we have 148 * |v|^2 = -2a*(qr[minor][minor]), so 149 * alpha = -<x,v>/(a*qr[minor][minor]) 150 */ 151 for (int col = minor+1; col < n; col++) { 152 final double[] qrtCol = qrt[col]; 153 double alpha = 0; 154 for (int row = minor; row < m; row++) { 155 alpha -= qrtCol[row] * qrtMinor[row]; 156 } 157 alpha /= a * qrtMinor[minor]; 158 159 // Subtract the column vector alpha*v from x. 160 for (int row = minor; row < m; row++) { 161 qrtCol[row] -= alpha * qrtMinor[row]; 162 } 163 } 164 } 165 } 166 } 167 168 /** 169 * Returns the matrix R of the decomposition. 170 * <p>R is an upper-triangular matrix</p> 171 * @return the R matrix 172 */ 173 public RealMatrix getR() { 174 175 if (cachedR == null) { 176 177 // R is supposed to be m x n 178 final int n = qrt.length; 179 final int m = qrt[0].length; 180 double[][] ra = new double[m][n]; 181 // copy the diagonal from rDiag and the upper triangle of qr 182 for (int row = FastMath.min(m, n) - 1; row >= 0; row--) { 183 ra[row][row] = rDiag[row]; 184 for (int col = row + 1; col < n; col++) { 185 ra[row][col] = qrt[col][row]; 186 } 187 } 188 cachedR = MatrixUtils.createRealMatrix(ra); 189 } 190 191 // return the cached matrix 192 return cachedR; 193 } 194 195 /** 196 * Returns the matrix Q of the decomposition. 197 * <p>Q is an orthogonal matrix</p> 198 * @return the Q matrix 199 */ 200 public RealMatrix getQ() { 201 if (cachedQ == null) { 202 cachedQ = getQT().transpose(); 203 } 204 return cachedQ; 205 } 206 207 /** 208 * Returns the transpose of the matrix Q of the decomposition. 209 * <p>Q is an orthogonal matrix</p> 210 * @return the transpose of the Q matrix, Q<sup>T</sup> 211 */ 212 public RealMatrix getQT() { 213 if (cachedQT == null) { 214 215 // QT is supposed to be m x m 216 final int n = qrt.length; 217 final int m = qrt[0].length; 218 double[][] qta = new double[m][m]; 219 220 /* 221 * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then 222 * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in 223 * succession to the result 224 */ 225 for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) { 226 qta[minor][minor] = 1.0d; 227 } 228 229 for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){ 230 final double[] qrtMinor = qrt[minor]; 231 qta[minor][minor] = 1.0d; 232 if (qrtMinor[minor] != 0.0) { 233 for (int col = minor; col < m; col++) { 234 double alpha = 0; 235 for (int row = minor; row < m; row++) { 236 alpha -= qta[col][row] * qrtMinor[row]; 237 } 238 alpha /= rDiag[minor] * qrtMinor[minor]; 239 240 for (int row = minor; row < m; row++) { 241 qta[col][row] += -alpha * qrtMinor[row]; 242 } 243 } 244 } 245 } 246 cachedQT = MatrixUtils.createRealMatrix(qta); 247 } 248 249 // return the cached matrix 250 return cachedQT; 251 } 252 253 /** 254 * Returns the Householder reflector vectors. 255 * <p>H is a lower trapezoidal matrix whose columns represent 256 * each successive Householder reflector vector. This matrix is used 257 * to compute Q.</p> 258 * @return a matrix containing the Householder reflector vectors 259 */ 260 public RealMatrix getH() { 261 if (cachedH == null) { 262 263 final int n = qrt.length; 264 final int m = qrt[0].length; 265 double[][] ha = new double[m][n]; 266 for (int i = 0; i < m; ++i) { 267 for (int j = 0; j < FastMath.min(i + 1, n); ++j) { 268 ha[i][j] = qrt[j][i] / -rDiag[j]; 269 } 270 } 271 cachedH = MatrixUtils.createRealMatrix(ha); 272 } 273 274 // return the cached matrix 275 return cachedH; 276 } 277 278 /** 279 * Get a solver for finding the A × X = B solution in least square sense. 280 * @return a solver 281 */ 282 public DecompositionSolver getSolver() { 283 return new Solver(qrt, rDiag, threshold); 284 } 285 286 /** Specialized solver. */ 287 private static class Solver implements DecompositionSolver { 288 /** 289 * A packed TRANSPOSED representation of the QR decomposition. 290 * <p>The elements BELOW the diagonal are the elements of the UPPER triangular 291 * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors 292 * from which an explicit form of Q can be recomputed if desired.</p> 293 */ 294 private final double[][] qrt; 295 /** The diagonal elements of R. */ 296 private final double[] rDiag; 297 /** Singularity threshold. */ 298 private final double threshold; 299 300 /** 301 * Build a solver from decomposed matrix. 302 * 303 * @param qrt Packed TRANSPOSED representation of the QR decomposition. 304 * @param rDiag Diagonal elements of R. 305 * @param threshold Singularity threshold. 306 */ 307 private Solver(final double[][] qrt, 308 final double[] rDiag, 309 final double threshold) { 310 this.qrt = qrt; 311 this.rDiag = rDiag; 312 this.threshold = threshold; 313 } 314 315 /** {@inheritDoc} */ 316 public boolean isNonSingular() { 317 for (double diag : rDiag) { 318 if (FastMath.abs(diag) <= threshold) { 319 return false; 320 } 321 } 322 return true; 323 } 324 325 /** {@inheritDoc} */ 326 public RealVector solve(RealVector b) { 327 final int n = qrt.length; 328 final int m = qrt[0].length; 329 if (b.getDimension() != m) { 330 throw new DimensionMismatchException(b.getDimension(), m); 331 } 332 if (!isNonSingular()) { 333 throw new SingularMatrixException(); 334 } 335 336 final double[] x = new double[n]; 337 final double[] y = b.toArray(); 338 339 // apply Householder transforms to solve Q.y = b 340 for (int minor = 0; minor < FastMath.min(m, n); minor++) { 341 342 final double[] qrtMinor = qrt[minor]; 343 double dotProduct = 0; 344 for (int row = minor; row < m; row++) { 345 dotProduct += y[row] * qrtMinor[row]; 346 } 347 dotProduct /= rDiag[minor] * qrtMinor[minor]; 348 349 for (int row = minor; row < m; row++) { 350 y[row] += dotProduct * qrtMinor[row]; 351 } 352 } 353 354 // solve triangular system R.x = y 355 for (int row = rDiag.length - 1; row >= 0; --row) { 356 y[row] /= rDiag[row]; 357 final double yRow = y[row]; 358 final double[] qrtRow = qrt[row]; 359 x[row] = yRow; 360 for (int i = 0; i < row; i++) { 361 y[i] -= yRow * qrtRow[i]; 362 } 363 } 364 365 return new ArrayRealVector(x, false); 366 } 367 368 /** {@inheritDoc} */ 369 public RealMatrix solve(RealMatrix b) { 370 final int n = qrt.length; 371 final int m = qrt[0].length; 372 if (b.getRowDimension() != m) { 373 throw new DimensionMismatchException(b.getRowDimension(), m); 374 } 375 if (!isNonSingular()) { 376 throw new SingularMatrixException(); 377 } 378 379 final int columns = b.getColumnDimension(); 380 final int blockSize = BlockRealMatrix.BLOCK_SIZE; 381 final int cBlocks = (columns + blockSize - 1) / blockSize; 382 final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns); 383 final double[][] y = new double[b.getRowDimension()][blockSize]; 384 final double[] alpha = new double[blockSize]; 385 386 for (int kBlock = 0; kBlock < cBlocks; ++kBlock) { 387 final int kStart = kBlock * blockSize; 388 final int kEnd = FastMath.min(kStart + blockSize, columns); 389 final int kWidth = kEnd - kStart; 390 391 // get the right hand side vector 392 b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y); 393 394 // apply Householder transforms to solve Q.y = b 395 for (int minor = 0; minor < FastMath.min(m, n); minor++) { 396 final double[] qrtMinor = qrt[minor]; 397 final double factor = 1.0 / (rDiag[minor] * qrtMinor[minor]); 398 399 Arrays.fill(alpha, 0, kWidth, 0.0); 400 for (int row = minor; row < m; ++row) { 401 final double d = qrtMinor[row]; 402 final double[] yRow = y[row]; 403 for (int k = 0; k < kWidth; ++k) { 404 alpha[k] += d * yRow[k]; 405 } 406 } 407 for (int k = 0; k < kWidth; ++k) { 408 alpha[k] *= factor; 409 } 410 411 for (int row = minor; row < m; ++row) { 412 final double d = qrtMinor[row]; 413 final double[] yRow = y[row]; 414 for (int k = 0; k < kWidth; ++k) { 415 yRow[k] += alpha[k] * d; 416 } 417 } 418 } 419 420 // solve triangular system R.x = y 421 for (int j = rDiag.length - 1; j >= 0; --j) { 422 final int jBlock = j / blockSize; 423 final int jStart = jBlock * blockSize; 424 final double factor = 1.0 / rDiag[j]; 425 final double[] yJ = y[j]; 426 final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock]; 427 int index = (j - jStart) * kWidth; 428 for (int k = 0; k < kWidth; ++k) { 429 yJ[k] *= factor; 430 xBlock[index++] = yJ[k]; 431 } 432 433 final double[] qrtJ = qrt[j]; 434 for (int i = 0; i < j; ++i) { 435 final double rIJ = qrtJ[i]; 436 final double[] yI = y[i]; 437 for (int k = 0; k < kWidth; ++k) { 438 yI[k] -= yJ[k] * rIJ; 439 } 440 } 441 } 442 } 443 444 return new BlockRealMatrix(n, columns, xBlocks, false); 445 } 446 447 /** {@inheritDoc} */ 448 public RealMatrix getInverse() { 449 return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length)); 450 } 451 } 452 }