Como resolver sistemas lineares no Python

Depois de falarmos sobre manipulação de matrizes e vetores, vamos continuar aqui com a série sobre álgebra linear no Python falando de sistemas lineares:

2020-9-24. Eu escrevi um texto semelhante para Julia.

Desta vez o tema é resolução de sistemas lineares, com exemplos de código em Python. Vamos achar solução para sistemas lineares (quando uma solução existe) e vamos ver o que podemos fazer quando não há solução única.

import numpy as np

Uma das utilidades de matrizes e vetores é a possibilidade de representar múltiplas relações lineares de maneira compacta. Vamos começar então com essas relações.

Definição de sistemas lineares

Um sistema linear é uma coleção de equações lineares:

a0,0x0+a0,1x1++a0,n1xn1+a0,nxn=b0a1,0x0+a1,1x1++a1,n1xn1+a1,nxn=b1am,0x0+am,1x1++am,n1xn1+am,nxn=bm\begin{aligned} a_{0,0} x_0 + a_{0,1} x_1 + \ldots + a_{0,n-1} x_{n-1} + a_{0,n} x_n &= b_0 \\ a_{1,0} x_0 + a_{1,1} x_1 + \ldots + a_{1,n-1} x_{n-1} + a_{1,n} x_n &= b_1 \\ &\:\vdots \\ a_{m,0} x_0 + a_{m,1} x_1 + \ldots + a_{m,n-1} x_{n-1} + a_{m,n} x_n &= b_m \end{aligned}

É possível escrever o sistema acima em notação matricial:

[a0,0a0,1a0,n1a0,na1,0a1,1a1,n1a1,nam,0am,1am,n1am,n][x0x1xn]=[b0b1bm]\begin{bmatrix}a_{0,0} & a_{0,1} & \ldots & a_{0,n-1} & a_{0,n}\\ a_{1,0} & a_{1,1} & \ldots & a_{1,n-1} & a_{1,n}\\ \vdots & & \ddots & & \vdots\\ a_{m,0} & a_{m,1} & \ldots & a_{m,n-1} & a_{m,n}\\ \end{bmatrix} \begin{bmatrix}x_0 \\ x_1 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix}b_0 \\ b_1 \\ \vdots \\ b_m \end{bmatrix}

Observe que a relação como escrita acima pode ser entendia de forma simplificada e abreviada como Ax=bA x = b. Em outras palavras, como xx é a incógnita, a matriz AA e o vetor bb caracterizam o sistema linear de maneira completa.

Solução de sistemas lineares: a ideia básica

A método tradicional de se resolver um sistema desses é por eliminação de Gauss, que consiste em operar linha por linha de maneira a reduzir o sistema a uma forma em que a solução é óbvia. Cada operação elementar é tal que a matriz AA e o vetor bb caracterizam um novo sistema que admite a mesma solução.

Talvez você já tenha feito isso; as três operações permitidas são:

  1. Somar kk vezes a linha ii à linha jj;

  2. Multiplicar a linha ii pelo escalar kk;

  3. Trocar a posição das linhas ii e jj.

Cada uma dessas operações é equivalente à multiplicação (pela esquerda) por uma matriz elementar. A primeira operação acima, por exemplo, que envolve multiplicar por kk a linha ii e adicionar o resultado à linha jj, consiste na multiplicação por uma matriz igual à identidade, exceto pela posição Eij=kE_{ij} = k:

E1 = np.eye(3)
E1[0, 2] = -3  # três vezes a linha dois somada à linha zero
E1
array([[ 1.,  0., -3.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

Imagine que tenhamos um sistema consistindo das equações abaixo:

x0+3x2=b03x0+2x1+x2=b1x0+x2=b2\begin{aligned} x_0 + 3 x_2 &= b_0 \\ 3 x_0 + 2 x_1 + x_2 &= b_1 \\ x_0 + x_2&= b_2 \end{aligned}

Então nossa matriz AA e o resultado da operação acima nela é:

A = np.array([[1, 0, 3],
              [3, 2, 1],
              [1, 0, 1]])
A
array([[1, 0, 3],
       [3, 2, 1],
       [1, 0, 1]])
E1 @ A
array([[-2.,  0.,  0.],
       [ 3.,  2.,  1.],
       [ 1.,  0.,  1.]])

A lógica da segunda operação é a mesma: multiplicar a linha ii por kk envolve multiplicar por uma matriz quase igual à identidade, exceto por Eii=kE_{ii} = k:

E2 = np.eye(3)
# a linha zero será multiplicada por um terço negativo
# (equivalentemente: dividida por menos três)
E2[0, 0] = -1/3
E2
array([[-0.33333333,  0.        ,  0.        ],
       [ 0.        ,  1.        ,  0.        ],
       [ 0.        ,  0.        ,  1.        ]])
E2 @ A
array([[-0.33333333,  0.        , -1.        ],
       [ 3.        ,  2.        ,  1.        ],
       [ 1.        ,  0.        ,  1.        ]])

Por último, permutar as linhas ii e jj envolve multiplicar pela matriz identidade com as respectivas linhas permutadas:

E3 = np.eye(3)
E3[[1, 2]] = E3[[2, 1]]  # troca as linhas um e dois
E3
array([[1., 0., 0.],
       [0., 0., 1.],
       [0., 1., 0.]])
E3 @ A
array([[1., 0., 3.],
       [1., 0., 1.],
       [3., 2., 1.]])

Podemos usar as operações acima em sequência para achar a matriz inversa de AA ao aplicarmos à matriz aumentada [AI][A|I]. Vamos aplicando operações de forma que obtenhamos uma matriz aumentada do tipo [IB][I|B], caso no qual BB será inversa de AA (isto é, B=A1B = A^{-1}):

AI = np.block([A, np.eye(3)])
AI
array([[1., 0., 3., 1., 0., 0.],
       [3., 2., 1., 0., 1., 0.],
       [1., 0., 1., 0., 0., 1.]])
# Subtrair a linha zero da linha dois
E1 = np.eye(3)
E1[2, 0] = -1
AI = E1 @ AI

# Subtrair três vezes a linha zero da linha um
E1 = np.eye(3)
E1[1, 0] = -3
AI = E1 @ AI

# Subtrair quarto vezes a linha dois da linha um
E1 = np.eye(3)
E1[1, 2] = -4
AI = E1 @ AI

# Somar três meios vezes a linha dois da linha zero
E1 = np.eye(3)
E1[0, 2] = 1.5
AI = E1 @ AI

AI
array([[ 1. ,  0. ,  0. , -0.5,  0. ,  1.5],
       [ 0. ,  2. ,  0. ,  1. ,  1. , -4. ],
       [ 0. ,  0. , -2. , -1. ,  0. ,  1. ]])

Resta agora dividir as linhas um e dois por dois e menos dois, respectivamente:

# Dividir a linha um por dois
E2 = np.eye(3)
E2[1, 1] = 0.5
AI = E2 @ AI

# Dividir a linha dois por menos dois
E2 = np.eye(3)
E2[2, 2] = -0.5
AI = E2 @ AI

AI
array([[ 1. ,  0. ,  0. , -0.5,  0. ,  1.5],
       [ 0. ,  1. ,  0. ,  0.5,  0.5, -2. ],
       [ 0. ,  0. ,  1. ,  0.5,  0. , -0.5]])

A inversa de AA se encontra a partir da coluna três da matriz aumentada:

Ainv = AI[:, 3:]
Ainv
array([[-0.5,  0. ,  1.5],
       [ 0.5,  0.5, -2. ],
       [ 0.5,  0. , -0.5]])

Observe que o resultado acima concorda perfeitamente com a saída da função numpy.linalg.inv, que inverte automaticamente uma matriz (quando a inversa existe):

np.linalg.inv(A)
array([[-0.5,  0. ,  1.5],
       [ 0.5,  0.5, -2. ],
       [ 0.5, -0. , -0.5]])

Por fim, o código abaixo confirma que a matriz acima encontrada é de fato inversa de AA ao comparar elemento por elemento do resultado da multiplicação A1AA^{-1} A com a matriz identidade:

np.allclose(Ainv @ A, np.eye(3))
True

Dado um sistema Ax=bA x = b, poderíamos ter realizado as mesmas operações acima à matriz aumentada [Ab][A|b]. Neste caso, teríamos obtido [Ix][I|x^*], onde xx^* é o vetor-solução desejado. A vantagem de se obter a inversa é que se pode obter a solução para qualquer bb:

x=A1bx^* = A^{-1} b
b = np.array([1., 2., 3.])
x = Ainv @ b
x
array([ 4. , -4.5, -1. ])

Ou, poderíamos ter obtido a solução da maneira mais fácil, usando a função numpy.linalg.solve:

np.linalg.solve(A, b)
array([ 4. , -4.5, -1. ])

O procedimento automático

Como visto acima, é em geral desnecessário inverter uma matriz para resolver um sistema linear. A função numpy.linalg.solve resolve sistemas lineares com solução única (aqueles em que a matriz AA é quadrada e seu determinane é diferente de zero) sem calcular a inversa. Por exemplo, o sistema Ax=bA x = b abaixo

[1234][x0x1]=[10]\begin{bmatrix}1 & 2 \\ 3 & 4\end{bmatrix} \begin{bmatrix}x_0 \\ x_1\end{bmatrix} = \begin{bmatrix}1 \\ 0\end{bmatrix}

tem solução x=(x0,x1)=(2,3/2)x = (x_0, x_1) = (-2, 3/2):

A = np.array([[1., 2.],
              [3., 4.]])
b = np.array([1., 0.])
x = np.linalg.solve(A, b)
x
array([-2. ,  1.5])

Para checar que o vetor obtido acima é de fato solução, basta checarmos a equação original, ou seja, que AxAx é igual a bb:

A @ x
array([1., 0.])

Poderíamos ter também feito a checagem diretamente no código com numpy.allclose, que, já vimos antes, compara matrizes componente a componente:

np.allclose(A @ x, b)
True

Não vou mostrar exemplos, mas vou fazer menção de que a mesma função numpy.linalg.solve também resolve sistemas lineares do tipo

AX=BA X = B

onde BB e a incógnita XX são matrizes. Isto é equivalente a resolver simultaneamente vários sistemas lineares do tipo Ax=bA x = b, onde cada xx e bb são colunas das matrizes XX e BB, respectivamente.

Então, nem a justificativa de se obter a inversa da matriz AA para resolver vários sistemas lineares com vetores bb diferentes serve: você provavelmente nunca vai precisar inverter uma matriz. Além disso, inverter matrizes custa em geral mais caro computacionalmente:

%%timeit
x = np.linalg.solve(A, b)
9.72 µs ± 78.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%%timeit
x = np.linalg.inv(A) @ b
12.4 µs ± 34.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Quadrados mínimos

Caso um sistema linear Ax=bA x = b não tenha solução, ainda assim é possível nos contentarmos com uma soluçao aproximada ótima, ou seja, uma que minimize o quadrado da norma entre os vetores AxAx e bb:

x=argminxAxb2x^* =\underset{x}{\operatorname{argmin}} |A x - b|^2

O problema acima é resolvido pela função numpy.linalg.lstsq. Observe que, quando uma solução para Ax=bAx = b existe, o mínimo xx^* acima é a própria solução:

x, _, _, _ = np.linalg.lstsq(A, b)
x
array([-2. ,  1.5])

(A função retorna quatro objetos, o código acima ignora os últimos três, uma vez que o primeiro contém a solução que desejamos.)

Um sistema sobredeterminado é exemplificado abaixo, onde a matriz AA possui mais linhas que colunas (ou seja, mais equações que incógnitas):

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

Sistemas desse tipo não possuem solução em geral. De fato, numpy.linalg.solve dará um erro do tipo LinAlgError devido à matriz AA não ser quadrada (tente fazer!). Mas podemos obter uma solução de quadrados mínimos:

b = np.array([1., 2., 2.])
x, residuals, _, _ = np.linalg.lstsq(A, b)
x
array([-0.66666667,  0.91666667])
A @ x
array([1.16666667, 1.66666667, 2.16666667])

Observe como a multiplicação de AA por xx não retorna bb (na verdade, eu escolhi bb propositalmente de maneira a não haver solução para o sistema). O segundo objeto retornado pela função numpy.linalg.lstsq nos dá a norma quadrada da diferença entre AxAx e bb:

residuals
array([0.16666667])

Podemos verificar o resíduo acima usando numpy.sum:

np.sum((A @ x - b)**2)
0.16666666666666652

É isso aí, talvez em outro momento eu fale sobre quadrados mínimos com mais profundidade, explicando como a numpy.linalg.lsqtsq faz o que faz. No entanto, no próximo texto, vamos ver alguma coisa sobre autovalores e autovetores com Python.