A = rand(1:4, 3, 3)
3×3 Matrix{Int64}: 4 1 2 1 1 4 1 2 1
Agora, vamos construir um vetor com 3 coordenadas todas iguais a 1:
v = fill(1.0, (3))
3-element Vector{Float64}: 1.0 1.0 1.0
b = A*v
3-element Vector{Float64}: 7.0 6.0 4.0
A matriz transpose(A)
é a transposta de A
:
transpose(A)
3×3 transpose(::Matrix{Int64}) with eltype Int64: 4 1 1 1 1 2 2 4 1
A matriz A'
é a adjunta de A
, ou seja, a transposta conjugada de A
, que coincide com a transposta quando A
é real.
A'
3×3 adjoint(::Matrix{Int64}) with eltype Int64: 4 1 1 1 1 2 2 4 1
Em Julia, podemos escrever simplesmente A'A
em vez de (A')*A
:
A'A
3×3 Matrix{Int64}: 18 7 13 7 6 8 13 8 21
(A')*A
3×3 Matrix{Int64}: 18 7 13 7 6 8 13 8 21
O problema $Ax = b$ pode ser resolvido usando a função \
.
x = A\b
3-element Vector{Float64}: 1.0 1.0 1.0
using LinearAlgebra
n = 1000
A = randn(n,n);
Frequentemente, é possível inferir a estrutura de uma matriz:
Asim = A + A'
issymmetric(Asim)
true
Todavia, muitas vezes os arredondamentos numéricos atrapalham:
Asim_ruido = copy(Asim)
ruido = 5eps()
println("ruido = $ruido")
Asim_ruido[1,2] += ruido
issymmetric(Asim_ruido)
ruido = 1.1102230246251565e-15
false
O comando acima inseriu um ruído da ordem de 10-15 no elemento (1,2) da matriz Asim
.
A linguagem Julia permite que declaremos explicitamente a estrutura de uma matriz usando, por exemplo, Diagonal
, Triangular
, Symmetric
, Hermitian
, Tridiagonal
e SymTridiagonal
.
Asim_explicita = Symmetric(Asim_ruido)
issymmetric(Asim_explicita)
true
Vamos comparar quanto tempo leva para a Julia calcular os autovalores de Asim
, Asim_ruido
e Asim_explicita
.
@time eigvals(Asim);
0.473789 seconds (697.98 k allocations: 43.036 MiB, 15.23% gc time, 85.25% compilation time)
@time eigvals(Asim_ruido);
0.608596 seconds (13 allocations: 7.920 MiB)
@time eigvals(Asim_explicita);
0.088237 seconds (13.74 k allocations: 8.711 MiB, 17.63% compilation time)
Observamos que o cálculo com Asim_explicita
(onde usamos Symmetric()
) foi muitas vezes mais rápido do que o cálculo com Asim_ruido
.
Usando os tipos Tridiagonal
e SymTridiagonal
, podemos armazenar matrizes de forma eficiente.
Isso torna possível fazer cálculos com matrizes muito grandes.
Por exemplo, não seria possível resolver em um laptop o seguinte problema se armazenássemos a matriz usando o tipo Array
.
n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)
0.752085 seconds (880.41 k allocations: 226.673 MiB, 2.62% gc time, 26.47% compilation time)
6.569112206378978
O cálculo acima demorou menos de 1 segundo. Explicamos a seguir o significado dos comandos.
O comando abaixo cria uma matriz simétrica 6x6 especificando 6 números aleatórios na diagonal e 5 números aleatórios na diagonal secundária.
B = SymTridiagonal(randn(6), randn(6-1))
6×6 SymTridiagonal{Float64, Vector{Float64}}: 0.148463 -1.29177 ⋅ ⋅ ⋅ ⋅ -1.29177 -0.190034 -2.01804 ⋅ ⋅ ⋅ ⋅ -2.01804 -0.557058 1.007 ⋅ ⋅ ⋅ ⋅ 1.007 0.642556 -0.677129 ⋅ ⋅ ⋅ ⋅ -0.677129 -1.56873 0.446176 ⋅ ⋅ ⋅ ⋅ 0.446176 1.09846
Portanto, a matriz A
acima tem estrutura semelhante mas é de ordem 1 milhão x 1 milhão:
A
1000000×1000000 SymTridiagonal{Float64, Vector{Float64}}: 0.0255158 0.702558 ⋅ … ⋅ ⋅ ⋅ 0.702558 -1.05006 -0.635067 ⋅ ⋅ ⋅ ⋅ -0.635067 -1.22556 ⋅ ⋅ ⋅ ⋅ ⋅ 0.510417 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋮ ⋱ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ -2.3868 ⋅ ⋅ ⋅ ⋅ ⋅ -0.0818654 -0.669316 ⋅ ⋅ ⋅ ⋅ -0.669316 1.27233 0.809435 ⋅ ⋅ ⋅ ⋅ 0.809435 -0.0141998
O comando randn(6)
cria uma lista de 6 números aleatórios.
randn(6)
6-element Vector{Float64}: 0.856298177874976 -1.08001047780381 -1.0546922643723062 -0.47569586118087903 -2.3399979251248864 -0.05306667490614464
O comando eigmax(A)
calcula o maior autovalor de A
.
Com os comandos a seguir, podemos ler as documentações sobre as funções (lembre da dica de traduzir usando https://translate.google.com):
?eigmax
search: eigmax
eigmax(A; permute::Bool=true, scale::Bool=true)
Return the largest eigenvalue of A
. The option permute=true
permutes the matrix to become closer to upper triangular, and scale=true
scales the matrix by its diagonal elements to make rows and columns more equal in norm. Note that if the eigenvalues of A
are complex, this method will fail, since complex numbers cannot be sorted.
jldoctest
julia> A = [0 im; -im 0]
2×2 Matrix{Complex{Int64}}:
0+0im 0+1im
0-1im 0+0im
julia> eigmax(A)
1.0
julia> A = [0 im; -1 0]
2×2 Matrix{Complex{Int64}}:
0+0im 0+1im
-1+0im 0+0im
julia> eigmax(A)
ERROR: DomainError with Complex{Int64}[0+0im 0+1im; -1+0im 0+0im]:
`A` cannot have complex eigenvalues.
Stacktrace:
[...]
?randn
search: randn lowrankdowndate lowrankdowndate! RankDeficientException rand
randn([rng=GLOBAL_RNG], [T=Float64], [dims...])
Generate a normally-distributed random number of type T
with mean 0 and standard deviation 1. Optionally generate an array of normally-distributed random numbers. The Base
module currently provides an implementation for the types Float16
, Float32
, and Float64
(the default), and their Complex
counterparts. When the type argument is complex, the values are drawn from the circularly symmetric complex normal distribution of variance 1 (corresponding to real and imaginary part having independent normal distribution with mean zero and variance 1/2
).
jldoctest
julia> using Random
julia> rng = MersenneTwister(1234);
julia> randn(rng, ComplexF64)
0.6133070881429037 - 0.6376291670853887im
julia> randn(rng, ComplexF32, (2, 3))
2×3 Matrix{ComplexF32}:
-0.349649-0.638457im 0.376756-0.192146im -0.396334-0.0136413im
0.611224+1.56403im 0.355204-0.365563im 0.0905552+1.31012im