Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- package Jama;
- /** Cholesky Decomposition.
- <P>
- For a symmetric, positive definite matrix A, the Cholesky decomposition
- is an lower triangular matrix L so that A = L*L'.
- <P>
- If the matrix is not symmetric or positive definite, the constructor
- returns a partial decomposition and sets an internal flag that may
- be queried by the isSPD() method.
- */
- public class CholeskyDecomposition implements java.io.Serializable {
- /* ------------------------
- Class variables
- * ------------------------ */
- /** Array for internal storage of decomposition.
- @serial internal array storage.
- */
- private double[][] L;
- /** Row and column dimension (square matrix).
- @serial matrix dimension.
- */
- private int n;
- /** Symmetric and positive definite flag.
- @serial is symmetric and positive definite flag.
- */
- private boolean isspd;
- /* ------------------------
- Constructor
- * ------------------------ */
- /** Cholesky algorithm for symmetric and positive definite matrix.
- Structure to access L and isspd flag.
- @param Arg Square, symmetric matrix.
- */
- public CholeskyDecomposition (Matrix Arg) {
- // Initialize.
- double[][] A = Arg.getArray();
- n = Arg.getRowDimension();
- L = new double[n][n];
- isspd = (Arg.getColumnDimension() == n);
- // Main loop.
- for (int j = 0; j < n; j++) {
- double[] Lrowj = L[j];
- double d = 0.0;
- for (int k = 0; k < j; k++) {
- double[] Lrowk = L[k];
- double s = 0.0;
- for (int i = 0; i < k; i++) {
- s += Lrowk[i]*Lrowj[i];
- }
- Lrowj[k] = s = (A[j][k] - s)/L[k][k];
- d = d + s*s;
- isspd = isspd & (A[k][j] == A[j][k]);
- }
- d = A[j][j] - d;
- isspd = isspd & (d > 0.0);
- L[j][j] = Math.sqrt(Math.max(d,0.0));
- for (int k = j+1; k < n; k++) {
- L[j][k] = 0.0;
- }
- }
- }
- /* ------------------------
- Temporary, experimental code.
- * ------------------------ *\
- \** Right Triangular Cholesky Decomposition.
- <P>
- For a symmetric, positive definite matrix A, the Right Cholesky
- decomposition is an upper triangular matrix R so that A = R'*R.
- This constructor computes R with the Fortran inspired column oriented
- algorithm used in LINPACK and MATLAB. In Java, we suspect a row oriented,
- lower triangular decomposition is faster. We have temporarily included
- this constructor here until timing experiments confirm this suspicion.
- *\
- \** Array for internal storage of right triangular decomposition. **\
- private transient double[][] R;
- \** Cholesky algorithm for symmetric and positive definite matrix.
- @param A Square, symmetric matrix.
- @param rightflag Actual value ignored.
- @return Structure to access R and isspd flag.
- *\
- public CholeskyDecomposition (Matrix Arg, int rightflag) {
- // Initialize.
- double[][] A = Arg.getArray();
- n = Arg.getColumnDimension();
- R = new double[n][n];
- isspd = (Arg.getColumnDimension() == n);
- // Main loop.
- for (int j = 0; j < n; j++) {
- double d = 0.0;
- for (int k = 0; k < j; k++) {
- double s = A[k][j];
- for (int i = 0; i < k; i++) {
- s = s - R[i][k]*R[i][j];
- }
- R[k][j] = s = s/R[k][k];
- d = d + s*s;
- isspd = isspd & (A[k][j] == A[j][k]);
- }
- d = A[j][j] - d;
- isspd = isspd & (d > 0.0);
- R[j][j] = Math.sqrt(Math.max(d,0.0));
- for (int k = j+1; k < n; k++) {
- R[k][j] = 0.0;
- }
- }
- }
- \** Return upper triangular factor.
- @return R
- *\
- public Matrix getR () {
- return new Matrix(R,n,n);
- }
- \* ------------------------
- End of temporary code.
- * ------------------------ */
- /* ------------------------
- Public Methods
- * ------------------------ */
- /** Is the matrix symmetric and positive definite?
- @return true if A is symmetric and positive definite.
- */
- public boolean isSPD () {
- return isspd;
- }
- /** Return triangular factor.
- @return L
- */
- public Matrix getL () {
- return new Matrix(L,n,n);
- }
- /** Solve A*X = B
- @param B A Matrix with as many rows as A and any number of columns.
- @return X so that L*L'*X = B
- @exception IllegalArgumentException Matrix row dimensions must agree.
- @exception RuntimeException Matrix is not symmetric positive definite.
- */
- public Matrix solve (Matrix B) {
- if (B.getRowDimension() != n) {
- throw new IllegalArgumentException("Matrix row dimensions must agree.");
- }
- if (!isspd) {
- throw new RuntimeException("Matrix is not symmetric positive definite.");
- }
- // Copy right hand side.
- double[][] X = B.getArrayCopy();
- int nx = B.getColumnDimension();
- // Solve L*Y = B;
- for (int k = 0; k < n; k++) {
- for (int j = 0; j < nx; j++) {
- for (int i = 0; i < k ; i++) {
- X[k][j] -= X[i][j]*L[k][i];
- }
- X[k][j] /= L[k][k];
- }
- }
- // Solve L'*X = Y;
- for (int k = n-1; k >= 0; k--) {
- for (int j = 0; j < nx; j++) {
- for (int i = k+1; i < n ; i++) {
- X[k][j] -= X[i][j]*L[i][k];
- }
- X[k][j] /= L[k][k];
- }
- }
- return new Matrix(X,n,nx);
- }
- private static final long serialVersionUID = 1;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement