FIT ČVUT Course Pages
BI-JUL — Programování v jazyku Julia
Jdi na navigaci předmětu

06: Standardní knihovna (Lineární algebra, Statistika,...)

Tento notebook je výukovým materiálem v předmětu BI-JUL.21 vyučovaném v zimním semestru akademického roku 2024/2025 Tomášem Kalvodou. Tvorba těchto materiálů byla podpořena NVS FIT.

Hlavní stránkou předmětu, kde jsou i další notebooky a zajímavé informace, je jeho Course Pages stránka.

versioninfo()

1 Standardní knihovna

Postupem času vykrystalizovalo několik balíčku, resp. modulů, které jsou distribuovány společně s Julia a tvoří tzv. Standardní knihovnu. Pro kompletní přehled a popis viz oficiální dokumentaci, část Standard Library.

V této lekci projdeme jen část. Naší snahou je vybrat ty nejzajímavější, které mohou být uživatelům nejčastěji k užitku. Během výkladu jejich funkcionalitu také využijeme při řešení různých úkolů. Následující výklad postupuje v abecedním pořádku dle názvu balíčku/modulu. Některé partie jsme už nakousli dříve, ale raději je zde na tomto místě přehledně zmiňuji. Na cvičení projdeme podrobněji jen některé.

Studentům doporučuji alespoň zběžně balíčky ze Standardní knihovny proklikat a získat tak povědomí o možnostech, které jsou k dispozici.

1.1 Dates

Zejména při zpracování různých reálných dat a měření je potřeba pracovat s časem a daty (ve smyslu časovém).

Dva základní konstruktory takovýchto údajů poskytované modulem Dates jsou

  • Date: pouze datum (rok, měsíc, den)
  • DateTime: datum i časový okamžík

Nejprve ovšem musíme modul importovat.

using Dates

Při tvorbě datumu stačí vždy zadat pouze část údajů:

Date(2025)
2025-01-01
Date(2025, 10)
2025-10-01
d = Date(2025, 10, 30)
2025-10-30
typeof(d)
Date

K různým částem daného údaje můžeme přistupovat pomocí k tomu určených metod.

day(d)
30
month(d)
10
year(d)
2025
Date(2024, 18, 345)
ArgumentError: Month: 18 out of range (1:12)

Stacktrace:
 [1] Date(y::Int64, m::Int64, d::Int64)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/types.jl:277
 [2] top-level scope
   @ In[10]:1
 [3] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
Date(2024, 11, 345)
ArgumentError: Day: 345 out of range (1:30)

Stacktrace:
 [1] Date(y::Int64, m::Int64, d::Int64)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/types.jl:277
 [2] top-level scope
   @ In[11]:1
 [3] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

Letos tento den nebyl.

Date(2025, 2, 29)
ArgumentError: Day: 29 out of range (1:28)

Stacktrace:
 [1] Date(y::Int64, m::Int64, d::Int64)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/types.jl:277
 [2] top-level scope
   @ In[17]:1
 [3] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

Ale v roce -4 ano!

Date(-4, 2, 29)
-0004-02-29

Při parsování časových údajů se může velmi hodit mít možnost specifikovat, v jakém formátu datum je. K tomu slouží DateFormat:

df = DateFormat("d. m. y")
dateformat"d. m. y"

Ekvivalentně lze tento formát vytvořit pomocí řetězce tohoto "typu":

dateformat"d. m. y." == DateFormat("d. m. y.")
true

Poté můžeme zadávat data i v tomto formátu:

Date("2024-10-22")
2024-10-22
Date("12. 3. 1999")
ArgumentError: Unable to parse date time. Expected directive Delim(-) at char 3

Stacktrace:
 [1] macro expansion
   @ /usr/share/julia/stdlib/v1.12/Dates/src/parse.jl:104 [inlined]
 [2] tryparsenext_core(str::String, pos::Int64, len::Int64, df::DateFormat{Symbol("yyyy-mm-dd"), Tuple{Dates.DatePart{'y'}, Dates.Delim{Char, 1}, Dates.DatePart{'m'}, Dates.Delim{Char, 1}, Dates.DatePart{'d'}}}, raise::Bool)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/parse.jl:38
 [3] macro expansion
   @ /usr/share/julia/stdlib/v1.12/Dates/src/parse.jl:150 [inlined]
 [4] tryparsenext_internal
   @ /usr/share/julia/stdlib/v1.12/Dates/src/parse.jl:125 [inlined]
 [5] parse(::Type{Date}, str::String, df::DateFormat{Symbol("yyyy-mm-dd"), Tuple{Dates.DatePart{'y'}, Dates.Delim{Char, 1}, Dates.DatePart{'m'}, Dates.Delim{Char, 1}, Dates.DatePart{'d'}}})
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/parse.jl:304
 [6] Date
   @ /usr/share/julia/stdlib/v1.12/Dates/src/io.jl:626 [inlined]
 [7] Date(d::String)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/io.jl:626
 [8] top-level scope
   @ In[22]:1
 [9] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
Date("12. 3. 1999", df)
1999-03-12

Časové údaje lze přirozeně porovnávat a operovat s nimi.

d1 = Date(2025, 1, 1)
d2 = Date(2025, 1, 7)
2025-01-07
d1 < d2
true
d2 - d1
6 days
typeof(d2 - d1)
Day
d3 = Date(55, 1, 1)
0055-01-01
d2 - d3
719534 days
d3 - d2
-719534 days
d2 - d2
0 days
Day(12)
12 days

Například ale součet už dobře definovaný není.

d2 + d1
MethodError: no method matching +(::Date, ::Date)
The function `+` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:642
  +(::Missing, ::Dates.AbstractTime)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/arithmetic.jl:91
  +(::Time, ::Date)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/arithmetic.jl:22
  ...


Stacktrace:
 [1] top-level scope
   @ In[38]:1
 [2] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
d2
2025-01-07
d2 + 10
MethodError: no method matching +(::Date, ::Int64)
The function `+` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:642
  +(::BigInt, ::Union{Int16, Int32, Int64, Int8})
   @ Base gmp.jl:556
  +(::Missing, ::Number)
   @ Base missing.jl:123
  ...


Stacktrace:
 [1] top-level scope
   @ In[40]:1
 [2] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
d2 + Day(10)
2025-01-17
d2 + Month(1)
2025-02-07
d2 - Year(10)
2015-01-07
d2 - Day(123)
2024-09-06

Naprosto analogicky se chová i metoda DateTime.

jul = DateTime(2025, 10, 30, 11, 25, 55)
2025-10-30T11:25:55
hour(jul)
11
minute(jul)
25
second(jul)
55

Výpočty s daty lze provádět poměrně snadno:

jul + Day(14) - Hour(10)
2025-11-13T01:25:55

Co změna času?

