Jdi na navigaci předmětu

05: Vícerozměrná pole a maticové typy

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 2025/2026 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()
Julia Version 1.13.0
Commit d1c37793dd2 (2026-09-09 19:00 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, skylake)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)

1. Vícerozměrná pole (Array)

Vícerozměrné pole si lze představovat jako "kolekci" objektů uložených ve vícerozměrné mřížce (matice/tabulka, kvádr,...).

I pole mají svůj parametrický typ, který určuje jaké objekty dokáží uchovávat a jaký mají rozměr.

Nejobecnějším případem je přirozeně Array{Any}, resp. Array{Any, n}, kde nn označuje dimenzionalitu (rozměry) pole. Pole těchto typů pojmou libovolné objekty, tj. objekty typu Any. Obecně objekty typu Array{T} představují vícerozměrná pole s prvky typu T.

Array{Any} <: Array
true
Array{Any, 2} <: Array{Any}
true

Jak už bylo zmíněno v třetí lekci, k prvkům, nebo částem, pole přistupujeme pomocí hranatých závorek. Nejrpve se ale pusťme do konstrukce prvních polí, aby pak bylo co indexovat.


1.1 Vytváření a indexování polí

Neinicializované pole s rozměry dims a prvky typu T vytvoříme pomocí příkazu Array{T}(undef, dims...). dims může být tuple, nebo popořadě velikosti jednotlivých "dimenzí". Ukažme si to na příkladech:

# 1d array of Int64s with 5 elements
Array{Int64}(undef, 5)
5-element Vector{Int64}:
 140487656736016
 140487656736496
 140487656736976
 140487656737456
               0

Ihned poznamenejme, že Vector{T} jen šikovná zkratka pro jednodimenzionální pole s prvky typu T.

Vector{Any} == Array{Any, 1}
true

Nyní se přesuňme k dvojrozměrným polím.

# 2d 10x5 array of Float64s
Array{Float64}(undef, 10, 5)
10×5 Matrix{Float64}:
 6.94101e-310    6.94101e-310  6.94101e-310  6.94101e-310  6.94101e-310
 6.94101e-310    6.94101e-310  6.94101e-310  5.0e-324      6.94101e-310
 6.94101e-310    6.94101e-310  0.0           6.94101e-310  6.94101e-310
 6.94101e-310    6.94101e-310  6.94101e-310  6.94101e-310  6.94101e-310
 6.94102e-310    6.94101e-310  6.94101e-310  6.94101e-310  6.94101e-310
 6.94101e-310    5.0e-324      6.94101e-310  6.94101e-310  6.94101e-310
 6.94101e-310  NaN             6.94101e-310  6.94101e-310  1.0e-323
 6.94101e-310    6.94101e-310  6.94101e-310  6.94101e-310  1.5e-323
 6.94101e-310    6.94101e-310  6.94101e-310  6.94101e-310  1.0e-323
 6.94101e-310    6.94101e-310  6.94101e-310  6.94101e-310  1.0e-323

Podobně jako výše, Matrix{T} je pouze konvenční krycí název pro dvourozměrné pole.

Matrix{Any} == Array{Any, 2}
true

Od tří dimenzí dále už výpis pole samozřejmě není úplně přehledný, je ale možný.

# 3d 4x4x4 array of strings
Array{String}(undef, (4, 4, 4))
4×4×4 Array{String, 3}:
[:, :, 1] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef

[:, :, 2] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef

[:, :, 3] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef

[:, :, 4] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
Array{Any}(undef, (4, 4, 4))
4×4×4 Array{Any, 3}:
[:, :, 1] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef

[:, :, 2] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef

[:, :, 3] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef

[:, :, 4] =
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef

Tato pole nejsou inicializována, jejich hodnoty jsou "náhodné", podle toho co zbylo na místě původní paměti. Vzhledem k explicitnímu zmínění slovíčka undef v předchozích příkazech by se mohlo zdát, že můžeme prvky polí tímto způsobem předvyplnit požadovanou hodnotou, ale není tomu tak:

Array{Int64}(42, 5)
MethodError: no method matching (Array{Int64})(::Int64, ::Int64)
The type `Array{Int64}` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (Array{T})(::UndefInitializer, ::Int64, ::Int64, ::Int64) where T
   @ Core boot.jl:668
  (Array{T})(::UndefInitializer, ::Int64, ::Int64) where T
   @ Core boot.jl:667
  (Array{T})(::UndefInitializer, ::Int64) where T
   @ Core boot.jl:666
  ...


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

Viz metodu fill níže.


Vytváření polí s požadovanými prvky (Maticové literály)

Jednorozměrná (vektory) a dvourozměrná (matice) pole můžeme zadat explicitně pomocí hranatých závorek, před kterými můžeme případně uvést požadovaný typ prvků. K oddělení prvků máme bohaté možnosti: používáme čárky (,), mezery, konce řádků, nebo středníky (;). Tímto způsobem můžeme kontrolovat rozměry pole.

[1, 2, 3]
[1, 2, 3.]
transpose([1,2,3])
[1.0; 2.0; 3.0]
Float64[1, 2, 3]

Ale pozor, následující pole je už dvourozměrné:

[1 2 3]

Dvourozměrná pole zadáváme po řádcích oddělených středníky, nebo koncemi řádků:

[1 0; 0 1]
Float64[
    1 2 3
    4 5 6
]

Pozor, čárky neprojdou.

[1 0, 0 1]

Tento způsob se samozřejmě hodí jen pro malá pole. Častěji potřebujeme pracovat například s maticemi s velkým počtem prvků a takovéto explicitní vypisování by bylo komplikované. Typicky třeba známe předpis pro ijij-tý prvek. Pak můžeme použít analog Pythonovského list comprehension.

Například Float64 5x5 verzi Hilbertovy matice, které má na ijij-tém políčku hodnotu 1i+j−1\frac{1}{i+j-1} můžeme definovat takto:

H = [ 1 / (i + j - 1) for i = 1:5, j = 1:5 ]

Zde první index je řádkový a druhý sloupcový, jak bude patrné z následujícího příkladu

[ (i, j) for i = 1:4, j = 1:5 ]
[ (i, j) for i = 1:4, j = 1:4, k = 1:2]
[ (i, j) for i = 1:4, j = 1:4, k = 1:2, ℓ = 1:2]

Před hranatou závorkou můžeme opět vynutit typ prvků:

[ j for j = 1:10 ]
10-element Vector{Int64}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

Dejte pozor na ne/přítomnost hranatých závorek, například následující buňka nevyústí v matici:

vv = [ [ j+k for j=1:3 ] for k=1:3 ]
hcat(vv...)
vcat(vv...)

Lze přidávat i podmínky jako v Pythonu:

[ j for j = 1:10 if iseven(j) ]

Indexování polí

Prvky pole jsou standardně indexovány od 11. K prvkům pole přistupujeme pomocí hranatých závorek. K indexování můžeme použít i rozsah (range) a získat tak část pole (slice, řez). Stejným způsobem můžeme prvky, nebo celé bloky matice, modifikovat pomocí přiřazení. Ukažme si to na následující Töplitzově matici (konstantní hodnoty na diagonále i všech vedlejších diagonálách):

T = [ 5 - i + j for i = 1:4, j = 1:4 ]
4×4 Matrix{Int64}:
 5  6  7  8
 4  5  6  7
 3  4  5  6
 2  3  4  5

Hodnoty v rozích matice:

println("Levý horní:  ", T[1, 1])
println("Pravý horní: ", T[1, 4])
println("Levý dolní:  ", T[4, 1])
println("Pravý dolní: ", T[4, 4])
Levý horní:  5
Pravý horní: 8
Levý dolní:  2
Pravý dolní: 5
T[1, 5]
BoundsError: attempt to access 4×4 Matrix{Int64} at index [1, 5]

Stacktrace:
 [1] throw_boundserror(A::Matrix{Int64}, I::Tuple{Int64, Int64})
   @ Base ./essentials.jl:14
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] getindex(::Matrix{Int64}, ::Int64, ::Int64)
   @ Base ./array.jl:962
 [4] top-level scope
   @ In[20]:1
 [5] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
