package util.linalg;
/**
* Calculates A = L*Lt where A is a symmetric
* positive definite matrix and L is a lower triangular matrix.
* @author Andrew Guillory gtg008g@mail.gatech.edu
* @version 1.0
*/
public class CholeskyFactorization {
/**
* The lower triangular matrix
*/
private LowerTriangularMatrix l;
/**
* The transpose of the matrix
*/
private UpperTriangularMatrix lt;
/**
* Create a cholesky factorization of the given matrix
* @param a the matrix to factor
*/
public CholeskyFactorization(Matrix a) {
l = new LowerTriangularMatrix(a);
decompose();
}
/**
* Factors the matrix l inplace.
*/
private void decompose() {
// loop through the diagonal of the matrix
for (int j = 0; j < l.n(); j++) {
// sqrt the diagonal
l.set(j,j, Math.sqrt(l.get(j,j)));
// divide the subdiagonal column by
// the diagonal
for (int i = j + 1; i < l.m(); i++) {
l.set(i,j, l.get(i,j) / l.get(j,j));
}
// symmetric rank 1 update
// subtract the crossproduct of the
// subdiagonal column from the remaining
// lower diagonal
for (int jj = j + 1; jj < l.n(); jj++) {
for (int ii = jj; ii < l.m(); ii++) {
l.set(ii,jj, l.get(ii,jj)
- l.get(ii,j) * l.get(jj,j));
}
}
}
// set the l transpose
lt = (UpperTriangularMatrix) l.transpose();
}
/**
* Get the lower triangular matrix
* @return the matrix
*/
public LowerTriangularMatrix getL() {
return l;
}
/**
* Get the transpose matrix
* @return the transpose
*/
public UpperTriangularMatrix getLt() {
return lt;
}
/**
* Calculate the determinant
* @return the determinant
*/
public double determinant() {
double d = l.determinant();
return d * d;
}
/**
* Solve the system of linear equations
* for the given vector. Find the column
* vector x such that A*x = b.
* @param b the column vector to solve for
* @return the solution
*/
public Vector solve(Vector b) {
// first solves L*y = b
Vector y = l.solve(b);
// now solves Lt*x = y
return lt.solve(y);
}
}