parameters: (M)
return_statement: return S
docstring: A numerically stable version of the Cholesky decomposition.
  
     Used in the GLE implementation. Since many of the matrices used in this
     algorithm have very large and very small numbers in at once, to handle a
     wide range of frequencies, a naive algorithm can end up having to calculate
     the square root of a negative number, which breaks the algorithm. This is
     due to numerical precision errors turning a very tiny positive eigenvalue
     into a tiny negative value.
  
     Instead of this, an LDU decomposition is used, and any small negative numbers
     in the diagonal D matrix are assumed to be due to numerical precision errors,
     and so are replaced with zero.
  
     Args:
        M: The matrix to be decomposed.
docstring_summary: A numerically stable version of the Cholesky decomposition.
docstring_tokens: A numerically stable version of the Cholesky decomposition .
function: def stab_cholesky(M):
     """ A numerically stable version of the Cholesky decomposition.
  
     Used in the GLE implementation. Since many of the matrices used in this
     algorithm have very large and very small numbers in at once, to handle a
     wide range of frequencies, a naive algorithm can end up having to calculate
     the square root of a negative number, which breaks the algorithm. This is
     due to numerical precision errors turning a very tiny positive eigenvalue
     into a tiny negative value.
  
     Instead of this, an LDU decomposition is used, and any small negative numbers
     in the diagonal D matrix are assumed to be due to numerical precision errors,
     and so are replaced with zero.
  
     Args:
        M: The matrix to be decomposed.
     """
  
     n = M.shape[1]
     D = np.zeros(n,float)
     L = np.zeros(M.shape,float)
     for i in range(n):
        L[i,i] = 1.
        for j in range(i):
           L[i,j] = M[i,j]
           for k in range(j):
              L[i,j] -= L[i,k]*L[j,k]*D[k]
           if (not D[j] == 0.0):
              L[i,j] = L[i,j]/D[j]
        D[i] = M[i,i]
        for k in range(i):
           D[i] -= L[i,k]*L[i,k]*D[k]
  
     S = np.zeros(M.shape,float)
     for i in range(n):
        if (D[i]>0):
           D[i] = math.sqrt(D[i])
        else:
           warning("Zeroing negative element in stab-cholesky decomposition: " + str(D[i]), verbosity.low)
           D[i] = 0
        for j in range(i+1):
           S[i,j] += L[i,j]*D[j]
     return S
function_tokens: def stab_cholesky ( M ) : n = M . shape [ 1 ] D = np . zeros ( n , float ) L = np . zeros ( M . shape , float ) for i in range ( n ) : L [ i , i ] = 1. for j in range ( i ) : L [ i , j ] = M [ i , j ] for k in range ( j ) : L [ i , j ] -= L [ i , k ] * L [ j , k ] * D [ k ] if ( not D [ j ] == 0.0 ) : L [ i , j ] = L [ i , j ] / D [ j ] D [ i ] = M [ i , i ] for k in range ( i ) : D [ i ] -= L [ i , k ] * L [ i , k ] * D [ k ] S = np . zeros ( M . shape , float ) for i in range ( n ) : if ( D [ i ] > 0 ) : D [ i ] = math . sqrt ( D [ i ] ) else : warning ( "Zeroing negative element in stab-cholesky decomposition: " + str ( D [ i ] ) , verbosity . low ) D [ i ] = 0 for j in range ( i + 1 ) : S [ i , j ] += L [ i , j ] * D [ j ] return S