Álgebra Linear Computacional II

Inicio apresentando o site de História da Análise Numérica e Computação Científica da SIAM. Depois apresento o NetLib, a BLAS e a LAPACK. Em seguida volto para a scipy.linalg para mostrar como utilizar as diversas rotinas para resolver sistemas lineares.envolvendo matrizes com características particulares (simétricas, de banda, etc). Depois mostro a fatoração LU da sicpy.linalg, com suas diversas versões para matrizes com características particulares. Fecho com uma discussão sobre a resolução de sistemas lineares via determinantes.

Transcript

1. Álgebra linear computacional II Prof. Paulo R. G. Bordoni UFRJ
2. Aula passada iniciamos nossa viagem às terras da scipy.linalg. Há muita coisa a ser vista! Passageiros com destino à scipy.linalg.
3. Surfista, procure na Internet por mapas da Álgebra Linear Computacional no mundo da Computação Científica. Vou buscar no Google, Mestre!
4. Comecei com scientific computation, Mestre:
5. Python e SciPy!
6. http://history.siam.org/ Já com history of scientific computation, achei:
7. Quando busquei por Numerical Linear Algebra, o 1º resultado foi da Wikipedia:
8. • BLAS: Basic Linear Algebra Subprograms, • LAPACK: Linear Algebra Package. Repetindo: BLAS e LAPACK, bibliotecas altamente otimizadas que implementam os algoritmos básicos de Álgebra Linear Computacional. Não deixem de olhar as listas indicadas de software de análise numérica e de bibliotecas numéricas.
9. Fui buscar a BLAS no Google!
10. Um repositório de software importantíssimo!
11. O que é a Netlib:
12. A BLAS na www.netlib.org
13. É isso aí, Loirinha:
14. O que é a BLAS:
15. A LAPACK na Wikipedia:
16. A continuação de LAPACK na Wikipedia:
17. O site oficial da LAPACK, claro: www.netlib.org.
18. 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, descobri o caminho das pedras. Muito mais que isto, percebi claramente a importância das bibliotecas BLAS e LAPACK !
19. 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!
20. No Tutorial da SciPy, clicando em “Resolvendo sistemas lineares”, obtemos:
21. Então vamos focar 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
22. As rotinas básicas no Guia de Referência do scipy.linalg para resolução de sistemas lineares.
23. As características particulares da matriz do sistema linear definirão nosso roteiro para estudar os métodos de resolução do sistema. Tipo da matriz A em Ax = b: 1. Triangular 2. Genérica 3. De banda 4. Simétrica e definida positiva (hermitiana)
24. 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
25. 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
26. = 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 det = 120 (= 1 × 3 × 5 × 2 × 4) = −3 0 0 0 −2 0 0 0 2 7 4 2 5 0 5 1 det = 0 (= −3 × 0 × 5 × 1) O determinante de uma matriz triangular é o produto dos termos da diagonal:
27. Se A é uma matriz ordem n, então det = det(). Um sistema linear = possui uma única solução quando, e apenas, quando det() ≠ 0. Os dois teoremas abaixo, são fatos sobejamente conhecidos sobre sistemas lineares e determinantes. Entretanto determinantes conduzem a conclusões erradas sobre sistemas lineares. Vejam a seguir:
28. Entretanto para = 0.01 temos det = 0.01100 = 10−200, um legítimo ZERO em qualquer computador! Então concluiríamos que o sistema não tem solução. 0 ⋮ 0 ⋯ 0 ⋮ ⋱ 0 0 0 0 1 2 ⋮ 100 = 1 1 ⋮ 1 A solução do sistema linear acima é, óbviamente 1 = 2 = ⋯ = 100 = 1/
29. 11 1 = 1 21 1 31 1 41 1 +22 2 +32 2 +42 2 +33 3 +43 3 = 2 = 3 +44 4 = 4 É muito simples, Mestra. Se os ≠ 0, da 1ª eq. tiramos 1 = 1 /11 . Substituindo esse valor na 2ª eq. temos 2 = (2 − 21 1 )/22 . Depois 3 = 3 − (31 1 + 32 2 ) /33 . Finalmente, 4 = 4 − (41 1 + 42 2 + 43 3 ) /44 Só para saborear, vamos resolver o sisteminha 4x4 abaixo, cuja matriz é triangular inferior.
30. Um programa simplérrimo! Marquei a parte que faz a substituição progressiva.
31. Surfista, observe no código que o número de multiplicações cresce de 1 a n como uma PA de razão 1. Logo, o total de multiplicações é ( + 1)/2. Um algoritmo 2 .
32. Esta rotina é da SciPy para resolver sistemas triangulares.
33. O detalhamento dos seus parâmetros.
34. Eis meu programa:
35. Examinando as informações sobre o comando solve( ).
36. Este programa recebe uma matriz A e um vetor b. Em seguida, chama solve(A,b) para resolver o sistema linear Ax = b.
37. Outra forma de obter informações do solve( )
38. Uma forma de calcular a inversa de uma matriz A é resolvendo a equação matricial = usando solve( )
39. 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):
40. = × × 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ó:
41. = × × 0 0 0 × × × 0 0 × 0 0 × × 0 × × 0 × × × × × × = ∎ × × × × × × × × × × × × × × × × ∎ ∎ ∎ ∎ ∎ ∎ ∎ Mestre, prefiro pensar em “deitar” as diagonais.
42. Esta é rotina específica da scipy.linalg para resolver sistemas lineares cuja matriz é de banda.
43. Novamente, os detalhes sobre os parâmetros
44. Um programa usando a solve_banded( )
45. Uma execução do programa.
46. Um exemplo maior! Nesse caso A possui 72 = 100-28 zeros.
47. Quando é igual à sua transposta: = . Simétrica Quando satisfaz: ≠ 0 ⟹ = ,=1 > 0 Definida positiva Quando sua transposta é a inversa: = −1 ⟺ = = Ortogonal Continuando a visita ao zoológico matricial:
48. 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. • Atenção: a 5. não permite concluir que a matriz produto é uma matriz simétrica. Algumas propriedades da transposição:
49. De fato, o produto de matrizes simétricas não é uma matriz simétrica! E, além disso, o produto 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
50. 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 = 2 1 0 obtemos = 0.
51. 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!
52. = 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 :
53. 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
54. As matrizes ortogonais (e as unitárias, no caso complexo) desempenham um papel fundamental em Álgebra Linear Computacional. Aguardem!
55. Esta rotina é apropriada para matrizes hermitianas (simétricas no caso real) e positivo-definidas.
56. Detalhamento dos parâmetros.
57. O exemplo de armazenamento do Tutorial da scipy.linalg. Armazena-se apenas a diagonal e a banda superior (ou inferior). Em vez de 36 = 6 × 6 elementos, apenas 15 = 3 × 6 − 3.
58. Esta 1ª parte do programa recebe uma matriz A (formato de banda), um termo independente b e devolve a solução x usando a rotina para matrizes simétricas, positivo-definidas e de banda da scipy.linalg.
59. Esta 2ª parte apresenta a matriz A no formato padrão n x n e confere a solução, , verificando se, de fato, = .
60. Uma execução do programa. Se a matriz A não for definida positiva, o programa acusa. Vejam a seguir.
61. O 4º elemento da diagonal principal é o responsável por tornar a matriz A não positiva definida. A linalg acusa o fato! ERRO
62. Que péssimo exemplo Mestre! Por que não incluiu um tratamento de erro, para evitar que seu programa pipocasse assim? Esta 2ª parte apresenta a matriz A no formato padrão n x n e confere a solução, , verificando se, de fato, = .
63. Há muita coisa a ser vista em nossa viagem à scipy.linalg. A paisagem é magnífica!
64. Á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.
65. A famosa fatoração LU.
66. Mestre, que vantagem Maria leva com essa tal de fatoração LU? Antes do Mestre te responder, jovem, proponho que discuta com seus colegas sobre essa cultura de levar vantagem em tudo. Grande Filósofo!
67. Respondendo tua pergunta, Loirinha: Tendo a fatoração = , resolver o sistema linear = é equivalente a resolver = . Então resolvemos primeiro o sistema = e depois o = . Note que esses dois são sistemas triangulares, facílimos de resolver! Conforme veremos depois, o maior trabalho reside em achar a fatoração.
68. Para completar a resposta do Mestre:
69. Sim, de = , segue que det() = det det det(). Veremos que L só possui 1’s na diagonal e que P é ortogonal. Assim, det = ± det = ±11 22 ⋯ A fatoração = também permite calcular det().
70. Rotinas existentes na Referência da scipy.linalg para a fatoração LU.
71. 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.
72. 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.
73. Conferindo a fatoração LU:
74. Um 1º exemplo em que não necessidade de trocas de linhas,
75. Já neste exemplo é evidente a necessidade de permutar a L0 com a L2 e depois a L2 com a L3. Atente para as matrizes A e P.
76. Uma outra forma de chamar a fatoração LU da SciPy. Como no 1º exemplo, não houve a necessidade de pivotamento.
77. Chamando a 2ª forma da fatoração LU da SciPy. Agora com pivotamento.
78. 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.
79. Os parâmetros retornados. Esta rotina usa a LAPACK.
80. 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.
81. Resolvendo um sistema linear Ax = b via fatoração LU com pivotamento.
82. Uma brevíssima explicação do que é a fatoração de Cholesky.
83. É essa a ideia Surfista. Além disso, vale a volta; veja na próxima transparência! Confirmando, se A é simétrica e positiva definida está garantida a existência de uma matriz triangular superior tal que = . É como se U fosse a raiz quadrada de A, Mestre!
84. Suponha que uma matriz não-singular A possa ser fatorada na forma = , com U triangular superior. Então A será simétrica e definida positiva. É difícil provar isto Mestra?
85. E também que A é definida positiva, isto é ≠ 0 ⟹ > 0. De fato assumindo ≠ 0 e = teremos ≠ 0, devido à não-singularidade de U. Então = () = = = 2 > 0. Bem, minha querida, de = segue que: • A é simétrica, pois = ()= ()= = • U é não-singular, pois det()2 = det() = det() ≠ 0.
86. As rotinas do SciPy para a fatoração de Cholesky.
87. Surfista e Loirinha, deixaremos com vocês a tarefa de construir alguns programas usando essas rotinas de Cholesky. É uma dica p’rá prova colega!
88. A explicação para passagem de parâmetros.
89. Esta rotina só efetua a fatoração. Deve ser usada com a próxima.
90. Esta rotina resolve o sistema após fatorado pela rotina cho_factor( ).
91. Esta é específica para matrizes de banda. Só fatora.
92. Esta resolve o sistema linear, após a fatoração pela cholesky_banded( ).
93. Tchau, até a próxima aula!
94. Uma leitura adicional sobre determinantes e sistemas lineares.
95. Cramer apresentou sua fórmula para resolver sistemas lineares em 1750. Antes disso, Cardano, no final do século XVI, trabalhou com determinantes 2x2. Leibniz usou determinantes em 1693 para resolver sistemas lineares.
96. Um determinante de 2ª ordem, é denotado e definido por: Notem que o cálculo exige 2 multiplicações e 1 subtração. No computador serão 3 flops.
97. ≠ A regra de Cramer para resolver um sistema 2x2 como este é (prove, Surfista): O cálculo exige 3x2 multiplicações, 3 adições e 2 divisões. Ao todo 11 flops.
98. = 2 3 1 −1 2 1 1 3 = −5 5 = −1 = 3 1 −1 3 2 1 1 3 = 10 5 = 2 Aplicando a regra de Cramer para resolver o sistema 2x2 + = + = − obtemos:
99. Um determinante de 3ª ordem, é definido por: Desenvolvendo os determinantes de 2ª ordem, obtemos:
100. A regra de Cramer para resolver um sistema 3x3 como o acima é dada por:
101. Portanto para resolver um sistema 3x3 pela regra de Cramer são necessárias: (3/2)4! multiplicações, 4!-4 adições e 3 divisões. Notem o fatorial no número de flops. Foram necessárias (3/2)3! multiplicações e 3!-1 adições/subtrações para cada determinante 3x3.
102. Este é um sistema linear constituído por n equações lineares à n incógnitas. Sendo breve: um sistema n x n
103. Ax = b Aliás, Lewis Carrol conta num livro, que a notação para matrizes que usamos hoje foi proposta por Leibniz. Na forma matricial esse sistema se reescreve:
104. O determinante det A de uma matriz A de ordem n é definido recursivamente, como na próxima transparência.
105. Para n > 1, o desenvolvimento pela − é linha do det é dado por det = 1 1 + 2 2 + ⋯ + onde = (−1)+ e é o determinante de ordem n-1 obtido a partir da matriz A excluindo-se sua − é linha e − é coluna. Há um desenvolvimento semelhante a partir da − é coluna. Para n = 1, det A = a 11 .
106. det(ℴ 4 ) α det(ℴ 3 ) β det(ℴ 3 ) γ det(ℴ 3 ) δ det(ℴ 3 ) + − + − + − α det(ℴ 2 ) β det(ℴ 2 ) γ det(ℴ 2 ) det(ℴ 3 ) + − + − α det(ℴ 1 ) β det(ℴ 1 ) det(ℴ 2 ) + − 2 mult. 1 adiç. 3 mult. 2 adiç. 4 mult. 3 adiç. O cálculo de um determinante de ordem 4 envolve: 4 × 3 × 2 × 1 = 4! multiplicações 3 × 2 × 1 = 3! adições/subtrações
107. A regra de Cramer para resolver um sistema linear n x n como este é: Se = det ≠ 0 a solução é dada por 1 = 1 , 2 = 2 , ⋯ , = onde é o determinante da matriz A com a coluna k subtituída pelo termo independente b.
108. Portanto, para resolver um sistema linear n x n precisaremos efetuar o cálculo de n+1 determinantes de ordem n e n divisões. Assim, o total de operações envolvidas é: • + 1 ! multiplicações • da ordem de n! adições • n divisões
109. Bobagem! Para o Titan – Cray XK7 é mole! O tempo gasto, por qualquer computador, para resolver sistemas 30 x 30 é inimaginável !
110. A tabela mostra o tempo gasto pelo Cray XK7, para resolver um sistema n x n: n tempo 20 48.41 min. 25 7.37 séculos 30 1.076 idade univ. Surfista querido, para um sistema 30x30 ele demora um pouco mais que a idade do universo!
111. O programinha que fez os cálculos.
112. Definitivamente: é inviável resolver sistemas lineares pela regra de Crammer! Jamais esqueça de contar o número de operações envolvidas num algoritmo!