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$

Avatar
Aritanan Gruber
Assistant Professor

“See, if y’all haven’t the same feeling for this, I really don’t give a damn. If you ain’t feeling it, then dammit this ain’t for you!"
(desconheço a autoria; agradeço a indicação)