2.1.5. Sistemas Lineares — Métodos Diretos — Decomposição QR

1. Decomposição QR

Se assumirmos que é possível decompor uma matriz real \(A \) em um produto de matrizes \(QR \) , onde \(Q \) seja ortogonal e \(R \) é uma matriz triangular-superior de diagonal não-nula, então quando \(A \) for uma matriz não-singular, esta decomposição é única e sempre existe.

Para obter esta decomposição, aplicamos em \(A \) a sequência \(P_k \) 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 \(P_j \) é ortogonal, temos que \(A = Q R \) , onde

\( Q = \left[ ~ P_{n-1} ~ P_{n-2} ~ \ldots ~ P_1 ~ \right]^T \)

Como cada \(P_j \) é unicamente determinado pela condição de não-negatividade, se o elemento da diagonal \(a_{jj}^{(j)} \) é não-nulo quando \(P_j \) é aplicado, então a decomposição é única para \(A \) não-singular.

Como o produto de matrizes triangulares é triangular, assim como o produto de matrizes ortogonais é ortogonal, quando temos que \(E_k \) é triangular-inferior ou é um fator ortogonal da \(k-esima \) transformação da original \(A_1 \) . Disso, a convergência do processo \(QR \) é determinado pelo comportamento da sequência \(E_k \) , desde que \(A_{k+1} = E_{k}^{-1} A_1 E_k \) .

 

1.1. Resolvendo sistemas lineares com decomposição QR

A matriz \(A \) de dimensões \((m,n) \) pode ser decomposta em termos de duas matrizes \(Q \) e \(R \) como

\( A = Q \left[ \begin{array}{c}
R \\
0
\end{array} \right] \)

com a matriz \(R \) sendo quadrada e triangular-superior \(R \) de dimensão \((n,n) \) .

Para o problema que tem \(A \) como matriz de coeficientes na forma

\( A x = b \)

podemos utilizar a decomposição QR para resolvê-lo com o sistema triangular \(R \) na forma
\( R x = Q^T b \)

Isto é, a solução
\( x = \left( A^T A \right)^{-1} A^T b \)

vêm de
\( x = \left( R^T Q^T Q R \right)^{-1} R^T Q^T b =
(R^T R)^{-1} R^T Q^T b = = R^{-1} Q^T b \)

portanto
\( R x = Q^T b \)

Separando este problemas, podemos solucionar \(x \) em dois passos:

  1. Passo 1: Calcular \(y = Q^T b \)
  2. Passo 2: Resolver o sistema \(R x = y \)

 

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/.

 

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/Householder_transformation

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