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

Como escrever um algoritmo para o processo de ortogonalização de Gram-Schmidt no Python, sem usar um código disponível para tal?

0 votos
212 visitas

2 Respostas

+1 voto
respondida Fev 16 por Rodrigo Fernandes (61 pontos)  
selecionada Fev 22 por Rodrigo Fernandes
 
Melhor resposta

Introdução
Em primeio lugar, é preciso definir o que é o Processo de Gram-Schmidt, para facilitar o passo-a-passo e sua compreensão.
O mesmo consiste produzir uma base ortogonal/normal para um espaço vetorial a partir de um conjunto de vetores linearmente independentes.
Assim, os vetores finais geram o mesmo subespaço vetorial inicial mas sendo ortogonais entre si e com norma igual à 1 (se normalizados).

De forma geral, essas são as etapas a serem seguidas para o \( \mathbb{R}^{n}\):
1ª etapa: escolher os vetores desejados e a quantidade escolhida deve ser igual a \( n\);
2ª etapa: testar se os \( n\) vetores escolhidos são linearmente independentes;
3ª etapa: ortogonalizar;
4ª etapa: normalizar os vetores;

Durante o código, as etapas serão aprofundadas e melhor explicadas.

Demonstração

# importando alguns pacotes usados de forma auxiliar
import numpy as np
import sympy as sy

A demonstração passo a passo será feito no \( \mathbb{R}^{3}\)

# 1ª etapa
v1 = np.array([[1],[1],[1]]) 
v2 = np.array([[0],[1],[1]]) 
v3 = np.array([[0],[0],[1]])

Ou seja, gerando 3 vetores para o \( \mathbb{R}^{3}\)

# 2ª etapa
matriz = np.column_stack((v1,v2,v3)) # gerando matriz com os vetores escolhidos na etapa 1
lambdas, V =  np.linalg.eig(matriz.T) # gerando os autovalores da matriz
_, colunas = sy.Matrix(matriz).T.rref() 
print("Quais colunas/vetores são LI?",colunas) # esse código basicamente nos informa quais colunas, da matriz transposta, são LI

Dessa forma, se todas as colunas/vetores são LI, o processo pode continuar.

# 3ª etapa
w1 = v1 # Aqui, começa o processo de fato. Um vetor, arbitrariamente escolhido, se mantêm igual
w2 = v2 - (v2.T.dot(w1)/(w1.T.dot(w1))*w1) # v2 - (projeção de v2 em w1). essa projeção é feita via: (produto interno de v2 e w1 / produto interno de w1 consigo mesmo) #TODO: explicar melhor projeção (imagem)
w3 = v3 - (v3.T.dot(w2)/(w2.T.dot(w2))*w2) - (v3.T.dot(w1)/(w1.T.dot(w1))*w1)  # v3 - (projeção de v3 em w2) - (projeção de v3 em w1)

Com essa etapa, conseguimos "rotacionar" os vetores com intuito de ficarem ortogonais entre si, usando suas projeções vetoriais.

print("Produto interno entre w1 e w2:", np.round(np.dot(w1.T,w2),4))
print("Produto interno entre w1 e w3:", np.round(np.dot(w1.T,w3),4))
print("Produto interno entre w1 e w3:", np.round(np.dot(w2.T,w3),4)) 

Como o produto interno é igual à zero, provamos que os vetores são ortogonais.

# 4ª etapa
w1 = w1/np.linalg.norm(w1)
w2 = w2/np.linalg.norm(w2)
w3 = w3/np.linalg.norm(w3)

Ao dividir os vetores pelas suas normas, conseguimos normalizar os mesmos. Assim, suas normas novas serão iguais à 1, como provamos com os códigos abaixo.

print("Norma final de w1:", np.round(np.linalg.norm(w1),4))
print("Norma final de w2:", np.round(np.linalg.norm(w2),4))
print("Norma final de w3:", np.round(np.linalg.norm(w3),4)) 

Exemplo 1 (\( \mathbb{R}^{2}\)):

print("############ R2 ############")
# 1ª etapa
v1 = np.array([[2],[3]]) 
v2 = np.array([[-3],[5]]) 

# 2ª etapa
matriz = np.column_stack((v1,v2)) 
lambdas, V =  np.linalg.eig(matriz.T) 
_, colunas = sy.Matrix(matriz).T.rref() 
print("Quais colunas/vetores são LI?",colunas)

# 3ª etapa
w1 = v1 
w2 = v2 - (v2.T.dot(w1)/(w1.T.dot(w1))*w1)
print("Produto interno entre w1 e w2:", np.round(np.dot(w1.T,w2),4))

# 4ª etapa
w1 = w1/np.linalg.norm(w1)
w2 = w2/np.linalg.norm(w2)
print("Norma final de w1:", np.round(np.linalg.norm(w1),4))
print("Norma final de w2:", np.round(np.linalg.norm(w2),4))

