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.12.0
Commit b907bd0600f (2025-10-07 15:42 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-18.1.7 (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 $n$ 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}:
   2203335000064
               1
 139690155076704
               0
 139690074177544

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

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

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.90161e-310  6.90161e-310  6.90161e-310  6.9016e-310   6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.9016e-310   6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310
 6.9016e-310   6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.9016e-310   6.90161e-310
 6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310  6.90161e-310

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

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

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))

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)

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 $ij$ -tý prvek. Pak můžeme použít analog Pythonovského list comprehension.

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

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

5×5 Matrix{Float64}:
 1.0       0.5       0.333333  0.25      0.2
 0.5       0.333333  0.25      0.2       0.166667
 0.333333  0.25      0.2       0.166667  0.142857
 0.25      0.2       0.166667  0.142857  0.125
 0.2       0.166667  0.142857  0.125     0.111111

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 ]

4×5 Matrix{Tuple{Int64, Int64}}:
 (1, 1)  (1, 2)  (1, 3)  (1, 4)  (1, 5)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)  (2, 5)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)  (3, 5)
 (4, 1)  (4, 2)  (4, 3)  (4, 4)  (4, 5)

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

4×4×2 Array{Tuple{Int64, Int64}, 3}:
[:, :, 1] =
 (1, 1)  (1, 2)  (1, 3)  (1, 4)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)
 (4, 1)  (4, 2)  (4, 3)  (4, 4)

[:, :, 2] =
 (1, 1)  (1, 2)  (1, 3)  (1, 4)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)
 (4, 1)  (4, 2)  (4, 3)  (4, 4)

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

4×4×2×2 Array{Tuple{Int64, Int64}, 4}:
[:, :, 1, 1] =
 (1, 1)  (1, 2)  (1, 3)  (1, 4)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)
 (4, 1)  (4, 2)  (4, 3)  (4, 4)

[:, :, 2, 1] =
 (1, 1)  (1, 2)  (1, 3)  (1, 4)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)
 (4, 1)  (4, 2)  (4, 3)  (4, 4)

[:, :, 1, 2] =
 (1, 1)  (1, 2)  (1, 3)  (1, 4)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)
 (4, 1)  (4, 2)  (4, 3)  (4, 4)

[:, :, 2, 2] =
 (1, 1)  (1, 2)  (1, 3)  (1, 4)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)
 (4, 1)  (4, 2)  (4, 3)  (4, 4)

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

UInt32[ j for j = 1:10 ]

10-element Vector{UInt32}:
 0x00000001
 0x00000002
 0x00000003
 0x00000004
 0x00000005
 0x00000006
 0x00000007
 0x00000008
 0x00000009
 0x0000000a

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 ]

3-element Vector{Vector{Int64}}:
 [2, 3, 4]
 [3, 4, 5]
 [4, 5, 6]

hcat(vv...)

3×3 Matrix{Int64}:
 2  3  4
 3  4  5
 4  5  6

vcat(vv...)

9-element Vector{Int64}:
 2
 3
 4
 3
 4
 5
 4
 5
 6

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

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

5-element Vector{Int64}:
  2
  4
  6
  8
 10

Indexování polí

Prvky pole jsou standardně indexovány od $1$ . 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:15
 [2] checkbounds
   @ ./abstractarray.jl:699 [inlined]
 [3] getindex(::Matrix{Int64}, ::Int64, ::Int64)
   @ Base ./array.jl:928
 [4] top-level scope
   @ In[19]: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:15
 [2] checkbounds
   @ ./abstractarray.jl:699 [inlined]
 [3] getindex(::Matrix{Int64}, ::Int64, ::Int64)
   @ Base ./array.jl:928
 [4] top-level scope
   @ In[20]: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

1:2

typeof(1:2)

UnitRange{Int64}

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

4×4 Matrix{Int64}:
 5  6  7  -2
 4  5  6   7
 3  4  1   0
 2  3  0   1

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:15
 [2] checkbounds
   @ ./abstractarray.jl:699 [inlined]
 [3] _getindex
   @ ./multidimensional.jl:955 [inlined]
 [4] getindex(::Matrix{Int64}, ::Function, ::StepRange{Int64, Int64})
   @ Base ./abstractarray.jl:1342
 [5] top-level scope
   @ In[34]:1
 [6] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

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

4-element Vector{Int64}:
 1
 2
 3
 4

x = pop!(a)

println(x)
a

3-element Vector{Int64}:
 1
 2
 3

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

append!(a, [5, 6])

a

5-element Vector{Int64}:
 1
 2
 3
 5
 6

append!(a, 42)

6-element Vector{Int64}:
  1
  2
  3
  5
  6
 42

append!(a, [42])

7-element Vector{Int64}:
  1
  2
  3
  5
  6
 42
 42

@which append!(a, 42)

append!(a::AbstractVector, iter) in Base at array.jl:1357

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)

5-element Vector{Int64}:
 42
 42
 42
 42
 42

# 2d 3x4 array of 0.5s
fill(0.5, 3, 4)

3×4 Matrix{Float64}:
 0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5
 0.5  0.5  0.5  0.5

fill(1.23, (2, 2))

2×2 Matrix{Float64}:
 1.23  1.23
 1.23  1.23

fill("Ahoj!", 2, 2)

