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.