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