2×2 Matrix{String}:
 "Ahoj!"  "Ahoj!"
 "Ahoj!"  "Ahoj!"

B = fill(undef, 2, 2)

2×2 Matrix{UndefInitializer}:
 UndefInitializer()  UndefInitializer()
 UndefInitializer()  UndefInitializer()

fill!(B, "Hello!")

MethodError: Cannot `convert` an object of type String to an object of type UndefInitializer
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  convert(::Type{T}, ::T) where T
   @ Base Base_compiler.jl:133


Stacktrace:
 [1] fill!(dest::Matrix{UndefInitializer}, x::String)
   @ Base ./array.jl:328
 [2] top-level scope
   @ In[49]:1
 [3] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

A = Array{String}(undef, 2, 2)
println(A)

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

[#undef #undef; #undef #undef]
["Hello!" "Hello!"; "Hello!" "Hello!"]

`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)

2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

ones(3)

3-element Vector{Float64}:
 1.0
 1.0
 1.0

ones(2, 2)

2×2 Matrix{Float64}:
 1.0  1.0
 1.0  1.0

ones(2, 2, 2)

2×2×2 Array{Float64, 3}:
[:, :, 1] =
 1.0  1.0
 1.0  1.0

[:, :, 2] =
 1.0  1.0
 1.0  1.0

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

zeros(Int64, 2, 2)

2×2 Matrix{Int64}:
 0  0
 0  0

ones(Int64, 3)

3-element Vector{Int64}:
 1
 1
 1

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

zero(Float64)

0.0

zero(0.1 + 2.1im)

0.0 + 0.0im

one(Int64)

one(38465)

ones(4, 4)

4×4 Matrix{Float64}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

one(Matrix{Int64}(undef, 2, 2))

2×2 Matrix{Int64}:
 1  0
 0  1

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

2-element Vector{Any}:
 42
   [1, 2]

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

a    = 11
b[1] = 9999

println(C)
println(D)

Any[42, [9999, 2]]
Any[42, [1, 2]]

`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)$ ), první typový argument je opět nepovinný a jeho výchozí hodnota je Float64.

rand(Int64, 4)

4-element Vector{Int64}:
  2959604866692875380
 -7875423552341897625
 -4425770935045202627
  2298360307598609262

rand(5, 5)

5×5 Matrix{Float64}:
 0.871191    0.876734  0.353404  0.37807    0.388049
 0.00121004  0.811214  0.110037  0.666142   0.770953
 0.500522    0.507722  0.433474  0.0771999  0.721393
 0.304884    0.521998  0.904505  0.9822     0.295554
 0.0813203   0.349198  0.137855  0.480715   0.902155

rand([1, 2, "Ahoj!"])

"Ahoj!"

rand(1:10)

rand([1  2 "Ahoj!"; π 7 Any], 2, 2)

2×2 Matrix{Any}:
 7     "Ahoj!"
  Any  Any

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 $0$ a standardní odchylkou $1$ , má dobrý smysl jen pro podtypy AbstractFloat).

randn(10)

10-element Vector{Float64}:
  0.8735074780757865
 -0.0769044165416915
  0.7233061367247997
  0.14524616937488646
  1.511012717318563
 -0.2936472883531717
  0.020498420063912552
  0.2149605134596732
  0.14169559183586325
 -0.5552972309311734

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

randn(Int64, 2, 2)

MethodError: no method matching randn(::Random.TaskLocalRNG, ::Type{Int64})
The function `randn` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  randn(::Random.AbstractRNG, ::Type{T}, ::NTuple{N, Int64} where N) where T
   @ Random /usr/share/julia/stdlib/v1.12/Random/src/normal.jl:269
  randn(::Random.AbstractRNG, ::Type{T}, ::Integer, ::Integer...) where T
   @ Random /usr/share/julia/stdlib/v1.12/Random/src/normal.jl:272
  randn(::Random.AbstractRNG, ::Type{Complex{T}}) where T<:AbstractFloat
   @ Random /usr/share/julia/stdlib/v1.12/Random/src/normal.jl:115
  ...


Stacktrace:
 [1] randn!(rng::Random.TaskLocalRNG, A::Matrix{Int64})
   @ Random /usr/share/julia/stdlib/v1.12/Random/src/normal.jl:230
 [2] randn
   @ /usr/share/julia/stdlib/v1.12/Random/src/normal.jl:272 [inlined]
 [3] randn(::Type{Int64}, ::Int64, ::Int64)
   @ Random /usr/share/julia/stdlib/v1.12/Random/src/normal.jl:274
 [4] top-level scope
   @ In[89]:1
 [5] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

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 $n$ té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)

4×3×2 Array{Float64, 3}:
[:, :, 1] =
 0.636342  0.127265  0.661351
 0.696332  0.400037  0.973281
 0.716652  0.828413  0.379522
 0.625074  0.899572  0.0763026

[:, :, 2] =
 0.726569  0.663072  0.0610617
 0.384091  0.571117  0.463332
 0.317993  0.7373    0.601152
 0.512442  0.934417  0.70682

eltype(A)

Float64

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)

Array{Float64, 3}

length(A)

4 * 3 * 2

length([[1, 2], [3, 4]])

ndims(A)

ndims([1, 2])

ndims([1 2; 3 4])

size(A)

(4, 3, 2)

println(size(A, 1))
println(size(A, 2))
println(size(A, 3))

4
3
2

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

size(A)[2]

axes(A)

(Base.OneTo(4), Base.OneTo(3), Base.OneTo(2))

axes(A, 2)

Base.OneTo(3)

N.B.:

Base.OneTo(4) == 1:4

true

for j = 1:size(A, 2)
    println(A[1, j, 1])
end

0.6363421000914545
0.12726452175658942
0.661350950650962

for j = axes(A, 2)
    println(A[1, j, 1])
end

0.6363421000914545
0.12726452175658942
0.661350950650962

1.2 Cvičení

Vytvořte matici A rozměru $10 \times 10$ mající na antidiagonálách postupně konstanty $1$ , $2$ , $3$ ,..., $19$ .
Vytvořte matici B rozměru $5 \times 10$ , která má všechny prvky nulové, pouze na "okraji" (obdélník) má jedničky.
Vytvořte čtvercovou matici C rozměru $20 \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.
Vytvořte čtvercovou matici D rozměru $10 \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

10×10 Matrix{Int64}:
  1   2   3   4   5   6   7   8   9  10
  2   3   4   5   6   7   8   9  10  11
  3   4   5   6   7   8   9  10  11  12
  4   5   6   7   8   9  10  11  12  13
  5   6   7   8   9  10  11  12  13  14
  6   7   8   9  10  11  12  13  14  15
  7   8   9  10  11  12  13  14  15  16
  8   9  10  11  12  13  14  15  16  17
  9  10  11  12  13  14  15  16  17  18
 10  11  12  13  14  15  16  17  18  19

B = [ i == 1 || j == 1 || i == 5 || j == 10 ? 1 : 0 for i = 1:5, j = 1:10]
B

5×10 Matrix{Int64}:
 1  1  1  1  1  1  1  1  1  1
 1  0  0  0  0  0  0  0  0  1
 1  0  0  0  0  0  0  0  0  1
 1  0  0  0  0  0  0  0  0  1
 1  1  1  1  1  1  1  1  1  1

B = ones(Int64, 5, 10)
B[2:4, 2:9] = zeros(Int64, 3, 8)
B

5×10 Matrix{Int64}:
 1  1  1  1  1  1  1  1  1  1
 1  0  0  0  0  0  0  0  0  1
 1  0  0  0  0  0  0  0  0  1
 1  0  0  0  0  0  0  0  0  1
 1  1  1  1  1  1  1  1  1  1

π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

10×10 Matrix{Float64}:
 0.476129  0.313752  0.174876   0.476975  0.875266   1.0  1.0  1.0  1.0  1.0
 0.867639  0.482137  0.0886191  0.760947  0.423206   1.0  1.0  1.0  1.0  1.0
 0.748895  0.023938  0.149141   0.566933  0.0705881  1.0  1.0  1.0  1.0  1.0
 0.332001  0.132999  0.207124   0.804796  0.654393   1.0  1.0  1.0  1.0  1.0
 0.417478  0.959089  0.0132395  0.513029  0.747776   1.0  1.0  1.0  1.0  1.0
 0.0       0.0       0.0        0.0       0.0        3.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  4.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  5.0

4×4 Matrix{Int64}:
 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

10×10 Matrix{Int64}:
 1  1   1   1    1    1   1   1  1  1
 1  2   3   4    5    6   7   8  9  0
 1  3   6  10   15   21  28  36  0  0
 1  4  10  20   35   56  84   0  0  0
 1  5  15  35   70  126   0   0  0  0
 1  6  21  56  126    0   0   0  0  0
 1  7  28  84    0    0   0   0  0  0
 1  8  36   0    0    0   0   0  0  0
 1  9   0   0    0    0   0   0  0  0
 1  0   0   0    0    0   0   0  0  0

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]

2×2 Matrix{Int64}:
 1  2
 3  4

B = [5 6 7; 8 9 10]

2×3 Matrix{Int64}:
 5  6   7
 8  9  10

C = [1 1; -1 -1]

2×2 Matrix{Int64}:
  1   1
 -1  -1

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

A + A

2×2 Matrix{Int64}:
 2  4
 6  8

A - C

2×2 Matrix{Int64}:
 0  1
 4  5

"Opačná matice":

-A

2×2 Matrix{Int64}:
 -1  -2
 -3  -4

A + B

DimensionMismatch: a has size (2, 2), b has size (2, 3), mismatch at dim 2

Stacktrace:
 [1] throw_promote_shape_mismatch(a::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, b::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, i::Int64)
   @ Base ./indices.jl:135
 [2] promote_shape
   @ ./indices.jl:196 [inlined]
 [3] promote_shape
   @ ./indices.jl:188 [inlined]
 [4] +(A::Matrix{Int64}, Bs::Matrix{Int64})
   @ Base ./arraymath.jl:14
 [5] top-level scope
   @ In[39]:1
 [6] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

[1, 2, 3] + [4, 5, 6]

3-element Vector{Int64}:
 5
 7
 9

[1, 2, 3] - [4, 5, 6]

3-element Vector{Int64}:
 -3
 -3
 -3

[1, 2, 3] + [1 2 3]

DimensionMismatch: a has size (1, 3), b has size (3,), mismatch at dim 1

