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          * Return the Householder vectors\r
 101          * @deprecated\r
 102          * @return Lower trapezoidal matrix whose columns define the reflections\r
 103          */\r
 104         private Matrix getH()\r
 105         {\r
 106                 Matrix X = new Matrix(_m, _n);\r
 107                 double[][] H = X.getArray();\r
 108                 for (int i = 0; i < _m; i++) {\r
 109                         for (int j = 0; j < _n; j++) {\r
 110                                 if (i >= j) {\r
 111                                         H[i][j] = _QR[i][j];\r
 112                                 } else {\r
 113                                         H[i][j] = 0.0;\r
 114                                 }\r
 115                         }\r
 116                 }\r
 117                 return X;\r
 118         }\r
 119 \r
 120         /**\r
 121          * Return the upper triangular factor\r
 122          * @deprecated\r
 123          * @return R\r
 124          */\r
 125         private Matrix getR()\r
 126         {\r
 127                 Matrix X = new Matrix(_n, _n);\r
 128                 double[][] R = X.getArray();\r
 129                 for (int i = 0; i < _n; i++) {\r
 130                         for (int j = 0; j < _n; j++) {\r
 131                                 if (i < j) {\r
 132                                         R[i][j] = _QR[i][j];\r
 133                                 } else if (i == j) {\r
 134                                         R[i][j] = _Rdiag[i];\r
 135                                 } else {\r
 136                                         R[i][j] = 0.0;\r
 137                                 }\r
 138                         }\r
 139                 }\r
 140                 return X;\r
 141         }\r
 142 \r
 143         /**\r
 144          * Generate and return the (economy-sized) orthogonal factor\r
 145          * @deprecated\r
 146          * @return Q\r
 147          */\r
 148         private Matrix getQ()\r
 149         {\r
 150                 Matrix X = new Matrix(_m, _n);\r
 151                 double[][] Q = X.getArray();\r
 152                 for (int k = _n - 1; k >= 0; k--) {\r
 153                         for (int i = 0; i < _m; i++) {\r
 154                                 Q[i][k] = 0.0;\r
 155                         }\r
 156                         Q[k][k] = 1.0;\r
 157                         for (int j = k; j < _n; j++) {\r
 158                                 if (_QR[k][k] != 0) {\r
 159                                         double s = 0.0;\r
 160                                         for (int i = k; i < _m; i++) {\r
 161                                                 s += _QR[i][k] * Q[i][j];\r
 162                                         }\r
 163                                         s = -s / _QR[k][k];\r
 164                                         for (int i = k; i < _m; i++) {\r
 165                                                 Q[i][j] += s * _QR[i][k];\r
 166                                         }\r
 167                                 }\r
 168                         }\r
 169                 }\r
 170                 return X;\r
 171         }\r
 172 \r
 173         /**\r
 174          * Least squares solution of A*X = B\r
 175          * @param B   A Matrix with as many rows as A and any number of columns\r
 176          * @return X that minimizes the two norm of Q*R*X-B\r
 177          * @exception IllegalArgumentException if matrix dimensions don't agree\r
 178          * @exception RuntimeException         if Matrix is rank deficient.\r
 179          */\r
 180         public Matrix solve(Matrix B)\r
 181         {\r
 182                 if (B.getNumRows() != _m) {\r
 183                         throw new IllegalArgumentException("Matrix row dimensions must agree.");\r
 184                 }\r
 185                 if (!isFullRank()) {\r
 186                         throw new RuntimeException("Matrix is rank deficient.");\r
 187                 }\r
 188 \r
 189                 // Copy right hand side\r
 190                 int nx = B.getNumColumns();\r
 191                 double[][] X = B.getArrayCopy();\r
 192 \r
 193                 // Compute Y = transpose(Q)*B\r
 194                 for (int k = 0; k < _n; k++) {\r
 195                         for (int j = 0; j < nx; j++) {\r
 196                                 double s = 0.0;\r
 197                                 for (int i = k; i < _m; i++) {\r
 198                                         s += _QR[i][k] * X[i][j];\r
 199                                 }\r
 200                                 s = -s / _QR[k][k];\r
 201                                 for (int i = k; i < _m; i++) {\r
 202                                         X[i][j] += s * _QR[i][k];\r
 203                                 }\r
 204                         }\r
 205                 }\r
 206                 // Solve R*X = Y;\r
 207                 for (int k = _n - 1; k >= 0; k--) {\r
 208                         for (int j = 0; j < nx; j++) {\r
 209                                 X[k][j] /= _Rdiag[k];\r
 210                         }\r
 211                         for (int i = 0; i < k; i++) {\r
 212                                 for (int j = 0; j < nx; j++) {\r
 213                                         X[i][j] -= X[k][j] * _QR[i][k];\r
 214                                 }\r
 215                         }\r
 216                 }\r
 217                 return (new Matrix(X, _n, nx).getMatrix(0, _n - 1, 0, nx - 1));\r
 218         }\r
 219 }\r