BCM0505-15 - Processamento da Informação

Aula A06 – 15/05: Matrizes

• Definição, soma, multiplicação por escalar, transposição e teste de simetria.
• Produtos simples, de Hadamard e de Krönecker.
• Potenciação e exponenciação.
• Matrizes em NumPy.

Definição

Formalmente, uma matriz $A$ de números reais é uma função $A:\mathcal{L}\times\mathcal{C}\to\mathbb{R}$, em que $\mathcal{L}$ e $\mathcal{C}$ são dois conjuntos discretos (e no nosso caso, finitos) de índices; o primeiro indexa “linhas” e o segundo, “colunas.” De outra forma, uma matriz real é uma função de duas variáveis discretas a um valor real.

É comum considerarmos $\mathcal{L}=\{0,1,\ldots,m-1\}$ e $\mathcal{C}=\{0,1,\ldots,n-1\}$ e dizermos que a matriz tem $m$ linhas e $n$ colunas. Aproveitamos e simplificamos a notação $$A:\{0,1,\ldots,m-1\}\times\{0,1,\ldots,n-1\}\to\mathbb{R} \qquad\text{para}\qquad A\in\mathbb{R}^{m\times n}.$$ A estrutura bidimensional de $A$ permite que ela seja arrumada graficamente de forma tabular, como ilustra o exemplo que segue. $$A=\left(\begin{array}{rrrr} -1 & 0 & 0 & 3\\ 2 & 5 & 8 & 0\\ 3 & 0 & 1 & -4 \end{array}\right)$$ Neste caso, $A$ tem $3$ linhas e $4$ colunas.

Uma matriz pode ser representada em Python como listas de listas: a lista externa armazena as linhas e cada lista interna, as colunas. Assim, a matriz $A$ do exemplo acima pode ser declarada em Python como

A = [
      [-1, 0, 0, 3], 
      [ 2, 5, 8, 0], 
      [ 3, 0, 1,-4]
    ]

ou, de forma mais sucinta, como

A = [[-1,0,0,3], [2,5,8,0], [3,0,1,-4]]

O elemento $A_{ij}$ é acessado via A[i][j] Exemplo: $A_{12} = 8 =\,$A[1][2].
Nesta representação, o número de linhas de $A$ é dado por $m=\,$len (A) e o de colunas, por $n=\,$len (A[0]).

Observe que poderíamos ter escolhido uma forma transposta de representação, ou ainda qualquer outra que fizesse sentido e que Python não oferece nenhum mecanismo especial de proteção neste caso: é nossa responsabilidade utilizar as “listas de listas” de forma correta para representar matrizes. Existem variações. Mencionaremos uma delas ao final desta aula.

Nota: Anotações de tipos para matrizes levam em consideração o fato de serem representadas como listas de listas. Logo, uma matriz de pontos flutuantes é anotada por [[float]].

---

Operações Elementares

Problema A6.1: Dada uma matriz quadrada $A\in\mathbb{R}^{n\times n}$, verifique se $A$ é a matriz identidade, isto é, se $$A[i][j] = \begin{cases} 1 & \text{se}~i=j,\\ 0 & \text{em caso contrário,} \end{cases}$$ para todo $0\leq i,j<n$.

Exemplo de matrizes identidades $1\times 1$, $2\times 2$ e $3\times 3$: $$ I_1 = (1), \qquad I_2 = \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right), \qquad I_3 = \left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right).$$

def matrix_is_identity (A: [[float]]) -> bool:
  n = len (A)
  for i in range (n):
    for j in range (n):
      if A[i][j] != (1 if i == j else 0):
        return False
  return True

Problema A6.2: Dada uma matriz quadrada $A\in\mathbb{R}^{n\times n}$, verifique se $A$ é simétrica, isto é, se $$A[i][j] = A[j][i] \qquad\text{para todo}\qquad 0\leq i,j<n.$$

Exemplo de uma matriz $3\times 3$ simétrica: $$A = \left(\begin{array}{ccc} 1 & 0 & 4 \\ 0 & 5 & 2 \\ 4 & 2 & 3 \end{array}\right)$$

def matrix_is_symmetric (A: [[float]]) -> bool:
  n = len (A)
  for i in range (n):
    for j in range (i):
      if A[i][j] != A[j][i]:
        return False
  return True

