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$;A
eX
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 amatrix_product
,A
eX
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!