Noções básicas de álgebra linear em

Sobre Noções básicas de álgebra linear

O que é a "Álgebra Linear"?

Há muitas definições técnicas, na Wikipédia, em manuais e em sites excelentes como o 3Blue1Brown, mas para os nossos objetivos podemos ser mais informais.

Os responsáveis pela manutenção do Julia (e do Python) adoram este tipo de coisas. A administração do Exercism, nem por isso.

Este conceito vai ficar limitado aos aspetos mais simples de um tema enorme, mas atenção: este é, inevitavelmente, um conceito bastante matemático.

O que é um vetor?

Depende de quem perguntas, e de como queres visualizar algo bastante abstrato.

• Em Julia, Vector{T} é apenas um sinónimo de Array{T, 1}, e o eltype T pode ser qualquer coisa.
• Em (grande parte da) Física, um vetor é uma seta com um length e uma direction no espaço N-dimensional, mas sem posição fixa.
• Em Álgebra Linear, a seta dos físicos fica fixa no espaço, com a cauda na origem. Isto torna o corpo da seta redundante, por isso podemos representar o vetor como um point no espaço N-dimensional, em que os N elementos representam a distância à origem ao longo de cada um dos N eixos.

Para já, ignora vetores de strings ou de carateres. Neste conceito, vamos trabalhar com tipos numéricos: Int, Float ou Complex.

O que é uma matriz?

Mais uma vez, vais encontrar respostas variadas.

• Um array retangular 2-D de números (sem linhas de comprimentos irregulares). Uma matriz quadrada é um caso especial comum.
• Uma linear combination de vetores coluna, empilhados lado a lado.
• Uma transformation que pode ser aplicada a um vetor, tal como uma function é aplicada a outros tipos de entrada.

Alguns tipos de matriz quadrada são suficientemente comuns para terem nomes especiais.

Matriz diagonal: Todas as entradas não nulas estão na main diagonal (do canto superior esquerdo ao canto inferior direito).

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

Matriz identidade: Uma matriz diagonal com apenas uns na diagonal (por razões que se vão tornar mais claras adiante). Frequentemente abreviada como I.

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Triangular superior: Valores não nulos na diagonal e acima dela, zeros abaixo. Triangular inferior: essa já adivinhas.

Transpor uma matriz

Se quiseres mesmo trocar as linhas pelas colunas, a função permutedims() faz isso e devolve-te uma nova matriz.

Isto implica copiar, o que é lento e consome muita memória quando trabalhas com matrizes grandes.

Para efeitos de Álgebra Linear, a função transpose() é mais útil, porque cria rapidamente um wrapper preguiçoso em torno da matriz original.

Ainda mais útil é a função adjoint(), que também troca o sinal da parte imaginária de quaisquer números complexos (se te perguntas porque é tão importante, uma resposta rápida é a Mecânica Quântica). Esta operação é suficientemente comum para podermos simplesmente acrescentar um apóstrofo ' ao nome da variável para criar a adjunta.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

Multiplicação

O resto deste documento inclui funcionalidades do módulo LinearAlgebra. Uma linha using LinearAlgebra traz esse módulo para o namespace, mas não vamos repetir isso nos exemplos (demasiado ruído visual).

Vetores, produto elemento a elemento

Isto foi abordado no conceito Operações com vetores. O operador é .*, que atua sobre os vetores de entrada aos pares, dando um resultado com o mesmo tamanho e tipo das entradas.

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

Vetores, produto escalar

Esta operação extremamente comum é escrita nos manuais como u ⋅ v e designa-se produto escalar.

Um produto escalar é equivalente à soma do produto elemento a elemento.

Dois vetores podem ser multiplicados com o operador habitual *, mas só se o vetor da esquerda for uma adjoint: convenientemente escrito u' * v. Isto ficará mais claro na secção mais à frente sobre multiplicação de matrizes.

O ponto elevado (ao centro) está disponível no Julia (introduzido com \cdot e a tecla de tabulação) como açúcar sintático para a função dot(). Não é preciso especificar a adjunta, pois esse detalhe é tratado automaticamente.

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

Para vetores com valores complexos, o vetor da esquerda tem de ser o conjugado (sinal invertido na parte imaginária). A função dot faz isso automaticamente, a sintaxe u' fá-lo explicitamente, mas sum(u .* v) falha se u e v forem complexos.

Vetores, produto vetorial

Um pouco de nostalgia para quem já fez um curso de Eletricidade e Magnetismo! E também para engenheiros habituados a calcular vetores de binário ou de momento angular.

Enquanto o produto escalar converte dois vetores num scalar, o produto vetorial converte dois 3-vetores num terceiro 3-vetor.

Na representação geométrica do espaço vetorial, o novo vetor é perpendicular ao plano que contém os dois vetores de entrada. Se as duas entradas forem paralelas (a menos de um sinal), não definem um plano, por isso o resultado será [0, 0, 0].

A ordem importa: u × v == -(v × u). Esta é a famosa Regra da Mão Direita, que já deixou muitos de nós a olhar para o polegar e dois dedos enquanto os rodávamos no espaço (muitas vezes com uma expressão de perplexidade).

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

Repara que esta operação está restrita a vetores de comprimento 3 (equivalente ao espaço euclidiano com eixos ortogonais x, y, z). Os matemáticos definem a operação de bom grado para outras dimensões, tendo como resultado uma salada de palavras incompreensível (tensores de ordem superior, produto cunha, produto exterior, multivectores...). A maioria de nós foge ou muda rapidamente de assunto nesse momento!

Normas

Quão "grande" é um vetor?

A norm é uma tentativa de captar isso, reduzindo o vetor a um escalar adequado.

