Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Partindo de elementos gerados da Fatorização QR para retornar à matriz principal

0 votos
11 visitas
perguntada Nov 10 em Ciência da Computação por Renan Abrantes (16 pontos)  

As bibliotecas para computar a Fatoração QR geram um retorno para uma matriz A: uma matriz R triângulo superior e os vetores de Householder. Escreva uma rotina que pegue a matriz as matrizes de saida, forme a matriz Q a partir dos vetores de Householder, e teste se Q é realmente ortogonal e se o produto de Q por R recupera a matriz de entrada A.

Exercício de computação 3.11 do livro Scientific Computing: An Introductory Survey, Michael T. Heath, Capítulo3: Linear Least Squares

Compartilhe

1 Resposta

0 votos
respondida Nov 10 por Renan Abrantes (16 pontos)  
editado Nov 10 por Renan Abrantes

A fatoração QR realiza uma decomposição de uma matriz A em uma matriz Q e uma matriz R triangular superior de forma que: A imagem será apresentada aqui. onde Q é uma matriz m x m ortogonal.

Este tipo de fatoração é utilizada para resolver o problema dos mínimos quadrados, para alcançar essa fatoração é utilizada a transformação de Householder. Esta transformação busca um vetor ortogonal que zera todas as entradas do vetor de entrada exceto a primeira.

A implementação do algoritmo foi realizado em python, vamos iniciar pela função da transformação de Householder:

def householder(A):
    n = len(A) # Linhas de A
    # Inicializa as variáveis
    R = A
    Q = [[0.0] * n for i in range(n)]

    # Transformação de Householder 
    for k in range(n-1):  # Não é realizado na matrix 1x1

        I = [[float(i == j) for i in range(n)] for j in range(n)] # Matrix identidade

        # Cria vetores utilizados na transformação
        x = [row[k] for row in R[k:]]
        e = [row[k] for row in I[k:]]
        alpha = -cmp(round(x[0]),0) * norm(x) # Função sinal  vezes a norma euclidiana de X

        u = list(map(lambda p,q: p + alpha * q, x, e)) # Realiza a equação de: elemento de x + alpha*elemento de e
        norm_u = norm(u)
        v = list(map(lambda p: p/norm_u, u)) # Faz a normalização do vetor u

        # Encontra os valores de menor complementar para as partições da lista e armazena em uma matriz 
        Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in range(n-k)] for j in range(n-k) ]

        # Preenche os elementos de menor da matriz Q com elementos da matriz identidade
        Q_t = [[ Q_i(Q_min,i,j,k) for i in range(n)] for j in range(n)]

        if k == 0: # Matriz completa
            Q = Q_t
            R = mult(Q_t,A)
        else: # Partições da matriz
            Q = mult(Q_t,Q) # Produto de matrizes Q com Q transposta, resultando em uma matriz simétrica
            R = mult(Q_t,R) 
    return trans_matrix(Q), R # é Realizado a transposição de Q novamente

Esta função tem como entrada uma matriz A e executa uma Decomposição QR baseada na transformação de Householder da matriz A. A função retorna Q, uma matriz ortogonal e R, um matriz triangular superior de modo que A = QR.

Para execução desta função foram utilizadas as funções auxiliares abaixo.

Multiplicação de matrizes na forma X*Y:

def gerar_matriz (n_linhas, n_colunas): # Gera uma matriz nula
    return [[0]*n_colunas for _ in range(n_linhas)]

def mult(X,Y): # Multiplicação de matrizes
    out = gerar_matriz(len(X),len(Y[0])) # Matriz resultado
    for i in range(len(X)): # Linhas de X
        for j in range(len(Y[0])): # Colunas de Y
            for k in range(len(Y)): # Linhas de Y
                out[i][j] += X[i][k] * Y[k][j]
    return out

Transposição de uma matriz M:

def trans_matrix(M): # Realiza a matriz transporta
    return [[ M[i][j] for i in range(len(M))] for j in range(len(M))]

Realizar a norma euclidiana (comprimento do vetor):

def norm(x): # Realiza a norma euclidiana do vetor X
    return pow(sum([x_i**2.0 for x_i in x]),0.5)

Montar a matriz Q transposta preenchendo a parte superior esquerda com elementos da matriz identidade e elementos da matriz com valores menor complementar dos vetores:

def Q_i(Q_min, i, j, k):
    if i < k or j < k:
        return float(i == j)
    else:
        return Q_min[i-k][j-k]

A função sinal (Assume valores -1, 0 e 1):

def cmp(a, b): # Antiga função cmp do Python 2
    return (a > b) - (a < b) 

