Álgebra linear computacional II

Resolução de sistemas lineares com as rotinas da scipy.linalg (matriz genérica,, de banda, e simétrica e definida positiva) e detalhes do método de eliminação, de Gauss, com argumentos para efetuar pivotamento parcial por coluna. Programação tradicional para obter a solução, sem pivotamento

Transcript

1. Álgebra linear computacional Prof. Paulo R. G. Bordoni UFRJ
2. Pois é Surfista, a história mostra que, desde 1947, a comunidade de Computação Científica, tem trabalhado arduamente no desenvolvimento da Álgebra Linear Computacional. Em particular na resolução de sistemas lineares. Mestre, achei o caminho das pedras. Muito mais que isto, percebi claramente a importância das bibliotecas BLAS e LAPACK !
3. Mestre, no tocante a resolução de sistemas lineares, proponho adotarmos a tática de operacionalizar primeiro. Depois vamos analisar os aspectos teóricos. Não gosto muito dessa ideia, mas como acredito na sua experiência didática, concordo!
4. Voltando então às rotinas básicas do Tutorial, foquemos nossa atenção no tópico “Resolvendo sistemas lineares”. Rotinas básicas Achando a inversa Resolvendo sistemas lineares Achando o determinante Calculando normas Resolvendo problemas lineares de mínimos quadrados e pseudo-inversas Inversa generalizada
5. Surfista, a orientação do Tutorial da scipy.linalg é usar a solve( ) [e não a inversa, inv( )] para obter a solução de um sistema linear.
6. O que é estabilidade numérica, Mestre? Se é mais rápido é porquê envolve menos operações, mas como contá-las?
7. O Surfista está certo, Mestra. Eis aí uma razão para não operacionalizarmos às cegas. Resmungão... Loirinha querida, mais adiante estudaremos erros, propagação de erros e estabilidade numérica.
8. Recortei as rotinas básicas na Referência do scipy.linalg para resolução de sistemas lineares. Além da solve( ), existem rotinas apropriadas para matrizes com característica especiais.
9. Aliás existe toda uma fauna de matrizes = [ ] cujas definições devemos recordar: Quando todos os seus elementos fora da diagonal são nulos: ≠ ⟹ = 0 Diagonal Quando todos os elementos acima da diagonal são nulos: < ⟹ = 0 Triangular inferior Quando todos os elementos abaixo da diagonal são nulos: > ⟹ = 0 Triangular superior
10. Exemplos: −2 0 0 0 3 0 0 0 6 Diagonal 3x3 3 0 0 0 0 3 0 0 1 8 4 2 2 0 5 1 Triangular inferior 4x4 1 4 2 5 2 0 3 1 0 7 0 0 0 0 0 0 5 0 0 9 2 0 3 1 4 Triangular superior 5x5
11. Continuando no zoológico: Quando é igual à sua transposta: = . Simétrica Quando satisfaz: ≠ 0 ⟹ = > ,=1 0 Definida positiva Quando sua transposta é a inversa: = −1 ⟺ = = Ortogonal
12. 1. + = + 2. ()= 3. = 4. ( )= 5. = Observações: • De 1 e 2 decorre que matrizes simétricas constituem um subespaço do espaço de todas as matrizes. • A 5 significa que a matriz produto é uma matriz simétrica. Algumas propriedades da transposição:
13. O produto de matrizes simétricas não é uma matriz simétrica!. E, além disso, também não é comutativo. Vejam os contra-exemplos! 3 2 2 1 0 1 1 −2 = 2 −1 1 0 0 1 1 −2 3 2 2 1 = 2 1 −1 0
14. A matriz A = 1 0 0 0 −2 0 0 0 3 não é definida-positiva, pois se = temos: = 1 0 0 0 −2 0 0 0 3 = = −2 3 = 2 − 22 + 32 e escolhendo = 1 2 0 obtemos = 0.
15. A matriz = 5 −1 2 −1 2 2 2 2 3 é definida-positiva, pois para = temos: = − 2 + + 2 2 + 2 + 2, uma soma de quadrados, que só se anula para = = 0 0 0 . Surfista, confira essa afirmação da Mestra!
16. = sen cos 0 −cos sen 0 0 0 1 sen − cos 0 cos sen 0 0 0 1 = 1 0 0 0 1 0 0 0 1 Loirinha, agora é tua vez de conferir a afirmação da Mestra! A matriz = sen − cos 0 cos sen 0 0 0 1 é ortogonal pois :
17. Examinando as informações sobre o comando solve( ).
18. Este programa recebe uma matriz A e um vetor b. Em seguida, chama solve(A,b) para resolver o sistema linear A*x = b.
19. Outra forma de obter informações do solve( )
20. Esta rotina é apropriada para resolver sistemas triangulares. Ela é mais rápida!
21. O detalhamento dos seus parâmetros.
22. Eis uma forma de calcular a inversa de uma matriz triangular superior.
23. Mestra esse algoritmo para achar a inversa é válido para qualquer tipo de matriz, não? Sim minha filha e é o mais rápido! Mas você tem que usar o solve( ) para um sistema cuja matriz não possui características particulares.
24. Uma matriz = [ ] possui banda inferior k quando = 0 para − > . E possui banda superior k, quando = 0 para − > . = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × Eis uma matriz A com banda inferior 2 (em laranja) e banda superior 1 (em verde):
25. = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ Se A é uma matriz 100x100 teríamos 10.000 elementos para armazenar. Nesse armazenamento de banda precisaremos apenas de 4x100-4 ∎ ∎ ∎ ∎ A ação para armazenar uma matriz de banda como essa é “esquecer os zeros e deixar cair”! Vejam só:
26. Esta é rotina específica da scipy.linalg para resolver sistemas lineares cuja matriz é de banda.
27. Novamente, os detalhes sobre os parâmetros
28. Um programa usando a solve_banded( )
29. Uma execução do programa.
30. Uma 2ª execução do programa.
31. Um zoológico bem complexo: A adjunta de = [ ], anotada AH, é definida por = [ ]. É a conjugada da transposta. Adjunta Quando satisfaz: = . É a correspondente de uma matriz simétrica, com elementos complexos, Hermitiana, ou auto-adjunta Quando sua adjunta é a inversa: = −1 ⟺ = = Unitária
32. E esta é apropriada para matrizes hermitianas e positivo-definidas.
33. Detalhamento dos parâmetos.
34. Exemplo do Tutorial da scipy.linalg
35. Álgebra linear Apresentação Rotinas básicas Decomposições Funções matriciais Matrizes especiais Agora vamos explorar as fatorações (Decompositions) existentes na scipy.linalg.
36. A famosa fatoração LU.
37. A famosa fatoração LU – continuação.
38. Rotinas existentes na Referência da scipy.linalg para a fatoração LU.
39. Esta é uma das rotinas do scipy.linalg para realizar a fatoração PLU. Ela só realiza a fatoração; não apresenta a solução.
40. A rotina pode efetuar o pivoteamento por linha (ou não). Ela retorna: • as três matrizes P, L e U • ou apenas duas, a PL que é a triangular inferior permutada e a U.
41. Esta é uma rotina alternativa à anterior. Ela recebe a matriz A e retorna uma matriz contendo a fatoração LU e um vetor contendo os índices de pivoteamento.
42. Os parâmetros retornados.
43. Esta rotina recebe a fatorada LU e o vetor piv fornecidos pela LU_factor( ). Recebe também o termo independente b e retorna a solução x.
44. Uma brevíssima explicação do que é a fatoração de Cholesky.
45. As rotinas do SciPy para a fatoração de Cholesky.
46. A explicação para passagem de parâmetros.
47. Esta rotina só efetua a fatoração. Deve ser usada com a próxima.
48. Esta rotina resolve o sistema após fatorado pela rotina cho_factor( ).
49. Esta é específica para matrizes de banda. Só fatora.
50. Esta resolve o sistema linear, após a fatoração pela cholesky_banded( ).
51. Carl Friedrich Gauss 30/04/1777 – 23/02/1855 Vamos encerrar este conjunto de transparências detalhando o famoso Método de Eliminação, de Gauss. Na verdade ele não foi criado por Gauss. Ele já aparece no cap. 8 do livro chinês “Os nove capítulos sobre a arte matemática”, composto por gerações de mestres ao longo dos séculos X a II A.C.
52. Fui buscar estas informações na Wikipedia!
53. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Vamos agora descrever como resolver um sistema linear como o abaixo, pelo método de eliminação. Primeiro só descreveremos o processo em linhas gerais, usando as equações.
54. Em cada etapa utilizaremos operações elementares sobre as linhas. Elas são três: Somar ou subtrair linhas Multiplicar uma linha por um número real Trocar duas linhas de posição Nosso objetivo será atingido após várias etapas consecutivas.
55. 11 1 + 12 2 + ⋯ + 1 = 1 21 1 + 22 2 + ⋯ + 2 = 2 ⋯ 1 1 + 2 2 + ⋯ + = Na 1ª etapa, eliminaremos os termos envolvendo 1 nas equações após a 1ª, efetuando operações elementares nas linhas. 11 1 + 12 2 + ⋯ + 1 = 1 22 2 + ⋯ + 2 = 2 ⋯ 2 2 + ⋯ + = Os coeficientes e os termos independentes mudam após essas operações – por isto o ~ .
56. 11 1 + 12 2 + ⋯ + 1 = 1 22 2 + ⋯ + 2 = 2 ⋯ 2 2 + ⋯ + = Repetimos o processo para o subsistema sem a 1ª linha: Um sistema com n-1 equações a n-1 incógnitas.
57. Repetindo n-1 vezes o processo chegamos a um sistema triangular superior. 11 1 + 12 2 + ⋯ + 1 = 1 = 22 2 + ⋯ + 2 = 2 ⋯ Todos os coeficientes e termos independentes são diferentes, exceto na 1ª linha.
58. O último passo consiste em resolver o sistema triangular. Algo extremamente simples. O processo consiste, portanto, em eliminarmos incógnitas nas equações. Daí o nome: processo de eliminação. Alguns textos falam em escalonamento.
59. Sim! Como é feita, de fato, a eliminação? Mestres, vocês descreveram o processo em linhas gerais! Só papo...
60. 11 12 ⋯ 1 21 22 ⋯ 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋯ 1 2 ⋮ = 1 2 ⋮ A x = b Pois é... Agora vamos refazer tudo usando a forma matricial do sistema linear:
61. Acrescentando o vetor b à última coluna da matriz A, obtemos uma matriz referida na literatura como “matriz aumentada”. | = 11 12 ⋯ 1 1 21 22 ⋯ 2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋮ ⋯
62. | = 11 12 ⋯ 1 1 21 22 ⋯ 2 2 ⋮ 1 ⋮ 2 ⋱ ⋮ ⋮ ⋯ Linha 2 O processo de eliminação consiste em efetuar as operações elementares sobre as linhas da matriz aumentada. Linha 1 Linha n
63. As operações elementares são: Linha 1 Linha 1 Linha 2 Linha 2 Linha 1 21 /11 × − Linha 3 Linha 3 Linha 1 31 /11 × − Linha n Linha n Linha 1 1 /11 × − Multiplicadores
64. O final do processo é o mesmo: uma matriz triangular superior. Repetimos apenas para evidenciar que não precisamos envolver as incógnitas no processo – só os vetores linha da matriz aumentada. | = 11 12 ⋯ 1 1 22 ⋯ 2 2 ⋱ ⋮ ⋮ Observem que todos os elementos da matriz triangular aumentada são, eventualmente, diferentes dos iniciais.
65. Quando pensamos nas equações, fica óbvio que as operações elementares não alteram a solução. Isto não é tão óbvio quando olhamos para a matriz aumentada.
66. 1 /11 31 /11 Multiplicadores Observem que os multiplicadores foram escolhidos de forma a anular todos os 1 abaixo do 11 . O 11 é intitulado “pivot”. Essas operações exigem que o pivot seja não-nulo: 11 ≠ 0. Caso 11 = 0 trocamos a linha 1 com outra onde 1 ≠ 0.
67. Um parêntese em cinza e negro!
68. Um outro motivo é oriundo da maneira como números reais são tratados no computador. A representação IEEE 754 de ponto flutuante, no padrão single, será ilustrada na próxima transparência. Evitar divisão por zero não é a única razão para efetuar trocas de linhas!
69. 1 bit para o sinal 4 bytes = 32 bits 8 bits para o expoente 23 bits para a fração O padrão single na representação IEEE 754 possui 32 bits, distribuídos da seguinte forma:
70. No single, a fração é constituída por 23 bits: 1 2 ⋯ 23 com = 0 1 . 23 bits para a fração E o valor da fração é dado por: = 1 2−1 + 2 2−2 + ⋯ + 23 2−23, sendo o restante dos bits desprezado!
71. Acreditem: A maioria esmagadora dos números reais na representação IEEE 754 de ponto flutuante, no padrão single, possui um erro na 24 casa binária: 2−24. O Mestre acabou de afirmar que a precisão de um single é de 23 bits.
72. Multiplicar por 2 desloca esse erro p casas binárias para a esquerda. Perdem-se p bits na precisão! Impressionante Sherlock, nunca tinha ouvido falar disto!
73. Retornando ao método da eliminação, de Gauss, e voltando a colocar os holofotes no pivot.
74. 1 /11 31 /11 Multiplicadores Assumam que o pivot 11 = 10−3. Então ao dividirmos por 11 estaremos multiplicando por 1.000 (≅ 1024 = 210) os elementos 21 , 31 , … , 1 . 11 ≅ 2−10
75. 1 /11 31 /11 Multiplicadores Em outras palavras, Mestre, os multiplicadores 21 /11 , 31 /11 , … , 1 /11 perderam 10 dos seus 23 preciosos bits de precisão. 11 ≅ 2−10
76. Que horror Mestra. É por isso que os pivots tem que ser (homens) grandes! Sim minha Loirinha querida! Despreze pivots pequenos, eles estragam o jogo. Sempre vá em busca dos maiores possíveis.
77. Uma forma de minimizar a propagação de erro consiste em localizar o 1 de maior valor absoluto e trocar a linha onde ele está com a linha 1. Esta é a estratégia conhecida na literatura como pivoteamento parcial por coluna. Na realidade estamos apenas apontando a “ponta do iceberg”. Analisaremos tudo isto em detalhes, mais adiante!
78. Antes que você reclame, Surfista, vamos implementar um programa para o método de eliminação.
79. Agora o famoso método de eliminação de Gauss. Porque eliminaram o Gauss, Mestre? Coitado!
80. Primeiro veremos a programação no estilo clássico para o método de eliminação, de Gauss
81. Para facilitar o entendimento, vamos mostrar 1º a execução do programa. Depois discutiremos seu código. Estes são a matriz A e o vetor b de nosso exemplo.
82. pivot Eis a 1ª etapa do processo de eliminação. As operações elementares visam colocar zeros na 1ª coluna, abaixo do pivot.
83. Mestre, seu programa está errado! A matriz A atualizada deveria ter sido esta aqui, com zeros abaixo do pivot. Mil perdões Loirinha! Esqueci de dizer que usei o espaço dos zeros para guardar os multiplicadores.
84. A 2ª etapa do processo de eliminação. As operações elementares visam colocar zeros na 2ª coluna, abaixo do pivot. Novamente, guardei os multiplicadores no lugar dos zeros.
85. A etapa final do processo de eliminação. A solução x foi obtida por retro-substituição.
86. Uma foto da parte final do processo: • A matriz A foi transformada na matriz triangular superior que eu marquei em azul. • E os multiplicadores estão em verde, abaixo da diagonal.
87. Esse é o início do código. Nada excepcional. Apenas definições básicas para iniciar.
88. Aqui são definidas duas funções que vão mostrar o que está acontecendo. Didática ...
89. Esta é a fase da eliminação. O loop em k define as etapas do processo de eliminação. É definido o pivot e, se ele é não-nulo, criamos uma matriz-coluna, Mult, para guardar os multiplicadores..
90. No loop em j atualizamos cada elemento da linha i efetuando a operação elementar sobre ele. O loop em i percorre as linhas abaixo do pivot e define um multiplicador para cada linha.
91. Ainda dentro do loop em i efetuamos a operação elementar sobre o elemento da linha i do termo independente b. E guardamos esse multiplicador na linha i da coluna do pivot da matriz A e também na linha i do vetor Mult.
92. Agora a fase da retro-substituição. • Da última equação achamos diretamente o valor de x n . • Substituímos esse valor na penúltima equação e achamos o valor de x n-1 . • Assim por diante, de trás para frente, até chegar à 1ª equação, de onde achamos x 1 .
93. Tchau, até a próxima aula!