Resolução de equações lineares em

Sobre Resolução de equações lineares

Lembras-te de quando aprendeste álgebra pela primeira vez, no secundário?

Normalmente, davam-nos um sistema de equações e pediam-nos para resolver em ordem a x e y.

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

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

Olha outra vez para o lado esquerdo do exemplo anterior. Parece uma matriz [4 -3; 2 7] a multiplicar um vetor desconhecido [x, y] para dar o vetor [-2, 16].

Em notação de Álgebra Linear: A x = b, seguindo a convenção habitual 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 no 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

Considera antes estas equações:

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

A segunda equação é apenas o dobro da primeira, não acrescentando informação nova. Quaisquer valores em que x = 0.75y serão uma solução.

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

As linhas da matriz não são linearmente independentes neste caso, e o posto da matriz é inferior ao número de linhas ou 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() devolve um único ponto, escolhido para ser normalizado para um vetor unitário, mas qualquer múltiplo escalar deste vetor também pertence ao espaço nulo.

Resolver A x = b

E se tiveres 1000 equações com 1000 incógnitas? Se isso soa disparatado, lembra-te de que os transdutores digitais sem fios são hoje baratos e versáteis (medem deformação, velocidade do vento, aceleração em 3 eixos, o que for...). Ter 1000 deles a monitorizar uma estrutura moderna como uma ponte suspensa é perfeitamente razoável. Tem de haver um sistema de controlo a interpretar o fluxo de dados, pronto para sinalizar um alerta se as coisas se tornarem alarmantes.

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

Um primeiro pensamento, para quem não está familiarizado com isto, é "dividir tudo por A" de alguma forma, para o passar para o lado direito.

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

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 matrizes grandes. Não é o ideal, se a tua ponte colapsar durante o cálculo!

Felizmente, outros algoritmos são drasticamente mais rápidos (neste caso, a eliminação de Gauss).

A Julia (a copiar o Matlab) usa simplesmente 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 pormenores a ter em atenção.

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

O comando rank() dá o número de linhas/colunas linearmente independentes, por isso o objetivo é que seja igual ao número de variáveis desconhecidas.

julia> rank(A)
2

Matrizes retangulares

Em criança, provavelmente ensinaram-te que precisas de um sistema de N equações para resolver N incógnitas.

Como acontece com a maioria das coisas na matemática, a realidade é um pouco mais matizada (e não apenas o pormenor da independência linear mencionado na secção anterior).

Poucas linhas

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 no espaço N-dimensional, mas N-1 dir-nos-á que a solução tem de estar algures numa reta.

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

A tua reta de soluções passa por uma região do espaço de parâmetros com que um engenheiro de renome se deva preocupar? Ou que pudesse justificar o encerramento da estrutura ao público (as pontes suspensas 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 dá-nos uma solução! Extraoficialmente, esta é aparentemente 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 dá 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

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

Demasiadas linhas

A situação oposta surge quando há mais equações do que incógnitas, mas, de alguma forma, continuam a ser linearmente independentes.

Na engenharia do mundo real, isto é perfeitamente 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 por mínimos quadrados aos dados com ruído.

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

Pode usar-se de forma muito semelhante à inversa de uma matriz quadrada (não singular), dando uma estimativa por 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

Valores próprios e vetores próprios

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

Esta secção é sobre soluções de A x = λ x, em que λ é um escalar.

Este é um conceito distintamente 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 alterar a sua direção (exceto que um λ negativo o inverte).

Isto soa a nicho, mas acaba por ser ridiculamente útil!

Terminologia: Os valores válidos de λ são os valores próprios de A, e os valores correspondentes de x são os vetores próprios de A. Infelizmente, temos de viver com palavras que mudam de língua a meio.

Normalmente ensina-se aos estudantes como calcular valores/vetores próprios para matrizes 2×2 à mão, mas usar um computador é muito mais fácil (embora ainda bastante lento, e escalando 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 valores próprios, embora nem sempre distintos. Pensa neles como as raízes de um polinómio de grau n (chamado polinómio característico), que podem ser repetidas, e são frequentemente complexas mesmo para uma matriz de valores reais.

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

Aplicações

Explicar porque os vetores próprios são importantes é, idealmente, um trabalho para um manual de 500 páginas, e não para alguns parágrafos curtos. Este tema permeia imensa matemática aplicada moderna.

Ao nível mais alto, os vetores próprios representam os eixos (direções) "mais importantes" de um conjunto de dados. Os valores próprios indicam a "importância relativa" de cada eixo (sujeita a alguns pressupostos sobre a normalização das entradas).

O que isto realmente significa depende da aplicação.

Análise de Componentes Principais

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

Precisamos de fazer redução de dimensionalidade, e a PCA é uma forma de trazer ordem ao caos:

• Calcula a matriz de covariância dos dados.
• Agora temos uma matriz quadrada, por isso a seguir calcula os valores e vetores próprios.
• Ordena os vetores próprios por ordem decrescente dos seus valores próprios (mais precisamente, do valor absoluto |λ|).
• Os k melhores vetores próprios são agora as tuas principal components, em que k é significativamente menor que o número original de características.
• Projeta o conjunto de dados completo nestes k eixos e começa a procurar padrões interessantes.

Nesta forma, a PCA tem sido usada por grupos tão diversos como cientistas médicos que tentam melhorar propriedades moleculares no desenho de fármacos, e ativistas políticos que tentam compreender as preferências dos eleitores a partir de dados de sondagens. Provavelmente também por grupos de marketing que tentam vender-te 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 da Julia, mas (fora do Exercism) o pacote MultivariateStats contém o que precisas.

Processamento de imagens

Levando a ideia da PCA um passo mais longe, as imagens digitais não passam de matrizes de valores de píxeis, e podemos calcular componentes principais para elas.

Isto tem sido usado há décadas na compressão de imagens, em que a PCA serve de guia para saber o que é mais importante preservar ao reduzir o tamanho do ficheiro.

Cada vez mais, a PCA é uma parte vital de classificadores de imagens, como o reconhecimento facial. Da próxima vez que passares por um aeroporto, o Big Brother não está só a observar: também usa Álgebra Linear para compreender o que vê!

Engenharia Mecânica

De forma menos polémica, os eixos principais de um componente mecânico são os vetores próprios do seu tensor de momento de inércia I (uma matriz com um nome ligeiramente diferente).

Cada roda do teu carro tem provavelmente um ou mais pesos de equilíbrio, ajustados para anular os elementos fora da diagonal neste tensor I. O mecânico da oficina evita certamente fazer as contas (ao contrário dos projetistas de aeronaves e dos engenheiros de foguetões), mas "oscilação" não passa de uma palavra do dia a dia para descrever os termos cruzados no tensor, que, neste caso, tornariam a tua viagem menos confortável e aumentariam o desgaste mecânico dos rolamentos.