Para verificar as condições do problema, se a multiplicação da matriz Q e R retorna a matriz A e se Q é ortogonal, foi criado um algoritmo de teste:

import numpy as np
from random import randint
import time 

sucessos = 0
falhas = 0
erro = 0.01

inicio = time.time()
for i in range(10**4):
    m = randint(2,10) # Gera alatoriamente o tamanho da matriz
    A = np.random.rand(m,m) # Gera aleatoriamente valores entre 0 e 1
    Q, R = householder(A)
    Ver = verificar_igualdade(A,mult(Q,R),erro)# Verifica se Q * R retorna a matriz A
    Ver2 = verificar_igualdade(np.linalg.inv(Q),trans_matrix(Q),erro) # Verifica se é ortogonal
    if Ver and Ver2:
        sucessos += 1
    else:
        falhas += 1
fim = time.time()

A biblioteca numpy foi utilizado apenas para faciliar a vizualização da matriz, gerar matrizes aleatórias e inverter a matriz. Neste teste realizado foi testado 10000 matrizes geradas aleatoriamente e de tamanho aleatório (variando de 2 a 10), um valor de erro foi atribuido para desprezar os erros decorrentes de arredondamento. Para verificação se Q é ortogonal foi verificado a igualdade da matriz transposta com a matriz insversa. Por fim, foi verificado o tempo de execução do teste.

A visualização dos dados podem ser observado a partir do seguinte código:

print(str((sucessos/(sucessos + falhas))*100) + '%' + " corresponderam aos critérios")
print()
print("Levou: "+ str(fim - inicio) + " segundos para testar")
print()
print("Última matriz analisada")
print('A')
print(np.matrix(arredonda(A)))
print()
print('Q')
print(np.matrix(arredonda(Q)))
print()
print('R')
print(np.matrix(arredonda(R)))
print()
print("Resultado de Q*R")
print()
print(np.matrix(arredonda(mult(Q,R))))
print()
print("Inversa de Q")
print()
print(np.matrix(arredonda(np.linalg.inv(Q))))
print()
print("Transposta de Q")
print()
print(np.matrix(trans_matrix(Q)))

O resultado obtido a partir da execução do teste pode ser visto abaixo:

100.0% correspondeu aos critérios

Levou: 45.72584366798401 segundos para testar

Última matriz analisada
A
[[0.68 0.28 0.47 0.95 0.49 0.71 0.3 0.59 0.23 0.93]
[0.11 0.82 0.94 0.25 0.74 0.32 0.43 0.74 0.5 0.14]
[0.87 0.28 0.21 0.61 0.65 0.79 0.52 0.64 0.96 0.41]
[0.83 0.05 0.65 0.09 0.03 0.84 0.4 0.37 0.21 0.09]
[0.32 0.09 0.92 0.1 0.92 0.46 0.66 0.68 0.46 0.73]
[0.43 0.65 0.56 0.78 0.55 0.08 0.45 0.62 0.15 0.11]
[0.29 0.61 0.56 0.49 0.04 0.47 0.97 0.25 0.34 0.14]
[0.73 0.71 0.19 0.76 0.08 0.92 0.65 0.01 0.02 1. ]
[0.08 0.34 0.96 0.76 0.82 0.58 0.14 0.9 0.12 0.21]
[0.31 0.01 0.46 0.3 0.96 0.47 0.43 0.66 0.68 0.96]]

Q
[[ 0.4 -0.08 0.43 0.04 -0.05 0.08 -0.33 -0.31 -0.54 0.37]
[ 0.06 0.65 -0.13 0. -0.03 -0.1 -0.18 0.51 -0.47 -0.2 ]
[ 0.51 -0.17 0.34 -0.15 -0.11 0.24 0.11 0.63 0.31 0.04]
[ 0.48 -0.35 -0.42 0.43 -0.29 -0.19 -0.07 -0.05 -0.09 -0.39]
[ 0.19 -0.07 0.2 0.17 0.86 -0.3 -0.06 0.03 0.04 -0.23]
[ 0.25 0.35 -0.01 -0.03 0.11 0.67 0.05 -0.39 0.11 -0.43]
[ 0.17 0.38 0.11 0.52 -0.06 -0.15 0.63 -0.1 0.08 0.32]
[ 0.42 0.27 -0.47 -0.35 0.12 -0.21 -0.22 -0.15 0.33 0.41]
[ 0.04 0.26 0.48 -0.15 -0.35 -0.49 -0.16 -0.24 0.31 -0.36]
[ 0.18 -0.13 -0.05 -0.58 0.03 -0.2 0.61 -0.12 -0.39 -0.18]]