T[-1, 1]
BoundsError: attempt to access 4×4 Matrix{Int64} at index [-1, 1]

Stacktrace:
 [1] throw_boundserror(A::Matrix{Int64}, I::Tuple{Int64, Int64})
   @ Base ./essentials.jl:14
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] getindex(::Matrix{Int64}, ::Int64, ::Int64)
   @ Base ./array.jl:962
 [4] top-level scope
   @ In[21]:1
 [5] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

Změníme hodnotu v pravo nahoře:

T[1, 4] = -2
T
4×4 Matrix{Int64}:
 5  6  7  -2
 4  5  6   7
 3  4  5   6
 2  3  4   5

2x3 podmatice v levém horním rohu:

1:2
typeof(1:2)
T
4×4 Matrix{Int64}:
 5  6  7  -2
 4  5  6   7
 3  4  5   6
 2  3  4   5
T[1:2, 1:3]
2×3 Matrix{Int64}:
 5  6  7
 4  5  6

Změna 2x2 bloku vpravo dole na jednotkovou matici 2x2:

T[3:4, 3:4] = [1 0; 0 1]
T
4×4 Matrix{Int64}:
 5  6  7  -2
 4  5  6   7
 3  4  1   0
 2  3  0   1

Pokud chceme celý řádek/sloupec, pak pro tyto účely můžeme použít pouze :. Rozsah "až do konce" můžeme vyjádřit pomocí end:

T[:, 2:end]
4×3 Matrix{Int64}:
 6  7  -2
 5  6   7
 4  1   0
 3  0   1
T[:, 2:-1]
4×0 Matrix{Int64}
T[:,1:2:end]
4×2 Matrix{Int64}:
 5  7
 4  6
 3  1
 2  0
T
T[:, 2:-1:1]
4×2 Matrix{Int64}:
 6  5
 5  4
 4  3
 3  2
T[:, 2:-1:-1]
BoundsError: attempt to access 4×4 Matrix{Int64} at index [1:4, 2:-1:-1]

Stacktrace:
 [1] throw_boundserror(A::Matrix{Int64}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, StepRange{Int64, Int64}})
   @ Base ./essentials.jl:14
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] _getindex
   @ ./multidimensional.jl:983 [inlined]
 [4] getindex(::Matrix{Int64}, ::Function, ::StepRange{Int64, Int64})
   @ Base ./abstractarray.jl:1363
 [5] top-level scope
   @ In[30]:1
 [6] eval(m::Module, e::Any)
   @ Core ./boot.jl:489
T
4×4 Matrix{Int64}:
 5  6  7  -2
 4  5  6   7
 3  4  1   0
 2  3  0   1
T[:, [3, 1, 2]]
4×3 Matrix{Int64}:
 7  5  6
 6  4  5
 1  3  4
 0  2  3

Nyní probereme metody, které se hodí při vytváření nových "obecných" polí. V případě matic máme více možností, kterým se budeme věnovat v příslušné lekci.


Konkatenace polí (vcat a hcat)

V metamatickém zápisu často používáme blokový zápis matic. Jeho analogem je konkatenace.

Mějme tři pole:

A = [1, 1, 1]
B = [2 2]
C = [3 3; 3 3; 3 3]

Pokud má operace rozměrově smysl, můžeme pole spojovat horizontálně:

[A C A]
hcat(A, C, A)

Nebo vertikálně:

[B; C; B]
vcat(B, C, B)

A oba přístupy lze dohromady i kombinovat právě v něčem, co připomíná blokový zápis:

[ B B; C C]
[
    B B
    C C
]

push! a pop!

Do jednorozměrného pole lze přidávat prvky pomocí push! a ubírat je pomocí pop! (odzadu):

a = [1, 2, 3]

push!(a, 4)

a
x = pop!(a)

println(x)
a

Dále máme i metodu append!, která ale spojuje více jednorzměrných polí dohromady.

append!(a, [5, 6])

a
append!(a, 42)
append!(a, [42])
@which append!(a, 42)

Existují i další související metody jako pushfirst!, popfirst!, popat!, insert! a další.


fill a fill!

