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 como o produto de uma matriz triangular (chamada de Fator de Cholesky) pela sua transposta , na forma:
onde será triangular superior.
Uma vez que se identifica que uma matriz é hermitiana, 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} } $$
Todavia, podemos evitar o cálculo da raiz quadrada na última expressão, caso encontremos uma decomposição 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
\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 e , 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/.