Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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}
share|cite|improve this question
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. – J. M. 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.

share|cite|improve this answer
Of course, Cholesky won't work if the symmetric matrix is not positive definite... – J. M. 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
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

share|cite|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.