K uniformní inicializaci prvků pole slouží dvě metody:

  • fill(x, dims...): vytvoří nové pole s prvky x a rozměry dims.
  • fill!(A, x): zaplní již existující pole A prvky x
# 1d array of 5 elemens equal to 42
fill(42, 5)
# 2d 3x4 array of 0.5s
fill(0.5, 3, 4)
fill(1.23, (2, 2))
fill("Ahoj!", 2, 2)
B = fill(undef, 2, 2)
fill!(B, "Hello!")
A = Array{String}(undef, 2, 2)
println(A)

fill!(A, "Hello!")
println(A)

zeros a ones

Často chceme předvyplnit matice jedničkami, nebo nulami. Přesně k tomuto účelu slouží metody zeros(T, dims...) a ones(T, dims...). První argument, tedy typ, je nepovinný, výchozí hodnota je Float.

zeros(2, 2)
ones(3)
ones(2, 2)
ones(2, 2, 2)

Kdybychom chtěli stejné matice, ovšem s Int64 prvky, stačí tento požadavek explicitně zmínit:

zeros(Int64, 2, 2)
ones(Int64, 3)

Pro vytvoření pouze nuly nebo jedničky v daném typu slouží analogické metody zero a one.

zero(Float64)
zero(0.1 + 2.1im)
one(Int64)
one(38465)
ones(4, 4)
one(Matrix{Int64}(undef, 2, 2))

Ale viz I v další lekci.


copy a deepcopy

Tyto metody asi nepotřebují vysvětlení. Vytvoří kopii, resp. hlubokou kopii, daného pole. Názorná ukázka:

A = Array{Any}(undef, 2)
a = 42
b = [1, 2]

A[1] = a
A[2] = b

A
C = copy(A)
D = deepcopy(A)

a    = 11
b[1] = 9999

println(C)
println(D)
b

rand a randn

Náhodně generovaným hodnotám, resp. pravděpodobnosti a statistice, se budeme podrobně věnovat později během semestru. Zde si zmíníme jenom dvě metody pro vytváření "náhodných" matic.

První z nich, rand(T, dims...) vytvoří pole dimenzí dims s prvky typu T rovnoměrně rozloženými (u podtypů AbstractFloat v intervalu [0,1)[0, 1)), první typový argument je opět nepovinný a jeho výchozí hodnota je Float64.

rand(Int64, 4)
rand(5, 5)
rand([1, 2, "Ahoj!"])
rand(1:10)
rand([1  2 "Ahoj!"; π 7 Any], 2, 2)

Příkaz randn(T, dims...) se od předchozího liší tím, že generuje prvky podle normálního rozdělení (se střední hodnotou 00 a standardní odchylkou 11, má dobrý smysl jen pro podtypy AbstractFloat).

randn(10)

Následující pokus přirozeně skončí nezdarem:

randn(Int64, 2, 2)

Pokud rand předáme pole, pak náhodně vybere jeden jeho prvek:

for j = 1:10
    println(rand([1, 2, 3, 4]))
end
for j = 1:10
    println(rand(1:4))
end
for j = 1:10
    println(rand([-1 0; 1 2]))
end

To by jako exkurze do světa náhody pro tento okamžik stačilo. Podrobněji se do něj vrátíme hned příště.


range

Pro vytváření jednorozměrných polí se často hodí metoda range(start; stop=stop, length=n), pomocí níž lze vytvořit pole s rovnoměrně rozmístěnými prvky:

range(1, stop=5)
range(1, stop=6, length=10)

Další užitečné metody

Pro práci s poli se mohou hodit následující metody:

  • eltype(A): typ prvků pole
  • length(A): počet prvků pole
  • ndims(A): počet rozměrů ("dimenzí") pole
  • size(A): tuple s velikostí jednotlivých rozměrů pole
  • size(A, n): velikost nntého rozměru pole
  • axes(A): tuple s hodnotami indexů v jednotlivých rozměrech (o přeindexování polí později)

Ukázka:

A = rand(4, 3, 2)
eltype(A)
length(A)

Vlastní implementace:

typ_prvku(matice::Matrix{N}) where {N} = N 
typ_prvku(ones(2, 2))
typ_prvku([1 1; "String" 1.0])
typeof(A)
length(A)
4 * 3 * 2
length([[1, 2], [3, 4]])
ndims(A)
ndims([1, 2])
ndims([1 2; 3 4])
size(A)
println(size(A, 1))
println(size(A, 2))
println(size(A, 3))

Tj. nemusíme dělat ošklivé:

size(A)[2]
axes(A)
axes(A, 2)

N.B.:

Base.OneTo(4) == 1:4
for j = 1:size(A, 2)
    println(A[1, j, 1])
end
for j = axes(A, 2)
    println(A[1, j, 1])
end

1.2 Cvičení

  1. Vytvořte matici A rozměru 10×1010 \times 10 mající na antidiagonálách postupně konstanty 11, 22, 33,..., 1919.
  2. Vytvořte matici B rozměru 5×105 \times 10, která má všechny prvky nulové, pouze na "okraji" (obdélník) má jedničky.
  3. Vytvořte čtvercovou matici C rozměru 20×2020 \times 20, která je složena ze čtyř stejně velkých bloků: vlevo nahoře je blok náhodně vygenerovaných strojových čísel, vpravo nahoře matice naplněná jedničkami, vlevo dole nulová matice a vpravo dole diagonální matice s počátečními ciframi čísla π\pi na diagonále.
  4. Vytvořte čtvercovou matici D rozměru 10×1010 \times 10, která bude v levém horním trojúhelníku mít Pascalův trojúhelník. Prvky vypočtete Pascalovou "konstrukcí", ne počítáním kombinačních čísel.
A = [ i + j + 1 for i = 0:9, j = 0:9 ]
A
B = [ i == 1 || j == 1 || i == 5 || j == 10 ? 1 : 0 for i = 1:5, j = 1:10]
B
B = ones(Int64, 5, 10)
B[2:4, 2:9] = zeros(Int64, 3, 8)
B
πvo = [3, 1, 4, 1, 5]
Pi = [ j == k ? πvo[j] : 0 for j = 1:5, k = 1:5 ]
C = [ rand(5, 5) ones(5, 5); zeros(5, 5) Pi ]
C
D = [ 
      1 1 1 1
      1 2 3 0 
      1 3 0 0
      1 0 0 0
    ]
# Toto jistě není nejjednodušší řešení!

D = zeros(Int64, 10, 10)
D[1, :] = ones(10)
D[:, 1] = ones(10)
for j = 2:10 # cisluje antidiagonalu
    for k = 2:(j-1) # iterace po antidiagonale
        D[1 + j - k, k] = D[1 + j - k - 1, k] + D[1 + j - k, k - 1] 
    end
end

D

2. Operace s poli/maticemi a jejich prvky

S poli a jejich prvky lze provádět celou řadu nejen matematických operací. Většinou jsme motivování operacemi s vektory a maticemi (na rozdíl od jiných programovacích jazyků).


2.1 Sčítání + a odčítání -

Pole lze přirozeně po složkách sčítat (+), odčítat (-). Tyto operátory fungují jak unárně, tak binárně. Pole musí být stejného rozměru.

A = [1 2; 3 4]
B = [5 6 7; 8 9 10]
C = [1 1; -1 -1]

Součet a rozdíl po složkách, pokud si odpovídají rozměry:

A + A
A - C

"Opačná matice":

-A
A + B
[1, 2, 3] + [4, 5, 6]
[1, 2, 3] - [4, 5, 6]
[1, 2, 3] + [1 2 3]
[1, 2, 3] + [1 2 3]'

Chování vůči typům:

1 + 2.4
C + rand(2, 2)

Násobení *

Hvězdička má význam "násobení" a to má více významů. Pole můžeme násobit číslem (skalárem), hvězdičku v tomto případě lze i vynechat:

A
4 * A
-2A
A * 2

Dále má hvězdička význam maticového násobení v pravém slova smyslu, tj. matice musí být správného rozměru:

A * A
A * B
B * A

Inverze inv

Maticovou inverzi lze počítat pouze pro regulární matice.

