package agents.Jama; /** Cholesky Decomposition.

For a symmetric, positive definite matrix A, the Cholesky decomposition is an lower triangular matrix L so that A = L*L'.

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.

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; }