tim/prune/function/estimate/jama/QRDecomposition.java

   1 package tim.prune.function.estimate.jama;\r
   2 \r
   3 /**\r
   4  * QR Decomposition.\r
   5  *\r
   6  * For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n\r
   7  * orthogonal matrix Q and an n-by-n upper triangular matrix R so that A = Q*R.\r
   8  *\r
   9  * The QR decomposition always exists, even if the matrix does not have full\r
  10  * rank, so the constructor will never fail. The primary use of the QR\r
  11  * decomposition is in the least squares solution of nonsquare systems of\r
  12  * simultaneous linear equations. This will fail if isFullRank() returns false.\r
  13  *\r
  14  * Original authors The MathWorks, Inc. and the National Institute of Standards and Technology\r
  15  * The original public domain code has now been modified and reduced,\r
  16  * and is placed under GPL2 with the rest of the GpsPrune code.\r
  17  */\r
  18 public class QRDecomposition\r
  19 {\r
  20 \r
  21         /** Array for internal storage of decomposition */\r
  22         private double[][] _QR;\r
  23 \r
  24         /** Row and column dimensions */\r
  25         private int _m, _n;\r
  26 \r
  27         /** Array for internal storage of diagonal of R */\r
  28         private double[] _Rdiag;\r
  29 \r
  30 \r
  31         /**\r
  32          * QR Decomposition, computed by Householder reflections.\r
  33          *\r
  34          * @param inA   Rectangular matrix\r
  35          * @return Structure to access R and the Householder vectors and compute Q.\r
  36          */\r
  37         public QRDecomposition(Matrix inA)\r
  38         {\r
  39                 // Initialize.\r
  40                 _QR = inA.getArrayCopy();\r
  41                 _m = inA.getNumRows();\r
  42                 _n = inA.getNumColumns();\r
  43                 _Rdiag = new double[_n];\r
  44 \r
  45                 // Main loop.\r
  46                 for (int k = 0; k < _n; k++)\r
  47                 {\r
  48                         // Compute 2-norm of k-th column without under/overflow.\r
  49                         double nrm = 0;\r
  50                         for (int i = k; i < _m; i++) {\r
  51                                 nrm = Maths.pythag(nrm, _QR[i][k]);\r
  52                         }\r
  53 \r
  54                         if (nrm != 0.0)\r
  55                         {\r
  56                                 // Form k-th Householder vector.\r
  57                                 if (_QR[k][k] < 0) {\r
  58                                         nrm = -nrm;\r
  59                                 }\r
  60                                 for (int i = k; i < _m; i++) {\r
  61                                         _QR[i][k] /= nrm;\r
  62                                 }\r
  63                                 _QR[k][k] += 1.0;\r
  64 \r
  65                                 // Apply transformation to remaining columns.\r
  66                                 for (int j = k + 1; j < _n; j++)\r
  67                                 {\r
  68                                         double s = 0.0;\r
  69                                         for (int i = k; i < _m; i++) {\r
  70                                                 s += _QR[i][k] * _QR[i][j];\r
  71                                         }\r
  72                                         s = -s / _QR[k][k];\r
  73                                         for (int i = k; i < _m; i++) {\r
  74                                                 _QR[i][j] += s * _QR[i][k];\r
  75                                         }\r
  76                                 }\r
  77                         }\r
  78                         _Rdiag[k] = -nrm;\r
  79                 }\r
  80         }\r
  81 \r
  82         /*\r
  83          * ------------------------ Public Methods ------------------------\r
  84          */\r
  85 \r
  86         /**\r
  87          * Is the matrix full rank?\r
  88          * @return true if R, and hence A, has full rank.\r
  89          */\r
  90         public boolean isFullRank()\r
  91         {\r
  92                 for (int j = 0; j < _n; j++) {\r
  93                         if (_Rdiag[j] == 0)\r
  94                                 return false;\r
  95                 }\r
  96                 return true;\r
  97         }\r
  98 \r
  99 \r
 100         /**\r
 101          * Least squares solution of A*X = B\r
 102          * @param B   A Matrix with as many rows as A and any number of columns\r
 103          * @return X that minimizes the two norm of Q*R*X-B\r
 104          * @exception IllegalArgumentException if matrix dimensions don't agree\r
 105          * @exception RuntimeException         if Matrix is rank deficient.\r
 106          */\r
 107         public Matrix solve(Matrix B)\r
 108         {\r
 109                 if (B.getNumRows() != _m) {\r
 110                         throw new IllegalArgumentException("Matrix row dimensions must agree.");\r
 111                 }\r
 112                 if (!isFullRank()) {\r
 113                         throw new RuntimeException("Matrix is rank deficient.");\r
 114                 }\r
 115 \r
 116                 // Copy right hand side\r
 117                 int nx = B.getNumColumns();\r
 118                 double[][] X = B.getArrayCopy();\r
 119 \r
 120                 // Compute Y = transpose(Q)*B\r
 121                 for (int k = 0; k < _n; k++) {\r
 122                         for (int j = 0; j < nx; j++) {\r
 123                                 double s = 0.0;\r
 124                                 for (int i = k; i < _m; i++) {\r
 125                                         s += _QR[i][k] * X[i][j];\r
 126                                 }\r
 127                                 s = -s / _QR[k][k];\r
 128                                 for (int i = k; i < _m; i++) {\r
 129                                         X[i][j] += s * _QR[i][k];\r
 130                                 }\r
 131                         }\r
 132                 }\r
 133                 // Solve R*X = Y;\r
 134                 for (int k = _n - 1; k >= 0; k--) {\r
 135                         for (int j = 0; j < nx; j++) {\r
 136                                 X[k][j] /= _Rdiag[k];\r
 137                         }\r
 138                         for (int i = 0; i < k; i++) {\r
 139                                 for (int j = 0; j < nx; j++) {\r
 140                                         X[i][j] -= X[k][j] * _QR[i][k];\r
 141                                 }\r
 142                         }\r
 143                 }\r
 144                 return (new Matrix(X, _n, nx).getMatrix(0, _n - 1, 0, nx - 1));\r
 145         }\r
 146 }\r