+++ /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