Existe toda uma família de normas definidas, mas de longe a mais comum é a norma 2, que é √(v ⋅ v).

Esta operação de raiz da média dos quadrados é a distância pitagórica à origem (no espaço N-dimensional). Se visualizarmos o vetor como uma seta, com a cauda na origem, a norma 2 é o comprimento da seta.

Qualquer norma p pode ser calculada fornecendo p como segundo argumento. A norma 1 é por vezes útil: é simplesmente a soma dos valores absolutos das entradas (por isso é muito rápida e fácil de calcular).

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

Multiplicação de matrizes

Por vezes parece que todos os matemáticos aplicados do mundo passaram grande parte dos últimos 80 anos a transformar cada cálculo numa série de multiplicações de matrizes.

Esta é uma operação em que os computadores são muito bons:

• É altamente repetitiva.
• Pode ser paralelizada com eficiência.
• Foi desenvolvido muito hardware especializado para a tornar mais rápida, desde o Cray-1, de 80 milhões de dólares, na década de 1970, até à GPU que provavelmente está integrada no teu portátil.

Os detalhes são bastante simples, embora não muito intuitivos à primeira vista.

Considera uma matriz A a multiplicar um vetor v para dar um resultado w (por convenção da Álgebra Linear, usamos letras maiúsculas para matrizes e minúsculas para vetores).

Faz-se o produto escalar da primeira linha de A com v para obter o primeiro elemento de w, a segunda linha dá o segundo elemento, e assim por diante.

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

A multiplicação matriz * matriz estende isto às colunas da matriz da direita.

Para C = A * B, podemos considerar que A multiplica cada coluna de B para dar a coluna correspondente de C: uma série de multiplicações matriz-vetor.

De forma equivalente, podemos dizer que a primeira linha de A faz produto escalar com cada coluna de B para dar a primeira linha de C, a segunda linha dá a segunda linha, e assim por diante.

Há várias representações mentais deste tipo, e os estudantes interessados podem ver uma aula completa do MIT a discuti-las.

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53

julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

Como mostra o exemplo acima, a multiplicação de matrizes não é comutativa: não há uma relação simples entre A*B e B*A.

A multiplicação de matrizes é provavelmente difícil de visualizar para quem está a começar, só de ler as palavras. O YouTube tem imensos vídeos que a demonstram graficamente, por isso procura por "multiplicação de matrizes" e escolhe um com o estilo, o nível de detalhe e o idioma que preferires.

Dimensões

O produto escalar de dois vetores depende de eles terem o mesmo comprimento.

Por extensão, na multiplicação de matrizes o número de colunas da matriz da esquerda tem de corresponder ao número de linhas da matriz da direita.

Se expressarmos os tamanhos como tuplos (nrows, ncols), como os que size(A) devolve, temos (a, b) * (b, c) -> (a, c). As dimensões "internas", aqui b e b, são compatíveis para se calcular o produto escalar. As dimensões "externas", aqui a e c, determinam as dimensões do resultado.

Um exemplo com matrizes retangulares:

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

Multiplicar pares de vetores ao estilo matricial tem duas possibilidades.

Por convenção, usaríamos u' * v como equivalente ao produto escalar. As dimensões são (1, 3) * (3, 1), e o Julia simplifica o resultado (1, 1) para um escalar (ao contrário, por exemplo, do R).

Em alternativa, poderíamos usar u * v', com dimensões (3, 1) * (1, 3) -> (3, 3), para multiplicar todos os pares possíveis de elementos e expandir os vetores para uma matriz.

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * v é por vezes chamado de produto interno.

O (menos comum) u * v' é, correspondentemente, o produto externo. Está relacionado com o produto tensorial, embora isso esteja muito para lá do nosso âmbito.

Rotações

Um tipo particularmente comum de multiplicação de matrizes envolve matrizes de rotação.

Em 2D, existe uma matriz relativamente simples para rodar um vetor no sentido anti-horário em θ radianos.

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0

As rotações em 3D precisam de uma matriz mais complexa, com dois ângulos. Isto é difícil de mostrar num documento Markdown sem incorporar muito LaTeX, por isso consulta a Wikipédia para ver a fórmula.

Para gráficos 3D, como o OpenGL e os seus sucessores, a convenção é trabalhar com Coordenadas Homogéneas.

Cada vértice é representado por um 4-vetor (ou, de forma equivalente, uma coluna de matriz): [x, y, z, 1.0] para um ponto em (x, y, z).

Isto permite matrizes de transformação mais elaboradas. A rotação continua em A[1:3, 1:3], as translações são [Δx, Δy, Δz] em A[1:3, 4], a escala está na diagonal, e existem outras possibilidades para inclinação, perspetiva, etc.

Desempenho

Nos departamentos de informática de todo o mundo há muitas teses de doutoramento com um ou mais capítulos sobre pequenas melhorias incrementais nos algoritmos de multiplicação de matrizes. Isto é mesmo importante, e várias organizações estão dispostas a financiar a investigação em seu próprio benefício.

O desempenho é um tema vasto, no qual não podemos entrar em profundidade, mas a título de ilustração podemos tentar multiplicar matrizes aleatórias de vários tamanhos.

O código abaixo é uma estimativa rápida e grosseira (usa o BenchmarkTools.jl para uma abordagem melhor).

O sistema usado foi um PC pequeno de 430 dólares (EUA): processador Ryzen 9, 32 GB de RAM, Linux Mint 22.1, Julia 1.11.6.

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

Grosso modo, um par de matrizes de um milhão de elementos demorou alguns milissegundos, e matrizes de 100 milhões de elementos demoraram alguns segundos. Acrescentar mais algumas ordens de grandeza vai exigir melhor hardware...