3
$\begingroup$

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
  • 3
    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.2012-08-12
  • 0
    Use Lapack. $ $2012-11-26

2 Answers 2

1

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.

  • 0
    Of course, Cholesky won't work if the symmetric matrix is not positive definite...2012-08-24
  • 0
    @J.M. Yes, that's correct. Sorry. I need to add a condition that the matrix is positive definite.2012-08-24
1

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