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

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

O que é "Álgebra Linear"?

Há muitas definições técnicas, na Wikipedia, em livros didáticos e em sites excelentes como o 3Blue1Brown, mas para os nossos propósitos podemos ser mais informais.

Os mantenedores do Julia (e do Python) adoram esse tipo de coisa. A administração do Exercism, nem tanto.

Este conceito vai se restringir aos aspectos mais simples de um assunto enorme, mas fique avisado: este é, inevitavelmente, um conceito bastante matemático.

O que é um vetor?

Depende de quem você pergunta e de como você quer visualizar algo bastante abstrato.

• No Julia, Vector{T} é apenas um apelido para Array{T, 1}, e o eltype T pode ser qualquer coisa.
• Em (boa parte da) Física, um vetor é uma flecha com length e direction em um espaço N-dimensional, mas sem posição fixa.
• Na Álgebra Linear, a flecha dos físicos fica fixa no espaço, com a cauda na origem. Isso torna a haste da flecha redundante, então podemos representar o vetor como um point no espaço N-dimensional, com os N elementos representando a distância até a origem ao longo de cada um dos N eixos.

Para os nossos propósitos aqui, ignore vetores de strings ou caracteres. Vamos trabalhar com tipos numéricos neste conceito: Int, Float ou Complex.

O que é uma matriz?

De novo, você vai encontrar respostas variadas.

• Um array 2-D retangular de números (sem bordas 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, assim como uma function é aplicada a outros tipos de entrada.

Alguns tipos de matriz quadrada são comuns o bastante para ter nomes especiais.

Matriz diagonal: Todas as entradas não nulas ficam na main diagonal (do canto superior esquerdo ao 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 motivos que ficarão mais claros depois). Muitas vezes 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 sobre e acima da diagonal, zeros abaixo. Triangular inferior você já imagina.

Transpor uma matriz

Se você quiser mesmo trocar linhas por colunas, a função permutedims() faz isso e retorna uma nova matriz.

Isso envolve cópia, o que é lento e consome muita memória com matrizes grandes.

Para os propósitos da Álgebra Linear, a função transpose() é mais útil, pois cria rapidamente um wrapper lazy em torno da matriz original.

Ainda mais útil é a função adjoint(), que também inverte o sinal da parte imaginária em números complexos (se você se pergunta por que isso é tão importante, uma resposta rápida é a Mecânica Quântica). Essa operação é comum o bastante para que possamos simplesmente adicionar um apóstrofo ' ao nome da variável para criar o adjoint.

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 restante deste documento inclui funcionalidades do módulo LinearAlgebra. Uma linha using LinearAlgebra o traz para o namespace, mas não vamos repetir isso nos exemplos (poluição visual demais).

Vetores, produto elemento a elemento

Isso foi discutido no conceito Operações com vetores. O operador é .*, que atua sobre os vetores de entrada par a par, gerando uma saída do mesmo tamanho e tipo das entradas.

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

Vetores, produto escalar

Essa operação extremamente comum é escrita nos livros didáticos como u ⋅ v e chamada de produto escalar.

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

Dois vetores podem ser multiplicados com o operador * usual, mas somente se o vetor da esquerda for um adjoint: convenientemente escrito u' * v. Isso deve ficar mais claro na seção sobre multiplicação de matrizes mais adiante.

O ponto elevado (no centro) está disponível no Julia (digitado com \cdot e depois tab) como açúcar sintático para a função dot(). Não é preciso especificar o adjoint, 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 precisa ser o conjugado (sinal invertido na parte imaginária). A função dot faz isso automaticamente, a sintaxe u' faz isso explicitamente, mas sum(u .* v) vai falhar 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 acostumados a calcular torque ou vetores de momento angular.

Enquanto o produto escalar converte dois vetores em um scalar, o produto vetorial converte dois vetores de comprimento 3 em um terceiro vetor de comprimento 3.

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), elas não definem um plano, então a saída será [0, 0, 0].

A ordem importa: u × v == -(v × u). Essa é a famosa regra da mão direita, que já deixou muitos de nós encarando o polegar e dois dedos enquanto os giramos no espaço (muitas vezes com uma expressão de dúvida).

