The cholesky-decomposition can be made "robust" - just keep track of zeros and negative signs in the diagonal. Here is a code-snippet in my (proprietary,sorry) MatMate-language, which should be easily translatable into a C- or Basic - or something-else routine:
TMP = X // X is the symmetric matrix L = Null(X) // Null-matrix of same size as X; // shall become the cholesky-factor sg = L[*,1] // vector which keeps track of the signs in the diagonal c = 0 // index for the current column/row // repeat the following up to the number of rows/colums of X c=c+1 d = w(TMP[c,c]) // the value of the diagonal element TMP[c,c] // if zero or lower/equal required precision finish procedure s = sign(d) sg[c] = s // store the sign of the diagonal element L[*,c] = TMP[*,c]/sqrt(s*d) // fill the c'th column TMP = TMP - s * L[*,c] * L[*,c]' // reduce tmp; // the apostroph means transposition
Check result:
chk = L * mkdiag(sg) * L' // check: re-compose cholesky-factors // with signed-diagonal err = X - chk // should be zero
We can even include, that the algorithm terminates when d becomes zero: the cholesky-factor L reflects then exactly the rank of X