Como manipular matrizes e vetores no Python

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

2020-9-23. Se você estiver interessado em um conteúdo semelhante, mas para a linguagem Julia, veja isso.

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:

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/Octave, 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:

b = np.array([1., 0.])
b
array([1., 0.])

O tipo de um array é sempre numpy.ndarray:

type(b)
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:

b.dtype
dtype('float64')

Vamos verificar que este é um array unidimensional:

b.ndim
1

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

b.size
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:

b.shape
(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:

b[1]
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:

A = np.array([[1., 2.],
              [3., 4.]])
A
array([[1., 2.],
       [3., 4.]])
type(A)
numpy.ndarray

Esse é um array bidimensional:

A.ndim
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):

A[0, 1]
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:

A[-1, -1]
4.0

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

A.size
4

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

A.shape
(2, 2)

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

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

Observe que a coluna obtida é um array unidimensional:

col.ndim
1

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

col.shape
(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):

col_2d = col.reshape(2, 1)
col_2d
array([[2.],
       [4.]])

Nosso novo array possui duas dimensões:

col_2d.ndim
2

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

col_2d.shape
(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):

B = np.block([[A, np.zeros([2, 2])],
              [np.zeros([2, 2]), np.eye(2)]])
B
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×22\times2 central de B:

sub_B = B[1:3, 1:3]
sub_B
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:

sub_B.ndim
2
sub_B.shape
(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:

A
array([[1., 2.],
       [3., 4.]])
2 * A
array([[2., 4.],
       [6., 8.]])

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

A * A
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:

u = np.array([1., -1.])
v = np.array([1.,  1.])
u + v
array([2., 0.])
2. * u
array([ 2., -2.])

Multiplicação matricial

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

A @ A
array([[ 7., 10.],
       [15., 22.]])

Vamos calcular 2AI+AB2 A - I + A B para

A=[1234]B=[3140]A = \begin{bmatrix}1 & 2\\3 & 4\end{bmatrix}\quad B = \begin{bmatrix}3 & 1\\4 & 0\end{bmatrix}

com II 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):

I = np.eye(2)
I
array([[1., 0.],
       [0., 1.]])

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

A
array([[1., 2.],
       [3., 4.]])
B = np.array([[3., 1.],
              [4., 0.]])
B
array([[3., 1.],
       [4., 0.]])
2. * A - I + A @ B
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 uu e vv definidos anteriormente são ortogonais:

u
array([ 1., -1.])
v
array([1., 1.])
u @ v
0.0

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

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

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

w.conj()
array([ 1.-2.j, -1.-0.j])

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

w.conj() @ u
(2-2j)

Norma de vetores

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

np.linalg.norm(u)
1.4142135623730951

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

np.sqrt(u @ u)
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):

np.linalg.matrix_power(A, 3)
array([[ 37.,  54.],
       [ 81., 118.]])

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

A @ A @ A
array([[ 37.,  54.],
       [ 81., 118.]])

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

np.linalg.matrix_power(A, 400)
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:

%%timeit
np.linalg.matrix_power(A, 50)
15.2 µs ± 118 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%%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
78 µs ± 672 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
2020-9-23. Achou rápido? Para os mesmos casos acima, Julia é pelo menos 13×13 \times mais rápida.

Transposição de matrizes

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

A.T
array([[1., 3.],
       [2., 4.]])

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

A
array([[1., 2.],
       [3., 4.]])

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

A @ A.T
array([[ 5., 11.],
       [11., 25.]])

Inversão de matrizes

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

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

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

A1A=AA1=I,A^{-1} A = A A^{-1} = I \text{,}

onde II é a matriz identidade

I=[1001].I = \begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}\text{.}
Ai @ A
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 101610^{-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):

np.allclose(Ai @ A, np.eye(2))
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:

np.linalg.matrix_rank(A)
2
A.ndim
2

Traço de matrizes

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

A
array([[1., 2.],
       [3., 4.]])

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

np.trace(A)
5.0

Determinantes de matrizes

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

det_A = np.linalg.det(A)
det_A
-2.0000000000000004

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

det_Ai = np.linalg.det(Ai)
det_Ai
-0.49999999999999967

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

np.allclose(det_A * det_Ai, 1.)
True

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