2.1.6 Sistemas Lineares — Métodos Diretos — Decomposição LDL

1. A Fatoração LDL

A Fatoração de Cholesky expressa uma matriz simétrica \(A \) como o produto de uma matriz triangular \(L \) (chamada de Fator de Cholesky) pela sua transposta \(L^T \) , na forma:
(eq1)

\( A = LL^T \)

onde \(L \) será triangular superior.

Uma vez que se identifica que uma matriz \(A \) é hermitiana, podemos utilizar a Equação 1 para gerar \(L \) e \(L^T \) a partir de \(A \) . Ou seja, ao multiplicarmos e igualarmos seus termos, obtemos as seguintes relações:
(eq2)


\( l_{ki} = \frac{ a_{ki} – \sum_{j=1}^{i-1} l_{ij}~l_{kj} }{ l_{ii} } \)

onde \(i = 1, 2, \ldots, k-1 \) . Além disso, o elemento da diagonal será
obtido por
(eq3)

\( l_{kk} = \sqrt{a_{kk} – \sum_{k-1}^{j=1} l_{kj}^{2} } \)

Todavia, podemos evitar o cálculo da raiz quadrada na última expressão, caso encontremos uma decomposição \(A = L D L^T \) tal que


\( A = \left[
\begin{array}{cc}
\alpha & v^T \\
v & C
\end{array}
\right] =
\left[
\begin{array}{cc}
1 & 0 \\
\frac{v}{\alpha} & I
\end{array}
\right]
\left[
\begin{array}{cc}
\alpha & 0 \\
0 & C – \frac{v v^T}{\alpha}
\end{array}
\right]
\left[
\begin{array}{cc}
1 & \frac{v^T}{\alpha} \\
0 & I
\end{array}
\right] \)

Por exemplo, a expressão acima para um problema \(3 \times 3 \)

\(\tiny A = \left[
\begin{array}{ccc}
1 & 0 & 0 \\
L_{21} & 1 & 0 \\
L_{31} & L_{32} & 1
\end{array}
\right]
\left[
\begin{array}{ccc}
D_1 & 0 & 0 \\
0 & D_2 & 0 \\
0 & 0 & D_3
\end{array}
\right]
\left[
\begin{array}{ccc}
1 & L_{21} & L_{31}\\
0 & 1 & L_{32} \\
0 & 0 & 1
\end{array}
\right] =
\left[
\begin{array}{ccc}
D_1 & ~ & ~ \\
L_{21} D_1 & L_{21}^{2}D_1 + D_2 & ~ \\
L_{31} D_1 & L_{31} L_{21} D_1 + L_{32} D_2 & L_{31}^2 D_1 + L_{32}^{2} D_2 + D_3
\end{array}
\right] \)

Portanto, podemos fazer a decomposição de Cholesky sem a necessidade da raiz quadrada através das relações


\( l_{ij} = \dfrac{\left( a_{ij} – \sum_{k=1}^{j-1} l_{ik} l_{jk} d_k \right)}{d_j} \)

e

\( d_i = a_{ii} – \sum_{k=1}^{i-1} l_{ik}^2 d_k \)

A partir destas funções, teremos as matrizes decompostas \(L = L \) e \(U = L^T \), que podem ser resolvidas por substituição regressiva os demais passos utilizados para resolução via decomposições LU.

 

2. Implementação

def LDL(A, b):
    """
    Return a list with the solution of linear system.
 
    LDL(A, b)
 
    INPUT:
    * A: a matrix of coefficients
    * b: a list, relative a vector of independent terms
 
    return: a list with unknown vector 'x' of system A*x = b.
 
    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br
 
    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Sep 2010
    """
 
    def transpose(M):
        """ Computes the transpose matrix of M.
 
        Return: M^T.
        """
        m = xrange(len(M))
        n = xrange(len(M[0]))
        p = product(m, n)
        MT = [[0.0 for i in n] for j in m]
        for i in p:
            MT[i[0]][i[1]] = M[i[1]][i[0]]
        return MT
 
    def matmul(A, B):
        """ Perform a matrix multiplication for two square matrices.
        """
        m = xrange(len(A))
        if type(B[0]) == list:
            C = [[0.0 for i in m] for j in m]
            p = product(m, m, m)
            for i in p:
                C[i[0]][i[1]] = C[i[0]][i[1]] + A[i[0]][i[2]] * B[i[2]][i[1]]
        else:
            C = [0.0 for i in m]
            p = product(m, m)
            for i in p:
                C[i[0]] = C[i[0]] + A[i[0]][i[1]] * B[i[1]]
        return C
 
    def isSimetric(M):
        """ As Cholesky factorization explore the symmetric proprierty of
        A, this function check it.
 
        Return: True if A is symmetric. False if not.
        """
        n = xrange(len(M))
        for i in n:
            for j in n:
                if M[i][j] != M[j][i]:
                    return False
        return True
 
    def LDL_decomposition(A):
        """ Perform the LDL decomposition on matriz A, returning L
        (Cholesky Factor).
        """
        n = xrange(len(A))
        M = deepcopy(A)
        L = [[0.0 for i in n] for j in n]
        Lt = [[0.0 for i in n] for j in n]
        D = [1.0 for i in n]
 
        # Calcule L and Lt
        for i in n:
            for j in n:
                suma = sum([M[i][k] * M[j][k] * D[k] for k in xrange(j)])
                M[i][j] = (A[i][j] - suma) / D[j]
            # Update diagonal
            D[i] = A[i][i] - sum([M[i][k] * M[i][k] * D[k] for k in xrange(i)])
 
        # L
        for k in n:
            L[k][0:k] = M[k][0:k]
            L[k][k] = 1.0
 
        # Lt is transpose of L
        Lt = transpose(L)
 
        return (L, D, Lt)
 
    def solveLDL(A, b):
        n = len(b)
        rn = xrange(n)
        x = [0.0 for i in rn]
        y = [0.0 for i in rn]
        (L, D, Lt) = LDL_decomposition(A)
 
       # Step 1: Solve the L system: L * y = b
        y[0] = b[0] / L[0][0]
        for i in xrange(1, n):
            suma = 0.0
            for j in xrange(i):
                suma += L[i][j] * y[j]
            y[i] = (b[i] - suma) / L[i][i]
 
        # Step 2: Solve the U system: U * x = y (regressive replacement)
        x[n - 1] = y[n - 1] / Lt[n - 1][n - 1]
        for i in xrange(n - 1, -1, -1):
            suma = y[i]
            for j in xrange(i + 1, n):
                suma = suma - Lt[i][j] * x[j]
            x[i] = suma / Lt[i][i]
 
        return x
 
    # Function's body
    if isSimetric(A):
        return solveLDL(A, b)
    else:
        return False

Estas funções estão disponíveis em http://www.sawp.com.br/code/linear/LDL.py, assim como um exemplo de como utilizá-las.

 

3. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

 

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001), 423-424.

[3] http://en.wikipedia.org/wiki/Cholesky_decomposition#Avoiding_taking_square_roots

[2] Lars Elden, Numerical Linear Algebra and Applications in Data Mining (Preliminary version), http://www.mai.liu.se/~ laeld/, (2005).