Aulas atrás, estudamos o básico da
NumPy e da MatPlotLib.
Iniciaremos agora nosso estudo da
SciPy – Scientific Python.
Slide 3
Slide 3 text
E o que é, precisamente,
Computação Científica
Mestre?
É o que você aprenderá ao longo
deste curso, Loirinha.
Vou começar pelo aspecto
histórico!
Slide 4
Slide 4 text
http://history.siam.org/
A visão da Sociedade de
Matemática Aplicada e
Industrial – SIAM:
Slide 5
Slide 5 text
Procure por Scientific
Computing na Web,
Surfista.
Farei a busca com
o Google, Mestra.
Slide 6
Slide 6 text
Peguei o retorno do Google
e separei o que achei
interessante:
Slide 7
Slide 7 text
Comecei pela Wikipedia.
Eis o que achei:
Slide 8
Slide 8 text
Nosso Mestre já citou
Maple, Mathematica e MatLab.
Como ele, a Wikipedia também
indica Python e SciPy!
Slide 9
Slide 9 text
O 2º que separeifoi o
site da INTEL de HPC.
HPC é High Performance
Computation
Slide 10
Slide 10 text
Separei algumas chamadas
desse site que me
interessaram – estão
todas na mídia recente!
Slide 11
Slide 11 text
A SciPy é a 5ª da lista retornada pelo
Google.
E é reconhecida como uma biblioteca
fundamental para computação científica!
É com ela que trabalharemos
Galileu.
Atualíssima, veja a seguir:
Slide 12
Slide 12 text
A EuroSciPy 2016 foi na
Alemanha e a de 2015
foi na Inglaterra.
Slide 13
Slide 13 text
A SciPy Latin America
2016 foi em Florianópolis
de 16 à 20 de maio.
Slide 14
Slide 14 text
Este é o “site”
da SciPy.
Slide 15
Slide 15 text
Clicando no ícone da SciPy
Library, chegamos aqui.
Agora escolha
Documentation, Loirinha!
Slide 16
Slide 16 text
Surfista, em
Documentation escolha o
Guia de Referência da
SciPy.
Slide 17
Slide 17 text
Todos os tópicos
tratados no Tutorial, na
última versão da SciPy:
Slide 18
Slide 18 text
Em Introduction temos uma
breve descrição da SciPy.
Lembrem-se, os Tutoriais são
para novatos!
Slide 19
Slide 19 text
É importante assinalar as
convenções de importação
utilizadas por padrão.
Slide 20
Slide 20 text
E agora todos os tópicos
tratados na Referência,
na última versão da SciPy:
Assinalei os sub-
pacotes que
estudaremos.
Slide 21
Slide 21 text
Antes de mergulharmos
fundo nas matrizes, fale um
pouco sobre Entrada/Saída
Manuel.
Vou listá-los, mas não
entrarei em detalhes!
Slide 22
Slide 22 text
Outros métodos além dos
padronizados para lidar
com arquivos:
Slide 23
Slide 23 text
Métodos padronizados?
Quais?
Aqueles da Numpy referidos
logo no início, filha:
Slide 24
Slide 24 text
Os métodos da Numpy
para lidar com arquivos:
Slide 25
Slide 25 text
Continuação:
Slide 26
Slide 26 text
A parte final:
Slide 27
Slide 27 text
Passageiros com destino à
scipy.linalg, embarque imediato pela
plataforma 9 ¼.
Slide 28
Slide 28 text
Precisaremos pesquisar muito
para descobrir tudo sobre
Álgebra Linear Computacional na
biblioteca do castelo.
Slide 29
Slide 29 text
Vou buscar com
Google, Galileu!
Surfista, antes disso, para
chegarmos preparados, faça uma
busca na Internet por Álgebra
Linear Computacional.
Slide 30
Slide 30 text
Quando busquei por Numerical
Linear Algebra, o 1º resultado
foi da Wikipedia:
Slide 31
Slide 31 text
Não deixem de olhar as listas
indicadas de software de análise
numérica e de bibliotecas
numéricas.
Muita indicações
importantes nesta
dica!
Slide 32
Slide 32 text
As primeiras estações da ferrovia, antes
de Hogwarts, são paradas obrigatórias,
chamadas BLAS, LAPACK e NETLIB.
Temos que desembarcar para conhecê-las.
Slide 33
Slide 33 text
Antes quero mostrar a
Estação da Luz e a
Júlio Prestes, na
cidade de São Paulo.
Hoje também são a sede do
Museu da Língua Portuguesa e
da Orquestra Sinfônica do
Estado de São Paulo, OSESP.
Slide 34
Slide 34 text
No content
Slide 35
Slide 35 text
A estação Leopoldina, no Rio.
Se a Princesa visse o
descaso iria chorar!
Os responsáveis
públicos pelo abandono
deveriam ser presos!
Slide 36
Slide 36 text
A Estação da Central do
Brasil, também no Rio.
Dias depois do comício do
Jango, veio o golpe de 64!
Slide 37
Slide 37 text
As bibliotecas fundamentais
sobre Álgebra Linear
Computacional são a BLAS e
LAPACK.
Slide 38
Slide 38 text
• 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.
Slide 39
Slide 39 text
Fui buscar BLAS
no Google!
Assinalei o
endereço da netlib!
Slide 40
Slide 40 text
Este é o principal repositório
de software científico que a
comunidade científica mundial
dispõe:
Slide 41
Slide 41 text
Respondendo antes de
você perguntar,
Loirinha:
Slide 42
Slide 42 text
A BLAS na
www.netlib.org
Slide 43
Slide 43 text
FAQ é comigo!
Inclusive já marquei minha
1ª pergunta:
Slide 44
Slide 44 text
Eis sua resposta,
Loirinha:
Slide 45
Slide 45 text
O endereço para chegar
aqui:
http://www.netlib.org/blas/
No final, os links para
todas as rotinas da
BLAS!
Slide 46
Slide 46 text
Sim, Surfista.
Clicando num link, somos
encaminhados para a rotina
correspondente na LAPACK.
O código está em FORTRAN.
Slide 47
Slide 47 text
O código, criado pelo
poderosíssimo feiticeiro Jack
Dongarra em 1978. A última
atualização é de 2011.
Slide 48
Slide 48 text
A LAPACK na
Wikipedia:
Slide 49
Slide 49 text
A continuação da LAPACK
na Wikipedia:
Slide 50
Slide 50 text
O site oficial da
LAPACK, claro:
www.netlib.org.
Slide 51
Slide 51 text
Agora vamos à biblioteca do
Castelo de Hogwarts.
É lá que encontraremos o
Tutorial da scipy.linalg.
Slide 52
Slide 52 text
Marquei o motivo para
usar a linalg da SciPy.
Mas já tem linalg na NumPy.
Por que repetir?
Slide 53
Slide 53 text
Antes que você pergunte,
Surfista, vá ler o resto desta
discussão lá no Tutorial!
Slide 54
Slide 54 text
A estrutura do tutorial Linear
Algebra é a seguinte:
Álgebra linear
Apresentação
Rotinas básicas
Decomposições
Funções matriciais
Matrizes especiais
Slide 55
Slide 55 text
Os tópicos em rotinas
básicas são :
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
Slide 56
Slide 56 text
Vamos começar pelos
determinantes, que já vimos
na aula passada!
E indo à referencia do
scipy.linalg, nas rotinas
básicas, temos:
Slide 57
Slide 57 text
A resposta à tua
pergunta está ao pé da
página Loirinha:
por outro algoritmo.
Mas vocês já provaram que o cálculo de
determinantes é inviável, na prática.
Então, como o pacote scipy faz isto?
Slide 58
Slide 58 text
E o que é fatoração
LU?
Mais adiante exploraremos tudo o que há
na Linalg sobre fatorações.
Em particular sobre a fatoração LU.
Slide 59
Slide 59 text
Respondendo à tua pergunta Loirinha,
escrevemos como um produto de
duas matrizes triangulares e :
= ∙
A só tem 1’ na diagonal. Assim
det = det ∙
= det ∙ det = det()
Slide 60
Slide 60 text
Neste programa entramos com
uma matriz A, via teclado, e em
seguida chamamos a função
det(A).
Slide 61
Slide 61 text
Voltando às rotinas básicas da
Scipy, vejamos o que há no
tópico “Calculando normas”.
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
Slide 62
Slide 62 text
A SciPy possibilita calcular
diversas normas.
Tanto normas de vetores como
de matrizes.
Slide 63
Slide 63 text
Por enquanto, vamos examinar apenas as
normas vetoriais.
Veremos as normas matriciais mais adiante.
Para elas precisaremos de um pouco mais
de teoria.
É importantíssimo observar que
todas aquelas cujo parâmetro
ord é negativo não são normas.
Slide 64
Slide 64 text
Surfista, apresente alguns
contra-exemplos que justifiquem
a última afirmação do Mestre!
Quando ord = inf temos a norma do máximo.
Para ord = 1 temos a norma da soma e para
ord = 2 temos norma euclidiana.
Elas são anotadas ∞
, 1
e 2
respectivamente.
As outras são anotadas
, ( = ).
Slide 65
Slide 65 text
Já provamos que 2
= 1
2 + 2
2 + ⋯ +
2
define uma norma. A verificação que
1
= 1
+ 2
+ ⋯ +
é uma norma é óbvia.
Da mesma forma, a comprovar que
∞
= max{ 1
, 2
, ⋯ ,
} define uma
norma também é fácil.
Slide 66
Slide 66 text
Para prová-la precisamos da
desigualdade de Hölder:
, ≤ ′ com ′ =
−1
.
A dificuldade para provar que
=
|1
| + |2
| + ⋯ + |
|
é uma norma está na desigualdade triangular:
+
≤
+
.
Ela é conhecida na literatura como desigualdade
de Minkowski.
Slide 67
Slide 67 text
Observem que a desigualdade de
Hölder desempenha um papel
análogo ao da desigualdade de
Cauchy-Schwarz.
Não vamos provar nenhuma das duas.
Porém, se você estiver interessado,
Surfista, já sabe as palavras-chave para
buscar as demonstrações.
Slide 68
Slide 68 text
Na referência do SciPy
encontramos os detalhes dos
parâmetros para usar a função
norma:
Slide 69
Slide 69 text
Por enquanto só veremos
normas de vetores.
Atenção para a referência!
Slide 70
Slide 70 text
Um exemplo de utilização da
função norm( ) para vetores.
Slide 71
Slide 71 text
Mestre, parece que os valores
numéricos das normas estão
decrescendo ao valor da ∞
.
É verdade minha filha.
Prove que para ∀ ∈ ℝ temos
lim
→∞
= ∞
Slide 72
Slide 72 text
Para cada norma vetorial em ∙ : ℝ → ℝ podemos
associar um subconjunto do próprio ℝ.
Estou me referindo às bolas abertas de raio > 0
com centro num ponto ∈ ℝ
∙ , ,
= ∈ ℝ . . − <
Slide 73
Slide 73 text
Para a norma euclidiana em ℝ2, é o conjunto dos
pontos em ℝ2 satisfazendo:
∙ 2, ,
= = (, ) ∈ ℝ2 . . − 2
< ,
sendo = , o centro
E como 2
= 2 + 2, temos
− 2
= ( − )2 + ( − )2.
Slide 74
Slide 74 text
Que é o conjunto dos pontos interiores
à circunferência 2 + 2 = 1, de centro
na origem e raio 1.
Portanto a bola unitária na ∙ 2
, com centro na
origem, corresponde ao conjunto definido pelos
pontos , ∈ ℝ2 cujas coordenadas satisfazem
a desigualdade
2 + 2 < 1
Slide 75
Slide 75 text
Fiz um programa que desenha a
circunferência e bola unitária em ℝ2.
Vejam o resultado:
Slide 76
Slide 76 text
O programa que desenha
as bolas unitárias na ∙
:
Slide 77
Slide 77 text
Bolas na ∙
,
para diversos
Slide 78
Slide 78 text
Surfista, não vá jogar bola
com 1
ou ∞
.
Você pode machucar seu pé!
É, são bolas que escapam
da nossa intuição
euclidiana.
Slide 79
Slide 79 text
Note que as bolas “crescem’
da ∙ 1
para a ∙ ∞
Slide 80
Slide 80 text
Agora vamos explorar o que há na
scipy.linalg para resolução de
sistemas lineares.
Slide 81
Slide 81 text
Vamos ao tópico
“Resolvendo sistemas
lineares” de Rotinas
básicas.
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
Slide 82
Slide 82 text
No Tutorial da SciPy, clicando
em “Resolvendo sistemas
lineares”, obtemos:
Slide 83
Slide 83 text
As rotinas básicas
no Guia de Referência do
scipy.linalg para
resolução de sistemas lineares.
Slide 84
Slide 84 text
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)
Slide 85
Slide 85 text
Recordando as
definições:
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
=
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:
Slide 88
Slide 88 text
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:
Slide 89
Slide 89 text
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 é, obviamente
1
= 2
= ⋯ = 100
= 1/
Mestra, fiz até este
programa:
Veja a
substituição
Slide 93
Slide 93 text
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 .
Slide 94
Slide 94 text
Esta é rotina da SciPy
para resolver sistemas
triangulares.
Slide 95
Slide 95 text
O detalhamento
dos seus
parâmetros.
Slide 96
Slide 96 text
Eis meu
programa:
Slide 97
Slide 97 text
Examinando as
informações sobre o
comando solve( ).
Slide 98
Slide 98 text
Este programa recebe uma
matriz A e um vetor b.
Em seguida, chama solve(A,b)
para resolver o sistema linear Ax = b.
Slide 99
Slide 99 text
Outra forma de
obter
informações do
solve( )
Slide 100
Slide 100 text
Uma forma de calcular a
inversa de uma matriz A é
resolvendo a equação matricial
= usando solve( )
Slide 101
Slide 101 text
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):
Slide 102
Slide 102 text
=
× × 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ó:
Esta é rotina específica da
scipy.linalg para resolver
sistemas lineares cuja matriz
é de banda.
Slide 105
Slide 105 text
Novamente, os
detalhes sobre os
parâmetros
Slide 106
Slide 106 text
Um programinha
usando a
solve_banded( )
Slide 107
Slide 107 text
Um outro
programa usando
a solve_banded( )
Slide 108
Slide 108 text
Uma execução do
programa.
Slide 109
Slide 109 text
Um exemplo maior!
Nesse caso A
possui 72 = 100-28
zeros.
Slide 110
Slide 110 text
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:
Slide 111
Slide 111 text
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:
Slide 112
Slide 112 text
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
Slide 113
Slide 113 text
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.
Slide 114
Slide 114 text
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!
Slide 115
Slide 115 text
=
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 :
Slide 116
Slide 116 text
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
Slide 117
Slide 117 text
As matrizes ortogonais
(e as unitárias, no caso complexo)
desempenham um papel fundamental em
Álgebra Linear Computacional.
Aguardem!
Slide 118
Slide 118 text
Esta rotina é apropriada
para matrizes hermitianas
(simétricas no caso real) e
positivo-definidas.
Slide 119
Slide 119 text
Detalhamento dos
parâmetros.
Slide 120
Slide 120 text
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.
Slide 121
Slide 121 text
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.
O nome do programa é:
sist_banda_sim_def_pos.py
Slide 122
Slide 122 text
Esta 2ª parte apresenta a matriz
A no formato padrão n x n e
confere a solução, , verificando
se, de fato, = .
Slide 123
Slide 123 text
Uma execução do programa.
Se a matriz A não for definida
positiva, o programa acusa.
Vejam a seguir.
Slide 124
Slide 124 text
O 4º elemento da diagonal principal
é o responsável por tornar a matriz
A não positiva definida.
A linalg acusa o fato!
ERRO
Slide 125
Slide 125 text
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, = .