6.13. Álgebra linear

A álgebra linear numa câmara é trabalho com matrizes pequenas: uma rotação 3x3 que combina uma amostra IMU na referência do mundo, uma matriz de calibração que retifica uma lente, a atualização da covariância de estado de um filtro de Kalman, um ajuste polinomial cujas equações normais resultam num pequeno sistema linear. Os submódulos numpy.linalg e scipy.linalg cobrem exatamente essa escala.

As funções residem em dois módulos. As operações de produto matricial estão no nível superior de numpy; as decomposições e a inversa da matriz estão em numpy.linalg; os solucionadores dedicados de sistemas lineares estão em scipy.linalg. Uma multiplicação de matrizes 2x2, por exemplo:

from ulab import numpy as np

A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
#        [43.0, 50.0]])

6.13.1. O que está disponível

  • dot() – produto matricial ou vetorial.

  • cross() – produto vetorial 3-D.

  • trace() – soma da diagonal.

  • inv() – inversa da matriz.

  • det() – determinante.

  • cholesky() – decomposição de Cholesky (entrada simétrica definida positiva).

  • eig() – valores próprios e vetores próprios de uma matriz simétrica real.

  • norm() – norma-2 de um vetor ou matriz.

  • qr() – decomposição QR com mode='reduced' (predefinição) ou mode='complete'.

  • solve_triangular() – resolve A @ x = b quando A é triangular.

  • cho_solve() – resolve A @ x = b dado um fator de Cholesky de A.

6.13.2. dot, cross, trace

dot() é a forma como a multiplicação de matrizes é escrita na câmara. O operador @ que o numpy de secretária fornece para a mesma tarefa não está implementado em ndarray, por isso cada produto matriz-vetor, matriz-matriz e produto interno de vetores passa pela chamada dot()

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.dot(a, b)             # 32.0 (scalar product)
np.cross(a, b)           # array([-3.0, 6.0, -3.0])

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a)             # matrix-vector product
np.dot(m, m)             # matrix-matrix product
np.trace(m)              # 1 + 5 + 9 = 15.0

O resultado de dot() é sempre do dtype float. cross() é o produto vetorial de dois vetores 3-D, trace() é a soma da diagonal principal de uma matriz quadrada.

6.13.3. inv e det

m = np.array([[1, 2, 3, 4],
              [4, 5, 6, 4],
              [7, 9, 9, 4],
              [3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))

A inversa é calculada por eliminação de Gauss-Jordan, pelo que inv() levanta ValueError quando a matriz é singular (uma entrada diagonal torna-se zero durante a eliminação). O custo em RAM é aproximadamente o dobro do tamanho da entrada.

O determinante reutiliza a mesma eliminação – o tempo de execução é essencialmente o mesmo que o da inversa.

Quando o objetivo da aplicação é resolver um sistema linear, não inverta e multiplique – prefira os solucionadores dedicados abaixo. Ambos são mais rápidos e melhor comportados numericamente.

6.13.4. cholesky

Para uma matriz simétrica definida positiva A, cholesky() devolve um L triangular inferior tal que A = L @ L.T

a = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]])
L = np.linalg.cholesky(a)

Se a entrada não for definida positiva ou não for simétrica, ValueError é levantado.

O fator de Cholesky representa metade do trabalho de uma fatorização LU e é o ponto de partida certo para qualquer problema cuja matriz se sabe ser simétrica definida positiva (atualizações de covariância, equações normais de um ajuste por mínimos quadrados).

6.13.5. eig

eig() funciona apenas em matrizes simétricas reais. Matrizes não simétricas levantam ValueError. Devolve um par ordenado (eigenvalues, eigenvectors)

a = np.array([[1, 2, 1, 4],
              [2, 5, 3, 5],
              [1, 3, 6, 1],
              [4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)

Notas:

  • Os valores próprios são devolvidos sem uma ordem específica. Aplique sort() (e a mesma permutação aos vetores próprios através de argsort()) quando a ordem importar.

  • Um vetor próprio é único apenas a menos de um escalar não nulo, pelo que o sinal dos vetores próprios individuais não está univocamente definido. Duas execuções corretas podem produzir vetores de sinal oposto; isso é inofensivo.

6.13.6. norm

A norma Euclidiana (de Frobenius) de um vetor ou matriz:

v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

np.linalg.norm(v)        # 7.416...
np.linalg.norm(m)        # 16.881...

A palavra-chave opcional axis= calcula a norma ao longo de um único eixo em vez de sobre todo o array.

6.13.7. qr

qr() fatora uma matriz retangular A (forma (M, N)) num Q ortonormal e num R triangular superior tais que A == Q @ R

A = np.arange(6).reshape((3, 2))

q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)

q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)

A decomposição é implementada através de rotações de Givens sucessivas. A escolha certa para problemas de mínimos quadrados em que a matriz não é simétrica.

6.13.8. Resolução de sistemas

Os dois solucionadores dedicados em ulab.scipy.linalg são ambos mais rápidos e mais precisos do que np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – resolve a @ x = b assumindo que a é triangular:

    A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

  • cho_solve(L, b)() – dado um fator de Cholesky L, resolve A @ x = b onde A = L @ L.T

    L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)

Recorra a estes em vez de inverter sempre que a estrutura de A o permitir – poupam trabalho de eliminação e a inversa explícita.

6.13.9. Resolver um pequeno sistema linear

A = np.array([[3, 0, 1, 1],
              [0, 1, 0, 2],
              [1, 0, 1, 1],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])

x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x))         # should equal b

O mesmo problema pode ser expresso fatorando A e chamando o solucionador adequado – mais rápido e mais preciso quando A tem a estrutura certa.

Para a referência completa ao nível dos argumentos, consulte numpy.linalg — Rotinas de álgebra linear e scipy.linalg — Rotinas de álgebra linear.