Exercício 1: Por que é suficiente que o laço $j$ esteja indexado de $0$ a $i-1$? Nota: a necessidade é óbvia.

Problema A6.3: Dadas duas matrizes $A,B\in\mathbb{R}^{m\times n}$, devolva a soma $C=A+B$, definida como $$C[i][j] = A[i][j] + B[i][j], \qquad\text{para todo}\qquad 0\leq i< m~\text{e}~0\leq j< n.$$

Neste ponto é instrutivo tentar definir a matriz $C$ como

C = [[0] * n] * m

que seria uma extensão natural da forma como definimos vetores. No entanto, a forma como o mecanismo de referências a objetos em Python funciona produz algo indesejável:

>>> m, n = 3, 4
>>> C = [[0] * n] * m
>>> print (C)
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
>>> C[0][1] = 2
>>> print (C)
[[0, 2, 0, 0], [0, 2, 0, 0], [0, 2, 0, 0]]

Observe que o comando C[0][1] = 2 modificou toda a coluna $1$, e não apenas a entrada/célula $(0,1)$. Isto ocorre porque [0] * n cria um vetor com $n$ entradas, digamos $X$, e [X] * m cria um outro vetor com $m$ entradas em que cada entrada aponta para o mesmo vetor $X$. O comportamento acima, do ponto de vista dos projetistas da linguagem, é correto; nós é que queremos algo diferente! Felizmente, existem várias formas de conseguirmos o que queremos.

A primeira variante consiste em utilizar a função de listas append:

>>> m, n = 3, 4
>>> C = []
>>> for i in range (m):
...   C.append ([0] * n)
>>> print (C)
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
>>> C[0][1] = 2
>>> print (C)
[[0, 2, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]

Agora, cada uma das $m$ linhas de $C$ é um vetor distinto com $n$ entradas (colunas). Um código para soma de matrizes com esta variação:

def matrix_sum (A: [[float]], B:[[float]]) -> [[float]]:
  m, n = len (A), len (A[0])

  # Cria matriz C cheia de zeros.
  C = []
  for i in range (m):
    C.append ([0] * n)

  # Realiza a soma A + B
  for i in range (m):
    for j in range (n):
      C[i][j] = A[i][j] + B[i][j]
  return C

A segunda variante consiste em utilizar list comprehensions:

>>> m, n = 3, 4
>>> C = [[0 for j in range (n)] for i in range (m)]
>>> print (C)
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
>>> C[0][1] = 2
>>> print (C)
[[0, 2, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]

Produz o mesmo efeito, mas é mais curta e possibilita uma expansão bacana: em vez de construir a lista $C$ e depois percorrê-la para preencher cada entrada individualmente, podemos fazer as duas tarefas de uma única vez:

def matrix_sum (A: [[float]], B:[[float]]) -> [[float]]:
  return [[A[i][j] + B[i][j] for j in range (len (A[0]))] for i in range (len (A))]

Claramente, a soma de duas matrizes com $m$ linhas e $n$ colunas consome tempo $O(mn)=O(n^2)$, quando $m=n$.

Problema A6.4: A transposta de uma matriz $A\in\mathbb{R}^{m\times n}$ é uma matriz $A^\top\in\mathbb{R}^{n\times m}$ tal que $$A[i][j] = A^\top[j][i] \qquad\text{para todo}\qquad 0\leq i,j<n.$$ Dada $A$, determine $A^\top$.

Exemplo da transposta de uma matriz $3\times 4$: $$A = \left(\begin{array}{cccc} 1 & 0 & 4 & 0\\ 0 & 5 & 2 & 1\\ 6 & 1 & 3 & 2\end{array}\right), \qquad A^\top = \left(\begin{array}{ccc} 1 & 0 & 6\\ 0 & 5 & 1\\ 4 & 2 & 3\\ 0 & 1 & 2\end{array}\right).$$

def matrix_transpose (A: [[float]]) -> [[float]]:
  return [[A[i][j] for i in range (len (A))] for j in range (len (A[0]))]

Caso $A$ seja quadrada ($m=n$), a transposição pode ser realizada na própria matriz.

def matrix_sq_transpose (A: [[float]]):
  n = len (A)
  for i in range (n-1):
    for j in range (i+1, n):
      A[i][j], A[j][i] = A[j][i], A[i][j]

Problema A6.5: Dado um inteiro $n>0$, devolva a matriz identidade $I_n\in\{0,1\}^{n\times n}$.

def matrix_identity (n: int) -> [[float]]:
  return [[(1. if i == j else 0.) for j in range (n)] for i in range (n)]

---

Produtos

Problema A6.6: Dadas uma matriz $A\in\mathbb{R}^{m\times n}$ e uma constante $c\in\mathbb{R}$, devolva o produto $c\cdot A$ na própria matriz $A$. Isto é, multiplique todas as entradas de $A$ por $c$.

def matrix_scalar (A: [[float]], c: float):
  for i in range (len (A)):
    for j in range (len (A[0])):
      A[i][j] *= c

O código acima funciona porque: (1) $A$ é uma referência a uma lista, e (2) listas são elementos mutáveis em Python.

De qualquer forma, para incrementar a usabilidade, é conveniente também devolver $A$ após a mudança.

def matrix_scalar (A: [[float]], c: float) -> [[float]]:
  for i in range (len (A)):
    for j in range (len (A[0])):
      A[i][j] *= c
  return A

Problema A6.7: Dadas duas matrizes $A,B\in\mathbb{R}^{m\times n}$ de mesmas dimensões, o produto de Hadamard $A\odot B$ entre $A$ e $B$ é dado por $$(A\odot B)[i][j] = A[i][j] \cdot B[i][j],$$ para todo $0\leq i< m$ e $0\leq j< n$.

Exemplo: $$\left(\begin{array}{ccc} 1 & 0 & 2\\ 3 & 1 & 2\end{array}\right) \odot \left(\begin{array}{ccc} 4 & 1 & 3\\ 5 & 1 & 0\end{array}\right) = \left(\begin{array}{ccc} 4 & 0 & 6\\ 15 & 1 & 0\end{array}\right)$$

def matrix_hadamard (A: [[float]], B:[[float]]) -> [[float]]:
  return [[A[i][j] * B[i][j] for j in range (len (A[0]))] for i in range (len (A))]

Problema A6.8: Dadas duas matrizes $A\in\mathbb{R}^{m\times n}$ e $B\in\mathbb{R}^{n\times\ell}$, devolva o produto (regular) $C=A\cdot B$, dado por $$C[i][j] = \sum_{k=0}^{n-1} A[i][k] \cdot B[k][j], \qquad\text{para todo}\qquad 0\leq i< m~\text{e}~0\leq j< \ell. \qquad (*)$$ $C$ é uma matriz com $m$ linhas e $\ell$ colunas.

Exemplo: $$\left(\begin{array}{ccc} 1 & 0 & 2\\ 3 & 1 & 2\end{array}\right) \cdot \left(\begin{array}{cc} 4 & 5\\ 1 & 1\\ 3 & 0 \end{array}\right) = \left(\begin{array}{cc} 10 & 5\\ 19 & 16 \end{array}\right)$$

O algoritmo fornecido em $(*)$ consome tempo $O(mn\ell)=O(n^3)$, quando $\ell=m=n$.

def matrix_product (A: [[float]], B:[[float]]) -> [[float]]:
  m, n, l = len (A), len (A[0]), len (B[0])

  C = [[0 for j in range (l)] for i in range (m)]
  for i in range (m):
    for j in range (l):
      for k in range (n):
        C[i][j] += A[i][k] * B[k][j]
  return C

Se quisermos garantir que as dimensões de $A$ e $B$ são compatíveis, podemos inserir um assert (n == len (B)) após a definição da variável $n$.

A título de curiosidade (e como um bom exercício), podemos escrever uma versão utilizando list comprehensions:

def matrix_product (A: [[float]], B:[[float]]) -> [[float]]:
  return [[sum (a*b for a, b in zip (A[i], matrix_transpose (B)[j])) for j in range (len (B))] for i in range (len (A))]

Observe o uso da função matrix_transpose que transpõe $B$. Logo, além do código ser bem menos legível, esta versão é mais lenta.
De qualquer forma, foi um bom exercício! Recomendo que você entenda direito o que está sendo feito na linha do return acima.

Problema A6.9: Dados dois vetores $x[0\,..m-1]$ e $y[0\,..n-1]$ com entradas em $\mathbb{R}$, o produto externo de $x$ e $y$ é dado por $$xy^\top = \left(\begin{array}{cccc} x_0y_0 & x_0y_1 & \cdots & x_0y_{n-1}\\ x_1y_0 & x_1y_1 & \cdots & x_1y_{n-1}\\ \vdots & \vdots & \ddots & \vdots\\ x_{m-1}y_0 & x_{m-1}y_1 & \cdots & x_{m-1}y_{n-1}\\ \end{array}\right).$$ Determine $xy^\top$.

Observe que o produto externo $xy^\top$ é uma matriz, ao passo que o produto interno $x^\top y$ é um escalar (claramente, este só funciona se $m=n$).

def external_product (x: [float], y:[float]) -> [[float]]:
  return [[x[i] * y[j] for j in range (len (y))] for i in range (len (x))]

Problema A6.10: Para duas matrizes $A\in\mathbb{R}^{k\times \ell}$ e $B\in\mathbb{R}^{m\times n}$, o produto de Krönecker ou tensorial entre $A$ e $B$, denotado por $A\otimes B$, é uma matriz de blocos que possui dimensões $km\times\ell n$ e é dada por: $$A\otimes B = \left(\begin{array}{ccc} a_{0,0}B & \cdots & a_{0,\ell-1}B\\ \vdots & \ddots & \vdots\\ a_{k-1,0}B & \cdots & a_{k-1,\ell-1}B \end{array}\right), \qquad\text{em que}\qquad a_{i,j}B = \left(\begin{array}{ccc} a_{i,j}b_{0,0} & \cdots & a_{i,j}b_{0,n-1}\\ \vdots & \ddots & \vdots\\ a_{i,j}b_{m-1,0} & \cdots & a_{i,j}b_{m-1,n-1} \end{array}\right).$$ Escreva uma função que recebe $A$ e $B$ e devolve $A\otimes B$.

Por exemplo, se $$A = \left(\begin{array}{cc} 1 & 2\\ 3 & 4 \end{array}\right) \quad\text{e}\quad B = \left(\begin{array}{cc} 0 & 5\\ 6 & 7 \end{array}\right), \qquad A\otimes B = \left(\begin{array}{cc} 1\left(\begin{array}{cc} 0 & 5\\ 6 & 7 \end{array}\right) & 2\left(\begin{array}{cc} 0 & 5\\ 6 & 7 \end{array}\right)\\ 3\left(\begin{array}{cc} 0 & 5\\ 6 & 7 \end{array}\right) & 4\left(\begin{array}{cc} 0 & 5\\ 6 & 7 \end{array}\right)\\ \end{array}\right) = \left(\begin{array}{cccc} 0 & 5 & 0 & 10\\ 6 & 7 & 12 & 14\\ 0 & 15 & 0 & 20\\ 18 & 21 & 24 & 28 \end{array}\right).$$

Não é difícil perceber que o produto de Krönecker é uma generalização para matrizes do produto externo de vetores.

Exercício 2: Escreva o código da função

def matrix_kronecker (A: [[float]], B:[[float]]) -> [[float]]:

que implementa o produto de Krönecker entre $A$ e $B$.

---

Potenciação e Exponenciação

Problema A6.11: Dada uma matriz quadrada $A\in\mathbb{R}^{n\times n}$ e um inteiro $k\geq 0$, determine $$A^k = \prod_{i=1}^k A,$$ em que $A^0 = I_n$, a matriz identidade $n\times n$.

Como temos uma solução para o Problema A6.8 (produto regular de duas matrizes), escrevemos rapidamente o código abaixo.

def matrix_power (A: [[float]], k: int) -> [[float]]:
  assert (k >= 0)
  X = [[(1 if i == j else 0) for j in range (len (A[0]))] for i in range (len (A))]
  while k > 0:
    X = matrix_product (A, X)
    k -= 1
  return X

Apesar de funcional, não é difícil perceber que ele realiza a multiplicação de $A$ pot $I_n$ quando $k=1$, uma operação de custo $O(n^3)$ — inaceitável! Felizmente, isso é fácil de se resolver.

def matrix_power (A: [[float]], k: int) -> [[float]]:
  assert (k >= 0)
  if   k == 0: return matrix_identity (len (A))
  elif k == 1: return A
  else:   
    X = A # <<<------- Cuidado!
    while k > 1:
      X = matrix_product (A, X)
      k -= 1
  return X

Qual a razão do comentário “Cuidado!” ao lado da atribuição X = A?
Nossas matrizes são listas de listas (objetos em Python). A variável A é uma referência para uma matriz. A atribuição X = A atribui a referência A a X. Ou seja, A e X passam a apontar para o mesmo objeto na memória; nenhuma cópia (clonagem) de matriz é feita. O código funciona porque:

• (1) a primeira chamada a matrix_product (A, X) é feita para calcularmos $A^2$; Ae X apontando para a mesma matriz produz o efeito desejado;
• (2) e mais importante, a função matrix_product cria e devolve uma nova matriz com o produto $A\cdot X$; logo, após a primeira chamada a matrix_product, A e X passam a apontar para matrizes diferentes.

Caso A e X apontassem sempre para a mesma matriz, o resultado de matriz_power seria $A^{2^k}$ e não $A^k$, como desejado. Colocamos o comentário “Cuidado!” ali para introduzir esta discussão e para mencionar que é preciso fazer uma análise semelhante à que fizemos sempre que for utilizar tal (tipo de) construção.

O código acima realiza $k-2$ produtos de matrizes, levando a um consumo de tempo $O(kn^3)$. Podemos reaproveitar algumas potências já feitas e reduzir $k$ para $\log k$.

def matrix_power (A: [[float]], k: int) -> [[float]]:
  assert (k >= 0)
  if   k == 0: return matrix_identity (len (A))
  elif k == 1: return A
  else:   
    X = matrix_power (A, k//2)
    if k % 2 == 0:
      return matrix_product (X, X)
    else:
      return matrix_product (A, matrix_product (X, X))

Tempo: $O(n^3\log k)$.

Problema A6.12: Escreva uma outra função que recebe uma matriz $A\in\mathbb{R}^{n\times n}$ e um inteiro $m > 1$ e devolve uma aproximação da exponencial de $A$, denotada por $e^A$ e definida como $$ e^A = I + A + \frac{A^2}{2} + \frac{A^3}{3!} + \cdots + \frac{A^k}{k!} + \cdots $$ que inclua os $m$ primeiros termos da série acima.

Naturalmente, a análise que fizemos da série de Taylor-Maclaurin para a função $e^x$ generaliza para a exponencial da matriz $A$.
Em outras palavras, de posse do $(k-1)$-ésimo termo da série, basta multiplicá-lo por $$\frac{A^k}{k!}\bigg/\frac{A^{k-1}}{(k-1)!} = \frac{1}{k}\cdot A$$ para gerar o $k$-ésimo termo.

def matrix_exp (A: [[float]], m: int) -> [[float]]:
  n, X = len (A), A
  eA = matrix_sum (matrix_identity (n), A)
  for k in range (2, m):
    X  = matrix_scalar (matrix_product (X, A), 1/k)
    eA = matrix_sum (eA, X)
  return eA   

---

Matrizes em NumPy

NumPy é uma biblioteca Python para computação científica (numérica). Um dos objetivos é oferecer os principais algoritmos básicos de álgebra linear e cálculo numérico a uma boa performance — os algoritmos são escritos em C e estão disponíveis para chamada em Python via NumPy.

NumPy disponibiliza a classe ndarray — que também pode ser referenciada por array — para representar vetores, matrizes e tensores.

>>> import numpy as np
>>> c = np.array ([0,1,2,3,4])
>>> c
array([0, 1, 2, 3, 4])
>>> A = np.array ([[1,2],[3,4]])
>>> A
array([[1, 2],
       [2, 4]])
>>> B = np.array ([[1,1],
...                [0,1]])       
>>> A + B # soma
array([[2, 3],
       [2, 5]])
>>> A * B # produto de Hadamard
array([[1, 2],
       [0, 4]])
>>> A @ B # produto regular, semelhante a A.dot(B)
array([[3, 6],
       [2, 4]])

Não temos a intenção de realizar um tratamento de NumPy neste curso; para tal, recomendamos o tutorial referido no link acima. De qualquer forma, observe que as declarações são feitas vias listas (de listas) da forma como fizemos, e que várias funcionalidades são encapsuladas na classe. NumPy é uma bela biblioteca e seu uso é recomendado em situações que requeiram um uso eficiente de vetores, matrizes ou tensores em Python — basicamente, aprenda-a e utilize-a depois que você aprender a programas direito em Python!