# 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

Repare que essa operação é restrita a vetores de comprimento 3 (equivalentes ao espaço euclidiano com eixos x, y, z ortogonais). Os matemáticos ficam felizes em definir a operação para outras dimensões, com uma salada de palavras incompreensível (tensores de ordem superior, produto wedge, produto exterior, multivectores...) como resultado. A maioria de nós foge ou muda de assunto rapidamente nesse ponto!

Normas

Quão "grande" é um vetor?

A norm é uma tentativa de capturar isso, reduzindo o vetor a um escalar apropriado.

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

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

Qualquer norma p pode ser calculada fornecendo p como segundo argumento. A norma 1 às vezes é útil: é simplesmente a soma dos valores absolutos das entradas (portanto 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

Às vezes parece que todo matemático aplicado do mundo passou boa parte dos últimos 80 anos transformando cada cálculo em uma série de multiplicações de matrizes.

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

• É altamente repetitiva.
• Pode ser paralelizada com eficiência.
• Muito hardware especializado foi desenvolvido para torná-la mais rápida, desde o Cray-1 de 80 milhões de dólares na década de 1970 até a GPU que provavelmente está integrada ao seu notebook.

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

Considere uma matriz A multiplicando um vetor v para dar uma saída w (por convenção da Álgebra Linear, usamos letras maiúsculas para matrizes e minúsculas para vetores).

A primeira linha de A é multiplicada escalarmente por 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]

Uma multiplicação matriz * matriz estende isso à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, poderíamos dizer que a primeira linha de A é multiplicada escalarmente por cada coluna de B para dar a primeira linha de C, a segunda linha dá a segunda linha, e assim por diante.

Existem várias representações mentais assim, e quem tiver interesse pode assistir a uma aula do MIT inteira discutindo elas.

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á começando, só pela leitura das palavras. O YouTube tem muitos vídeos demonstrando isso graficamente, então pesquise por "multiplicação de matrizes" e escolha um com o estilo, o nível de detalhe e o idioma que você preferir.

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 deve ser igual ao número de linhas da matriz da direita.

Expressando os tamanhos como tuplas (nrows, ncols), como a saída de size(A), temos (a, b) * (b, c) -> (a, c). As dimensões "internas", aqui b e b, são compatíveis para o cálculo de produtos escalares. As dimensões "externas", aqui a e c, determinam as dimensões da saída.

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, no estilo matricial, tem duas possibilidades.

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

Como 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 em 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 às vezes é chamado de produto interno.

O (menos comum) u * v' é, correspondentemente, o produto externo. Isso está relacionado ao produto tensorial, embora isso esteja muito além do nosso escopo.

Rotações

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

Em 2D, existe uma matriz relativamente simples para girar 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

Rotações em 3D precisam de uma matriz mais complexa, com dois ângulos. Isso é difícil de exibir em um documento Markdown sem incluir muito LaTeX, então consulte a Wikipedia para ver a fórmula.

Para gráficos 3-D, como OpenGL e seus sucessores, a convenção é trabalhar com coordenadas homogêneas.

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

Isso permite matrizes de transformação mais elaboradas. A rotação ainda fica em A[1:3, 1:3], as translações são [Δx, Δy, Δz] em A[1:3, 4], a escala fica na diagonal, e existem outras possibilidades para cisalhamento, perspectiva etc.

Desempenho

Nos departamentos de Ciência da Computação do mundo todo, há muitas teses de doutorado com um ou mais capítulos sobre pequenas melhorias incrementais em algoritmos de multiplicação de matrizes. Isso é muito importante, e várias organizações estão dispostas a financiar a pesquisa em benefício próprio.

Desempenho é um assunto amplo, no qual não podemos nos aprofundar, 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 improvisada (use 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)

Aproximadamente, um par de matrizes de um milhão de elementos levou alguns milissegundos, e matrizes de 100 milhões de elementos levaram alguns segundos. Adicionar várias ordens de magnitude a mais vai exigir um hardware melhor...