R
[[ 1.71 0.94 1.28 1.45 1.15 1.79 1.37 1.28 1.03 1.42]
[-0. 1.17 0.9 0.77 0.58 0.24 0.66 0.61 0.17 0.15]
[ 0. 0. 0.46 0.61 0.82 0.14 0.05 0.77 0.42 0.23]
[ 0. 0. 0.23 -0.32 -0.61 -0.1 0.21 -0.21 -0.21 -0.74]
[-0. 0. 0.26 -0.17 0.47 -0.08 0.41 0.11 0.18 0.59]
[-0. 0. -0.71 0.01 -0.41 -0.67 -0.29 -0.34 -0.13 -0.52]
[-0. 0. 0.06 -0.06 0.21 0. 0.54 0.12 0.51 0.09]
[-0. -0. -0.13 -0.47 0.12 0. -0. 0.04 0.58 -0.31]
[-0. 0. -0.37 0.05 -0.4 -0. 0. -0.36 -0.24 -0.37]
[ 0. -0. -0.88 0.08 -0.82 -0. 0. -0.83 -0.27 0.28]]

Resultado de Q*R

[[0.68 0.28 0.47 0.95 0.49 0.7 0.3 0.59 0.23 0.92]
[0.1 0.82 0.94 0.24 0.74 0.31 0.42 0.74 0.5 0.13]
[0.87 0.28 0.2 0.6 0.65 0.78 0.52 0.64 0.96 0.41]
[0.82 0.04 0.64 0.07 0.02 0.82 0.39 0.37 0.2 0.08]
[0.32 0.1 0.93 0.11 0.93 0.47 0.67 0.69 0.47 0.75]
[0.43 0.64 0.57 0.77 0.55 0.08 0.44 0.62 0.14 0.11]
[0.29 0.6 0.56 0.49 0.05 0.46 0.96 0.25 0.34 0.13]
[0.72 0.71 0.19 0.75 0.08 0.92 0.65 0.01 0.01 0.99]
[0.07 0.34 0.95 0.76 0.83 0.57 0.13 0.9 0.12 0.21]
[0.31 0.02 0.46 0.29 0.96 0.47 0.44 0.66 0.68 0.96]]

Inversa de Q

[[ 0.4 0.07 0.51 0.48 0.18 0.26 0.18 0.43 0.04 0.18]
[-0.07 0.64 -0.17 -0.34 -0.07 0.35 0.39 0.27 0.25 -0.13]
[ 0.43 -0.13 0.34 -0.42 0.19 -0.01 0.11 -0.47 0.48 -0.05]
[ 0.04 0.01 -0.15 0.43 0.17 -0.03 0.53 -0.35 -0.15 -0.59]
[-0.06 -0.03 -0.11 -0.29 0.86 0.12 -0.06 0.12 -0.36 0.03]
[ 0.08 -0.1 0.24 -0.2 -0.3 0.67 -0.14 -0.21 -0.5 -0.2 ]
[-0.33 -0.18 0.11 -0.08 -0.06 0.06 0.62 -0.21 -0.17 0.6 ]
[-0.31 0.5 0.62 -0.05 0.04 -0.39 -0.1 -0.15 -0.24 -0.12]
[-0.55 -0.47 0.31 -0.09 0.05 0.11 0.08 0.33 0.31 -0.39]
[ 0.37 -0.2 0.04 -0.39 -0.23 -0.42 0.32 0.41 -0.36 -0.18]]

Transposta de Q

[[ 0.4 0.06 0.51 0.48 0.19 0.25 0.17 0.42 0.04 0.18]
[-0.08 0.65 -0.17 -0.35 -0.07 0.35 0.38 0.27 0.26 -0.13]
[ 0.43 -0.13 0.34 -0.42 0.2 -0.01 0.11 -0.47 0.48 -0.05]
[ 0.04 0. -0.15 0.43 0.17 -0.03 0.52 -0.35 -0.15 -0.58]
[-0.05 -0.03 -0.11 -0.29 0.86 0.11 -0.06 0.12 -0.35 0.03]
[ 0.08 -0.1 0.24 -0.19 -0.3 0.67 -0.15 -0.21 -0.49 -0.2 ]
[-0.33 -0.18 0.11 -0.07 -0.06 0.05 0.63 -0.22 -0.16 0.61]
[-0.31 0.51 0.63 -0.05 0.03 -0.39 -0.1 -0.15 -0.24 -0.12]
[-0.54 -0.47 0.31 -0.09 0.04 0.11 0.08 0.33 0.31 -0.39]
[ 0.37 -0.2 0.04 -0.39 -0.23 -0.43 0.32 0.41 -0.36 -0.18]]

Referência: https://www.quantstart.com/articles/QR-Decomposition-with-Python-and-NumPy/

...