Manipulando matrizes e vetores no Python

Este é o início de uma mini série sobre sobre álgebra linear com Python. Vamos falar sobre matrizes, vetores e como manipulá-los com a biblioteca NumPy.

Este texto fala da construção de vetores e matrizes, suas dimensões, operações e algumas funções matriciais, com exemplos de código em Python.

O pacote que vamos precisar é o numpy. Esse pacote é a base para operações matemáticas em Python pois, além de prover um objeto para representação de tensores de dimensão arbitrária chamado array (o que inclui matrizes e vetores), também fornece implementações eficientes de uma série de operações e funções.

Vamos importar esse módulo e adotar a convenção de chamá-lo no código de np:

In [1]:
import numpy as np

Uma série de operações importantes para álgebra linear se encontram no submódulo numpy.linalg, que já é importado por padrão com a linha acima.

Arrays

Representamos matrizes e vetores através do objeto array. Podemos pensar em arrays unidimensionais como listas de números e arrays bidimensionais como matrizes. É possível construir arrays de dimensões arbitrárias (representando tensores de dimensões três, quatro, etc.), mas não vamos falar disso aqui.

Quando selecionamos uma linha ou coluna de um array bidimensional, o resultado é um array unidimensional (isso difere do MATLAB, onde o resultado seria também bidimensional). Por isso é importante ficarmos atentos nos tamanhos e dimensões, já que pode ficar por vezes ambíguo quando estamos considerando arrays unidimensionais como vetores linha ou vetores coluna.

Vamos criar um array unidimensional, a partir de uma lista de elementos, e determinar alguns dos seus parâmetros:

In [2]:
b = np.array([1., 0.])
b
Out[2]:
array([1., 0.])

O tipo de um array é sempre numpy.ndarray:

In [3]:
type(b)
Out[3]:
numpy.ndarray

Já o tipo dos seus elementos é armazenado na propriedade dtype. No caso os elementos de b são de números de ponto flutuante:

In [4]:
b.dtype
Out[4]:
dtype('float64')

Vamos verificar que este é um array unidimensional:

In [5]:
b.ndim
Out[5]:
1

O número total de elementos do array é dois:

In [6]:
b.size
Out[6]:
2

O formato do array é a disposição desses elementos nas dimensões. No array b, os dois elementos estão obviamente dispostos na única dimensão:

In [7]:
b.shape
Out[7]:
(2,)

Os elementos de um array unidimensional são indexados por suas posições (começando em zero). O segundo elemento de b, portanto, é zero:

In [8]:
b[1]
Out[8]:
0.0

Vamos agora criar uma matriz (array bidimensional). Matrizes são definidas linha por linha de maneira similar aos arrays unidimensionais, com o fornecimento de uma lista de listas de elementos:

In [9]:
A = np.array([[1., 2.], [3., 4.]])
A
Out[9]:
array([[1., 2.],
       [3., 4.]])
In [10]:
type(A)
Out[10]:
numpy.ndarray

Esse é um array bidimensional:

In [11]:
A.ndim
Out[11]:
2

A primeira dimensão corresponde à direção vertical (contagem de linhas) e a segunda corresponde à direção horizontal (contagem de colunas). Por exemplo, o segundo elemento da primeira linha é (Python usa contagens começando em zero):

In [12]:
A[0, 1]
Out[12]:
2.0

Observe que a indexação por números negativos significa contar do final daquela direção. O último elemento da última coluna, portanto, é quatro:

In [13]:
A[-1, -1]
Out[13]:
4.0

O número total de elementos em A é quatro:

In [14]:
A.size
Out[14]:
4

Esses elementos estão dispostos dois a dois nas duas dimensões (ou seja, $2 \times 2 = 4$ elementos):

In [15]:
A.shape
Out[15]:
(2, 2)

Podemos fatiar a segunda coluna inteira da matriz A (ou seja, todos os índices de uma coluna):

In [16]:
col = A[:, 1]
col
Out[16]:
array([2., 4.])

Observe que a coluna obtida é um array unidimensional:

In [17]:
col.ndim
Out[17]:
1

Seu formato portanto é (2,) (ambos elementos na única dimensão):

In [18]:
col.shape
Out[18]:
(2,)

No entanto, é possível transformar esse array unidimensional em bidimensional através do método reshape, que recebe o formato desejado do array. No nosso caso, para obtermos um vetor coluna, desejamos um formato (2, 1):

In [19]:
col_2d = col.reshape(2, 1)
col_2d
Out[19]:
array([[2.],
       [4.]])

Nosso novo array possui duas dimensões:

In [20]:
col_2d.ndim
Out[20]:
2

O formato é (2, 1) como requisitado pelo método reshape:

In [21]:
col_2d.shape
Out[21]:
(2, 1)

Ambos objetos col e col_2d guardam as mesmas informações, mas são fundamentalmente diferentes.

Construção de matrizes em bloco

É comum precisar contruir matrizes bloco a bloco. Isto é obtido com a função numpy.block (abaixo usamos também numpy.eye, para construir uma matriz identidade, e numpy.zeros, para construir uma matriz de zeros, em ambos os casos da dimensões desejadas):

In [22]:
B = np.block([[A, np.zeros([2, 2])],
              [np.zeros([2, 2]), np.eye(2)]])
B
Out[22]:
array([[1., 2., 0., 0.],
       [3., 4., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

A função numpy.block facilita a construção de matrizes grandes e com muita complexidade.

A operação inversa, ou seja, se fatiar blocos de arrays, pode ser feita com : e índices de início (incluso) e fim (não incluso). Por exemplo, para se fatiar a submatriz $2\times2$ central de B:

In [23]:
subB = B[1:3, 1:3]
subB
Out[23]:
array([[4., 0.],
       [0., 1.]])

Os índices incluídos são [1, 3), o que não inclui o três. Assim como iniciar a contagem por zero, não incluir o índice final é padrão em Python.

O resultado é um array menor, mas de mesma dimensão:

In [24]:
subB.ndim
Out[24]:
2
In [25]:
subB.shape
Out[25]:
(2, 2)

Operações com matrizes e vetores

Operações aritméticas de matrizes e vetores

As operações de adição (+), subtração (-), divisão (/), multiplicação (*) e exponenciação (**) podem ser realizadas entre arrays e escalares:

In [26]:
A
Out[26]:
array([[1., 2.],
       [3., 4.]])
In [27]:
2 * A
Out[27]:
array([[2., 4.],
       [6., 8.]])

Quando realizadas entre arrays, essas operações são realizadas elemento a elemento:

In [28]:
A * A
Out[28]:
array([[ 1.,  4.],
       [ 9., 16.]])

Podemos também manipular arrays unidimensionais como esperado para vetores, ou seja, em termos de soma de vetores e multiplicação por escalar:

In [29]:
u = np.array([1., -1.])
v = np.array([1.,  1.])
u + v
Out[29]:
array([2., 0.])
In [30]:
2. * u
Out[30]:
array([ 2., -2.])

Multiplicação matricial

A multiplicação matricial possui um operador próprio (@) em Python:

In [31]:
A @ A
Out[31]:
array([[ 7., 10.],
       [15., 22.]])

Vamos calcular $2 A - I + A B$ para

$$A = \begin{bmatrix}1 & 2\\3 & 4\end{bmatrix}\quad B = \begin{bmatrix}3 & 1\\4 & 0\end{bmatrix}$$

com $I$ sendo a matriz identidade de dimensão dois (como vimos antes na parte de matrizes bloco, as matrizes identidade podem ser construídas com numpy.eye):

In [32]:
I = np.eye(2)
I
Out[32]:
array([[1., 0.],
       [0., 1.]])

A matriz A já foi definida anteriormente, basta definirmos a matriz B:

In [33]:
A
Out[33]:
array([[1., 2.],
       [3., 4.]])
In [34]:
B = np.array([[3., 1.], [4., 0.]])
B
Out[34]:
array([[3., 1.],
       [4., 0.]])
In [35]:
2. * A - I + A @ B
Out[35]:
array([[12.,  5.],
       [31., 10.]])

Produto interno entre vetores

O produto interno entre vetores é obtido com o operador de multiplicação matricial entre arrays unidimensionais (@). Ele nos revela que os os $u$ e $v$ definidos anteriormente são ortogonais:

In [36]:
u
Out[36]:
array([ 1., -1.])
In [37]:
v
Out[37]:
array([1., 1.])
In [38]:
u @ v
Out[38]:
0.0

O produto interno com vetores complexos é levemente diferente, uma vez que se usa o complexo conjugado do vetor à esquerda ($\sum_{i = 0}^n u^*_i v_i$). Primeiro, vamos definir um vetor complexo (a unidade complexa em Python é 1j):

In [39]:
w = np.array([1. + 2.j, -1.])
w
Out[39]:
array([ 1.+2.j, -1.+0.j])

Podemos objetor o complexo conjugado de um array unidimensional com o método conj, que basicamente inverte o sinal da parte complexa de todos os elementos:

In [40]:
w.conj()
Out[40]:
array([ 1.-2.j, -1.-0.j])

O produto interno de w com u é $2 - 2i$:

In [41]:
w.conj() @ u
Out[41]:
(2-2j)

Norma de vetores

A norma euclidiana pode ser calculada com uma função numpy.linalg.norm (presente no submódulo numpy.linalg):

In [42]:
np.linalg.norm(u)
Out[42]:
1.4142135623730951

Observe que a norma euclidiana é a raiz quadrada (numpy.sqrt) do produto interno de um vetor com ele próprio:

In [43]:
np.sqrt(u @ u)
Out[43]:
1.4142135623730951

Potenciação matricial

Não existe um símbolo para potenciação matricial (lembre-se: ** significa potenciação elemento a elemento), mas existe uma função para esta operação (numpy.linalg.matrix_power):

In [44]:
u
Out[44]:
array([ 1., -1.])
In [45]:
np.linalg.matrix_power(A, 3)
Out[45]:
array([[ 37.,  54.],
       [ 81., 118.]])

Naturalmente, a potenciação é equivalente a multiplicações sequenciais:

In [46]:
A @ A @ A
Out[46]:
array([[ 37.,  54.],
       [ 81., 118.]])

No entanto, não tente fazer $A^{400}$ multiplicação por multiplicação!

In [47]:
np.linalg.matrix_power(A, 400)
Out[47]:
array([[2.76493466e+291, 4.02969073e+291],
       [6.04453610e+291, 8.80947076e+291]])

Aliás, internamente a potenciação não é feita multiplicação por multiplicação (em um loop, por exemplo) por outro motivo que não a conveniência: eficiência. Para entender melhor, compare os tempos de execução abaixo:

In [48]:
%%timeit
np.linalg.matrix_power(A, 50)
15.9 µs ± 123 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [49]:
%%timeit
A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A \
@ A @ A @ A @ A @ A @ A @ A @ A @ A @ A
77.2 µs ± 435 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Transposição de matrizes

Podemos transpor uma matriz com a propriedade .T dos arrays:

In [50]:
A.T
Out[50]:
array([[1., 3.],
       [2., 4.]])

Observe que os números fora da diagonal trocaram de posição com relação à matriz original:

In [51]:
A
Out[51]:
array([[1., 2.],
       [3., 4.]])

É interessante observar também que, para qualquer matriz $A$, a operação $A A^T$ produz uma matriz simétrica:

In [52]:
A @ A.T
Out[52]:
array([[ 5., 11.],
       [11., 25.]])

Inversão de matrizes

Matrizes são invertidas com numpy.linalg.inv:

In [53]:
Ai = np.linalg.inv(A)
Ai
Out[53]:
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

A matriz inversa $A^{-1}$ é tal que

$$A^{-1} A = A A^{-1} = I \text{,}$$

onde $I$ é a matriz identidade

$$I = \begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}\text{.}$$
In [54]:
Ai @ A
Out[54]:
array([[1.00000000e+00, 0.00000000e+00],
       [1.11022302e-16, 1.00000000e+00]])

A menos de erros numéricos, essa parece ser a resposta correta (valores da ordem de $10^{-16}$ são praticamente zero). Vamos checar isso por comparação da multiplicação acima com a matriz identidade (usando a função numpy.allclose, que compara arrays a menos de erros numéricos):

In [55]:
np.allclose(Ai @ A, np.eye(2))
Out[55]:
True

Lembre-se: nem toda matriz possui inversa! Uma maneira de identificar se uma matriz possui inversa é com o cálculo do rank da mesma (numpy.linalg.matrix_rank), que deve ser igual à dimensão do array:

In [56]:
np.linalg.matrix_rank(A)
Out[56]:
2
In [57]:
A.ndim
Out[57]:
2

Traço de matrizes

O traço de uma matriz é a soma de seus valores diagonais. Ele pode ser calculado com numpy.trace:

In [58]:
A
Out[58]:
array([[1., 2.],
       [3., 4.]])

Para o caso da matriz acima, o traço deve ser $1 + 4 = 5$:

In [59]:
np.trace(A)
Out[59]:
5.0

Determinantes de matrizes

Determinantes são calculados com numpy.linalg.det:

In [60]:
det_A = np.linalg.det(A)
det_A
Out[60]:
-2.0000000000000004

A matriz inversa tem a propriedade de que $\det(A^{-1}) = \det(A)^{-1}$:

In [61]:
det_Ai = np.linalg.det(Ai)
det_Ai
Out[61]:
-0.49999999999999967

De fato, podemos verificar $\det(A) \det(A^{-1}) = 1$ a menos de erros numéricos com numpy.allclose:

In [62]:
np.allclose(det_A * det_Ai, 1.)
Out[62]:
True

É isso aí, no próximo texto vamos falar um pouco da resolução de sistemas de equações lineares.

Deixe uma resposta