2.1.3 Sistemas Lineares — Métodos Diretos — Decomposição de Cholesky

1. A Fatoração de Cholesky

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* \) , na forma:
(eq1)


\( A = LL* \)

onde \(L \) será triangular superior.

Sabendo-se que \(A \) é simétrica quando seus elementos obedecem à uma formação \(a_{ij} = a_{ji} \) para toda linha \(i \) e coluna \(j \) . Desta forma, a matriz \(A \) será idêntica à sua transposta: \(A = A^{T} \) . Matrizes com este tipo de formação fornece algumas vantagens computacionais, pois favorecem a eficiência na solução de sistemas lineares, quando usados métodos adequados. Além disso, essa identidade permite a criação de algoritmos simples de criptografia e verificação de dados.

Uma vez que se identifica que uma matriz \(A \) é simétrica, podemos utilizar a Equação 1 para gerar \(L \) e \(L* \) 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} } \)

As expressões acima são obtidas como uma especificação da decomposição LU para o caso da matriz decomposta ser simétrica. Tais equações são desenvolvidas em A First Course in Numerical Analysis [1].

2. Implementação

def cholesky(A, b):
    """
    Return a list with the solution of linear system.
 
    cholesky(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/
 
    Apr 2010
    """
 
    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 = len(M)
        for i in xrange(n):
            for j in xrange(n):
                if M[i][j] != M[j][i]:
                    return False
        return True
 
    def cholesky_decomposition(A):
        """ Perform the Cholesky decomposition on matriz A, returning L
        (Cholesky Factor).
        """
        n = len(A)
        M = deepcopy(A)
        L = [[0.0 for i in xrange(n)] for j in xrange(n)]
        Lt = [[0.0 for i in xrange(n)] for j in xrange(n)]
 
        for k in xrange(0, n):
            for i in xrange(0, k):
                soma = 0.0
                for j in xrange(0, i):
                    soma += M[i][j] * M[k][j]
                M[k][i] = (M[k][i] - soma)/M[i][i]
 
            # Diagonal
            soma = 0.0
            for j in xrange(0, k):
                soma += M[k][j] * M[k][j]
            M[k][k] = sqrt(M[k][k] - soma)
 
        # This is to show the LU form of Cholesky
        for k in xrange(0, n):
            L[k][0:k+1] = M[k][0:k+1]
 
        # Lt is transpose of L
        for i in xrange(0, n):
            for j in xrange(0, n):
                Lt[i][j] = L[j][i]
        return (L, Lt)
 
    def solveLU(A, b):
        n = len(b)
        y = [0.0 for i in xrange(n)]
        x = [0.0 for i in xrange(n)]
        (L, U) = cholesky_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]/U[n-1][n-1]
        for i in xrange(n-1, -1, -1):
            suma = y[i]
            for j in xrange(i+1, n):
                suma = suma - U[i][j] * x[j]
            x[i] = suma/U[i][i]
        return x
 
    # Function's body
    if isSimetric(A):
        return solveLU(A, b)
    else:
        return False

A partir deste código, podemos verificar que a única diferença que ele possui em relação ao algoritmo apresentado no artigo sobre decomposição LU está em utilizar a função cholesky_decomposition ao invés do método de Crout. Isto porque a fatoração de Cholesky, assim como muitos algoritmos de solução direta, é um método de decomposição LU. Ainda pelo código, é possível notar que o fator de Cholesky é a própria matriz \(L \), enquanto que sua transposta é \(U \).

Esta função está disponível em http://www.sawp.com.br/code/linear/cholesky.py, assim como um exemplo de como utilizá-la.

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.

[2] Neide Bertoldi Franco, Cálculo Numérico, Pearson Prentice Hall, (2006), ISBN 85-7605-087-0, pag. 141-145.