SimpleCholesky.java

package mikera.matrixx.decompose.impl.chol;

import mikera.matrixx.AMatrix;
import mikera.matrixx.Matrix;
import mikera.matrixx.Matrixx;
import mikera.matrixx.decompose.ICholeskyResult;
import mikera.util.Maths;

/**
 * Class implementing a standard Cholesky decomposition
 * 
 *    A = L.L*
 *    
 * Where: A  is a symmetric (actually Hermitian), positive-definite matrix 
 * and:   L  is a lower triangular matrix
 *        L* is the conjugate transpose of L (which is equal to its transpose, since A is real in Vectorz)
 * 
 * See: http://en.wikipedia.org/wiki/Cholesky_decomposition
 * 
 * @author Mike
 *
 */
public class SimpleCholesky {

	/**
	 * Decompose a matrix according the the Cholesky decomposition A = L.L*
	 * 
	 * @param a Any symmetric, positive definite matrix
	 * @return The decomposition result
	 */
	public static final ICholeskyResult decompose(AMatrix a) {
		return decompose(a.toMatrix());
	}
	
	public static final ICholeskyResult decompose(Matrix a) {
		if (!a.isSquare()) throw new IllegalArgumentException("Matrix must be square for Cholesky decomposition");
		int n=a.rowCount();
		
		Matrix u=Matrix.create(n,n);
		for (int i=0; i<n;i++) {
			double squareSum=0.0;

			for (int j=0; j<i; j++) {
				double crossSum=0.0;
				
				for (int k=0; k<j; k++) {
					crossSum+=u.get(i,k)*u.get(j,k);
				}
				
				final double aij=a.get(i,j);
				double uij=(aij-crossSum);
				
				double ujj=u.get(j,j);
				if (ujj==0.0) return null;
				uij=uij/ujj;
				u.set(i,j,uij);
				squareSum+=uij*uij;
			}	
			
			double aii =a.get(i,i);
			double uii=Maths.sqrt(aii-squareSum);
			u.set(i,i,uii);
		}
		
		// TODO: should be null return for a failed decomposition?
		
		AMatrix L = Matrixx.extractLowerTriangular(u);
		return new CholeskyResult(L);
	}
}