2.1.5. Sistemas Lineares — Métodos Diretos — Decomposição QR
1. Decomposição QR
Se assumirmos que é possível decompor uma matriz real em um produto de matrizes , onde seja ortogonal e é uma matriz triangular-superior de diagonal não-nula, então quando for uma matriz não-singular, esta decomposição é única e sempre existe.
Para obter esta decomposição, aplicamos em a sequência de transformações de Holseholder, garantindo que os elementos da diagonal de cada estágio (Rdiag) sejam não-negativos. Isto é possível, uma vez que cada transformação é determinada pela mudança de sinal. Assim, teremos que
$$ P_{n-1} P_{n-2} \cdots P_{1} A = R $$
Como cada é ortogonal, temos que , onde
$$ Q = \left[ ~ P_{n-1} ~ P_{n-2} ~ \ldots ~ P_1 ~ \right]^T $$
Como cada é unicamente determinado pela condição de não-negatividade, se o elemento da diagonal é não-nulo quando é aplicado, então a decomposição é única para não-singular.
Como o produto de matrizes triangulares é triangular, assim como o produto de matrizes ortogonais é ortogonal, quando temos que é triangular-inferior ou é um fator ortogonal da transformação da original . Disso, a convergência do processo é determinado pelo comportamento da sequência , desde que . </p>
1.1. Resolvendo sistemas lineares com decomposição QR </p>
A matriz de dimensões pode ser decomposta em termos de duas matrizes e como
R \\
0
\end{array} \right] $$
com a matriz sendo quadrada e triangular-superior de dimensão .
Para o problema que tem como matriz de coeficientes na forma
podemos utilizar a decomposição QR para resolvê-lo com o sistema triangular na forma
Isto é, a solução
vêm de
(R^T R)^{-1} R^T Q^T b = = R^{-1} Q^T b $$
portanto
Separando este problemas, podemos solucionar em dois passos:
- Passo 1: Calcular
- Passo 2: Resolver o sistema
2. Implementação
def QR(A):
"""
Return the (Q,R) from QR Decomposition.
QR(A)
INPUT:
* A: matrix to decompose
return: the tuple (Q, R)
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/
OBS: This is a implementation of a QR with Householder transformation and
the code is adapted from the NIST Template Numerical Toolkit (TNT)
by Roldan Pozo. see: http://math.nist.gov/tnt/
Sep 2010
"""
(m, n) = (len(A), len(A[0]))
QR = deepcopy(A)
nrange = xrange(n)
mrange = xrange(m)
minor = min(m, n)
Rdiag = [0.0 for i in nrange]
Q = [[0.0 for i in nrange] for j in mrange]
R = [[0.0 for i in nrange] for j in nrange]
for least in xrange(minor):
norm = 0.0
for row in xrange(least, m):
norm += QR[row][least] * QR[row][least]
a = sqrt(norm)
if QR[least][least] > 0:
a = -a
Rdiag[least] = a
if a != 0.0:
QR[least][least] -= a
for col in xrange(least + 1, n):
phi = 0.0
for row in xrange(least, m):
phi -= QR[row][col] * QR[row][least]
phi /= a * QR[least][least]
for row in xrange(least, m):
QR[row][col] -= phi * QR[row][least]
# get R
for row in xrange(minor - 1, -1, -1):
R[row][row] = Rdiag[row]
for col in xrange(row + 1, n):
R[row][col] = QR[row][col]
# get Q
for least in xrange(m - 1, minor - 1, -1):
Q[least][least] = 1.0
for least in xrange(minor - 1, -1, -1):
Q[least][least] = 1.0
if QR[least][least] != 0.0:
for col in xrange(least, m):
phi = 0.0
for row in xrange(least, m):
phi -= Q[row][col] * QR[row][least]
phi /= Rdiag[least] * QR[least][least]
for row in xrange(least, m):
Q[row][col] -= phi * QR[row][least]
return (Q, R)
2.1. Utilizando a decomposição para resolver sistemas lineares
def solveQR(A, b):
"""
Return a list with the solution of linear system.
solveQR(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 solve(A, b):
""" Solution by QR decomposition R x = Q^T b.
"""
n = len(b)
x = [0.0 for i in xrange(n)]
(Q, R) = QR(A)
# Step 1: Solve the system: y = Q^T b
y = matmul(transpose(Q), b)
# Step 2: Solve the R system: Rx = y (regressive replacement)
x[n - 1] = y[n - 1] / R[n - 1][n - 1]
x[n - 2] = (y[n - 2] - R[n - 2][n - 1] * x[n - 1]) / (R[n - 2][n - 2])
for j in xrange(n - 3, -1, -1):
sum = 0.0
for k in xrange(j + 1, n):
sum += R[j][k] * x[k]
x[j] = (y[j] - sum) / R[j][j]
return x
# Function's body
return solve(A, b)
Estas funções estão disponíveis em http://www.sawp.com.br/code/linear/QR.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/.