Resolução de equações lineares em

Sobre Resolução de equações lineares

Lembra quando você aprendeu álgebra pela primeira vez, no ensino médio?

Normalmente, recebíamos um sistema de equações e tínhamos que resolver para x e y.

4x - 3y = -2
2x + 7y = 16

Duas equações, duas incógnitas: um pouco de substituição e logo você percebe que x = 1, y = 2.

Olhe de novo para o lado esquerdo do exemplo anterior. Ele se parece com uma matriz [4 -3; 2 7] multiplicando um vetor desconhecido [x, y] para dar o vetor [-2, 16].

Na notação de Álgebra Linear: A x = b, seguindo a convenção usual de letras maiúsculas para matrizes e minúsculas para vetores.

A x = 0 e o espaço nulo

Um caso especial importante é quando todas as nossas equações têm zero do lado direito.

4x - 3y = 0
2x + 7y = 0

No exemplo acima, a única solução é quando x = y = 0, o que normalmente não é interessante.

Soluções não triviais para A x = 0 só existem quando A tem determinante zero.

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

Considere, em vez disso, estas equações:

4x - 3y = 0
8x + 6y = 0

A segunda equação é apenas o dobro da primeira, sem acrescentar nenhuma informação nova. Qualquer valor em que x = 0.75y será uma solução.

Esses valores de (x, y) ficam ao longo de uma reta no espaço 2D. Há um número infinito de soluções, e elas formam o espaço nulo da matriz.

As linhas da matriz não são linearmente independentes nesse caso, e o posto da matriz é menor que o número de linhas ou de colunas.

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

A função nullspace() retorna um único ponto, escolhido para ser normalizado como vetor unitário, mas qualquer múltiplo escalar desse vetor também está no espaço nulo.

Resolvendo A x = b

Agora suponha que você tenha 1000 equações com 1000 incógnitas? Se isso parece bobagem, lembre-se de que os transdutores digitais sem fio hoje são baratos e versáteis (medem deformação, velocidade do vento, aceleração em 3 eixos, o que for...). Ter 1000 deles monitorando uma estrutura moderna como uma ponte pênsil é perfeitamente razoável. É preciso haver um sistema de controle interpretando o fluxo de dados, pronto para sinalizar um alerta se as coisas ficarem preocupantes.

De novo, A x = b, em que temos A (do projeto de engenharia) e b (valores medidos pelos transdutores), mas precisamos encontrar x.

Um primeiro pensamento, para quem não está familiarizado, é de algum modo "dividir tudo por A" para passá-lo para o lado direito.

Na verdade, A tem uma inversa (a maioria das matrizes quadradas tem, mas não todas: o determinante precisa ser diferente de zero), e o cálculo funciona (devagar!).

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

Infelizmente, calcular a inversa é lento para matrizes pequenas e glacial para as grandes. Nada ideal, se sua ponte desabar durante o cálculo!

Felizmente, outros algoritmos são dramaticamente mais rápidos (neste caso, a eliminação gaussiana).

O Julia (copiando o Matlab) simplesmente usa uma barra invertida para o solucionador (tecnicamente, "divisão à esquerda").

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

Este é um exemplo trivial, mas há alguns detalhes com que tomar cuidado.

• Os coeficientes precisam ficar alinhados em colunas, então insira zeros conforme necessário para os termos que faltam.
• As equações (e, portanto, as linhas da matriz) precisam ser linearmente independentes. Se a linha 3 for apenas a soma das linhas 1 e 2, ela não acrescenta nenhuma informação nova e o problema fica subdeterminado.

O comando rank() dará o número de linhas/colunas linearmente independentes, então procure fazer com que isso seja igual ao número de variáveis desconhecidas.

julia> rank(A)
2

Matrizes retangulares

Em algum momento, quando você era criança, provavelmente lhe ensinaram que você precisa de um sistema de N equações para resolver N incógnitas.

Como na maioria das coisas em matemática, a realidade é um pouco mais sutil (não só o detalhe da independência linear, mencionado na seção anterior).

Linhas "de menos"

Com N-1 equações (linhas na matriz A), o problema fica subdeterminado. Desistimos em desespero?

Depende!

As N-1 equações ainda contêm muita informação. Em termos geométricos, precisamos de N para determinar a solução como um ponto em um espaço de N dimensões, mas N-1 nos dirá que a solução deve ficar em algum lugar de uma reta.

Isso nos deixa com um número infinito de soluções, mas também com uma quantidade infinita de espaço de parâmetros sem solução alguma.

Sua reta de soluções passa por uma região do espaço de parâmetros com que um engenheiro de reputação deveria se preocupar? Ou que poderia justificar fechar a estrutura ao acesso do público (pontes pênsis podem ficar animadas com ventos laterais fortes)? Talvez valha a pena fazer mais alguns cálculos para determinar isso!

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

O solucionador com barra invertida nos dá uma solução! Não oficialmente, aparentemente esta é a solução com a menor norma (a mais próxima da origem, geometricamente).

Para obter as outras soluções, precisamos do nullspace de A2.

Somar qualquer múltiplo escalar do espaço nulo à solução original dará outra solução válida.

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