Stacktrace:
 [1] throw_promote_shape_mismatch(a::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, b::Tuple{Base.OneTo{Int64}}, i::Int64)
   @ Base ./indices.jl:135
 [2] promote_shape
   @ ./indices.jl:196 [inlined]
 [3] promote_shape
   @ ./indices.jl:193 [inlined]
 [4] promote_shape
   @ ./indices.jl:188 [inlined]
 [5] +(A::Vector{Int64}, Bs::Matrix{Int64})
   @ Base ./arraymath.jl:14
 [6] top-level scope
   @ In[42]:1
 [7] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

[1, 2, 3] + [1 2 3]'

3×1 Matrix{Int64}:
 2
 4
 6

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

1 + 2.4

3.4

C + rand(2, 2)

2×2 Matrix{Float64}:
  1.23138    1.95815
 -0.648761  -0.38124

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:

2×2 Matrix{Int64}:
 1  2
 3  4

4 * A

2×2 Matrix{Int64}:
  4   8
 12  16

-2A

2×2 Matrix{Int64}:
 -2  -4
 -6  -8

A * 2

2×2 Matrix{Int64}:
 2  4
 6  8

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

2×2 Matrix{Int64}:
  7  10
 15  22

A * B

2×3 Matrix{Int64}:
 21  24  27
 47  54  61

B * A

DimensionMismatch: incompatible dimensions for matrix multiplication: tried to multiply a matrix of size (2, 3) with a matrix of size (2, 2). The second dimension of the first matrix: 3, does not match the first dimension of the second matrix: 2.

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] _generic_matmatmul!(C::Matrix{Int64}, A::Matrix{Int64}, B::Matrix{Int64}, alpha::Bool, beta::Bool)
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1020
  [5] generic_matmatmul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1011 [inlined]
  [6] generic_matmatmul_wrapper!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:343 [inlined]
  [7] _mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:328 [inlined]
  [8] mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:297 [inlined]
  [9] mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:265 [inlined]
 [10] mul
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:118 [inlined]
 [11] *(A::Matrix{Int64}, B::Matrix{Int64})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:114
 [12] top-level scope
    @ In[52]:1
 [13] eval(m::Module, e::Any)
    @ Core ./boot.jl:489

Inverze `inv`

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

inv([1 2; 3 4])

2×2 Matrix{Float64}:
 -2.0   1.0
  1.5  -0.5

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

2×2 Matrix{Float64}:
 1.0          0.0
 8.88178e-16  1.0

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

2×2 Matrix{Float64}:
 1.0          0.0
 1.11022e-16  1.0

inv([1 1; 1 1])

LinearAlgebra.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[56]:1
 [11] eval(m::Module, e::Any)
    @ Core ./boot.jl:489

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

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[57]:1
 [4] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

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])

2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

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

ratA    = Rational{Int64}[1 1; 1 -1]
ratA

2×2 Matrix{Rational{Int64}}:
 1   1
 1  -1

ratAinv = inv(ratA)

2×2 Matrix{Rational{Int64}}:
 1//2   1//2
 1//2  -1//2

ratA * ratAinv

2×2 Matrix{Rational{Int64}}:
 1  0
 0  1

ratAinv * ratA

2×2 Matrix{Rational{Int64}}:
 1  0
 0  1

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

