Normas matriciais, condicionamento e fatorações

...

Transcript

1. Normas, condicionamento e fatorações Prof. Paulo R. G. Bordoni UFRJ
2. Vamos começar analisando o conceito de norma, para matrizes.
3. Uma norma num espaço vetorial é uma função ∙ ∶ → ℝ satisfazendo: Para , ∈ e ∈ ℝ: I. ≥ 0, II. = 0 se e somente se = 0, III. = , IV. + ≤ + .
4. Portanto é possível definir uma norma para matrizes em ℳ× . Basta garantir as quatro propriedades relacionadas pelo Galileu. Já conferimos que o conjunto ℳ× das matrizes × com a adição de matrizes e a multiplicação por fator de escala é um espaço vetorial.
5. A norma de Frobenius é definida por: = ෍ ,=1 2 = () Surfista, prove que é uma norma. Operacionalmente é como se fosse a norma euclidiana para um vetor de 2 componentes.
6. Na realidade nosso foco estará na aplicação linear associada à uma matriz ∈ ℳ× . Pensar na aplicação linear : ℝ → ℝ associada à matriz envolverá três normas: uma em ℝ outra em ℝ e uma terceira em ℳ× .
7. Vamos provar que: se é uma norma vetorial, então a função ∶ ℳ× → ℝ definida por ⟼ = sup ≠0 é uma norma. É uma norma matricial denominada norma matricial induzida (pelas normas em ℝ e ℝ).
8. Sim, basta considerar o supremo sobre a bola unitária, isto é: = sup =1 . Antes, vamos mostrar que não precisamos tomar o supremo sobre todos os vetores ≠ 0 de ℝ.
9. Observe, Surfista, que para todo ∈ ℝ, temos ∙ = = () Considerando o caso particular em que = Τ 1 , para ∈ ℝ e ≠ 0, temos ∙= 1 = Τ e, obviamente, = 1 para = Τ .
10. Portanto = sup ≠0 = sup =1 E isso vale para qualquer uma das três bolas unitárias, nas ∙ 1 , ∙ ∞ e ∙ 2 .
11. Fantástico, Surfista. Pela interpretação geométrica fornecida pela decomposição SVD segue daí que 2 = sup 2=1 2 = 1 , o maior valor singular de . 1 = 2
12. É o que nossos alunos farão! Senhores, essa afirmação que = sup ≠0 = sup =1 é uma norma é tão importante que merece ser provada!
13. A condição I da definição de norma é imediata: = sup =1 ≥ ≥ 0. A condição III também: = sup =1 = sup =1 = .
14. A parte ⇐ da exigência II da definição de norma é trivial. Se = 0, claramente = 0. Mostre a parte ⇒, Cabelos de Fogo! Se = 0 então sup ≠0 = 0 e portanto = 0, ∀ ∈ ℝ. Logo = 0, ∀ ∈ ℝ e assim = 0.
15. Só falta a condição IV, a desigualdade triangular. Antes de prová-la, vamos provar que ≤ , ∀ ∈ ℝ, uma desigualdade importante em si mesma.
16. Ora, Cabelos de Fogo, sua afirmação decorre da definição de supremo: = sup ≠0 ≥ , ∀ ∈ ℝ, ≠ 0 Claro, pois daí temos ≤ , ∀ ∈ ℝ.
17. Para vetores , ∈ ℝ temos + ≤ + , portanto ∀ ∈ ℝ, segue daí que + = + ≤ + . Então, para na bola unitária, da desigualdade que acabamos de provar, temos + ≤ + ≤ +
18. Notem que + ≤ + , ∀ ∈ ℝ com = 1. Consequentemente + = sup =1 + ≤ + . Brilhante Cabelos de Fogo! Assim é válida a exigência IV, a desigualdade triangular: + ≤ + , ∀, ∈ ℳ× .
19. Agora vamos examinar as normas induzidas 1 e ∞ de uma matriz .
20. Sim Surfista, basta aplicarmos aos vetores dos vértices , , , tanto da 1 (1) como da ∞ (1). Mestra, essas duas bolas unitárias são quadrados. Obter o efeito de uma matriz sobre elas é imediato, não?
21. Por exemplo, a imagem por = 1.1 0.6 0.4 1.8 da ∞ (1) é:
22. E a imagem pela mesma matriz = 1.1 0.6 0.4 1.8 da 1 (1) é:
23. Por inspeção visual, o vetor de tamanho máximo da imagem por A da ∞ (1) é ∗ . Temos ∗ ∞ = max{1.7, 2.2} = 2.2. Portanto ∞ = 2.2.
24. Já, examinando a imagem por da 1 (1), vemos que o vetor de tamanho máximo é ∗ . Temos ∗ 1 = 0.6 + 1.8 = 2.4. Portanto 1 = 2.4.
25. Entretanto, não precisamos gerar as imagens por das bolas ∞ (1) e 1 (1) para calcular as ∞ e 1 . Basta calcularmos a matriz nos vértices , , , das bolas. Como, Mestre?
26. Ora Loirinha, porque ∞ = max ∞ , ∈ ∞ 1 = = max ∞ , = , , , . E, com a taquigrafia óbvia: = , , , = 1 1 , −1 1 , −1 −1 , 1 −1 ⟹ = 11 12 21 22 1 2 = 1 2 = = 11 + 12 21 + 22 , −11 + 12 −21 + 22 , −11 − 12 −21 − 22 , 11 − 12 21 − 22 .
27. Portanto ∞ = max ∞ , x = a, b, c, d = = max{ 11 + 12 , 21 + 22 } Claro, pois ൠ ()1 = 11 + 12 , −11 +12 , −11 −12 , 11 −12 ()2 = 21 + 22 , −21 +22 , −21 −22 , 21 −22 ⟹ ∞ = max ()1 , ()2 = = max{ 11 + 12 , 21 + 22 }
28. Da mesma forma, porque 1 = max 1 , ∈ 1 1 = = max 1 , = , , , . E com a taquigrafia, também óbvia: = , , , = +1 0 , 0 +1 , −1 0 , 0 −1 , ⟹ = 11 12 21 22 1 2 = ()1 ()2 = 11 21 , 12 22 , −11 −21 , −12 −22 .
29. Portanto 1 = max 1 , = , , , = = max{ 11 + 21 , 12 + 22 } Claro, pois neste caso: ቊ = , ⟹ 1 = 1 + 2 = 11 + 21 = , ⟹ 1 = 1 + 2 = 12 + 22
30. Repetindo, para = 11 12 21 22 : • 1 = max{ 11 + 21 , 12 + 22 } • ∞ = max{ 11 + 12 , 21 + 22 } Por esses dois resultados: • 1 é referida como norma (...) das colunas de e • ∞ como norma (...) das linhas de . (...) = (do máximo da soma dos valores absolutos)
31. Conferindo, para = 1.1 0.6 0.4 1.8 tenho: • 1 = 2.4 (colunas de ), • ∞ = 2.2 (linhas de ).
32. Apresento novamente a scipi.linalg.norm().
33. Agora vamos calcular normas matriciais.
34. Especificamente, vamos calcular as normas ∙ , ∙ 2 , ∙ 1 , ∙ ∞ de uma matriz e conferir a igualdade 2 = 1 .
35. Vejam:
36. Passaremos a analisar a sensibilidade da multiplicação matricial. Isto porque ela está diretamente associada à resolução de sistemas lineares.
37. No início do curso estudamos o número de condicionamento de um problema/função. Sim Mestre. Tanto o condicionamento absoluto , = () como o condicionamento relativo , = , ∙ Τ () .
38. Nosso interesse imediato é calcular os números de condicionamento de funções ∶ ℝ → ℝ dadas por ↦ ∈ ℝ, ∈ ℝ. O 1º passo da generalização é imediato: substituímos os módulos por normas. Assim eles serão os números dados por: , = () e , = , ∙ Τ ()
39. Para funções ∶ ℝ → ℝ a () é a matriz jacobiana, = Τ () O 2º passo envolve o entendimento do que é ().
40. Claro e a é a multiplicação de matriz por vetor, de forma que cada é definido por = (1 , 2 , ⋯ , ) = 1 1 +2 2 + ⋯ + Lembrem-se colegas, sistemas lineares são descritos por funções lineares : ℝ → ℝ: = 1 2 ⋮ ∈ ℝ ⟼ = 1 2 ⋮ ∈ ℝ, onde = .
41. Sendo prático, vamos imaginar um sistema 3 × 3 típico: ቐ 21 − 2 + 43 = 1 −1 + 52 + 33 = 3 61 − 22 + 3 = 0 Para esse sistema, temos: ൞ 1 (1 , 2 , 3 ) = 21 − 52 + 43 2 (1 , 2 , 3 ) = −1 + 52 + 3 3 (1 , 2 , 3 ) = 61 − 22 + 3
42. Portanto para a : ℝ3 → ℝ3 do sistema linear do Surfista teremos Τ () = 2 −1 4 −1 5 3 6 −2 1 . Claro, Loirinha! A matriz jacobiana dela é constante, e igual à própria matriz do sistema!
43. Portanto o número de condicionamento absoluto para sistemas lineares é = , = . Resumindo, se ∶ ℝ → ℝ é a multiplicação matricial, = , seu número de condicionamento absoluto é , = .
44. Claro que não colega! A solução de um sistema linear = é dada por = −1. Assim, a função é : ↦ = −1. Brilhante Cabelos de Fogo. O número de condicionamento para o problema da resolução de sistemas lineares é, portanto , = −1 .
45. Assim, o número de condicionamento absoluto ( = , ) de um sistema linear só está definido quando é inversível. E independe tanto do termo independente quanto da incógnita . Só depende da matriz .
46. O número de condicionamento relativo é dado por , = , Τ ∙ () . Bem, só falta entender o termo Τ () . Mestre, se = () = −1 é claro que = .
47. Sim Cabelos de Fogo e portanto Τ () = Τ ≤ . Assim o condicionamento relativo da resolução de sistemas lineares satisfaz , ≤ −1 ∙
48. Quando é ortogonal a igualdade em , ≤ −1 ∙ é atingida. Claro Mestra, pois nesse caso a : ↦ = −1 não altera o tamanho de , não importa a forma de medir tamanho.
49. Conclusão: o número de condicionamento relativo de um sistema linear = é = , = −1 ∙ Portanto os números de condicionamento absoluto e relativo para sistemas lineares = dependem apenas da matriz (e da −1): • = , = , • = , = ∙ −1 .
50. Assim suprimiremos a indicação do sistema = e de , escrevendo apenas () e ().
51. Quando não é inversível, é usual escrever = +∞ e = +∞ Quando é inversível de ∙ −1 = decorre que 1 = = ∙ −1 ≤ −1 . Portanto ≥ 1, para toda matriz inversível .
52. Como antes, quando os números de condicionamento forem pequenos o problema de resolução de sistemas lineares é bem condicionado; caso contrário é mal condicionado.
53. Na numpy.linalg temos a função cond( ):
54. Fiz este programinha, que recebe uma matriz A e devolve os números de condicionamento oferecidos pela numpy.linalg.con():
55. Surfista, observe que: • 2 = - o maior valor singular de A , • −1 2 = Τ 1 - o menor valor singular de A . Assim = Τ Em outras palavras = Τ , é a excentricidade do hiper-elipse que é a imagem por A da circunferência unitária.
56. Os teoremas 12.1 e 12.2 do livro de Trefethen & Bau, exibidos a seguir, mostram a importância dos números de condicionamento de uma matriz A.
58. Eis como estimar o número de casas decimais confiáveis na solução de um sistema linear:
59. A paisagem é magnífica; há muita coisa a ser vista!
60. Á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.
61. A famosa fatoração LU.
62. Rotinas existentes na Referência da scipy.linalg para a fatoração LU.
63. 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 = permite calcular det().
64. Mestre, que outras vantagens Maria Loirinha 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!
65. 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 já vimos, o maior trabalho reside em fatorar (i.é, efetuar o processo de eliminação).
66. Os detalhes para usar a linalg.lu( ):
67. O que linalg.lu( ) retorna:
68. A rotina pode efetuar o pivotamento 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.
69. Conferindo a fatoração LU:
70. Um 1º exemplo em que não necessidade de trocas de linhas,
71. 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.
72. Uma outra forma de chamar a fatoração LU da SciPy. Como no 1º exemplo, não houve a necessidade de pivotamento.
73. Chamando a 2ª forma da fatoração LU da SciPy. Agora com pivotamento.
74. Esta é lu_factor( ) que efetua a fatoração . Ela recebe a matriz A e retorna uma matriz contendo a fatoração LU e um vetor P contendo os índices de pivoteamento.
75. Os parâmetros retornados. Esta rotina usa a LAPACK.
76. 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.
77. Resolvendo um sistema linear Ax = b via fatoração LU com pivotamento.
78. Ela é boa para resolver sistemas Ax=b com muitos b’s, mas a mesma A.
79. Vejam, entrei com uma matriz 4x4 e 2 termos independentes. Poderiam ser muitos: 10, 50, etc.
80. Uma brevíssima explicação do que é a fatoração de Cholesky.
81. É 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, pois = .
82. 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?
83. 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.
84. As rotinas do SciPy para a fatoração de Cholesky.
85. A explicação para passagem de parâmetros.
86. Surfistas e Loirinhas, fiz este programa usando a fatoração de Cholesky.
87. Uma execução. Deixo a cho_fator() e a cho_solve() para vocês.
88. Esta rotina só efetua a fatoração. Deve ser usada com a cho_solve().
89. Fácil, fácil. Só troquei o nome da rotina e acrescentei o parâmetro de retorno Eh_inf.
90. Fácil nada, você errou. O apressado come cru! Repeti os dados da Mestra e não conferiu.
91. Errou porque não prestou atenção no aviso: Vou corrigir!
92. Foi só aplicar a rotina triu() na Cho. Vejam o resultado a seguir:
93. Como você disse Surfista, fácil, fácil!
94. Esta rotina resolve o sistema após fatorado pela rotina cho_factor().
95. Vou resolver o sistema linear abaixo. A matriz do sistema é a do exemplo anterior. 5 2 1 0 2 7 3 1 1 0 3 1 6 −1 −1 2 = 1 2 3 4
96. O Mestre chamou a Cho_fator() e depois a Cho_solve(). Percebam que ele não precisou zerar a parte inferior da matriz Cho. Além disso, b entrou como vetor-linha e foi transformado em vetor-coluna.
97. A execução, com conferência a resposta.
98. Esta é específica para matrizes de banda. Só fatora.
99. Esta resolve o sistema linear, após a fatoração pela cholesky_banded().
100. Façam um programa para resolver um sistema linear = no qual , além de ser simétrica e positiva-definida, também é de banda. Vai cair na prova Mestra?
101. Tchau, até a próxima aula!