Da mesma forma, N-2 equações vão restringir as soluções (infinitas) a um plano, e os mesmos princípios se aplicam.

Linhas "de mais"

A situação oposta surge quando há mais equações do que incógnitas, mas, de algum modo, elas ainda são linearmente independentes.

Na engenharia do mundo real, isso é totalmente normal e visto como algo bom!

Os transdutores têm precisão limitada, o vetor b tem precisão limitada, e há ruído no cálculo. Agora a solução é um ajuste de mínimos quadrados aos dados ruidosos.

A técnica mais simples usa a pseudoinverse da matriz, que o Julia implementa como a função pinv().

Ela pode ser usada de forma muito parecida com a inversa de uma matriz quadrada (não singular), dando uma estimativa de mínimos quadrados das variáveis em vez de uma solução exata.

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

Autovalores e autovetores

A seção anterior descreveu soluções para A x = b, equivalentes à álgebra familiar.

Esta seção trata das soluções de A x = λ x, em que λ é um escalar.

Este é um conceito tipicamente de Álgebra Linear, mas podemos interpretá-lo geometricamente:

Para valores adequados de x e λ, a matriz quadrada A escala o vetor não nulo x por um fator λ em comprimento, sem mudar sua direção (exceto que um λ negativo a inverte).

Isso parece muito específico, mas acontece que é ridiculamente útil!

Terminologia: Os valores válidos de λ são os autovalores de A, e os valores correspondentes de x são os autovetores de A. Infelizmente, só nos resta conviver com palavras que trocam de idioma no meio do caminho.

Normalmente, ensina-se aos estudantes como calcular autovalores/autovetores de matrizes 2×2 manualmente, mas usar um computador é bem mais fácil (embora ainda seja razoavelmente lento, e escale como O(n^3) para uma matriz n×n no caso geral).

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

Em geral, uma matriz n×n terá n autovalores, embora nem sempre distintos entre si. Pense neles como as raízes de um polinômio de grau n (chamado de polinômio característico), que podem se repetir e muitas vezes são complexas, mesmo para uma matriz de valores reais.

Cada autovetor representa uma direção, e qualquer múltiplo escalar também é um autovetor válido. Por conveniência em cálculos posteriores, o Julia retorna vetores unitários com norma igual a 1.

Aplicações

Explicar por que os autovetores são importantes é, idealmente, tarefa para um livro de 500 páginas, e não para alguns parágrafos curtos. Esse tópico permeia tanta coisa da matemática aplicada moderna.

No nível mais alto, os autovetores representam os eixos (direções) "mais importantes" de um conjunto de dados. Os autovalores indicam a "importância relativa" de cada eixo (sujeitos a algumas suposições sobre a normalização das entradas).

O que isso de fato significa depende da aplicação.

Análise de componentes principais

Em um problema típico de ciência de dados, podemos ter 100 ou mais "características", armazenadas como colunas de dados. Quase inevitavelmente, haverá ruído, redundância e correlações indesejadas.

Precisamos reduzir a dimensionalidade, e a PCA é uma forma de trazer ordem ao caos:

• Calcule a matriz de covariância dos dados.
• Agora temos uma matriz quadrada, então calcule em seguida os autovalores e autovetores.
• Ordene os autovetores em ordem decrescente de seus autovalores (mais precisamente, do valor absoluto |λ|).
• Os k primeiros autovetores são agora seus principal components, em que k é bem menor que o número original de características.
• Projete o conjunto de dados completo sobre esses k eixos e comece a procurar padrões interessantes.

Nessa forma, a PCA tem sido usada por grupos tão diversos quanto cientistas médicos tentando melhorar propriedades moleculares no desenvolvimento de medicamentos e ativistas políticos tentando entender as preferências dos eleitores a partir de dados de pesquisas. Provavelmente também por equipes de marketing tentando lhe vender coisas, mas qualquer tecnologia pode ser usada para o bem ou para o mal.

A PCA não faz parte de uma instalação mínima do Julia, mas (fora do Exercism) o pacote MultivariateStats contém o que você precisa.

Processamento de imagens

Levando a ideia da PCA um passo adiante, imagens digitais são apenas matrizes de valores de pixels, e podemos calcular componentes principais para elas.

Isso tem sido usado há décadas na compressão de imagens, em que a PCA orienta o que é mais importante preservar ao reduzir o tamanho do arquivo.

Cada vez mais, a PCA é parte vital de classificadores de imagens, como o reconhecimento facial. Da próxima vez que você passar por um aeroporto, o Big Brother não está apenas observando: ele também usa Álgebra Linear para entender o que vê!

Engenharia mecânica

De forma menos controversa, os eixos principais de um componente mecânico são os autovetores de seu tensor de momento de inércia I (uma matriz com um nome um pouco diferente).

Cada roda do seu carro provavelmente tem um ou mais contrapesos, ajustados para zerar os elementos fora da diagonal nesse tensor I. O mecânico da oficina sem dúvida evita fazer as contas (ao contrário dos projetistas de aeronaves e engenheiros de foguetes), mas "tremida" é apenas uma palavra do dia a dia para descrever os termos cruzados do tensor, que, nesse caso, deixariam sua viagem menos confortável e aumentariam o desgaste mecânico dos rolamentos.