m = [ 1//1 2//1; 3//1 4//1 ]
invm = inv(m)

2×2 Matrix{Rational{Int64}}:
 -2      1
 3//2  -1//2

m * invm

2×2 Matrix{Rational{Int64}}:
 1  0
 0  1

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

8×8 Matrix{Rational{Int64}}:
  1    1//2  1//3   1//4   1//5   1//6   1//7   1//8
 1//2  1//3  1//4   1//5   1//6   1//7   1//8   1//9
 1//3  1//4  1//5   1//6   1//7   1//8   1//9   1//10
 1//4  1//5  1//6   1//7   1//8   1//9   1//10  1//11
 1//5  1//6  1//7   1//8   1//9   1//10  1//11  1//12
 1//6  1//7  1//8   1//9   1//10  1//11  1//12  1//13
 1//7  1//8  1//9   1//10  1//11  1//12  1//13  1//14
 1//8  1//9  1//10  1//11  1//12  1//13  1//14  1//15

inv(H)

8×8 Matrix{Rational{Int64}}:
      64      -2016       20160      -92400  …       192192      -51480
   -2016      84672     -952560     4656960       -10594584     2882880
   20160    -952560    11430720   -58212000       141261120   -38918880
  -92400    4656960   -58212000   304920000      -776936160   216216000
  221760  -11642400   149688000  -800415000      2118916800  -594594000
 -288288   15567552  -204324120  1109908800  …  -3030051024   856215360
  192192  -10594584   141261120  -776936160      2175421248  -618377760
  -51480    2882880   -38918880   216216000      -618377760   176679360

H * inv(H)

8×8 Matrix{Rational{Int64}}:
 1  0  0  0  0  0  0  0
 0  1  0  0  0  0  0  0
 0  0  1  0  0  0  0  0
 0  0  0  1  0  0  0  0
 0  0  0  0  1  0  0  0
 0  0  0  0  0  1  0  0
 0  0  0  0  0  0  1  0
 0  0  0  0  0  0  0  1

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

15×15 Matrix{Float64}:
 1.0        0.5        0.333333   …  0.0769231  0.0714286  0.0666667
 0.5        0.333333   0.25          0.0714286  0.0666667  0.0625
 0.333333   0.25       0.2           0.0666667  0.0625     0.0588235
 0.25       0.2        0.166667      0.0625     0.0588235  0.0555556
 0.2        0.166667   0.142857      0.0588235  0.0555556  0.0526316
 0.166667   0.142857   0.125      …  0.0555556  0.0526316  0.05
 0.142857   0.125      0.111111      0.0526316  0.05       0.047619
 0.125      0.111111   0.1           0.05       0.047619   0.0454545
 0.111111   0.1        0.0909091     0.047619   0.0454545  0.0434783
 0.1        0.0909091  0.0833333     0.0454545  0.0434783  0.0416667
 0.0909091  0.0833333  0.0769231  …  0.0434783  0.0416667  0.04
 0.0833333  0.0769231  0.0714286     0.0416667  0.04       0.0384615
 0.0769231  0.0714286  0.0666667     0.04       0.0384615  0.037037
 0.0714286  0.0666667  0.0625        0.0384615  0.037037   0.0357143
 0.0666667  0.0625     0.0588235     0.037037   0.0357143  0.0344828

inv(H) * H

15×15 Matrix{Float64}:
  1.0          -5.58794e-9    7.82398e-9  …  -1.22192e-9  -1.86265e-9
  3.59734e-7    0.999998     -3.93468e-6     -1.13755e-6  -4.76837e-7
  1.61656e-5    4.24385e-5    1.00008        -6.2141e-6    0.0
 -0.000196082  -0.000553131  -0.00133044      0.00013677  -0.000244141
 -0.00195368    0.00274658    0.0077043      -0.00278265   0.0
  0.088988     -0.0244141    -0.0708378   …   0.0149771   -0.015625
 -0.0282454     0.0546875     0.217265       -0.0334611    0.15625
 -0.393002     -0.84375      -0.806325        0.0161129   -0.4375
  0.350812      1.21875       0.951084        0.258436     0.75
  1.70837      -0.0625       -0.812464       -0.490984    -1.75
 -1.44207      -0.75         -1.08491     …   0.186639     1.40625
  0.152157     -1.10938      -0.237068       -0.879307    -1.65625
  0.0353458     0.617188      0.219147        0.966645     1.06233
  0.00351884   -0.0712891     0.0546616       0.777002    -0.273587
 -0.00547683    0.0078125    -0.0317454       0.0336043    1.02711

big"1" // 2 |> typeof

Rational{BigInt}

H = [ big"1" // (i + j - 1) for i = 1:15, j = 1:15 ]

15×15 Matrix{Rational{BigInt}}:
  1     1//2   1//3   1//4   1//5   …  1//11  1//12  1//13  1//14  1//15
 1//2   1//3   1//4   1//5   1//6      1//12  1//13  1//14  1//15  1//16
 1//3   1//4   1//5   1//6   1//7      1//13  1//14  1//15  1//16  1//17
 1//4   1//5   1//6   1//7   1//8      1//14  1//15  1//16  1//17  1//18
 1//5   1//6   1//7   1//8   1//9      1//15  1//16  1//17  1//18  1//19
 1//6   1//7   1//8   1//9   1//10  …  1//16  1//17  1//18  1//19  1//20
 1//7   1//8   1//9   1//10  1//11     1//17  1//18  1//19  1//20  1//21
 1//8   1//9   1//10  1//11  1//12     1//18  1//19  1//20  1//21  1//22
 1//9   1//10  1//11  1//12  1//13     1//19  1//20  1//21  1//22  1//23
 1//10  1//11  1//12  1//13  1//14     1//20  1//21  1//22  1//23  1//24
 1//11  1//12  1//13  1//14  1//15  …  1//21  1//22  1//23  1//24  1//25
 1//12  1//13  1//14  1//15  1//16     1//22  1//23  1//24  1//25  1//26
 1//13  1//14  1//15  1//16  1//17     1//23  1//24  1//25  1//26  1//27
 1//14  1//15  1//16  1//17  1//18     1//24  1//25  1//26  1//27  1//28
 1//15  1//16  1//17  1//18  1//19     1//25  1//26  1//27  1//28  1//29

inv(H)

15×15 Matrix{Rational{BigInt}}:
          225           -25200            928200  …            1163381400
       -25200          3763200        -155937600            -244310094000
       928200       -155937600        6892441920           12704124888000
    -16707600       2994001920     -137848838400         -287960164128000
    174594420     -32590958400     1543414672800         3563507031084000
  -1163962800     223480857600   -10803902709600  …    -27082653436238400
   5237832600   -1026615189600    50418212644800       135413267181192000
 -16461759600    3277719244800  -162984589447680      -464274058906944000
  36810323550   -7420961227680   372734643481200      1117159454244834000
 -58896517680   11993472691200  -607419419006400     -1903308699824532000
  66927861000  -13742520792000   700868560392000  …   2283970439789438400
 -52731042000   10903156992000  -559370893536000     -1887578875859040000
  27379579500   -5694952536000   293669719106400      1022438557756980000
  -8424486000    1761279206400   -91228758894000      -326696343898680000
   1163381400    -244310094000    12704124888000        46670906271240000

inv(H) * H

15×15 Matrix{Rational{BigInt}}:
 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
 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  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  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  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  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  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  0  0  1

Dělení `/`, resp. `\`

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

A / A

2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

A \ B

2×3 Matrix{Float64}:
 -2.0  -3.0  -4.0
  3.5   4.5   5.5

inv(A) * B

2×3 Matrix{Float64}:
 -2.0  -3.0  -4.0
  3.5   4.5   5.5

Umocňování `^`

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

A ^ 2

2×2 Matrix{Int64}:
  7  10
 15  22

A * A

2×2 Matrix{Int64}:
  7  10
 15  22

A ^ 6

2×2 Matrix{Int64}:
  5743   8370
 12555  18298

A ^ -1

2×2 Matrix{Float64}:
 -2.0   1.0
  1.5  -0.5

inv(A)

2×2 Matrix{Float64}:
 -2.0   1.0
  1.5  -0.5

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):

2×3 Matrix{Int64}:
 5  6   7
 8  9  10

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

2×3 Matrix{Bool}:
 0  1  0
 1  0  1

iseven.(B)

2×3 BitMatrix:
 0  1  0
 1  0  1

map(iseven, B)

2×3 Matrix{Bool}:
 0  1  0
 1  0  1

iseven(B)

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

Closest candidates are:
  iseven(::BigInt)
   @ Base gmp.jl:368
  iseven(::Missing)
   @ Base missing.jl:101
  iseven(::AbstractFloat)
   @ Base float.jl:957
  ...


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

2×3 Matrix{Int64}:
 5  6   7
 8  9  10

B .* B

2×3 Matrix{Int64}:
 25  36   49
 64  81  100

B * B

DimensionMismatch: incompatible dimensions for matrix multiplication: tried to multiply a matrix of size (2, 3) with a matrix of size (2, 3). The second dimension of the first matrix: 3, does not match the first dimension of the second matrix: 2.

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] _generic_matmatmul!(C::Matrix{Int64}, A::Matrix{Int64}, B::Matrix{Int64}, alpha::Bool, beta::Bool)
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1020
  [5] generic_matmatmul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:1011 [inlined]
  [6] generic_matmatmul_wrapper!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:343 [inlined]
  [7] _mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:328 [inlined]
  [8] mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:297 [inlined]
  [9] mul!
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:265 [inlined]
 [10] mul
    @ /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:118 [inlined]
 [11] *(A::Matrix{Int64}, B::Matrix{Int64})
    @ LinearAlgebra /usr/share/julia/stdlib/v1.12/LinearAlgebra/src/matmul.jl:114
 [12] top-level scope
    @ In[100]:1
 [13] eval(m::Module, e::Any)
    @ Core ./boot.jl:489

Operand ani nemusí být maticový:

4 .* B

2×3 Matrix{Int64}:
 20  24  28
 32  36  40

4 * B

2×3 Matrix{Int64}:
 20  24  28
 32  36  40

4 .+ B

2×3 Matrix{Int64}:
  9  10  11
 12  13  14

4 + B

MethodError: no method matching +(::Int64, ::Matrix{Int64})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
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
  +(::Array, ::Array...)
   @ Base arraymath.jl:12
  +(::Real, ::Complex{Bool})
   @ Base complex.jl:322
  ...


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

2×3 Matrix{Int64}:
 5  6   7
 8  9  10

B[:, 2:3] .= 0
B

2×3 Matrix{Int64}:
 5  0  0
 8  0  0

B .<< 1

2×3 Matrix{Int64}:
 10  0  0
 16  0  0

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

sin.(A)

2×2 Matrix{Float64}:
 0.841471   0.909297
 0.14112   -0.756802

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)

2×2 Matrix{Float64}:
 -0.465581  -0.148424
 -0.222637  -0.688218

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)

2×2 Matrix{Float64}:
 -0.467181  -0.150756
 -0.226135  -0.693316

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

#53 (generic function with 1 method)

2×2 Matrix{Int64}:
 1  2
 3  4

f.(A)

2×2 Matrix{Int64}:
  1   8
 27  64

map(f, A)

2×2 Matrix{Int64}:
  1   8
 27  64

2×2 Matrix{Int64}:
 1  2
 3  4

2×2 Matrix{Int64}:
  1   1
 -1  -1

min.(A, C)

2×2 Matrix{Int64}:
  1   1
 -1  -1

min.(A, 3)

2×2 Matrix{Int64}:
 1  2
 3  3

min(A, 3)

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

Closest candidates are:
  isless(::Missing, ::Any)
   @ Base missing.jl:87
  isless(::Real, ::AbstractFloat)
   @ Base operators.jl:223
  isless(::Real, ::Real)
   @ Base operators.jl:477
  ...


Stacktrace:
 [1] min(x::Matrix{Int64}, y::Int64)
   @ Base ./operators.jl:545
 [2] top-level scope
   @ In[126]:1
 [3] eval(m::Module, e::Any)
   @ Core ./boot.jl:489

min.(A, C, ones(Int64, 2, 2))

2×2 Matrix{Int64}:
  1   1
 -1  -1

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.

2×2 Matrix{Int64}:
 1  2
 3  4

C = [0 3; 4 -1]

2×2 Matrix{Int64}:
 0   3
 4  -1

A .< C

2×2 BitMatrix:
 0  1
 1  0

mask = A .< C

2×2 BitMatrix:
 0  1
 1  0

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

isless.(A, C)

2×2 BitMatrix:
 0  1
 1  0

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

A[mask]

2-element Vector{Int64}:
 3
 2

A[A .< C]

2-element Vector{Int64}:
 3
 2

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

2×2 Matrix{Int64}:
 0   3
 4  -1

C .<= 2

2×2 BitMatrix:
 1  0
 0  1

C[C .<= 5]

4-element Vector{Int64}:
  0
  4
  3
 -1

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

mask = iseven.(A)

2×2 BitMatrix:
 0  1
 0  1

A[mask]

2-element Vector{Int64}:
 2
 4

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

filter(iseven, A)

2-element Vector{Int64}:
 2
 4

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)

2×1000000 Matrix{Float64}:
 0.718102  0.678207  0.139518  0.217687  …  0.52297   0.04458   0.103312
 0.61242   0.390676  0.133131  0.958779     0.206085  0.375988  0.665905

@benchmark f1($mat)

BenchmarkTools.Trial: 477 samples with 1 evaluation per sample.
 Range (min … max):   6.612 ms … 16.704 ms  ┊ GC (min … max): 0.00% … 8.10%
 Time  (median):      9.718 ms              ┊ GC (median):    4.67%
 Time  (mean ± σ):   10.421 ms ±  2.278 ms  ┊ GC (mean ± σ):  4.75% ± 2.84%

           ██▅▄▄▃▂ ▁            ▁                              
  ▄▄▁▂▃▃▆▄▆███████████▄▆▄▄█▄▄▆▅▄█▆▆▆█▆█▅▅▄▆▄▄▄▆▄▃▄▇▇▄▆▅▄▄▄▃▂▃ ▄
  6.61 ms         Histogram: frequency by time        15.4 ms <

 Memory estimate: 15.26 MiB, allocs estimate: 6.

@benchmark f2($mat)

BenchmarkTools.Trial: 2789 samples with 1 evaluation per sample.
 Range (min … max):  1.150 ms …   5.784 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.341 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.781 ms ± 795.238 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃▆█▇▆▄▃▂▁▁               ▁  ▁▁▁▁ ▁  ▁▁▁     ▁        ▁      ▁
  ███████████▇▇▇▇▇▆▆▇▆▇▄▇▇▇█▇▇██████▇████▇██████████▇▇▇█▇▇▆▇▆ █
  1.15 ms      Histogram: log(frequency) by time      3.89 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark f3($mat)

BenchmarkTools.Trial: 379 samples with 1 evaluation per sample.
 Range (min … max):   7.932 ms … 19.438 ms  ┊ GC (min … max): 0.00% … 4.77%
 Time  (median):     12.641 ms              ┊ GC (median):    7.11%
 Time  (mean ± σ):   13.181 ms ±  2.652 ms  ┊ GC (mean ± σ):  7.19% ± 3.86%

          ▁▁ ▁  ▆▇█▆▅▁▁▇▅▃▁▃▅  ▁   ▃▁  ▁   ▁▁▃ ▂▁  ▄▁          
  ▇▃▁▁▃▄▅▄██▇█▇▇█████████████▁▇█▇▄▅███▆█▇▇████▆██▆▇███▅▇▆▆▄▆▄ ▅
  7.93 ms         Histogram: frequency by time        18.4 ms <

 Memory estimate: 22.89 MiB, allocs estimate: 9.

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

my_filter (generic function with 1 method)

@benchmark my_filter($mat, 0.5)

BenchmarkTools.Trial: 455 samples with 1 evaluation per sample.
 Range (min … max):   7.599 ms … 17.812 ms  ┊ GC (min … max): 0.00% … 3.31%
 Time  (median):     10.660 ms              ┊ GC (median):    4.19%
 Time  (mean ± σ):   10.994 ms ±  1.942 ms  ┊ GC (mean ± σ):  4.58% ± 4.59%

   ▁       ▂▄▃▁   ▁█▄▁▃▂▃ ▁    ▁ ▁ ▁▃                          
  ██▄▄▇▃▅▃▇████▆▆▇███████▃█▅▆███▇█▇██▇▇▆▆▆▄▇▇▆▄▆▄▄▃▄▃▃▃▃▄▁▃▄▃ ▄
  7.6 ms          Histogram: frequency by time        15.6 ms <

 Memory estimate: 11.69 MiB, allocs estimate: 7.

@benchmark filter((x) -> x < 0.5, $mat)

BenchmarkTools.Trial: 914 samples with 1 evaluation per sample.
 Range (min … max):  2.804 ms … 12.084 ms  ┊ GC (min … max):  0.00% … 4.74%
 Time  (median):     4.629 ms              ┊ GC (median):    11.69%
 Time  (mean ± σ):   5.453 ms ±  1.810 ms  ┊ GC (mean ± σ):  11.32% ± 5.69%

        ▆▆▁  █▆█▃                                             
  ▃▄▃▂▅▄███▅██████▆▅▅▄▆▄▄▃▄▄▃▃▃▄▄▃▅▃▅▅▆▆▅▄▅▅▆▅▄▅▄▄▄▃▄▂▁▃▁▂▂▂ ▄
  2.8 ms         Histogram: frequency by time        9.96 ms <

 Memory estimate: 26.71 MiB, allocs estimate: 5.

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

Pro reálnou čtvercovou matici $\mathbb{A}$ definujeme její operátorovou normu předpisem $\|\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 $\mathbb{R}^n$ .

Jinak řečeno, každá hodnota $\| \mathbb{A} x\|$ pro jednotkový vektor $x$ představuje spodní odhad $\|\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\in\mathbb{R}^n$ je $\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

mc_op_norm

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)

0.0

Diagonální matice.

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

3.9999972285050065

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

mc_op_norm([1. 4; 5 6], 10^6)

8.683348976417081

sqrt(39 + 5sqrt(53))

8.683348976426238

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

J(n) = Float64[ abs(j-k) == 1 ? 1 : 0 for j=1:n, k=1:n ]

J (generic function with 1 method)

J(5)

5×5 Matrix{Float64}:
 0.0  1.0  0.0  0.0  0.0
 1.0  0.0  1.0  0.0  0.0
 0.0  1.0  0.0  1.0  0.0
 0.0  0.0  1.0  0.0  1.0
 0.0  0.0  0.0  1.0  0.0

mc_op_norm(J(2), 10^6)

1.0

mc_op_norm(J(20), 10^6)

1.9002363092758066

mc_op_norm(J(20), 10^7)

Tossing coins furiously... 100%|█████████████████████████| Time: 0:00:04

1.9205840018257356

mc_op_norm(J(30), 10^9)

Tossing coins furiously...   6%|█▋                       |  ETA: 0:09:40

InterruptException:

Stacktrace:
  [1] mapreduce_impl(f::typeof(identity), op::typeof(Base.add_sum), A::Vector{Float64}, ifirst::Int64, ilast::Int64, blksize::Int64)
    @ Base ./reduce.jl:245
  [2] mapreduce_impl
    @ ./reduce.jl:269 [inlined]
  [3] _mapreduce
    @ ./reduce.jl:436 [inlined]
  [4] _mapreduce_dim
    @ ./reducedim.jl:334 [inlined]
  [5] mapreduce
    @ ./reducedim.jl:326 [inlined]
  [6] _sum
    @ ./reducedim.jl:984 [inlined]
  [7] _sum
    @ ./reducedim.jl:983 [inlined]
  [8] sum
    @ ./reducedim.jl:979 [inlined]
  [9] macro expansion
    @ ./In[154]:41 [inlined]
 [10] macro expansion
    @ ~/.julia/packages/ProgressMeter/N660J/src/ProgressMeter.jl:1026 [inlined]
 [11] mc_op_norm(A::Matrix{Float64}, trials::Int64)
    @ Main ./In[154]:39
 [12] top-level scope
    @ In[167]:1
 [13] eval(m::Module, e::Any)
    @ Core ./boot.jl:489

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

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

50×50 Matrix{Float64}:
 1.0        0.5        0.333333   …  0.0208333  0.0204082  0.02
 0.5        0.333333   0.25          0.0204082  0.02       0.0196078
 0.333333   0.25       0.2           0.02       0.0196078  0.0192308
 0.25       0.2        0.166667      0.0196078  0.0192308  0.0188679
 0.2        0.166667   0.142857      0.0192308  0.0188679  0.0185185
 0.166667   0.142857   0.125      …  0.0188679  0.0185185  0.0181818
 0.142857   0.125      0.111111      0.0185185  0.0181818  0.0178571
 0.125      0.111111   0.1           0.0181818  0.0178571  0.0175439
 0.111111   0.1        0.0909091     0.0178571  0.0175439  0.0172414
 0.1        0.0909091  0.0833333     0.0175439  0.0172414  0.0169492
 0.0909091  0.0833333  0.0769231  …  0.0172414  0.0169492  0.0166667
 0.0833333  0.0769231  0.0714286     0.0169492  0.0166667  0.0163934
 0.0769231  0.0714286  0.0666667     0.0166667  0.0163934  0.016129
 ⋮                                ⋱                        
 0.025641   0.025      0.0243902     0.0116279  0.0114943  0.0113636
 0.025      0.0243902  0.0238095     0.0114943  0.0113636  0.011236
 0.0243902  0.0238095  0.0232558  …  0.0113636  0.011236   0.0111111
 0.0238095  0.0232558  0.0227273     0.011236   0.0111111  0.010989
 0.0232558  0.0227273  0.0222222     0.0111111  0.010989   0.0108696
 0.0227273  0.0222222  0.0217391     0.010989   0.0108696  0.0107527
 0.0222222  0.0217391  0.0212766     0.0108696  0.0107527  0.0106383
 0.0217391  0.0212766  0.0208333  …  0.0107527  0.0106383  0.0105263
 0.0212766  0.0208333  0.0204082     0.0106383  0.0105263  0.0104167
 0.0208333  0.0204082  0.02          0.0105263  0.0104167  0.0103093
 0.0204082  0.02       0.0196078     0.0104167  0.0103093  0.0102041
 0.02       0.0196078  0.0192308     0.0103093  0.0102041  0.010101

mc_op_norm(H, 10_000_000)

Tossing coins furiously... 100%|█████████████████████████| Time: 0:00:09

1.8545276741011454

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 $2$ ; 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

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

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

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

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

Indexování polí

Konkatenace polí (vcat a hcat)

push! a pop!

fill a fill!

zeros a ones

copy a deepcopy

rand a randn

range