Resolvendo 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.

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.

In [1]:
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:

$$\begin{align} 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{align}$$

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

$$\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 $A x = b$. Em outras palavras, como $x$ é a incógnita, a matriz $A$ e o vetor $b$ 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 $A$ e o vetor $b$ 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 $k$ vezes a linha $i$ à linha $j$;
  2. Multiplicar a linha $i$ pelo escalar $k$;
  3. Trocar a posição das linhas $i$ e $j$.

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 $k$ a linha $i$ e adicionar o resultado à linha $j$, consiste na multiplicação por uma matriz igual à identidade, exceto pela posição $E_{ij} = k$:

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

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

$$\begin{align} x_0 + 3 x_2 &= b_0 \\ 3 x_0 + 2 x_1 + x_2 &= b_1 \\ x_0 + x_2&= b_2 \end{align}$$

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

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

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

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

Por último, permutar as linhas $i$ e $j$ envolve multiplicar pela matriz identidade com as respectivas linhas permutadas:

In [7]:
E3 = np.eye(3)
E3[[1, 2]] = E3[[2, 1]]  # troca as linhas um e dois
E3
Out[7]:
array([[1., 0., 0.],
       [0., 0., 1.],
       [0., 1., 0.]])
In [8]:
E3 @ A
Out[8]:
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 $A$ ao aplicarmos à matriz aumentada $[A|I]$. Vamos aplicando operações de forma que obtenhamos uma matriz aumentada do tipo $[I|B]$, caso no qual $B$ será inversa de $A$ (isto é, $B = A^{-1}$):

In [9]:
AI = np.block([A, np.eye(3)])
AI
Out[9]:
array([[1., 0., 3., 1., 0., 0.],
       [3., 2., 1., 0., 1., 0.],
       [1., 0., 1., 0., 0., 1.]])
In [10]:
# Subtrair 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
Out[10]:
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:

In [11]:
# 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
Out[11]:
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 $A$ se encontra a partir da coluna três da matriz aumentada:

In [12]:
Ainv = AI[:, 3:]
Ainv
Out[12]:
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):

In [13]:
np.linalg.inv(A)
Out[13]:
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 $A$ ao comparar elemento por elemento do resultado da multiplicação $A^{-1} A$ com a matriz identidade:

In [14]:
np.allclose(Ainv @ A, np.eye(3))
Out[14]:
True

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

$$x^* = A^{-1} b$$
In [15]:
b = np.array([1., 2., 3.])
x = Ainv @ b
x
Out[15]:
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:

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

O procedimento automático

Como visto acima, é em geral desecessá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 $A$ é quadrada e seu determinane é diferente de zero) sem calcular a inversa. Por exemplo, o sistema $A x = b$ abaixo

$$\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 = (x_0, x_1) = (-2, 3/2)$:

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

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

In [18]:
A @ x
Out[18]:
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:

In [19]:
np.allclose(A @ x, b)
Out[19]:
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

$$A X = B$$

onde $B$ e a incógnita $X$ são matrizes. Isto é equivalente a resolver simultaneamente vários sistemas lineares do tipo $A x = b$, onde cada $x$ e $b$ são colunas das matrizes $X$ e $B$, respectivamente.

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

In [20]:
%%timeit
x = np.linalg.solve(A, b)
9.03 µs ± 194 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [21]:
%%timeit
x = np.linalg.inv(A) @ b
13.5 µs ± 394 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Quadrados mínimos

Caso um sistema linear $A 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 $Ax$ e $b$:

$$x^* =\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 = b$ existe, o mínimo $x^*$ acima é a própria solução:

In [22]:
x, _, _, _ = np.linalg.lstsq(A, b)
x
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  """Entry point for launching an IPython kernel.
Out[22]:
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 $A$ possui mais linhas que colunas (ou seja, mais equações que incógnitas):

In [23]:
A = np.array([[1., 2.], [3., 4.], [5., 6.]])
A
Out[23]:
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 $A$ não ser quadrada (tente fazer!). Mas podemos obter uma solução de quadrados mínimos:

In [24]:
b = np.array([1., 2., 2.])
x, residuals, _, _ = np.linalg.lstsq(A, b)
x
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:2: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  
Out[24]:
array([-0.66666667,  0.91666667])
In [25]:
A @ x
Out[25]:
array([1.16666667, 1.66666667, 2.16666667])

Observe como a multiplicação de $A$ por $x$ não retorna $b$ (na verdade, eu escolhi $b$ 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 $Ax$ e $b$:

In [26]:
residuals
Out[26]:
array([0.16666667])

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

In [27]:
np.sum((A @ x - b)**2)
Out[27]:
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.

Deixe uma resposta