--- /dev/null
+package tim.prune.function.estimate.jama;\r
+\r
+/**\r
+ * QR Decomposition.\r
+ *\r
+ * For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n\r
+ * orthogonal matrix Q and an n-by-n upper triangular matrix R so that A = Q*R.\r
+ *\r
+ * The QR decomposition always exists, even if the matrix does not have full\r
+ * rank, so the constructor will never fail. The primary use of the QR\r
+ * decomposition is in the least squares solution of nonsquare systems of\r
+ * simultaneous linear equations. This will fail if isFullRank() returns false.\r
+ *\r
+ * Original authors The MathWorks, Inc. and the National Institute of Standards and Technology\r
+ * The original public domain code has now been modified and reduced,\r
+ * and is placed under GPL2 with the rest of the GpsPrune code.\r
+ */\r
+public class QRDecomposition\r
+{\r
+\r
+ /** Array for internal storage of decomposition */\r
+ private double[][] _QR;\r
+\r
+ /** Row and column dimensions */\r
+ private int _m, _n;\r
+\r
+ /** Array for internal storage of diagonal of R */\r
+ private double[] _Rdiag;\r
+\r
+\r
+ /**\r
+ * QR Decomposition, computed by Householder reflections.\r
+ *\r
+ * @param inA Rectangular matrix\r
+ * @return Structure to access R and the Householder vectors and compute Q.\r
+ */\r
+ public QRDecomposition(Matrix inA)\r
+ {\r
+ // Initialize.\r
+ _QR = inA.getArrayCopy();\r
+ _m = inA.getNumRows();\r
+ _n = inA.getNumColumns();\r
+ _Rdiag = new double[_n];\r
+\r
+ // Main loop.\r
+ for (int k = 0; k < _n; k++)\r
+ {\r
+ // Compute 2-norm of k-th column without under/overflow.\r
+ double nrm = 0;\r
+ for (int i = k; i < _m; i++) {\r
+ nrm = Maths.pythag(nrm, _QR[i][k]);\r
+ }\r
+\r
+ if (nrm != 0.0)\r
+ {\r
+ // Form k-th Householder vector.\r
+ if (_QR[k][k] < 0) {\r
+ nrm = -nrm;\r
+ }\r
+ for (int i = k; i < _m; i++) {\r
+ _QR[i][k] /= nrm;\r
+ }\r
+ _QR[k][k] += 1.0;\r
+\r
+ // Apply transformation to remaining columns.\r
+ for (int j = k + 1; j < _n; j++)\r
+ {\r
+ double s = 0.0;\r
+ for (int i = k; i < _m; i++) {\r
+ s += _QR[i][k] * _QR[i][j];\r
+ }\r
+ s = -s / _QR[k][k];\r
+ for (int i = k; i < _m; i++) {\r
+ _QR[i][j] += s * _QR[i][k];\r
+ }\r
+ }\r
+ }\r
+ _Rdiag[k] = -nrm;\r
+ }\r
+ }\r
+\r
+ /*\r
+ * ------------------------ Public Methods ------------------------\r
+ */\r
+\r
+ /**\r
+ * Is the matrix full rank?\r
+ * @return true if R, and hence A, has full rank.\r
+ */\r
+ public boolean isFullRank()\r
+ {\r
+ for (int j = 0; j < _n; j++) {\r
+ if (_Rdiag[j] == 0)\r
+ return false;\r
+ }\r
+ return true;\r
+ }\r
+\r
+\r
+ /**\r
+ * Least squares solution of A*X = B\r
+ * @param B A Matrix with as many rows as A and any number of columns\r
+ * @return X that minimizes the two norm of Q*R*X-B\r
+ * @exception IllegalArgumentException if matrix dimensions don't agree\r
+ * @exception RuntimeException if Matrix is rank deficient.\r
+ */\r
+ public Matrix solve(Matrix B)\r
+ {\r
+ if (B.getNumRows() != _m) {\r
+ throw new IllegalArgumentException("Matrix row dimensions must agree.");\r
+ }\r
+ if (!isFullRank()) {\r
+ throw new RuntimeException("Matrix is rank deficient.");\r
+ }\r
+\r
+ // Copy right hand side\r
+ int nx = B.getNumColumns();\r
+ double[][] X = B.getArrayCopy();\r
+\r
+ // Compute Y = transpose(Q)*B\r
+ for (int k = 0; k < _n; k++) {\r
+ for (int j = 0; j < nx; j++) {\r
+ double s = 0.0;\r
+ for (int i = k; i < _m; i++) {\r
+ s += _QR[i][k] * X[i][j];\r
+ }\r
+ s = -s / _QR[k][k];\r
+ for (int i = k; i < _m; i++) {\r
+ X[i][j] += s * _QR[i][k];\r
+ }\r
+ }\r
+ }\r
+ // Solve R*X = Y;\r
+ for (int k = _n - 1; k >= 0; k--) {\r
+ for (int j = 0; j < nx; j++) {\r
+ X[k][j] /= _Rdiag[k];\r
+ }\r
+ for (int i = 0; i < k; i++) {\r
+ for (int j = 0; j < nx; j++) {\r
+ X[i][j] -= X[k][j] * _QR[i][k];\r
+ }\r
+ }\r
+ }\r
+ return (new Matrix(X, _n, nx).getMatrix(0, _n - 1, 0, nx - 1));\r
+ }\r
+}\r
projects / GpsPrune.git / blobdiff