2.1.3 Sistemas Lineares — Métodos Diretos — Decomposição de Cholesky
1. A Fatoração de Cholesky </p>
A Fatoração de Cholesky expressa uma matriz simétrica como o produto,de uma matriz triangular (chamada de Fator de Cholesky) pela sua transposta , na forma:
$$ A = LL* $$
onde será triangular superior.
Sabendo-se que é simétrica quando seus elementos obedecem à uma formação para toda linha e coluna . Desta forma, a matriz será idêntica à sua transposta: . 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 é simétrica, podemos utilizar a Equação 1 para gerar e a partir de . Ou seja, ao multiplicarmos e igualarmos seus termos, obtemos as seguintes relações:
$$ l_{ki} = \frac{ a_{ki} – \sum_{j=1}^{i-1} l_{ij}~l_{kj} }{ l_{ii} } $$
onde . Além disso, o elemento da diagonal será obtido por
$$ 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]. </p>
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
</p>
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 , enquanto que sua transposta é .
Esta função está disponível em http://www.sawp.com.br/code/linear/cholesky.py, assim como um exemplo de como utilizá-la. </p> </p>
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/.