Como comparar matrizes maior ou igual numpy 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 = 42×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\times22×2 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 2A−I+AB2 A - I + A B2A−I+AB para

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

com III 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 uuu e vvv 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=0nui∗vi\sum_{i = 0}^n u^*_i v_i∑i=0n​ui∗​vi​). 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 é 2−2i2 - 2i2−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}A400 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 \ 78 µs ± 672 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:

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 AAA, a operação AATA A^TAAT 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 A−1A^{-1}A−1 é tal que

A−1A=AA−1=I,A^{-1} A = A A^{-1} = I \text{,}A−1A=AA−1=I,

onde III é a matriz identidade

I=[1001].I = \begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}\text{.}I=[10​01​]. 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 10−1610^{-16}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):

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 = 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⁡(A−1)=det⁡(A)−1\det(A^{-1}) = \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⁡(A−1)=1\det(A) \det(A^{-1}) = 1det(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.