Como manipular matrizes e vetores em Julia

Há um ano eu escrevi uma minissérie sobre álgebra linear com Python. Aqui vou fazer o mesmo com a liguagem Julia. Vamos falar basicamente sobre matrizes e vetores nativos da linguagem (não há necessidade de nenhum pacote extra):

Este texto fala da construção de vetores e matrizes, suas dimensões, operações e algumas funções matriciais, com exemplos de código em Julia.

Arrays

Representamos matrizes e vetores através do objeto array. Podemos pensar em arrays unidimensionais como listas de números e arrays bidimensionais como matrizes. É possível construir arrays de dimensões arbitrárias (representando tensores de dimensões três, quatro, etc.), mas não vamos falar disso aqui.

Uma diferença do Python é que em Julia, assim como no MATLAB/Octave, quando selecionamos uma linha ou coluna de um array bidimensional, o resultado também é bidimensional. Portanto, não há ambiguidade quando estamos considerando arrays unidimensionais como vetores-linha ou vetores-coluna.

Vamos criar arrays unidimensionais, a partir de uma lista de elementos, e determinar alguns dos seus parâmetros:

b = [1., 0.]
2-element Array{Float64,1}:
 1.0
 0.0

O tipo deste array pode ser determinado por typeof:

typeof(b)
Array{Float64,1}

O tipo dos elementos e as dimensões são determinadas pelo tipo do array neste caso. Ou seja, temos um vetor de números de ponto flutuante (Float64, iguais a double em C/C++). Também podemos verificar as dimensões usando ndims:

ndims(b)
1

O número total de elementos do array é dois:

length(b)
2

O formato do array é a disposição desses elementos nas dimensões. No array b, os dois elementos estão obviamente dispostos na única dimensão:

size(b)
(2,)

Os elementos de um array unidimensional são indexados por suas posições (começando em um, ao contrário do Python, C/C++, etc., mas seguindo a conveção do MATLAB/Octave e do Fortran). O segundo elemento de b, portanto, é zero:

b[2]
0.0

Uma nota antes de entrarmos em matrizes: b é um vetor-coluna. Estes são definidos com vírgulas ([1., 0.]). Ao usarmos apenas espaços como separadores, obtemos vetores-linha (se você já usou MATLAB/Octave, você está em casa):

c = [1. 0.]
1×2 Array{Float64,2}:
 1.0  0.0

Observe que vetores-linha são essencialmente matrizes 1×n1 \times n (array bidimensional). Já matrizes m×nm \times n são definidas linha por linha de maneira similar aos vetores-linha, com ponto e vírgula separando linhas (na verdade, o ponto e vírgula é opcional):