DateTime(2025, 10, 25, 23, 0, 0) + Hour(5)
2025-10-26T04:00:00

Zde není brána v potaz.

Pomocné metody

Mám k dispozici i další užitečné metody jako

jul
2025-10-30T11:25:55
dayofweek(jul)
4
dayname(jul)
"Thursday"
monthname(jul)
"October"
dayofmonth(jul)
30
dayofyear(jul)
303
# pořadí týdne v roce?
week(jul)
44

Lokalizace názvů (dny, měsíce) je možná, je ale nutné ji přidat do (viz Dates.DateLocale a dokumentaci):

Dates.LOCALES
Dict{String, Dates.DateLocale} with 1 entry:
  "english" => DateLocale(["January", "February", "March", "April", "May", "Jun…

Dates se chová naivně vůči časovým zónám, ignoruje je. Pokud potřebujete pracovat s časovými zónami, můžete použít balížek TimeZones.jl, který ho o tuto funkcinalitu rozšiřuje.

Aktuální čas:

today()
2025-10-30
now()
2025-10-30T11:32:20.615

Cvičení

  1. Kolik minut uplynulo od výbuchu jaderného reaktoru v Černobylu?
  2. Který den v týdnu se narodil Isaac Newton?
  3. Jak provádět převody mezi různými časovými intervaly? Např. kolik minut je v sedmi dnech (Day(7))?
#1
explosion = DateTime(1986, 4, 26, 1, 24, 40) - Hour(3)
t = now() - explosion
1246972077415 milliseconds
typeof(t)
Millisecond
Second(120)
Millisecond(120)
Nanosecond(40)
round(t, Minute)
20782868 minutes
round(t, Day)
14433 days
Year(1)
t.value
t.value / 60_000
#2
birthday = Date(1642, 12, 25)
1642-12-25
dayname(birthday)
"Thursday"
dayname(Date(2025, 12, 25))
"Thursday"
#3
Minute(Day(7))
10080 minutes

Pokusy.

Hour(12) + Hour(12) == Day(1)
true
Hour(12) + Hour(12)
24 hours
Hour(Day(1))
24 hours
Day(0.5)
InexactError: Int64(0.5)

Stacktrace:
 [1] Int64
   @ ./float.jl:923 [inlined]
 [2] convert
   @ ./number.jl:7 [inlined]
 [3] Day(v::Float64)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/types.jl:55
 [4] top-level scope
   @ In[81]:1
 [5] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
Day(1) / Hour(2)
12.0
Day(2) / 2
1 day
Hour(4) / 2
2 hours
2 * Day(2)
4 days
Day(2) * Day(2)
MethodError: no method matching *(::Day, ::Day)
The function `*` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:642
  *(::Real, ::Period)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/periods.jl:91
  *(::AbstractRange{<:Real}, ::Period)
   @ Dates /usr/share/julia/stdlib/v1.12/Dates/src/ranges.jl:68
  ...


Stacktrace:
 [1] top-level scope
   @ In[86]:1
 [2] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

1.2 Linear Algebra

V předchozí kapitole jsme se bavili o vícerozměrných polích. Dívali jsme se na ně "abstraktně", nebavili jsme se o tom, co ve skutečnosti mohou popisovat (matice lineárního operátoru, incidenční matice grafu, obrázek koťátka...).

Jak název napovídá, modul Linear Algebra exportuje metody pro práci s maticemi (případně tenzory). V této části lekce probereme ty nejzásadnější, které jsou většinou studentům dobře známy, nebo by jim mohly být v budoucnu užitečné.

Neprobereme zdaleka všechna zákoutí tohoto modulu. Například se zcela vyhneme sekci Low-level matrix operations, která se ovšem může velmi hodit při optimalizaci maticových výpočtů (vyhnutí se alokaci zbytečných mezivýsledků, provádění operací "na místě", atp.).

Nakonec tohoto úvodu ještě poznamenejme, že Julia na spoustu maticových operací ve strojových číslech přirozeně používá BLAS/LAPACK knihovny. Čili opět stojíme na ramemou FORTRANu...

using LinearAlgebra

Maticové typy

Typ Array{T, 2}, resp Matrix, odpovídá obecné matici bez jakékoliv struktury, kterou by bylo možné využít k uložení matice nebo při provádění operací s ní. V modulu LinearAlgebra máme k dispozici hned několik dalších užitečných maticových typů, zmiňme alespoň ty nejdůležitější a studentům většinou známé:

  • Symmetric: symetrická matice (AT=A\mathbb{A}^T = \mathbb{A}),
  • Hermitian: hermitovská matice (A†=A\mathbb{A}^\dagger = \mathbb{A}),
  • UpperTriangular, LowerTriangular: horní a dolní trojúhelníková matice,
  • Diagonal: diagonální matice,
  • Tridiagonal, SymTridiagonal: tridiagonální a symetrická tridiagonální matice.

Pokud víme, že naše matice má některou z těchto struktur, tak je výhodné tyto typy využít. Získáme tím například efektivnější algoritmy pro výpočet některých vlastností (maticové funkce, vlastní čísla...), spousta metod je pro tyto typy optimalizovaná.

Maticové konstruktory I, diagm

V Lineární algebře velmi často pracujeme s jednotkovou maticí. Konstanta I reprezentuje jednotkovou matici správné velikosti a typu.

I
UniformScaling{Bool}
true*I
2.0 * I
UniformScaling{Float64}
2.0*I
1 // 2 * I
UniformScaling{Rational{Int64}}
1//2*I
[1 1; 1 1] + I
2×2 Matrix{Int64}:
 2  1
 1  2
I * [1.0, 2.0, 3.0]
3-element Vector{Float64}:
 1.0
 2.0
 3.0
I
UniformScaling{Bool}
true*I
I.λ
true
moje_I = 3.3 * I
UniformScaling{Float64}
3.3*I
moje_I.λ
3.3
@which I * [1, 2, 3]
*(J::UniformScaling, A::AbstractVecOrMat) in LinearAlgebra at /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/uniformscaling.jl:268

Pokud chceme vytvořit jednotkovou matici konkrétního rozměru, stačí tento předat metodě I: 

I(4)
4×4 Diagonal{Bool, Vector{Bool}}:
 1  ⋅  ⋅  ⋅
 ⋅  1  ⋅  ⋅
 ⋅  ⋅  1  ⋅
 ⋅  ⋅  ⋅  1
I(4) + zeros(4, 4)
4×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
I(5) * zeros(4, 4)
DimensionMismatch: incompatible dimensions for matrix multiplication: tried to multiply a matrix of size (5, 5) with a matrix of size (4, 4). The second dimension of the first matrix: 5, does not match the first dimension of the second matrix: 4.

Stacktrace:
  [1] matmul_size_check_error(sizeA::Tuple{Int64, Int64}, sizeB::Tuple{Int64, Int64})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:431
  [2] matmul_size_check
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:418 [inlined]
  [3] matmul_size_check
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:441 [inlined]
  [4] _mul_diag!(out::Matrix{Float64}, A::Diagonal{Bool, Vector{Bool}}, B::Matrix{Float64}, alpha::Bool, beta::Bool)
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/diagonal.jl:504
  [5] _mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/diagonal.jl:519 [inlined]
  [6] mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:297 [inlined]
  [7] mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:265 [inlined]
  [8] mul
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:118 [inlined]
  [9] *(A::Diagonal{Bool, Vector{Bool}}, B::Matrix{Float64})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:114
 [10] top-level scope
    @ In[101]:1
 [11] eval(m::Module, e::Any)
    @ Core ./boot.jl:489

Typ Diagonal je speciální maticový typ, který v paměti ukládá pouze prvky na diagonále ("řídká" matice).

Diagonal([1, 2, 3])
3×3 Diagonal{Int64, Vector{Int64}}:
 1  ⋅  ⋅
 ⋅  2  ⋅
 ⋅  ⋅  3

Funkce diagm nám umožňuje konstruovat matice zadáváním hlavních a vedlejších diagonál. Vždy ovšem konstruuje hustou matici (tj. i nuly jsou v paměti uloženy). Tato funkce má hned několik metod, zmiňme alespoň ty nejužitečnější:

diagm([1, 2, 3, 4])
4×4 Matrix{Int64}:
 1  0  0  0
 0  2  0  0
 0  0  3  0
 0  0  0  4
diagm(1 => [1, 2, 3], -2 => [2, 3])
4×4 Matrix{Int64}:
 0  1  0  0
 0  0  2  0
 2  0  0  3
 0  3  0  0

Výsledná matice nemusí být nutně čtvercová.

diagm(3, 4, [1, 2, 3])
3×4 Matrix{Int64}:
 1  0  0  0
 0  2  0  0
 0  0  3  0
diagm(3, 4, 0 => [1, 2, 3], 1 => [-1, -2, -3])
3×4 Matrix{Int64}:
 1  -1   0   0
 0   2  -2   0
 0   0   3  -3

Dále máme k dispozici typ Diagonal. Porovnejte paměťovou náročnost pomocí Base.summarysize.

Maticové operace tr, det, inv, diag, rank, transpose, adjoint

Operace budeme demonstrovat na matici

A = [
    -1  1  1
     1 -1  1
     1  1 -1
]
3×3 Matrix{Int64}:
 -1   1   1
  1  -1   1
  1   1  -1

Determinant matice:

det(A)
4.0
det([1//2 0; 0 1//3])
1//6

Stopa matice, tr A=∑j=1nAjj\mathrm{tr}\, A = \sum_{j=1}^n A_{jj}:

tr(A)
-3

Inverze matice:

inv(A)
3×3 Matrix{Float64}:
 0.0   0.5  0.5
 0.5  -0.0  0.5
 0.5   0.5  0.0
A * inv(A)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Přirozeně lze používat i mocnění matice se záporným exponentem.

A^(-1)
3×3 Matrix{Float64}:
 0.0   0.5  0.5
 0.5  -0.0  0.5
 0.5   0.5  0.0

Pro singulární, nebo nečtvercovou, matici tato operace přirozeně selže.

B = [1 -1; -1 1]

inv(B)
SingularException(2)

Stacktrace:
  [1] checknonsingular
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/factorization.jl:69 [inlined]
  [2] _check_lu_success
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:84 [inlined]
  [3] #lu!#180
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:92 [inlined]
  [4] lu!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:90 [inlined]
  [5] lu!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:89 [inlined]
  [6] _lu
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:347 [inlined]
  [7] lu(::Matrix{Int64}; kwargs::@Kwargs{})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:341
  [8] lu
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/lu.jl:341 [inlined]
  [9] inv(A::Matrix{Int64})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/dense.jl:1087
 [10] top-level scope
    @ In[115]:3
 [11] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
inv(A[1:2, :])
DimensionMismatch: matrix is not square: dimensions are (2, 3)

Stacktrace:
 [1] checksquare
   @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/LinearAlgebra.jl:342 [inlined]
 [2] inv(A::Matrix{Int64})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/dense.jl:1079
 [3] top-level scope
   @ In[116]:1
 [4] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

Máme k dispozici i pseudoinverzi pinv.

Dále zmiňme metodu diag(A, k=0), které z matice A extrahuje kk-tou "diagonálu" (00 odpovídá skutečné diagonále, 11 nad diagonálou, −1-1 pod diagonálou, atd.). 

C = [
    1 2 3 4
    5 6 7 8
    0 0 0 0
    1 1 1 1
]
4×4 Matrix{Int64}:
 1  2  3  4
 5  6  7  8
 0  0  0  0
 1  1  1  1
diag(C)
4-element Vector{Int64}:
 1
 6
 0
 1
diag(C[:,1:2])
2-element Vector{Int64}:
 1
 6
diag(C, 1)
3-element Vector{Int64}:
 2
 7
 0
diag(C, -1)
3-element Vector{Int64}:
 5
 0
 1

Funguje i s nečtvercovou maticí.

C = [
    1 2 3
    4 5 6
    ]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
diag(C,  -2)
Int64[]

Hodnost matice lze počítat pomocí metody rank. Pozor, výpočet hodnosti v neexaktní aritmetice není zcela triviální a je zatížen numerickou chybou.

A
3×3 Matrix{Int64}:
 -1   1   1
  1  -1   1
  1   1  -1
rank(A)
3
B
2×2 Matrix{Int64}:
  1  -1
 -1   1
rank(B)
1
rank(C)
2

Metoda transpose implementuje línou transpozici matice. K dispozici máme i transpose!, která mění matici na které transpozici provádíme.

C
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
transpose(C)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

Vedle toho adjoint představuje sdružení matice. Také má verzi adoint! a navíc ji můžeme velmi zkráceně zapisovat pomocí '.

D = [1 2im; 3 4im]
2×2 Matrix{Complex{Int64}}:
 1+0im  0+2im
 3+0im  0+4im
adjoint(D)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1+0im  3+0im
 0-2im  0-4im
C'
3×2 adjoint(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6
D'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1+0im  3+0im
 0-2im  0-4im

Cvičení

  1. Experimentálně ověřte vztah det⁡(eA)=etr A\det(e^\mathbb{A}) = e^{\mathrm{tr}\, \mathbb{A}} platný pro čtvercovou matici A\mathbb{A}.
A = rand(7, 7)
7×7 Matrix{Float64}:
 0.82423    0.662486  0.693849   0.886411    0.530206  0.360918  0.691788
 0.114072   0.446791  0.179846   0.00105323  0.460695  0.450705  0.522562
 0.639288   0.709996  0.843532   0.435651    0.81695   0.14675   0.261625
 0.798997   0.721967  0.554162   0.43418     0.383171  0.291747  0.314134
 0.44236    0.050238  0.171674   0.46302     0.804846  0.792962  0.263198
 0.831153   0.554976  0.0147654  0.00304053  0.349219  0.515718  0.476831
 0.0629234  0.154115  0.869568   0.116271    0.116345  0.326194  0.0512608
exp(A)
7×7 Matrix{Float64}:
 6.66513  4.99594  5.12035  4.16409  5.17549  4.07935  4.13561
 1.92378  3.00925  1.78014  1.11732  2.28307  2.06631  1.86353
 4.69291  4.28009  5.4583   3.23615  4.91552  3.37682  3.17403
 4.47808  4.03373  3.8622   3.94166  3.96407  3.15437  3.01865
 3.72466  2.72588  2.73197  2.55296  4.78074  3.369    2.50863
 3.53144  2.9186   2.38588  1.8079   2.88933  3.72913  2.53363
 1.96754  1.87262  2.56549  1.33532  2.03104  1.70261  2.34571
exp.(A)
7×7 Matrix{Float64}:
 2.28012  1.93961  2.0014   2.4264   1.69928  1.43465  1.99728
 1.12083  1.56329  1.19703  1.00105  1.58518  1.56942  1.68634
 1.89513  2.03398  2.32456  1.54597  2.26359  1.15806  1.29904
 2.22331  2.05848  1.74048  1.5437   1.46693  1.33876  1.36907
 1.55638  1.05152  1.18729  1.58886  2.23635  2.20993  1.30108
 2.29596  1.7419   1.01487  1.00305  1.41796  1.67484  1.61096
 1.06495  1.16663  2.38588  1.1233   1.12338  1.38568  1.0526
det(exp(A))
50.428509108251
tr(A)
3.9205566720635003
exp(tr(A))
50.428509108251

Vektorové operace dot, cross

Metoda dot implementuje standardní skalární součin v Cn\mathbb{C}^n, který je antilineární v první složce (tj. komplexně se sdružuje v první složce).

v1 = [1, 2, 3]
v2 = [3, 2, 1]
v3 = [1im, 2im, 3im];
dot(v1, v2)
10

Obyčejná hvězdička takto nefunguje!

v1 * v1
MethodError: no method matching *(::Vector{Int64}, ::Vector{Int64})
The function `*` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:642
  *(::Adjoint{<:Number, <:AbstractVector}, ::AbstractVector{<:Number})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/adjtrans.jl:501
  *(::LinearAlgebra.AbstractQ, ::AbstractVector)
   @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/abstractq.jl:180
  ...


Stacktrace:
 [1] top-level scope
   @ In[149]:1
 [2] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
sum(v1 .* v2)
10
dot(v1, v3)
0 + 14im
dot(v3, v1)
0 - 14im

Lze použít i \cdot.

v3 ⋅ v1
0 - 14im

Velmi často, například v kvantově-mechanických výpočtech, je potřeba vyhodnocovat skalární součiny tvaru ⟨x,Ay⟩\langle x, A y\rangle, kde AA je jistá zadaná matice. K tomu účelu slouží metoda dot(x, A, y), která tento výpočet provede efektivně bez ukládání mezivýsledku maticového násobení.

A = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9
dot(v1, A, v2)
204

Dále máme díky metodě cross k dispozici známý vektorový součin dvou třísložkových vektorů.

cross(v1, v2)
3-element Vector{Int64}:
 -4
  8
 -4

Nebo pomocí makra \times.

v1 × v2
3-element Vector{Int64}:
 -4
  8
 -4
[0, 2, 0] × [3, 0, 0]
3-element Vector{Int64}:
  0
  0
 -6

Nelze:

[1,2,3,4] × [4,5,6,7]
DimensionMismatch: cross product is only defined for vectors of length 3

Stacktrace:
 [1] cross(a::Vector{Int64}, b::Vector{Int64})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/generic.jl:413
 [2] top-level scope
   @ In[161]:1
 [3] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

Vlastní čísla matic

Výpočty vlastních čísel patří k velmi důležitým maticovým algoritmům. Uživatel nejčastěji využije následující metody:

  • eigen: vrací vlastní čísla i vektory,
  • eigvals: vrací vlastní čísla,
  • eigvecs: vrací vlastní vektory.

Dále jsou k dispozici metody:

  • eigmax a eigmin: největší a nejmenší vlastní číslo,
  • svdvals: singulární hodnoty,
  • svd: singulární rozklad.

Jako ukázku použijeme například matici:

A = [1. 0. 0.; 0. 0. 1.; 0. -1. 0.]
3×3 Matrix{Float64}:
 1.0   0.0  0.0
 0.0   0.0  1.0
 0.0  -1.0  0.0

Vlastní čísla:

eigvals(A)
3-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
 1.0 + 0.0im

Vlastní vektory:

eigvecs(A)
3×3 Matrix{ComplexF64}:
      0.0-0.0im            0.0+0.0im       1.0+0.0im
 0.707107-0.0im       0.707107+0.0im       0.0+0.0im
      0.0-0.707107im       0.0+0.707107im  0.0+0.0im

Nebo dohromady:

e = eigen(A)
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
3-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
 1.0 + 0.0im
vectors:
3×3 Matrix{ComplexF64}:
      0.0-0.0im            0.0+0.0im       1.0+0.0im
 0.707107-0.0im       0.707107+0.0im       0.0+0.0im
      0.0-0.707107im       0.0+0.707107im  0.0+0.0im
typeof(e)
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
e.values
3-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
 1.0 + 0.0im
e.vectors
3×3 Matrix{ComplexF64}:
      0.0-0.0im            0.0+0.0im       1.0+0.0im
 0.707107-0.0im       0.707107+0.0im       0.0+0.0im
      0.0-0.707107im       0.0+0.707107im  0.0+0.0im

Vlastní vektory jsou v matici uloženy po sloupcích, to znamená, že platí následující maticová rovnost:

A * e.vectors ≈ e.vectors * diagm(e.values)
true

Respektive jednotlivě:

for k = 1:3
    println(A * e.vectors[:, k] ≈ e.values[k] * e.vectors[:, k])
end
true
true
true

eigen(I(4))
Eigen{Float64, Bool, Diagonal{Float64, Vector{Float64}}, Vector{Bool}}
values:
4-element Vector{Bool}:
 1
 1
 1
 1
vectors:
4×4 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅   1.0
eigen(diagm([2, 1, 2]))
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 1.0
 2.0
 2.0
vectors:
3×3 Matrix{Float64}:
 0.0  1.0  0.0
 1.0  0.0  0.0
 0.0  0.0  1.0
eigen([0 1;-1 0])
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
vectors:
2×2 Matrix{ComplexF64}:
 0.707107-0.0im       0.707107+0.0im
      0.0-0.707107im       0.0+0.707107im
eigen([
        1 1
        0 1
      ])
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 1.0
 1.0
vectors:
2×2 Matrix{Float64}:
 1.0  -1.0
 0.0   2.22045e-16

Cvičení

  1. Na symetrické reálné tridiagonální matici většího rozměru porovnejte rychlost výpočtu vlastních vektorů a vlastních čísel pokud použijete "hustou" matici a matici vhodného typu.
  2. Pomocí "matice společnice" a výpočtu vlastních čísel nalezněte kořeny polynomu ∑j=020(−1)jxj\sum_{j=0}^{20} (-1)^j x^j.
#1
n = 1_000
x = rand(n - 1)

A1 = diagm(0 => rand(n), 1 => x, -1 => x)
A2 = SymTridiagonal(A1);
A1
1000×1000 Matrix{Float64}:
 0.169353  0.878686  0.0       0.0       …  0.0       0.0       0.0
 0.878686  0.754788  0.265703  0.0          0.0       0.0       0.0
 0.0       0.265703  0.586549  0.596859     0.0       0.0       0.0
 0.0       0.0       0.596859  0.85906      0.0       0.0       0.0
 0.0       0.0       0.0       0.275606     0.0       0.0       0.0
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 ⋮                                       ⋱                      
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.530605  0.0       0.0
 0.0       0.0       0.0       0.0          0.873151  0.655321  0.0
 0.0       0.0       0.0       0.0          0.655321  0.349781  0.577912
 0.0       0.0       0.0       0.0          0.0       0.577912  0.573991
A2
1000×1000 SymTridiagonal{Float64, Vector{Float64}}:
 0.169353  0.878686   ⋅         ⋅        …   ⋅         ⋅         ⋅ 
 0.878686  0.754788  0.265703   ⋅            ⋅         ⋅         ⋅ 
  ⋅        0.265703  0.586549  0.596859      ⋅         ⋅         ⋅ 
  ⋅         ⋅        0.596859  0.85906       ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅        0.275606      ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
 ⋮                                       ⋱                      
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        …   ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅           0.530605   ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅           0.873151  0.655321   ⋅ 
  ⋅         ⋅         ⋅         ⋅           0.655321  0.349781  0.577912
  ⋅         ⋅         ⋅         ⋅            ⋅        0.577912  0.573991
Base.summarysize(A1)
8000048
Base.summarysize(A2)
16088
using BenchmarkTools
@benchmark eigen($A1)
BenchmarkTools.Trial: 16 samples with 1 evaluation per sample.
 Range (min … max):  278.248 ms … 330.677 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     321.944 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   317.486 ms ±  14.960 ms  ┊ GC (mean ± σ):  0.13% ± 0.19%

  ▁                      ▁  ▁▁                 ▁ ▁  █▁  ▁ ▁▁ ██  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁██▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁█▁▁██▁▁█▁██▁██ ▁
  278 ms           Histogram: frequency by time          331 ms <

 Memory estimate: 15.61 MiB, allocs estimate: 21.
@benchmark eigen($A2)
BenchmarkTools.Trial: 41 samples with 1 evaluation per sample.
 Range (min … max):  118.469 ms … 141.110 ms  ┊ GC (min … max): 0.00% … 0.88%
 Time  (median):     123.553 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   124.330 ms ±   4.352 ms  ┊ GC (mean ± σ):  0.19% ± 0.34%

    ▁   █  ▄▄ ▁█ ▁▄ ▁     ▁                                      
  ▆▁█▆▁▁█▆▆██▆██▆██▁█▆▁▁▆▁█▆▆▁▆▁▁▁▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆▁▁▁▁▁▁▁▁▁▁▁▆ ▁
  118 ms           Histogram: frequency by time          141 ms <

 Memory estimate: 7.89 MiB, allocs estimate: 38.
#2
cs = [ (-1)^k for k=0:19 ]
P  = hcat(diagm(20, 19, -1 => ones(19)), -cs)
20×20 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  1.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     1.0  0.0  0.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  1.0  0.0  0.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  1.0  0.0  0.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  1.0  0.0  0.0   1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  1.0  0.0  -1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  1.0   1.0
roots = eigvals(P)
20-element Vector{ComplexF64}:
  -0.9555728057861391 - 0.2947551744109035im
  -0.9555728057861391 + 0.2947551744109035im
  -0.8262387743159954 - 0.5633200580636218im
  -0.8262387743159954 + 0.5633200580636218im
  -0.6234898018587338 - 0.7818314824680289im
  -0.6234898018587338 + 0.7818314824680289im
  -0.3653410243663954 - 0.9308737486442047im
  -0.3653410243663954 + 0.9308737486442047im
 -0.07473009358642396 - 0.9972037971811801im
 -0.07473009358642396 + 0.9972037971811801im
  0.22252093395631406 - 0.9749279121818216im
  0.22252093395631406 + 0.9749279121818216im
   0.5000000000000004 - 0.8660254037844389im
   0.5000000000000004 + 0.8660254037844389im
    0.733051871829826 - 0.6801727377709196im
    0.733051871829826 + 0.6801727377709196im
   0.9009688679024193 - 0.43388373911755873im
   0.9009688679024193 + 0.43388373911755873im
   0.9888308262251289 - 0.14904226617617436im
   0.9888308262251289 + 0.14904226617617436im
p(x) = sum([ (-1)^k * x^k for k=0:20 ])
p (generic function with 1 method)
p(1)
1
p(2)
699051
p.(roots)
20-element Vector{ComplexF64}:
  3.9968028886505635e-15 + 1.237343560944737e-13im
  3.9968028886505635e-15 - 1.237343560944737e-13im
 -1.8318679906315083e-14 - 4.9960036108132044e-15im
 -1.8318679906315083e-14 + 4.9960036108132044e-15im
 -1.1102230246251565e-14 + 1.9761969838327786e-14im
 -1.1102230246251565e-14 - 1.9761969838327786e-14im
  -8.271161533457416e-15 - 6.328271240363392e-15im
  -8.271161533457416e-15 + 6.328271240363392e-15im
  3.4416913763379853e-15 - 2.55351295663786e-15im
  3.4416913763379853e-15 + 2.55351295663786e-15im
   2.148281552649678e-14 + 1.5432100042289676e-14im
   2.148281552649678e-14 - 1.5432100042289676e-14im
 -2.9976021664879227e-15 - 6.106226635438361e-15im
 -2.9976021664879227e-15 + 6.106226635438361e-15im
  -5.551115123125783e-16 + 3.9968028886505635e-15im
  -5.551115123125783e-16 - 3.9968028886505635e-15im
 -4.9960036108132044e-15 + 4.274358644806853e-15im
 -4.9960036108132044e-15 - 4.274358644806853e-15im
 -2.9976021664879227e-15 - 1.609823385706477e-15im
 -2.9976021664879227e-15 + 1.609823385706477e-15im
isapprox(p.(roots), zeros(20), atol=1e-10)
true

pp-norma vektoru/matice a normalizace

Metoda norm(A, p::Real=2) vypočte pp normu vektoru nebo matice. Ve výchozím stavu se volí p=2p=2, tedy Euklidovská norma. Lze zadat i hodnotu Inf pro nekonečno (tedy maximovou). Volba p=1p=1 pak odpovídá Manhattanské (taxicab) normě).

norm([1, 1, 1, 1]) # p = 2
2.0
norm([1, -1, 3], 1)
5.0
norm([-20, 0, 10], Inf)
20.0
norm([-20, 0, 10], 40)
20.00000000000046
norm([1im, 2])
2.23606797749979
sqrt(5)
2.23606797749979

Často je potřeba vektory tzv. normalizovat. Tj. vynásobit je takovou konstantou, aby jistá jejich norma byla jednotková. Přesně k tomuto účelu slouží metoda normalize, resp. normalize!, která tuto operaci provádí vzhledem k právě udvedené pp-normě.

normalize([-20, 0, 10], Inf)
3-element Vector{Float64}:
 -1.0
  0.0
  0.5
u = [1., 2. , 3., 4.]
println(norm(u, 3))
normalize!(u, 3)
println(u)
println(norm(u, 3))
4.641588833612778
[0.2154434690031884, 0.4308869380063768, 0.6463304070095652, 0.8617738760127536]
1.0000000000000002

Maticové funkce exp, log, sin, cos a sqrt

Exponenciální, logaritmická a trigonometrické funkce jsou pro čtvercové matice definovány také pomocí Taylorových polynomů. V případě logaritmu ne přo všechny matice.

Maticová metoda sqrt vrací matici, pro kterou platí sqrt(A) * sqrt(A) = A. Takováto matice přirozeně nemusí existovat pro všechny matice.

A = [1 -1; 1 1]
2×2 Matrix{Int64}:
 1  -1
 1   1
exp(A)
2×2 Matrix{Float64}:
 1.46869  -2.28736
 2.28736   1.46869
log(A)
2×2 Matrix{Float64}:
 0.346574  -0.785398
 0.785398   0.346574
sin(A)
2×2 Matrix{Float64}:
 1.29846   -0.634964
 0.634964   1.29846
cos(A)
2×2 Matrix{Float64}:
  0.83373   0.988898
 -0.988898  0.83373
sqrt(A)
2×2 Matrix{Float64}:
 1.09868  -0.45509
 0.45509   1.09868
sqrt(A) * sqrt(A)
2×2 Matrix{Float64}:
 1.0  -1.0
 1.0   1.0
A
2×2 Matrix{Int64}:
 1  -1
 1   1
A
I + A + A^2 / 2 + A^3 / 6 + A^4 / 24
2×2 Matrix{Float64}:
 1.5      -2.33333
 2.33333   1.5
exp(A)
2×2 Matrix{Float64}:
 1.46869  -2.28736
 2.28736   1.46869

Cvičení: Maticová exponenciála

Maticová exponenciála má celou řadu podobných vlastností jako obyčejná reálná (nebo komplexní) exponenciála. Například je řešením diferenciální rovnice y′=Ayy' = \mathbb{A} y, y(0)=y0∈Cny(0) = y_0 \in \mathbb{C}^n. Toto řešení má tvar y(t)=etAy0y(t) = e^{t\mathbb{A}} y_0.

Ukažme si to na příkladu matice

A=(0−110)\mathbb{A} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}

a počáteční podmínky y0=(1,2)Ty_0 = (1, 2)^T. Přibližně ověřte platnost diferenciální rovnice pro t=1t = 1.

Poznámka: V dalším notebooku budeme toto řešení i vizualizovat.

ts = LinRange(0., 5., 20)
A  = [0 -1; 1 0]
y0 = [1, 2]
function solution(mat, ival, times)
    data = zeros(length(times), length(ival))
    
    for j = 1:length(times)
        data[j, :] = exp(times[j] * A) * ival
    end
    
    data
end
ysI  = solution(A, [1, 2], ts)
ysII = solution(A, [-1, -1], ts);
using PyPlot
xlabel("\$y_1\$")
ylabel("\$y_2\$")
title("Trajektorie")
grid()
scatter(ysI[:,1], ysI[:, 2], color="r")
scatter(ysII[:,1], ysII[:, 2], color="b")

1.3 Logging

Tento modul zajišťuje makra

  • @debug: zprávy pro vývojáře,
  • @info: informační sdělení pro uživatele,
  • @warn: nefatální varování,
  • @error: hlášení zásadnějšího problému.

a další nástroje pro logování.

using Logging

Informace se vypisují podle nastavené úrovně logování (výchozí je Info; lze ji měnit i pomocí proměnné prostředí JULIA_DEBUG v které se předá jméno modulu, kde chceme úroveň Debug nastavit, pro REPL jde o Main; případně lze i vytvořit více loggerů, mohou zapisovat do souborů atd.).

x = 1

global_logger(SimpleLogger(stdout, Logging.Debug))  # <=
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"
Debug
┌ Debug: x = 1
└ @ Main In[219]:6
┌ Info: All done!
└ @ Main In[219]:7
┌ Warning: This is fishy!
└ @ Main In[219]:8
┌ Error: Run away!
└ @ Main In[219]:9

x = 1

global_logger(SimpleLogger(stdout, Logging.Info))  # výchozí hodnota
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"
Info
┌ Info: All done!
└ @ Main In[220]:7
┌ Warning: This is fishy!
└ @ Main In[220]:8
┌ Error: Run away!
└ @ Main In[220]:9

x = 1

global_logger(SimpleLogger(stdout, Logging.Warn))
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"
Warn
┌ Warning: This is fishy!
└ @ Main In[221]:8
┌ Error: Run away!
└ @ Main In[221]:9

x = 1

global_logger(SimpleLogger(stdout, Logging.Error))
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x"
@info "All done!"
@warn "This is fishy!"
@error "Run away!"

Cvičení

Vytvořte logger, který bude zapisovat všechny zprávy od úrovně debug do souboru log.txt.

stdout
IJulia.IJuliaStdio{Base.PipeEndpoint}(IOContext(Base.PipeEndpoint(RawFD(44) open, 0 bytes waiting)))
x = 1
y = "Cože?"

io = open("/home/kalvin/jul/tmp/log.txt", "w")

global_logger(SimpleLogger(io, Logging.Debug))
println(Logging.min_enabled_level(global_logger()))

@debug "x = $x" rand(2, 2) y
@info "All done!"
@warn "This is fishy!"
@error "Run away!"
Debug

close(io)

1.4 Pkg

Modul Pkg umožňuje spravovat balíčky. Už jsme se s ním vlastně setkali v Pkg módu. Následující kód je v podstatě ekvivalentí příkazu status v Pkg módu REPLu (po stisku ]).

import Pkg
Pkg.status()

Podrobněji do možností Pkg zabíhat na tomto místě nebudeme.

1.5 Random

Tento modul poskytuje podporu nejen pro generování (pseudo) náhodných čísel.

Například můžeme vybrat náhodný prvek matice:

A = [1 2 3; 4 5 6; 1 -1 1]
rand(A)
6
rand()
0.7484753105570906

Metoda rand nám také umožňuje generovat náhodné hodnoty s rovnoměrným rozdělením v intervalu [0,1)[0,1).

rand(10)
10-element Vector{Float64}:
 0.08141530874769598
 0.5560913883881384
 0.10694442071317978
 0.3731922499389829
 0.10982852682574384
 0.931439937142186
 0.7097203580185815
 0.3017112887502822
 0.24231900346544755
 0.748247708380128

Metoda randn podobně generuje hodnoty s náhodným rozdělením se střední hodnotou 00 a rozptylem 11.

randn(10)
10-element Vector{Float64}:
 -0.4722218179594087
  0.3184501267354136
 -2.718176111209226
  1.6762722912760657
  0.7612724403808242
 -0.35702555268003267
  0.8479271746102384
 -1.1080837424993843
 -1.2815958236847502
 -1.9512909504751033

Metoda randexp nám umožňuje generovat hodnoty s exponenciálním rozdělením s parametrem 11.

randexp
UndefVarError: `randexp` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name also exists in Random.

Stacktrace:
 [1] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
using Random
randexp(10)
10-element Vector{Float64}:
 0.8717275566138638
 0.5543842255335715
 0.46471192821883933
 0.1126213546402189
 0.22760147806171144
 0.5387085491009882
 1.5334899149482935
 0.1341687380512306
 0.7754869177324581
 0.9432610888038618

Můžem i generovat náhodný řetězec pomocí metody randstring.

randstring(10)
"Xzic0k955g"

Někdy je potřeba pracovat s permutacemi, resp. generovat náhodné permutace (zde chápány jako uspořádané nn-tice čísel 11 až nn). K tomu slouží metoda randperm:

randperm(10)
10-element Vector{Int64}:
  8
  6
  1
  4
  9
  2
  5
  7
 10
  3

Alternativně s může hodit metoda shuffle, která náhodně přeuspořádá prvky pole:

A
3×3 Matrix{Int64}:
 1   2  3
 4   5  6
 1  -1  1
shuffle(A)
3×3 Matrix{Int64}:
 -1  1  1
  1  5  2
  4  3  6
shuffle([1, 2, 3])
3-element Vector{Int64}:
 2
 1
 3

Všem těmto metodám lze případně předat "náhodný generátor", který se při generování na pozadí použije. pomocí něho lze i pseudonáhodnou posloupnost "nastartovat" (seed).

K dispozici máme dva generátory: MerseneTwister a RandomDevice. Což jsou konkrétní podtypy abstraktního typu AbstractRNG

MersenneTwister <: AbstractRNG
true
RandomDevice <: AbstractRNG
true
rng = MersenneTwister(1234)
rand(rng, 10)
10-element Vector{Float64}:
 0.5383210129299967
 0.9973545274591418
 0.027541876868637738
 0.07036318103103678
 0.11585668372488978
 0.9055435622074075
 0.00539224541327199
 0.3693301527613009
 0.627423891322153
 0.4968573924871824
rng = MersenneTwister(1234)
rand(rng, 10)
10-element Vector{Float64}:
 0.5383210129299967
 0.9973545274591418
 0.027541876868637738
 0.07036318103103678
 0.11585668372488978
 0.9055435622074075
 0.00539224541327199
 0.3693301527613009
 0.627423891322153
 0.4968573924871824
rng = MersenneTwister(1234)
shuffle(rng, [1, 2, 3, 4])
4-element Vector{Int64}:
 3
 4
 1
 2

RandomDevice nelze inicializovat pomocí seed. Náhodná data bere z hardware.

?RandomDevice()
rng = RandomDevice()
rand(rng, 10)
10-element Vector{Float64}:
 0.7903818357778811
 0.4932264282837786
 0.970130409041666
 0.5465212912937434
 0.8772235234006864
 0.4804446634573851
 0.6827763684102353
 0.3783874985722635
 0.7070075915286067
 0.5257503787129933

1.6 Sparse Arrays a Sparse Linear algebra

Tyto dva moduly umožňující efektivně pracovat s řídkými maticemi. Druhý modul poskytuje různé lineár algebraické algoritmy optimalizované pro řídké matice.

1.7 Statistics

Tento modul poskytuje metody pro výpočet známých statistických veličin jako jsou

  • standardní odchylka, std a stdm

  • variance, var a varm

  • korelace, cor

  • kovariance, cov

  • průměr, mean a mean!

  • medián, median a median!

  • kvantily, quantile a quantile!

    Následuje několik ukázek použití někerých z těchto metod.

v = rand(100)
100-element Vector{Float64}:
 0.12799280436067273
 0.36585645268759415
 0.8488470849684135
 0.7161572899851125
 0.9619699303436153
 0.5286043865256276
 0.8839744618439771
 0.836390491752679
 0.2952355289621884
 0.8005219778828653
 0.14027267874193505
 0.8988606763499496
 0.6985049355092245
 ⋮
 0.37117431232064624
 0.8143351518839034
 0.29935539841637837
 0.8751161821144668
 0.6089614408482413
 0.8079154587993485
 0.49748306931062125
 0.38558273975233304
 0.758031744576532
 0.5902859779518077
 0.255955516645207
 0.7835197185139483
using Statistics
mean(v)
0.5227650371243865
std(v)
0.3099197284144358
median(v)
0.5397435742171455

Pearsonův korelační koeficient dvou vektorů.

cor(v, v)
1.0
cor(v, -v)
-1.0
cor(v, v/2)
1.0
cor(v, rand(100))
-0.014343813471058888

Cvičení: PRNG

Logistické zobrazení je definováno na intervalu [0,1][0,1] předpisem x↦αx(1−x)x \mapsto \alpha x (1-x), kde α∈(0,4⟩\alpha \in (0,4\rangle je parametr. Pro α=4\alpha = 4 se chová "chaoticky", přesněji: sestrojíme-li rekurentní posloupnost xn+1=4xn(1−xn)x_{n+1} = 4x_n(1-x_n), n∈Nn\in\mathbb{N}, pak množina jejích hromadných bodů je celý interval ⟨0,1⟩\langle 0, 1 \rangle pro libovolnou počáteční podmínku (seed) x0∈(0,1)x_0 \in (0, 1).

Jak prozkoumáme v dalším notebooku, takto definovaná posloupnost bohužel preferuje okraje intervalu (0,1)(0,1), toho se lze zbavit škálováním y=1πarccos⁡(1−2x)y = \frac{1}{\pi}\arccos(1 - 2x).

Vygenerujte prvních několik členů těchto posloupností a spočtěte průměr, medián, standardní odchylku a Pearsonův korelační koeficient mezi sousedními členy.

function my_prng(n, seed)
    data = zeros(n)
    data[1] = seed
    
    for j = 2:n
        data[j] = 4data[j-1]*(1 - data[j-1])
    end
    
    return data
end

Volba prvního členu (seed) je důležitá.

my_prng(30, 0.7500001)
xs = my_prng(10^6, 0.765)
ys = [ acos(1 - 2x)/pi for x in xs ];
mean(xs), mean(ys)
median(xs), median(ys)
std(xs), std(ys)
cor(xs[1:(end-1)], xs[2:end])
cor(ys[1:(end-1)], ys[2:end])
zs = rand(10^6)
cor(zs[1:end-1], zs[2:end])

Pro odhad entropie použijeme balíček Discreet.jl.

using Discreet #  Nedělá přesně to co chceme!
estimate_entropy(ones(1_000_000))
estimate_entropy([ (-1)^n for n=1:1_000 ])
estimate_entropy(rand(1_000_000))
xs = my_prng(1_000, 0.765)
ys = [ acos(1 - 2x)/pi for x in xs ];
estimate_entropy(xs)
estimate_entropy(ys)
estimate_entropy(rand(1000))
estimate_entropy(LinRange(1, 10, 1000))

1.8 Test

Modul test poskytuje pomocná makra pro provádění unit testů. Tyto testy jsou v případě balíčků sdruženy v jednom adresáři. Můžeme s nimi pracovat ale i zcela samostatně.

using Test

@test

Prvním jednoduchým makrem je makro @test, které prostě testuje správnost (true) uvedeného výrazu.

@test true
Test Passed
@test false
Test Failed at In[24]:1
  Expression: false


There was an error during testing

Stacktrace:
 [1] record(ts::Test.FallbackTestSet, t::Union{Test.Error, Test.Fail})
   @ Test /usr/share/julia/stdlib/v1.12/Test/src/Test.jl:1030
 [2] do_test(result::Test.ExecutionResult, orig_expr::Any)
   @ Test /usr/share/julia/stdlib/v1.12/Test/src/Test.jl:713
 [3] top-level scope
   @ In[24]:1
 [4] macro expansion
   @ /usr/share/julia/stdlib/v1.12/Test/src/Test.jl:680 [inlined]
 [5] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

Čili například

@test isodd(15) == true
Test Passed

Dále toto makro interpretuje @test f(args...) key=val... jako @test f(args..., key=val...). To nám umožňuje použít metodu isapprox v infixové notaci pomocí ≈\approx i s uvedením absolutní tolerance:

@test π ≈ 3.0 atol=1.0
Test Passed
@test isapprox(π, 3.0, atol=1.0)
Test Passed

@test_throws

Jak název napovídá, toto makro testuje výjimky. Nejprve je uvedena očekávaná výjimka, poté testovaný výraz. 

@test_throws MethodError 1 + "10"
Test Passed
      Thrown: MethodError
@test_throws DimensionMismatch 1 + "10"
Test Failed at In[33]:1
  Expression: 1 + "10"
    Expected: DimensionMismatch
      Thrown: MethodError
      MethodError: no method matching +(::Int64, ::String)
      The function `+` exists, but no method is defined for this combination of argument types.
      
      Closest candidates are:
        +(::Any, ::Any, ::Any, ::Any...)
         @ Base operators.jl:642
        +(::Real, ::Complex{Bool})
         @ Base complex.jl:322
        +(::Real, ::Complex)
         @ Base complex.jl:334
        ...
      
      Stacktrace:
       [1] top-level scope
         @ In[33]:1
       [2] macro expansion
         @ /usr/share/julia/stdlib/v1.12/Test/src/Test.jl:774 [inlined]
       [3] macro expansion
         @ In[33]:1 [inlined]


There was an error during testing

Stacktrace:
 [1] record(ts::Test.FallbackTestSet, t::Union{Test.Error, Test.Fail})
   @ Test /usr/share/julia/stdlib/v1.12/Test/src/Test.jl:1030
 [2] do_test_throws(result::Test.ExecutionResult, orig_expr::Any, extype::Any)
   @ Test /usr/share/julia/stdlib/v1.12/Test/src/Test.jl:870
 [3] top-level scope
   @ In[33]:1
 [4] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

Sady testů (Test sets)

Pro přehlednost je vhodné testy sdružit do sad, ideálně podle souvisejících cílů testování. K tomu slouží makro @testset. Toto makro umožňuje testy parametrizovat zadaným parametrem či parametry.

Jednoduchým příkladem by tedy mohlo být například následující testování vlastností funkce exp:

@testset "Exponential" begin
    x = div(2, 2)
    
    @test exp(0)    == x
    @test exp(2.0)  ≈ exp(1.0) * exp(1.0)
    @test exp(-3.3) ≈ 1 / exp(3.3)
end;
Test Summary: | Pass  Total  Time
Exponential   |    3      3  0.0s

Parametrizovaný test by mohl třeba kontrolovat správnou délku pole náhodných čísel:

@testset "rand length ($i)" for i in 1:10
    @test length(rand(i)) == i
end;
Test Summary:   | Pass  Total  Time
rand length (1) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (2) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (3) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (4) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (5) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (6) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (7) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (8) |    1      1  0.0s
Test Summary:   | Pass  Total  Time
rand length (9) |    1      1  0.0s
Test Summary:    | Pass  Total  Time
rand length (10) |    1      1  0.0s

Pro další možnosti modulu Test odkazujeme čtenáře na jeho dokumentaci.

Cvičení

Sestavte hypotetické sady testů pro metodu * implemetující maticové násobení (výjimky při různých rozměrech, správnost výsledků násobení různě velikých matic, asociativita maticového násobení). Diskutujte vhodnost různých přístupů k otestování těchto vlastností.

Řešení některých příkladů

K maticové exponenciále:

A    = [0. -1.; 1. 0.]
y(t) = exp(t*A) * [1, 2]

eps = 1e-4

println((y(1 + eps) - y(1)) / eps)
println(A * y(1))

K PRNG:

function my_prng(n, seed)
    x   = seed
    seq = zeros(n)
    seq[1] = x
    
    for j = 2:n
        seq[j] = 4 * seq[j-1] * (1-seq[j-1])
    end
    
    seq
end

Reference

Všechny moduly standardní knihovny jsou popsány v oficiální dokumentaci, část Standard Library (bohužel nelze dát přímo odkaz na tuto část v HTML verzi dokumentace, jde pouze o rozklikávací menu v levé části stránky).