inv([1 2; 3 4])
[1 2; 3 4] * inv([1 2; 3 4])
inv([1 2; 3 4]) * [1 2; 3 4] 
inv([1 1; 1 1])
inv([1 2 3; 4 5 6])

Vidíte, že je-li typem prvků matice integer, výsledkem bude matice strojových čísel. Podobné chování jsme už viděli i u obyčejného dělení. I pro jednotkovou matici, kde k popsání inverze strojová čísla potřeba nejsou, k tomuto efektu dojde:

inv([1 0; 0 1])

Invertovat ale můžeme i přímo v aritmetice racionálních čísel.

ratA    = Rational{Int64}[1 1; 1 -1]
ratA
ratAinv = inv(ratA)
ratA * ratAinv
ratAinv * ratA

Zpět k příklad výše.

m = [ 1//1 2//1; 3//1 4//1 ]
invm = inv(m)
m * invm
H = [ 1 // (i + j - 1) for i = 1:8, j = 1:8 ]
inv(H)
H * inv(H)
H = [ 1 / (i + j - 1) for i = 1:15, j = 1:15 ]
inv(H) * H
big"1" // 2 |> typeof
H = [ big"1" // (i + j - 1) for i = 1:15, j = 1:15 ]
inv(H)
inv(H) * H

Dělení /, resp. \

Přesněji, mezi maticovými operandy jde o násobení maticovou inverzí zprava, resp. zleva.

A / A
A \ B
inv(A) * B

Umocňování ^

Matice lze umocňovat na celočíselné exponenty, se standardními omezeními.

A ^ 2
A * A
A ^ 6
A ^ -1
inv(A)

Provádění operací po složkách pomocí . (broadcasting)

Binární operace i metody lze provádět a aplikovat po složkách, resp. na prvcích matice, pomocí symbolu ..

Například můžeme vynásobit odpovídající si prvky obdélníkové matice (které by v pravém slova smyslu maticově násobit nešlo):

B
[ iseven(B[j, k]) for j in axes(B, 1), k in axes(B, 2) ]
iseven.(B)
map(iseven, B)
iseven(B)
B
B .* B
B * B

Operand ani nemusí být maticový:

4 .* B
4 * B
4 .+ B
4 + B
B
B[:, 2:3] .= 0
B
B .<< 1

Můžeme například aplikovat sinus na všechny prvky matice:

sin.(A)

Pozor! Následující kód je naprosto v pořádku, většinu elementárních funkcí lze zavést i pro operátory (čtvercové matice):

sin(A)
A - A^3 / 6 + A^5 / factorial(5) - A^7 / factorial(7)  + A ^ 9 / factorial(9) - A ^ 11 / factorial(11) +
    A ^ 13 / factorial(13) - A ^ 15 / factorial(15)

Tento přístup (broadcasting) funguje i pro funkce mající více či méně argumentů, nebo pro další binární operátory!

f = (x) -> x^3
A
f.(A)
map(f, A)
A
C
min.(A, C)
min.(A, 3)
min(A, 3)
min.(A, C, ones(Int64, 2, 2))

Maskování

Porovnáváním (například pomocí <) matic prvek po prvku můžeme vytvořit masku, pomocí které lze poté z matice extrahovat požadované prvky.

A
C = [0 3; 4 -1]
A .< C
mask = A .< C

V zápisu výše jde v podstatě o infixovou notaci, explicitněji:

isless.(A, C)

Indexováním pomocí takovéto matice Dostaneme pouze prvky se zkoumanou vlastností:

A[mask]
A[A .< C]

Jeden z operandů opět nutně nemusí být maticový:

C
C .<= 2
C[C .<= 5]

Nebo (všimněte si, že nedostaneme matici s Boolean prvky, jak bychom možná očekávali):

mask = iseven.(A)
A[mask]

Alternativně bychom k podobnému účelu mohli použít metodu filter:

filter(iseven, A)

Porovnání rychlosti

Explicitní for loop vs. broadcasting.

using BenchmarkTools

function f1(A::Matrix{T}) where { T <: Number }
    A[1, :] .*= A[1, :]
end

function f2(A::Matrix{T}) where { T <: Number }
    for j in axes(A, 2)
        A[1, j] *= A[1, j]
    end
end

function f3(A::Matrix{T}) where { T <: Number }
    A[1, :] = A[1, :] .* A[1, :]
end

mat = rand(2, 1_000_000)
@benchmark f1($mat)
@benchmark f2($mat)
@benchmark f3($mat)

Zde vidíme, že Julia má ráda rychlé výpočetní for smyčky, není nutná snaha o vektorizaci, jak je tomu třeba u Pythonu nebo Matlabu (kde pak díky tomu dané operace probíhají rychle mimo Python/Matlab díky C/FORTRAN).

Maskování vs. filter.

function my_filter(A::Matrix{T}, number::T) where { T <: Number }
    A[A .< number]
end
@benchmark my_filter($mat, 0.5)
@benchmark filter((x) -> x < 0.5, $mat)

Cvičení: Monte Carlo výpočet operátorové normy symetrické matice

Pro reálnou čtvercovou matici A\mathbb{A} definujeme její operátorovou normu předpisem ∥A∥=sup⁡{∥Ax∥∣x∈Rn, ∥x∥=1}.\|\mathbb{A}\| = \sup \big\{ \|\mathbb{A}x\| \mid x \in \mathbb{R}^n, \ \|x\| = 1 \big\}. Na pravé straně mají normy význam Euklidovké normy na Rn\mathbb{R}^n.

Jinak řečeno, každá hodnota ∥Ax∥\| \mathbb{A} x\| pro jednotkový vektor xx představuje spodní odhad ∥A∥\|\mathbb{A}\|.

Implementujte metodu, která se náhodnou volbou různých vektorů snaží dostat co nejlepší odhad této normy.

Zatím nepoužívejme nástroje z LinearAlgebra, vystačíme si se základními nástroji pro práci s poli. Připomeňme, že Euklidovská norma vektoru x∈Rnx\in\mathbb{R}^n je ∑j=1nxj2\sqrt{\sum_{j=1}^n x_j^2}.

using ProgressMeter

function rands2(n::Int64)
    x = rand(n) .- 0.5
    return x / sqrt(sum(x .* x))
end

"""
Náhodný vektor na jednotkové `n`-dimenzionální sféře.
"""
function rands(n::Int64)
    v = ones(n)
    f = () -> 2 * (rand() - 0.5)

    v[1] = f()
    v[2] = sqrt(1 - v[1]^2)
    
    for j = 2:n-2
        r = f()
        v[j], v[j+1] = v[j] * r, v[j] * sqrt(1 - r^2)
    end
    
    if n > 2
        r = f()
        v[n-1], v[n] = v[n-1] * r, v[n-1] * rand([-1, 1]) * sqrt(1 - r^2)
    end
    
    return v
end

"""
Spodní odhad operátorové normy pomocí Monte Carlo přístupu.
"""
function mc_op_norm(A::Matrix{Float64}, trials::Int64)
    size(A, 1) == size(A, 2) || error("`A` is not a square matrix!")
    
    max_norm = 0.0
    
    @showprogress 1 "Tossing coins furiously..." for j = 1:trials
        x = A * rands(size(A, 1))
        val = sum(x .* x)
        val > max_norm && (max_norm = val)
    end
    
    return sqrt(max_norm)
end

Test náhodného generátoru.

rands2(2)
x = rands(2)
println(x)
sum(x .* x)
x = rands(10)
println(x)
sum(x .* x)

Normy některých matic. Začněme nulovou maticí.

mc_op_norm(zeros(3, 3), 10)

Diagonální matice.

mc_op_norm([-4. 0; 0 3], 1_000_000)

A ještě jedna zajímavá matice:

mc_op_norm([1. 4; 5 6], 10^6)
sqrt(39 + 5sqrt(53))

Pro nn rostoucí do nekonečna se norma následujících matic blíží 22.

J(n) = Float64[ abs(j-k) == 1 ? 1 : 0 for j=1:n, k=1:n ]
J(5)
mc_op_norm(J(2), 10^6)
mc_op_norm(J(20), 10^6)
mc_op_norm(J(20), 10^7)
mc_op_norm(J(30), 10^9)

Pro rostoucí rozměry by tato hodnota měla konvergovat k 22.

H = [ 1 / (j + k - 1) for j = 1:50, k = 1:50 ]
mc_op_norm(H, 10_000_000)

Cvičení 2.3: Erastothenovo síto

Pomocí maskování (indexování pomocí BitArray) lze elegantně využít k implementaci jednoduchého Erastothenova síta.

Erastothenovo síto nám umožňuje nalézt všechna prvočísla menší nebo rovna zadané mezi (n). Základní verze bez optimalizací je následující

  • Jednička není prvočíslo.
  • Dalším číslem je dvojka, tedy prvočíslo, zahodíme všechny další násobky dvou.
  • Dalším nevyškrtnutým číslem je trojka, tedy prvočíslo, zahodíme všechny další násobky tří.
  • Čtyřku jsme už vyškrtli, takže se dostaneme k pětce, prvočíslu, a opět zahodíme všechny další násobky pěti.
  • atd.

I. Implementujte metodu erastothenes(n) počítající všechna prvnočísla menší nebo rovna n.

# N.B. BitArray: 
m = BitArray([1 0; 1 0])
println(eltype(m))
println(m[1,1] == true)
# Range se chová jako vektor
(1:5)[[true, true, false, false, true]]
function erastothenes(n::Integer)
    mask = ones(Bool, n)
    mask[1] = false
    
    for j = 2:ceil(Int64,sqrt(n))
        if mask[j]
            for k = 2j:j:n
                mask[k] = 0
            end
        end
    end
    
    return (1:n)[mask]
end
erastothenes(50)

II. Spočtěte všechna prvnočísla menší nebo rovna milion (OEIS:A000040).

ps = erastothenes(100_000_000)

III. S pomocí napočtených dat nalezněte všechna prvočíselná dvojčata v uvedeném rozsahu (tj. prvočísla lišící se o 22; OEIS:A077800). Implementujte metodu, která nalezne všechna dvojčata lišící se o zadané číslo. Tato metoda bude vracet pole tuplů (dvojčat).

function twins(primes::Vector{Int64}, gap::Int64)
    # ...
end
twins(ps, 2)
twins(ps, 4)

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

Možné řešení Erastothenova síta:

function erastothenes(n::Integer)
    mask = ones(Bool, n)
    mask[1] = 0

    for k = 2:n
        if mask[k]
            for j = 2k:k:n
                mask[j] = 0
            end
        end
    end
    
    (1:n)[mask]
end

Možné řešení twins metody.

function twins(primes::Vector{Int64}, gap::Int64)
    twins = []
    
    for j in 2:length(primes)
        if primes[j] - primes[j-1] == gap
            push!(twins, (ps[j-1], ps[j]))
        end
    end

    return twins
end

Monte Carlo.

using ProgressMeter

"""
Náhodný vektor na jednotkové `n`-dimenzionální sféře.
"""
function rands(n::Int64)
    v = ones(n)
    f = () -> 2 * (rand() - 0.5)

    v[1] = f()
    v[2] = sqrt(1 - v[1]^2)
    
    for j = 2:n-2
        r = f()
        v[j], v[j+1] = v[j] * r, v[j] * sqrt(1 - r^2)
    end
    
    if n > 2
        r = f()
        v[n-1], v[n] = v[n-1] * r, v[n-1] * rand([-1, 1]) * sqrt(1 - r^2)
    end
    
    return v
end

"""
Spodní odhad operátorové normy pomocí Monte Carlo přístupu.
"""
function mc_op_norm(A::Matrix{Float64}, trials::Int64)
    size(A, 1) == size(A, 2) || error("`A` is not a square matrix!")
    
    max_val = 0.0
    
    @showprogress 1 "Tossing coins furiously..." for j = 1:trials
        v = A * rands(size(A, 1))
        val = sum(v .* v)
        
        val > max_val && (max_val = val)
    end
    
    sqrt(max_val)
end

Reference

V oficiální dokumentaci je tato problematika shrnuta v kapitole Multi-dimensional Arrays