A = [1. 2.;
     3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
typeof(A)
Array{Float64,2}

Esse é um array bidimensional:

ndims(A)
2

Semelhante a outras linguagens, a primeira dimensão corresponde à direção vertical (contagem de linhas) e a segunda corresponde à direção horizontal (contagem de colunas). Por exemplo, o segundo elemento da primeira linha é (Julia usa contagens começando em um):

A[1, 2]
2.0

Contar do final de uma direção (ou seja, de trás para frente) requer o uso de end (isto difere do Python, que usa números negativos). O último elemento da última coluna, portanto, é quatro:

A[end, end]
4.0

O número total de elementos em A é quatro:

length(A)
4

Esses elementos estão dispostos dois a dois nas duas dimensões (ou seja, 2×2=42 \times 2 = 4 elementos):

size(A)
(2, 2)

Podemos fatiar a segunda coluna inteira da matriz A (ou seja, todos os índices de uma coluna):

col = A[:, 2]
2-element Array{Float64,1}:
 2.0
 4.0

Observe que a coluna obtida é um array unidimensional (vetores em Julia são vetores-coluna por padrão, assim como em MATLAB/Octave):

ndims(col)
1

Seu formato portanto é (2,) (ambos elementos na única dimensão):

size(col)
(2,)

Como transformar um vetor-coluna em vetor-linha? Por transposição e, assim como em MATLAB/Octave, podemos usar o apóstrofo (') para isso:

col_t = col'
1×2 LinearAlgebra.Adjoint{Float64,Array{Float64,1}}:
 2.0  4.0

Como já dito antes sobre vetores-linha em geral, nosso novo array possui duas dimensões:

ndims(col_t)
2

O formato, como esperado, é (1, 2):

size(col_t)
(1, 2)

Ambos objetos col e col_t guardam as mesmas informações, mas são fundamentalmente diferentes. Por exemplo, irão se comportar como esperado frente multiplicações, como veremos adiante.

Mas, antes disso, vamos ver alguma coisa sobre construções um pouco mais avançadas de matrizes.

Construção de matrizes em bloco

É comum precisar contruir matrizes bloco a bloco. Não há necessidade de nenhuma função especial para isso em Julia (abaixo usamos também zeros, para construir uma matriz de zeros, e I, do pacote padrão LinearAlgebra, para construir uma matriz identidade, e em ambos os casos das dimensões desejadas):

using LinearAlgebra

B = [A zeros(2, 2)
     zeros(2, 2) I]
4×4 Array{Float64,2}:
 1.0  2.0  0.0  0.0
 3.0  4.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

Observe que não precisamos especificar as dimensões de I. A flexibilidade da lingugagem se equipara à sintaxe do MATLAB/Octave neste e em muitos outros casos.

A operação inversa, ou seja, fatiar blocos de arrays, pode ser feita com : (assim como em Python) e índices de início e fim (ambos inclusos). Por exemplo, para se fatiar a submatriz 2×22 \times2 central de B:

subB = B[2:3, 2:3]
2×2 Array{Float64,2}:
 4.0  0.0
 0.0  1.0

Os índices incluídos são [2, 3], incluindo o três, ao contrário do que ocorreria em Python, que não incluiria o último índice (lá isso faz sentido, já que as contagens iniciam em zero, pois assim o tamanho dos resultados continua sendo o esperado).

O resultado é um array menor, mas de mesma dimensão:

ndims(subB)
2
size(subB)
(2, 2)

Operações com matrizes e vetores

Operações aritméticas de matrizes e vetores

As operações de adição (+), subtração (-), divisão (/), multiplicação (*) e exponenciação (**) podem ser realizadas entre arrays e escalares:

A
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
2 * A
2×2 Array{Float64,2}:
 2.0  4.0
 6.0  8.0

Podemos portanto manipular arrays unidimensionais como esperado para vetores, ou seja, em termos de soma de vetores e multiplicação por escalar:

u = [1., -1.]
v = [1.,  1.]
u + v
2-element Array{Float64,1}:
 2.0
 0.0
2.0 * u
2-element Array{Float64,1}:
  2.0
 -2.0

Quando realizadas entre arrays, no entanto, essas operações são como esperado em matemática (e não elemento a elemento como em Python). Este comportamento é mais próximo novamente do MATLAB/Octave:

A * A
2×2 Array{Float64,2}:
  7.0  10.0
 15.0  22.0

A operação acima foi, portanto, uma multiplicação matricial (veja detalhes a seguir). Para forçarmos operações elemento a elemento, usamos o .:

A .* A
2×2 Array{Float64,2}:
 1.0   4.0
 9.0  16.0

Multiplicação matricial

Como vimos, a multiplicação matricial é a operação padrão * (ao contrário do Python, por exemplo, que usa @):

A * A
2×2 Array{Float64,2}:
  7.0  10.0
 15.0  22.0

Vamos calcular 2AI+AB2 A - I + A B, onde

A=[1234]B=[3140]A = \begin{bmatrix}1 & 2\\3 & 4\end{bmatrix}\quad B = \begin{bmatrix}3 & 1\\4 & 0\end{bmatrix}

com II sendo a matriz identidade vista anteriormente na parte de matrizes bloco:

A matriz A já foi definida anteriormente, basta definirmos a matriz B:

A
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
B = [3. 1.
     4. 0.]
2×2 Array{Float64,2}:
 3.0  1.0
 4.0  0.0
2.0 * A - I + A * B
2×2 Array{Float64,2}:
 12.0   5.0
 31.0  10.0

Produto interno entre vetores

O produto interno entre vetores é simplesmente uma multiplicação matricial entre um vetor-linha e um vetor-coluna. Uma vez que representamos vetores como colunas por padrão, vamos precisar de uma transposição (ou adjunto, veja a seguir). Ele nos revela que os uu e vv definidos anteriormente são ortogonais:

u
2-element Array{Float64,1}:
  1.0
 -1.0
v
2-element Array{Float64,1}:
 1.0
 1.0
u' * v
0.0

Em situações como esta, em que o vetor a esquerda está sendo transposto, a multiplicação pode ser dada implicitamente pois é entendida do contexto:

u'v
0.0

O produto interno com vetores complexos é levemente diferente, uma vez que se usa o complexo conjugado do vetor à esquerda (i=0nuivi\sum_{i = 0}^n u^*_i v_i). Primeiro, vamos definir um vetor complexo (a unidade complexa em Julia é im):

w = [1. + 2im, -1.]
2-element Array{Complex{Float64},1}:
  1.0 + 2.0im
 -1.0 + 0.0im
w'u
2.0 - 2.0im

O valor acima está correto! Acontece que o operador ' significa adjunto para vetores e matrizes complexas (observe os negativos nos complexos):

w'
1×2 LinearAlgebra.Adjoint{Complex{Float64},Array{Complex{Float64},1}}:
 1.0-2.0im  -1.0-0.0im

A função dot do pacote LinearAlgebra também calcula produtos internos:

dot(w, u) ≈ w'u
true

De fato, essa função costuma ser mais rápida que a multiplicação acima. Compare os tempos de execução de ambos os métodos:

using BenchmarkTools  # instale usando "] add BenchmarkTools"

@benchmark dot(w, u)
BenchmarkTools.Trial: 
  memory estimate:  32 bytes
  allocs estimate:  1
  --------------
  minimum time:     24.498 ns (0.00% GC)
  median time:      31.728 ns (0.00% GC)
  mean time:        36.622 ns (10.36% GC)
  maximum time:     5.577 μs (99.03% GC)
  --------------
  samples:          10000
  evals/sample:     996
@benchmark w'u
BenchmarkTools.Trial: 
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     41.254 ns (0.00% GC)
  median time:      51.668 ns (0.00% GC)
  mean time:        61.376 ns (10.18% GC)
  maximum time:     7.851 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     989

Norma de vetores

A norma euclidiana pode ser calculada com uma função norm (presente no pacote LinearAlgebra):

norm(u)
1.4142135623730951

Observe que a norma euclidiana é a raiz quadrada (sqrt) do produto interno de um vetor com ele próprio:

sqrt(u' * u)
1.4142135623730951

O símbolo indica a mesma função:

√(u' * u)
1.4142135623730951

Em Julia, matemática parece matemática.

Potenciação matricial

Ao contrário do Python, onde não existe um símbolo para potenciação matricial (o operador de potenciação ** em Python é expecificamente potenciação elemento a elemento), em Julia as coisas "simplesmente funcionam" (observe que o operador de potenciação em Julia é ^:

A ^ 3
2×2 Array{Float64,2}:
 37.0   54.0
 81.0  118.0

Naturalmente, a potenciação é equivalente a multiplicações sequenciais:

A * A * A
2×2 Array{Float64,2}:
 37.0   54.0
 81.0  118.0

No entanto, não tente fazer A400A^{400} multiplicação por multiplicação!

A ^ 400
2×2 Array{Float64,2}:
 2.76493e291  4.02969e291
 6.04454e291  8.80947e291

Aliás, internamente a potenciação não é feita multiplicação por multiplicação (em um loop, por exemplo) por outro motivo que não a conveniência: eficiência. Para entender melhor, compare os tempos de execução abaixo:

@benchmark A ^ 50
BenchmarkTools.Trial: 
  memory estimate:  784 bytes
  allocs estimate:  7
  --------------
  minimum time:     273.820 ns (0.00% GC)
  median time:      394.421 ns (0.00% GC)
  mean time:        453.472 ns (11.97% GC)
  maximum time:     15.810 μs (97.02% GC)
  --------------
  samples:          10000
  evals/sample:     233
@benchmark begin
    (A * A * A * A * A * A * A * A * A * A *
     A * A * A * A * A * A * A * A * A * A *
     A * A * A * A * A * A * A * A * A * A *
     A * A * A * A * A * A * A * A * A * A *
     A * A * A * A * A * A * A * A * A * A)
end
BenchmarkTools.Trial: 
  memory estimate:  5.77 KiB
  allocs estimate:  51
  --------------
  minimum time:     3.525 μs (0.00% GC)
  median time:      4.425 μs (0.00% GC)
  mean time:        5.074 μs (8.78% GC)
  maximum time:     472.227 μs (98.66% GC)
  --------------
  samples:          10000
  evals/sample:     8

Observe que no segundo caso o tempo de execução médio é de cerca de 10×10\times maior.

Transposição de matrizes

Como indicado antes, podemos transpor uma matriz com o operador ', o que representa adjuntos (transpostos com complexos conjugados) para arrays de números complexos:

A'
2×2 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}:
 1.0  3.0
 2.0  4.0

Observe que os números fora da diagonal trocaram de posição com relação à matriz original:

A
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

É comum usarmos isto para produzir matrizes simétricas, usando AATA A^T, que produz uma matriz simétrica para qualquer AA:

A * A'
2×2 Array{Float64,2}:
  5.0  11.0
 11.0  25.0

Inversão de matrizes

Matrizes podem ser invertidas com inv:

Ai = inv(A)
2×2 Array{Float64,2}:
 -2.0   1.0
  1.5  -0.5

Obviamente, a matriz inversa A1A^{-1} é tal que

A1A=AA1=I,A^{-1} A = A A^{-1} = I \text{,}

onde II é a matriz identidade

I=[1001].I = \begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}\text{.}
Ai * A
2×2 Array{Float64,2}:
 1.0          0.0
 2.22045e-16  1.0

A menos de erros numéricos, essa parece ser a resposta correta (valores da ordem de 101610^{-16} são praticamente zero dadas as limitações na representação de pontos flutuantes). Vamos checar isso por comparação da multiplicação acima com a matriz identidade (usando a função isapprox através de seu alias , que compara arrays a menos de erros numéricos):

Ai * A ≈ I
true

Em Julia, matemática parece matemática.

Lembre-se: nem toda matriz possui inversa! Uma maneira de identificar se uma matriz possui inversa é com o cálculo do rank da mesma (rank, disponível de LinarAlgebra), que deve ser igual à dimensão do array (para o caso de matrizes quadradas, claro):

rank(A)
2
ndims(A)
2

Traço de matrizes

O traço de uma matriz é a soma de seus valores diagonais. Ele pode ser calculado com tr, disponível de LinearAlgebra:

A
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

Para o caso da matriz acima, o traço deve ser 1+4=51 + 4 = 5:

tr(A)
5.0

Determinantes de matrizes

Determinantes são calculados com det de LinearAlgebra:

det_A = det(A)
-2.0

A matriz inversa tem a propriedade de que det(A1)=det(A)1\det(A^{-1}) = \det(A)^{-1}:

det_Ai = det(Ai)
-0.49999999999999967

De fato, podemos verificar det(A)det(A1)=1\det(A) \det(A^{-1}) = 1 a menos de erros numéricos com :

det_A * det_Ai ≈ 1.0
true

É isso aí, no próximo texto vamos falar um pouco da resolução de sistemas de equações lineares.