Como resolver sistemas lineares em Julia

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

Assim como eu fiz com a série em Python, vou continuar com o tema da resolução de sistemas lineares, com exemplos de código agora em Julia. 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.

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=b1 am,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:

using LinearAlgebra  # define diagm

E1 = diagm(ones(3))
E1[1, 3] = -3  # três vezes a linha três somada à linha um
E1
3×3 Array{Float64,2}:
 1.0  0.0  -3.0
 0.0  1.0   0.0
 0.0  0.0   1.0

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 = [1 0 3
     3 2 1
     1 0 1]
3×3 Array{Int64,2}:
 1  0  3
 3  2  1
 1  0  1
E1 * A
3×3 Array{Float64,2}:
 -2.0  0.0  0.0
  3.0  2.0  1.0
  1.0  0.0  1.0

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 = diagm(ones(3))
# a linha um será multiplicada por um terço negativo
# (equivalentemente: dividida por menos três)
E2[1, 1] = -1/3
E2
3×3 Array{Float64,2}:
 -0.333333  0.0  0.0
  0.0       1.0  0.0
  0.0       0.0  1.0
E2 * A
3×3 Array{Float64,2}:
 -0.333333  0.0  -1.0
  3.0       2.0   1.0
  1.0       0.0   1.0

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

E3 = diagm(ones(3))
# troca as linhas dois e três
E3[2, :], E3[3, :] = E3[3, :], E3[2, :]
E3
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  0.0  1.0
 0.0  1.0  0.0
E3 * A
3×3 Array{Float64,2}:
 1.0  0.0  3.0
 1.0  0.0  1.0
 3.0  2.0  1.0

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á a inversa de AA (isto é, B=A1B = A^{-1}):

AI = [A I]
3×6 Array{Int64,2}:
 1  0  3  1  0  0
 3  2  1  0  1  0
 1  0  1  0  0  1
# Subtrair a linha um da linha três
E1 = diagm(ones(3))
E1[3, 1] = -1
AI = E1 * AI

# Subtrair três vezes a linha um da linha dois
E1 = diagm(ones(3))
E1[2, 1] = -3
AI = E1 * AI

# Subtrair quarto vezes a linha três da linha dois
E1 = diagm(ones(3))
E1[2, 3] = -4
AI = E1 * AI

# Somar três meios vezes a linha três da linha um
E1 = diagm(ones(3))
E1[1, 3] = 1.5
AI = E1 * AI

AI
3×6 Array{Float64,2}:
 1.0  0.0   0.0  -0.5  0.0   1.5
 0.0  2.0   0.0   1.0  1.0  -4.0
 0.0  0.0  -2.0  -1.0  0.0   1.0

Resta agora dividir as linhas dois e três por dois e menos dois, respectivamente:

# Dividir a linha dois por dois
E2 = diagm(ones(3))
E2[2, 2] = 0.5
AI = E2 * AI

# Dividir a linha três por menos dois
E2 = diagm(ones(3))
E2[3, 3] = -0.5
AI = E2 * AI

AI
3×6 Array{Float64,2}:
 1.0  0.0  0.0  -0.5  0.0   1.5
 0.0  1.0  0.0   0.5  0.5  -2.0
 0.0  0.0  1.0   0.5  0.0  -0.5

A inversa de AA se encontra a partir da quarta coluna da matriz aumentada:

Ainv = AI[:, 4:end]
3×3 Array{Float64,2}:
 -0.5  0.0   1.5
  0.5  0.5  -2.0
  0.5  0.0  -0.5

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

inv(A)
3×3 Array{Float64,2}:
 -0.5  0.0   1.5
  0.5  0.5  -2.0
  0.5  0.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:

Ainv * A ≈ I
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 (embora na prática seja ainda muito caro):

x=A1bx^* = A^{-1} b
b = [1., 2., 3.]
x = Ainv * b
3-element Array{Float64,1}:
  4.0
 -4.5
 -1.0

Ou, poderíamos ter obtido a solução da maneira mais fácil, usando a função backslash (sim, igual ao MATLAB/Octave):

A \ b
3-element Array{Float64,1}:
  4.0
 -4.5
 -1.0

O procedimento automático

Como visto acima, é em geral desnecessário inverter uma matriz para resolver um sistema linear. A função backslash resolve sistemas lineares com solução única (aqueles em que a matriz AA é quadrada e seu determinane é diferente de zero) sem calcular a inversa (ela não faz só isso, veremos mais para frente). 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 = [1. 2.
     3. 4.]
b = [1., 0.]
x = A \ b
2-element Array{Float64,1}:
 -1.9999999999999998
  1.4999999999999998

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
2-element Array{Float64,1}:
 0.9999999999999998
 0.0

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

A * x ≈ b
true

Não vou mostrar exemplos, mas vou fazer menção de que a mesma função backslash 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 (na vida? 🤔). Na prática, inverter matrizes custa em geral mais caro computacionalmente:

using BenchmarkTools

@benchmark x = A \ b
BenchmarkTools.Trial: 
  memory estimate:  336 bytes
  allocs estimate:  4
  --------------
  minimum time:     322.857 ns (0.00% GC)
  median time:      467.143 ns (0.00% GC)
  mean time:        501.969 ns (8.92% GC)
  maximum time:     31.018 μs (98.29% GC)
  --------------
  samples:          10000
  evals/sample:     210
@benchmark x = inv(A) * b
BenchmarkTools.Trial: 
  memory estimate:  1.48 KiB
  allocs estimate:  6
  --------------
  minimum time:     573.404 ns (0.00% GC)
  median time:      802.128 ns (0.00% GC)
  mean time:        895.248 ns (9.01% GC)
  maximum time:     35.123 μs (97.12% GC)
  --------------
  samples:          10000
  evals/sample:     94

No caso acima, a diferença é de metade do tempo, mas para matrizes maiores a diferença será maior ainda.

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 assim formulado é dito "de quadrados mínimos". Interessantemente, ele também é resolvido pela função backslash! Acontece que esta função possui múltiplos despachos (Bezanson (2017)):

Observe que, quando uma solução para Ax=bAx = b existe, o mínimo xx^* acima é a própria solução e, portanto, não há ambiguidades no retorno da função backslash.

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

A = [1. 2.
     3. 4.
     5. 6.]
3×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
 5.0  6.0

Sistemas desse tipo não possuem solução em geral. (De fato, em Python, numpy.linalg.solve daria um erro do tipo LinAlgError devido à matriz AA não ser quadrada). Mas uma solução de quadrados mínimos é em geral suficiente:

b = [1., 2., 2.]
x = A \ b
2-element Array{Float64,1}:
 -0.666666666666669
  0.9166666666666686
A * x
3-element Array{Float64,1}:
 1.1666666666666683
 1.6666666666666674
 2.166666666666667

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). Podemos verificar o resíduo da equação (7) usando sum:

sum((A * x - b) .^ 2)
0.16666666666666682

É isso aí! No próximo texto, vamos ver alguma coisa sobre autovalores e autovetores em Julia.

Referências