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.