# Fast algorithm for LU factorization

If A is a symmetric matrix, is there a fast algorithm for LU factorization? I know this algorithm for non-symmetric matrix.

    For k = 1,..,n
For i = k + 1,...,n
mult := a_{ik}/a_{kk}
a_{jk} := mult
For j = k+1,...,n
a_{ij} := a_{ij} - mult * a_{kj}
Endfor
Endfor
Endfor

-
Look up the $\mathbf L\mathbf D\mathbf L^\top$ and Cholesky decompositions. Of course, these method will only work if all the leading submatrices of your symmetric matrix are nonsingular; otherwise, (symmetric) pivoting is necessary. –  Guess who it is. Aug 12 '12 at 5:26
Use Lapack.  –  Inquest Nov 26 '12 at 16:23

See Wikipedia article, try using Cholesky decomposition.

Added Later: The matrix has to be positive definite. So this is just a partial answer to the question.

-
Of course, Cholesky won't work if the symmetric matrix is not positive definite... –  Guess who it is. Aug 24 '12 at 8:13
@J.M. Yes, that's correct. Sorry. I need to add a condition that the matrix is positive definite. –  HbCwiRoJDp Aug 24 '12 at 8:17

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

-