O processo pode ser realizado e, dessa forma, os vetores são ortogonais e com norma unitária.

Exemplo 2 (\( \mathbb{R}^{4}\)):

print("############ R4 ############")
# 1ª etapa
v1 = np.array([[1],[1],[0],[0]]) 
v2 = np.array([[0],[1],[1],[0]]) 
v3 = np.array([[0],[0],[1],[1]]) 
v4 = np.array([[1],[0],[0],[1]]) 

# 2ª etapa
matriz = np.column_stack((v1,v2,v3,v4)) 
lambdas, V =  np.linalg.eig(matriz.T) 
_, colunas = sy.Matrix(matriz).T.rref() 
print("Quais colunas/vetores são LI?",colunas)

Como apenas 3 vetores são LI, os quatro vetores escolhidos não são LI em \( \mathbb{R}^{4}\), logo, o processo não pode ser finalizado!

Observação: Um problema do código é sua extensão. Eu poderia, por exemplo, ter criado uma função para cada etapa do processo mas como minha programação em Python ainda é inicial e escrever funções não é trivial, acabei optando por escrever caso a caso. Também tive que arredondar pelo Python não trabalhar com frações.

+1 voto
respondida Fev 22 por Thiago Lappicy (6 pontos)  

Uma outra opção para resolver o problema é criar um código em R ao invés de Python para resolver isso. Para não ficar igual a outra resposta já descrita, será feito um código interativo para qualquer "n" dimensões em R. O programa proposto aqui, retorna as coordenadas dos vetores ortogonais.

As primeiras linhas do código são apenas para que o usuário insira o número de vetores, e depois insira cada vetor em seu devido local. Depois o programa já diz se existe, e quais são, as linhas que NÃO são LI (Linearmente Independentes). Em seguida é criada uma matriz chamada resposta para ser printada no final (no momento se encontra vazia).

Primeiro é feito um for loop para criarmos vetores W como feito na resposta de Rodrigo. Tomamos W1 arbitrariamente para ser o vetor "fixo" e analisamos a projeção dos outros n-1 vetores em cima dele.

Para esse cálculo, criamos um novo for loop para fazer o somatório referente ao cálculo de Wn. Por último preenchemos a matriz resposta criada acima com os valores de Wn já divididos pela norma do vetor (para terem a sua norma = 1). Para facilitar ao usuário, modificamos também o nome das colunas, para que essas sejam representando os valores de Wn (com n variando de 1 até o número de dimensões).

Segue abaixo o código em R com comentários:

GramSchmidt <- function() {

# Perguntas para aparecer no console do R de forma interativa
pergunta <- readline("Quantos vetores se tem? ")
n_vetor <<- as.numeric(pergunta)
matriz <<- matrix(data = NA, ncol = n_vetor, nrow = n_vetor)

for(i in 1:n_vetor) {
  pergunta <- readline(paste0("Escreva seu ", i, "o vetor: "))
  matriz[,i] <<- as.numeric(unlist(strsplit(pergunta, ",")))
  assign(paste0("V", i), matriz[,i])
}

# Aqui é analisado se o vetor tem linhas que são LI ou não
vetor_LI <<- qr(matriz)$rank

# Caso tenha alguma linha LD, ela será printada na tela!
if(vetor_LI < n_vetor) {
  return(paste0("Vetor LD: ",
                which(as.numeric(sprintf("%.5f", 
                                 qr(matriz)$qraux)) == 0)))}

# Cria-se uma matriz para ter os valores das respostas
resposta <<- matrix(data = NA, nrow = n_vetor, ncol = n_vetor)

# Aqui de fato se começa com os cálculos de álgebra linear
for(i in 1:n_vetor){

  if(i == 1) {W1 <<- V1}

  if(i > 1){
    somatorio <- matriz[,i]

  # Aqui é feita a projeção de um vetor i sobre outro i-j
    for(j in 1:(i-1)){
      somatorio <<- somatorio -
        as.vector((get(paste0("V", i)) %*% get(paste0("W", i-j))/
                   get(paste0("V", i-j)) %*% get(paste0("W", i-j)))) *
        get(paste0("W", i-j))
    }

    # É atribuido ao W o valor do somatório
    assign(paste0("W", i), somatorio)
  }
  # Preenche-se a matriz de resposta
  resposta[,i] <- get(paste0("W", i))/
                  sqrt(sum(get(paste0("W", i))^2)) 
}  
colnames(resposta) <- c(paste0("W", 1:n_vetor))

# É retornado na tela a resposta
return(resposta)  
}

Para rodar o código acima, basta escrever no R o seguinte comando:

if(interactive()) GramSchmidt()

Abaixo segue a resposta do programa para os exemplos em R2, R3 e R4 utilizados pelo Rodrigo na outra resposta.

Para R2:

A imagem será apresentada aqui.

Para R3:

A imagem será apresentada aqui.

Para R4:

A imagem será apresentada aqui.

...