BCM0505-15 - Processamento da Informação
Laboratório L06 – 26/05
Instruções
- Membros das turmas NA1 e NB{4,5} deverão submeter códigos em Python para todos os exercícios nesta página em um único arquivo
.py
via Tidia até 22/06 às 19h – sob a entrada E10.
Exercícios
-
(01) Matrizes e Sistemas Tridiagonais.
Uma matriz quadrada $T\in\mathbb{R}^{n\times n}$ é tridiagonal se $T[i][j]==0$ sempre que $|j-i|\not\in\{0,1\}$. Logo, uma matriz tridiagonal $T$ possui a seguinte estrutura: $$T = \left(\begin{array}{ccccccccc} t_{0,0} & t_{0,1} & 0 & 0 & 0 & \cdots & 0 & 0\\ t_{1,0} & t_{1,1} & t_{1,2} & 0 & 0 & \cdots & 0 & 0\\ 0 & t_{2,1} & t_{2,2} & t_{2,3} & 0 & \cdots & 0 & 0\\ 0 & 0 & t_{3,2} & t_{3,3} & t_{3,4} & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots\\ \vdots & \vdots & \vdots & \vdots & \ddots & \ddots & \ddots & \vdots\\ 0 & 0 & 0 & 0 & 0 & \cdots & t_{n-1,n-2} & t_{n-1,n-1}
\end{array}\right).$$ Claramente, uma matriz tridiagonal $T$ pode ser representada utilizando-se três vetores — o que garante uma economia de $n^2 - (3n-2)$ entradas.
Escreva uma função que recebe três vetores que representam uma matriz tridiagonal $T$ e um vetor $b\in\mathbb{R}^n$ e determina um vetor $x\in\mathbb{R}^n$ tal que $$Tx=b.$$ Sua função deve determinar $x$ em tempo $O(n)$, isto é, linear no tamanho da entrada.
Dica: especialize os algoritmos de eliminação gaussiana e retro-substituição.Solução:
Supondo que os três vetores recebidos são $$ p = [t_{1,0}, t_{2,1}, \ldots, t_{n-1,n-2}],\qquad q = [t_{0,0}, t_{1,1}, \ldots, t_{n-1,n-1}],\qquad r = [t_{0,1}, t_{1,2}, \ldots, t_{n-2,n-1}], $$ com $n-1$, $n$ e $n-1$ elementos, respectivamente,def linsys_tridiag (p: [float], q: [float], r: [float], b: [float]) -> [float]: < a ser escrita >
def gaussian (A: [[float]], b: [float]):
n = len (A)
for k in range (n-1):
l = max (range (k, n), key = lambda i: abs (A[i][k]))
if A[l][k] == 0: continue
A[l], A[k] = A[k], A[l]
b[l], b[k] = b[k], b[l]
pivot = A[k][k]
for i in range (k+1, n):
c = -A[i][k] / pivot
A[i][k] = 0
for j in range (k+1, n):
A[i][j] += c * A[k][j]
b[i] += c * b[k]
$\square$
-
(02) Traço de Matrizes.
O traço de uma matriz quadrada $A\in\mathbb{R}^{n\times n}$ é dado por: $$\mathsf{tr}(A) = \sum_{i=0}^{n-1} A[i][i].$$ Isto é, o traço, $\mathsf{tr}(A)$, de $A$ é igual à soma dos elementos na diagonal de $A$.
Escreva uma função que recebe duas matrizes quadradas $A,B\in\mathbb{R}^{n\times n}$ e devolve $\mathsf{tr}(AB)$. Note que não é preciso realizar o produto $AB$.Solução:
Seja $C=AB$. Temos que $C[i][i] = \sum_{k=0}^{n-1} A[i][k]\cdot B[k][i]$. Logo, $$ \mathsf{tr}(C) = \sum_{i=0}^{n-1} C[i][i] = \sum_{i=0}^{n-1} \sum_{k=0}^{n-1} A[i][k]\cdot B[k][i] $$def trace (A: [[float]], B: [[float]]) -> float: return sum (sum (A[i][k] * B[k][i] for k in range (len (A))) for i in range (len (A)))
$\square$
-
(03) Quadrados Mágicos.
Uma matriz quadrada $Q\in\mathbb{Z}^{n\times n}$ é um quadrado mágico se todo elemento em $\{1,2,\ldots,n^2\}$ ocorre em $Q$ uma única vez e cada linha, cada coluna, a diagonal principal e a diagonal secundária de $Q$ somam a um mesmo valor. Por exemplo, a matriz $Q$ de dimensões $3\times 3$ abaixo $$Q=\left(\begin{array}{ccc} 2 & 7 & 6 \\ 9 & 5 & 1 \\ 4 & 3 & 8 \end{array}\right)$$ é um quadrado mágico. De fato, cada inteiro em $\{1,2,\ldots,9\}$ ocorre uma única vez em $Q$ e as três linhas, três colunas e duas diagonais de $Q$ somam $15$.
Escreva uma função que recebe $Q$ e decide se $Q$ é um quadrado mágico.Solução:
def is_magic_square (Q: [[int]]) -> bool: # verifica se todos os números entre 1 e n*n são usados em Q x = [True] * (n := len (Q)) m = n*n for i in range (n): for j in range (n): if not (1 <= Q[i][j] <= m): return False else: x[Q[i][j]-1] = True if not all (x): return False # verifica se as diagonais principal e secundária possuem a mesma soma r = s = 0 for i in range (n): r += Q[i][i] s += Q[i][n-1-i] if r != s: return False # verifica se as linhas e colunas possuem a mesma soma das diagonais for i in range (n): s = t = 0 for j in range (n): s += Q[i][j] t += Q[j][i] if s != r or t != r: return False return True
$\square$
-
(04) Matrizes Simétricas Positivas Definidas.
Uma matriz simétrica $A\in\mathbb{R}^{n\times n}$ é positiva definida se $x^\top Ax>0$ para todo $x\in\mathbb{R}^n\setminus\{(0,0,\ldots,0)\}$.
Equivalentemente, $A$ é simétrica positiva definida (SPD) se todos os seus auto-valores são positivos.
Escreva uma função que recebe uma matriz $A$ e determina se $A$ é simétrica positiva definida.
Dica: utilize eliminação gaussiana.Solução:
Como $A$ é simétrica, todos os seus auto-valores — considerando possíveis multiplicidades — são reais (a prova deste fato é simples, tente fazê-la). Como $A$ é positiva definida, todos estes auto-valores são positivos. Logo, durante uma eliminação gaussiana em $A$, nenhum pivô é nulo. O código abaixo está devidamente adaptado.def is_SPD (A: [[float]]) -> bool: n = len (A) for k in range (n-1): if (pivot := A[k][k]) == 0: return False for i in range (k+1, n): c = -A[i][k] / pivot A[i][k] = 0 for j in range (k+1, n): A[i][j] += c * A[k][j] return True
$\square$
-
(05) Decomposição de Cholesky.
Uma matriz simétrica positiva definida $A\in\mathbb{R}^{n\times n}$ pode ser escrita como $$A = LL^\top,$$ em que $L\in\mathbb{R}^{n\times n}$ é triangular inferior. Exemplo: $$ \left(\begin{array}{rrr} 4 & 12 & -16 \\ 12 & 37 & -43 \\ -16 & -43 & 98 \end{array}\right) = \left(\begin{array}{rrr} 2 & 0 & 0 \\ 6 & 1 & 0 \\ -8 & 5 & 3 \end{array}\right) \left(\begin{array}{rrr} 2 & 6 & -8 \\ 0 & 1 & 5 \\ 0 & 0 & 3 \end{array}\right). $$
Escreva uma função que recebe $A$ e devolve $L$ como enunciadas.Solução:
Vamos usar uma matriz simétrica positiva definida (SPD) $3\times 3$ para descrever o algoritmo de Cholesky-Banachiewicz, que obtém a decomposição desejada.
Seja $A$ dada e suponha que $L$ é triangular inferior tal que $A=LL^\top$. Temos então, que $$ \left(\begin{array}{ccc} a_{0,0} & a_{0,1} & a_{0,2} \\ a_{1,0} & a_{1,1} & a_{1,2} \\ a_{2,0} & a_{2,1} & a_{2,2} \end{array}\right) = \left(\begin{array}{ccc} \ell_{0,0} & 0 & 0 \\ \ell_{1,0} & \ell_{1,1} & 0 \\ \ell_{2,0} & \ell_{2,1} & \ell_{2,2} \end{array}\right)
\left(\begin{array}{ccc} \ell_{0,0} & \ell_{1,0} & \ell_{2,0} \\ 0 & \ell_{1,1} & \ell_{2,1} \\ 0 & 0 & \ell_{2,2} \end{array}\right) = \left(\begin{array}{ccc} \ell_{0,0}^2 & \ell_{0,0} \ell_{1,0} & \ell_{0,0} \ell_{2,0} \\ \ell_{0,0} \ell_{1,0} & \ell_{1,0}^2 + \ell_{1,1}^2 & \ell_{1,0} \ell_{2,0} + \ell_{1,1} \ell_{2,1} \\ \ell_{0,0} \ell_{2,0} & \ell_{1,0} \ell_{2,0} + \ell_{1,1} \ell_{2,1} & \ell_{1,0}^2 + \ell_{1,1}^2 + \ell_{2,2}^2 \end{array}\right), $$ donde decorre que $$ L = \left(\begin{array}{ccc} \sqrt{a_{0,0}} & 0 & 0 \\ a_{1,0} / \ell_{0,0} & \sqrt{a_{1,1} - \ell_{1,0}^2} & 0 \\ a_{2,0} / \ell_{0,0} & (a_{2,1} - \ell_{1,0} \ell_{2,0}) / \ell_{1,1} & \sqrt{a_{2,2} - \ell_{2,0}^2 - \ell_{2,1}^2} \end{array}\right). $$ Escolhemos os elementos das diagonais como as raízes positivas; poderíamos ter escolhido as negativas já que $A$ é real.
O princípio é facilmente extensível para matrizes SPD $n\times n$, fornecendo que: $$ \begin{align*} \ell_{i,i} & = \pm\sqrt{a_{i,i} - \sum_{k=0}^{j-1} \ell_{j,k}^2},\\ \ell_{i,j} & = \frac{1}{\ell_{j,j}}\left(a_{i,j} - \sum_{k=0}^{j-1} \ell_{i,k}\ell_{j,k}\right) \qquad\text{para}~ i > j. \end{align*}
$$
Obtemos assim, o código abaixo.from math import sqrt def cholesky (A: [[float]]) -> [[float]]: n = len (A) L = [[0] for j in range (n)] for i in range (n)] for i in range (n): for j in range (i+1): d = A[i][j] - sum (L[i][k] * L[j][k] for k in range (j)) L[i][j] = sqrt (d) if i == j else d / L[j][j] return L
$\square$