6.13. Álgebra linear

Álgebra linear em uma câmera é trabalho com matrizes pequenas: uma rotação 3x3 que funde uma amostra da IMU ao referencial do mundo, uma matriz de calibração que retifica uma lente, a atualização de estado-covariância de um filtro de Kalman, um ajuste polinomial cujas equações normais resultam em uma pequena resolução 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 de matrizes estão no nível superior de numpy; as decomposições e a inversa de matriz estão em numpy.linalg; os solucionadores dedicados de sistemas lineares estão em scipy.linalg. Uma multiplicação de matrizes 2 por 2, 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 de matrizes ou vetores.

  • cross() – produto vetorial de vetores 3-D.

  • trace() – soma da diagonal.

  • inv() – inversa de matriz.

  • det() – determinante.

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

  • eig() – autovalores e autovetores de uma matriz real simétrica.

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

  • qr() – decomposição QR com mode='reduced' (padrã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() é como a multiplicação de matrizes é escrita na câmera. O operador @ que o numpy de desktop fornece para a mesma tarefa não é implementado em ndarray, então todo produto matriz-vetor, matriz-matriz e produto escalar 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 de 3 elementos, e 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, então inv() lança ValueError quando a matriz é singular (uma entrada da diagonal se torna zero durante a eliminação). O custo de RAM é aproximadamente o dobro do tamanho da entrada.

O determinante reutiliza a mesma eliminação – o tempo de execução é essencialmente o mesmo 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 numericamente mais bem-comportados.

6.13.4. cholesky

Para uma matriz simétrica positiva definida A, cholesky() retorna uma matriz triangular inferior L 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 positiva definida ou não for simétrica, ValueError é lançada.

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

6.13.5. eig

eig() funciona apenas em matrizes reais simétricas. Matrizes não simétricas lançam ValueError. Ela retorna uma 2-tupla (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 autovalores retornam sem ordem específica. Aplique sort() (e a mesma permutação aos autovetores via argsort()) quando uma ordem classificada for importante.

  • Um autovetor é único apenas a menos de um escalar não nulo, então o sinal de autovetores individuais não é definido de forma única. 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 o array inteiro.

6.13.7. qr

qr() fatora uma matriz retangular A (formato (M, N)) em uma Q ortonormal e uma 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 por meio de sucessivas rotações de Givens. É a escolha certa para problemas de mínimos quadrados em que a matriz não é simétrica.

6.13.8. Resolvendo sistemas

Os dois solucionadores dedicados em ulab.scipy.linalg são ao mesmo tempo 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 supondo que a seja 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 permitir – eles economizam tanto o trabalho de eliminação quanto a inversa explícita.

6.13.9. Resolvendo 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 apropriado – mais rápido e mais preciso quando A tem a estrutura adequada.

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