Pole a výpočty lineární algebry

Vektory a matice mají ve fyzice nepřebené množství užití. Proto digitální počítače a jazyky určené k jejich programování braly v úvahu potřeby takových výpočtů, jako je například řešení soustav rovnic.

Zde si na několik příkladech ilustrujeme elementární funkce z lineární algebry, jak se v nich pracuje s poli a jak používat funkce z knihoven.

Jako opakování ověříme, že pro tři konkrétní vektory platí

\[\vec a \times (~\vec b \times \vec c~) = \vec b ~(\vec a\cdot \vec c)-\vec c ~(\vec a\cdot \vec b)\]

Napíšeme si fuknci pro skalární součin definovaný vztahem

\[\vec a \cdot \vec b := \sum_i a_i b_i\]

a také funkci pro vektorový součin.

[4]:
import numpy as np

def scalar_product(a,b):
    "Počítá skalární součin dvou vektorů"
    dim = len(a)

    s = 0.0
    for i in range(dim):
        s = s + a[i]*b[i]

    return s

def vector_product(a,b):
    "Počítá vektorový součin dvou vektorů"
    assert len(a)==len(b)==3

    c = np.zeros(3)
    c[0] = a[1] * b[2] - a[2] * b[1]
    c[1] = a[2] * b[0] - a[0] * b[2]
    c[2] = a[0] * b[1] - a[1] * b[0]

    return c

a = np.array([ 1,2,0 ])
b = np.array([ 2,0,1 ])
c = np.array([ 0,1,2 ])

print( vector_product(a, vector_product(b,c)) )
print( b*scalar_product(a,c) - c*scalar_product(a,b) )
[ 4. -2. -2.]
[ 4. -2. -2.]

Problém: Povšimněte si, že obě funkce vracejí správnou hodnotu i tehdy, jsou-li jejich argumenty seznamy typu list nikoli pole numpy.ndarray. Přesto si rozmystlete, na které operaci by program ohlásil chybu, kdybychom proměnné a,b,c inicializovali jako typ list, tedy a = [ 1,2,3 ] atd. Vyzkoušejte, zda Váš odhad vyšel.

Co by se změnilo, kdybychom na čtvrém řádku měli s = 0, tedy celé číslo místo čísla reálného?

Samozřejmě, že v knihovně numpy existují funkce, které počítají totéž. Jmenují se numpy.dot a numpy.cross. Zejména ta první je velmi užitečná a tak ji zanedlouho o něco podrobněji zmíníme znovu.

Cvičení: Napište funkci, která na seznamu vektorů provede Gramm-Schmidtovu ortogonalizaci. (Tento postup akceptuje první vektor v seznamu a od každého dalšího odečítá vhodné násobky těch předchozích, již ortogonálních tak, aby byl každý vektor ortogonální k těm předchozím.)

  • Rozmyslete si, jak v programu uložíte počáteční a výsledný seznam vektorů.

  • Rozmyslete si, zda chcete mít funkci, která vrací seznam nebo proceduru, která modifikuje odkazem předávaný argument.

2D a 3D vektory

I u těchto krátkých vektorů sice numpy nepřináší zrychlení výpočtů, ale z hlediska psaní kódu je práce s nimi je pohodlější, protože se hodí, že vektory můžeme sčítat a násobit. Jako příklad si uvedeme klasickou demonstraci vývoje vertikální složky rychlosti ozářením volně padajícího objektu stroboskopem:

[21]:
# Program vykreslí polohy hmotného bodu
# podléhajícímu volnému pádu
# a to v ekvidistantních časových okamžicích

import numpy as np
import matplotlib.pyplot as plt

x0 = np.array([0, 0])
v0 = np.array([1, 9.81])
g  = np.array([0,-9.81])

seznam_casu = np.linspace(0, 2, 30)

xy = [ x0 + v0 * t + 0.5 * g * t**2 for t in seznam_casu ]
xy = np.array(xy)

plt.rcParams['axes.facecolor'] = 'black'
plt.axis('equal')
plt.plot( xy[:,0], xy[:,1],'o', color='yellow')
plt.show()
../_images/Prednaska_0090_Pole_a_linearni_algebra_4_0.png

Diskuse:

  • používáme 2D vektory

  • chceme vyrobit seznam takových vektorů představující soubor poloh volně padajícího bodu

  • bohužel pak nemůžeme napsat v0*seznam_casu, protože definovány jsou jen násobení stejně dlouhých vektorů nebo vektoru a skaláru

  • proto je zkonstruujeme jako xy = [f(t) for t in seznam_casu]

  • funkce plot vyžaduje zvlášť x-ové a zvlášť y-ové složky polí, proto používáme řezy, např. xy[:,0]

Cvičení: V kódu se používá seznamový výraz (list comprehension)

xy = [ x0 + v0 * t + 0.5 * g * t**2 for t in seznam_casu ]
xy = np.array(xy)

Přepište tuto pohodlnou konstrukci do podoby cyklu, abyste si vyzkoušeli, jaké to je programovat v běžných jazycích, které tuto vychytávku nemají.

  • vytvořte pole xy naplněné hodnotami \(\vec x_0\) (použijte np.full((len(seznam_casu),2), x0)), tedy pole xy bude obsahovat tolik 2D vektorů, kolik je počet prvků proměnné seznam_casu

  • pro všechny časy t z proměnné seznam_casu

    • k příslušnému prvku pole xy přičtěte \(\vec v_0 t+ \tfrac12 \vec g t^2\)

V uvedeném návodu je drobná zrada: - pokud budete přičítat příkazem xy[i] = xy[i] + v0 * t + .... program poběží, ale obrázek nebude dobře

  • pokud budete přičítat příkazem

    xy[i] += v0 * t + ....
    

    program ohlásí chybu Cannot cast ufunc 'add' output ...

Rozmyslete si, v čem chyba spočívá a opravte ji.

Matice a vektory

Jak víme, klíčovou operací je násobení matice vektorem. Knihovna numpy to umí skrze funkci numpy.dot, ale my se to musíme naučit sami:

[33]:
import numpy as np

def matrix_dot_vector(a, x):
    dim_x = len(x)
    dim_y = len(a)

    assert dim_x == len(a[0])

    y = np.zeros(dim_y)

    for i in range(dim_y):
        s = 0.0
        for j in range(dim_x):
            s = s + a[i,j]*x[j]
        y[i] = s

    return y

matice = np.array( [
    [1,2,3,4],
    [1,0,0,1]
])

matrix_dot_vector(matice,[-1,3,-3,1])
[33]:
array([0., 0.])

JIž víme, že knihovna numpy naučí operátory +, * atp. pracovat s poli. Základní podoby jsou

  • pole * 4, kdy se každý prvek pole vynásobí čtyřmi.

  • pole * pole2, kdy obě pole jsou stejně dlouhá a násobení probího prvek po prvku.

Proto výraz matice*matice2 předpokládá, že obě matice mají stejný tvar a vynásobí je prvek po prvku, tedy podle vzorečku \(c_{ij}=a_{ij} b_{ij}\).

Při práci s poli s více indexy, kde matice jsou nejjednodušším příkladem, je nutné přinejměnším vědět, že je pokud * vložíte mezi matici a vektor, nemusí program nahlásit chybu! Knihovna numpy to pochopí tak, že operaci * použije na každý řádek. Pozor, nejde o násobení \({\mathbb A}\cdot \vec x\) z lineární algebry a tedy dostanete stejný výsledek pro obě pořadí: matice*pole i pole*matice.

[38]:
print( 'matice = \n' , matice , "\n" )
print( 'matice*2 =\n', matice * 2 , "\n" )
print( 'matice * [1,2,3,4] =\n', matice * [1,2,3,4], "\n" )
print( '[1,2,3,4] * matice =\n', [1,2,3,4] * matice , "\n" )
print( 'matice * matice =\n', matice * matice , "\n" )
matice =
 [[1 2 3 4]
 [1 0 0 1]]

matice*2 =
 [[2 4 6 8]
 [2 0 0 2]]

matice * [1,2,3,4] =
 [[ 1  4  9 16]
 [ 1  0  0  4]]

[1,2,3,4] * matice =
 [[ 1  4  9 16]
 [ 1  0  0  4]]

matice * matice =
 [[ 1  4  9 16]
 [ 1  0  0  1]]

Knihovna numpy obsahuje velké množství funkcí, které realizují operace z lineární algebry. Z hledniska našeho předmětu - Nehrozí, že bychom vás u zkoušky trápili dotazem: Jak se jmenuje funkce knihovny numpy, která vypočte skalární součin dvou vektorů - Potřebujeme v přednášce praktickou ukázku toho, co to je metoda objektu. Následující řádky se tedy hodí, kdybychom se vás zeptali: Uveďte nějaký příklad proměnné o vhodném typu, kde lze dát příklad volání její metody, a význam takové operace.

Zejména je důležitá funkce numpy.dot. Realizuje operace jako jsou - skalární součin vektorů \(s = x_i y_i\) - součin matice a vektoru \(y_i = \sum_j a_{ij} x_j\) - součin matice a matice \(c_{ik} = \sum_j a_{ij} b_{jk}\)

Máme tři možnosti

  • funkce dot z knihovny numpy, která má dva argumenty

    • s = np.dot(vec_x, vec_y)

    • vec_y = np.dot(mat_A, vec_x)

    • mat_C = np.dot(mat_A, mat_B)

  • použít metodu dot prvního operandu, která jako argument dostane druhý operand

    • s = vec_x.dot(vec_y)

    • vec_y = mat_A.dot(vec_x)

    • mat_C = mat_A.dot(mat_B)

  • použít operátor @ (ve významu dot)

    • s = vec_x @ vec_y

    • vec_y = mat_A @ vec_x

    • mat_C = mat_A @ mat_B

Cvičení: Napište kód testující devět výše uvedených operací.

Gaussova-Jordanova eliminace

Řešení soustavy lineárních rovnic

\[{\mathbf A}\vec x = \vec b\]

si ukážeme za použití Gaussovy-Jordanovy eliminace. Ta pomocí sady ekvivalentních úprav redukuje soustavu rovnic do tvaru, kdy výsledná matice soustavy \(\mathbf A\) je jednotková matice a tedy výsledná pravá strana je hledaným řešením. Jak je obvyklé, spojíme \({\mathbf A}\) a \(\vec b\) do jedné obdélníkové matice \({\mathbf A'}\).

Ekvivalentní úpravy jsou tři: - přehazování rovnic, tedy řádků \({\mathbf A'}\), - násobení rovnic nenulovoým činitelem, - odečítání násobku jedné rovnice od druhé.

Tato metoda je sice pomalejší, než Gaussova metoda se zpětnou substitucí, ale pro naše potřeby dostatečně dobře ilustruje operace, které s maticemi při řešení soustav rovnic provádíme.

Nejprve si ukažme tyto kroky postupně s použitím řezů matic a výhod interaktivního režimu.

[ ]:
import numpy as np
# Nastavíme srozumitelný formát výpisu čísel
np.set_printoptions(floatmode="fixed",precision=4)

# Společná matice soustavy a vektoru pravé strany
matice = np.array( [
        [3, 2, 1, 1,        3],
        [0, 2, 3, -1,      -4],
        [1, 9, 0, 4,       -5],
        [2, 3, 1, 1.0,      0] ] )

# Výpis
matice
array([[ 3.0000,  2.0000,  1.0000,  1.0000,  3.0000],
       [ 0.0000,  2.0000,  3.0000, -1.0000, -4.0000],
       [ 1.0000,  9.0000,  0.0000,  4.0000, -5.0000],
       [ 2.0000,  3.0000,  1.0000,  1.0000,  0.0000]])

Abychom získali na diagonále jedničku, je třeba první řádek vydělit třemi. První řádek matice můžeme zapsat jako matice[0] nebo matice[0,:].

[ ]:
matice[0]
array([3.0000, 2.0000, 1.0000, 1.0000, 3.0000])
[ ]:
matice[0] /= 3

matice
array([[ 1.0000,  0.6667,  0.3333,  0.3333,  1.0000],
       [ 0.0000,  2.0000,  3.0000, -1.0000, -4.0000],
       [ 1.0000,  9.0000,  0.0000,  4.0000, -5.0000],
       [ 2.0000,  3.0000,  1.0000,  1.0000,  0.0000]])
[ ]:
matice[2] -= matice[0]
matice[3] -= 2*matice[0]

matice
array([[ 1.0000,  0.6667,  0.3333,  0.3333,  1.0000],
       [ 0.0000,  2.0000,  3.0000, -1.0000, -4.0000],
       [ 0.0000,  8.3333, -0.3333,  3.6667, -6.0000],
       [ 0.0000,  1.6667,  0.3333,  0.3333, -2.0000]])
[ ]:
# Druhý řádek
matice[1] /= matice[1,1]

matice
array([[ 1.0000,  0.6667,  0.3333,  0.3333,  1.0000],
       [ 0.0000,  1.0000,  1.5000, -0.5000, -2.0000],
       [ 0.0000,  8.3333, -0.3333,  3.6667, -6.0000],
       [ 0.0000,  1.6667,  0.3333,  0.3333, -2.0000]])
[ ]:
matice[0] -= matice[0,1] * matice[1]
matice[2] -= matice[2,1] * matice[1]
matice[3] -= matice[3,1] * matice[1]

matice
array([[  1.0000,   0.0000,  -0.6667,   0.6667,   2.3333],
       [  0.0000,   1.0000,   1.5000,  -0.5000,  -2.0000],
       [  0.0000,   0.0000, -12.8333,   7.8333,  10.6667],
       [  0.0000,   0.0000,  -2.1667,   1.1667,   1.3333]])
[ ]:
# Třetí řádek
matice[2] /= matice[2,2]

matice[0] -= matice[0,2] * matice[2]
matice[1] -= matice[1,2] * matice[2]
matice[3] -= matice[3,2] * matice[2]

matice
array([[ 1.0000,  0.0000,  0.0000,  0.2597,  1.7792],
       [ 0.0000,  1.0000,  0.0000,  0.4156, -0.7532],
       [-0.0000, -0.0000,  1.0000, -0.6104, -0.8312],
       [ 0.0000,  0.0000,  0.0000, -0.1558, -0.4675]])
[ ]:
# Čtvrtý řádek
matice[3] /= matice[3,3]

matice[0] -= matice[0,3] * matice[3]
matice[1] -= matice[1,3] * matice[3]
matice[2] -= matice[2,3] * matice[3]

# Obnovíme výchozí formátování
np.set_printoptions(floatmode="maxprec_equal",precision=8)   # set default

# Jsme hotovi
matice
array([[ 1.,  0.,  0.,  0.,  1.],
       [ 0.,  1.,  0.,  0., -2.],
       [-0., -0.,  1.,  0.,  1.],
       [-0., -0., -0.,  1.,  3.]])

Následující funkce GJelim(matice, prava_strana) provádí tyto kroky najednou. Obsahuje navíc kód, který prohazuje řádky. Lze snadno vidět, že by předchozí postup havaroval, kdyby se na diagonále objevila nula, protože tu násobením na jedničku nepředěláme. Protože jsou rovnice záměnné, můžeme přehodit řádky. Pro regulární matici nějaký musí dát nenulový diagonální prvek. Také nás nepřekvapí, že by mohlo vadit i číslo velmi blízké nule. Podrobná analýza ale ukazuje, že je dobré řádky prohazovat vždy, když lze pod diagonálou najít nějaký větší prvek. Protže by bylo možné ještě navíc prohazovat i sloupce, řáká se takovému postupu částečná pivotace.

Funkce je doplněna o výpis hodnot během výpočtu (print_A(...)), který demonstruje, jak eliminace postupuje.

Na začátku kódu funkce nejprve spojíme matici soustavy a pravou stranu do jediné obdélníkové matice. Nejde o plýtvání časem, stejně bychom si museli pořídit kopie matice i pravé strany, abychom v průběhu eliminace neměnili proměnné, které byly použity při volání.

[ ]:
import numpy as np

def print_A(A, sleep=0.5):
  "Výpis upravované matice a pravé strany"
  # pro přehlednost skrývám import nepodstatných knihoven
  from IPython.display import clear_output
  import time

  with np.printoptions(floatmode="fixed",precision=3):
    print(A)

  time.sleep(sleep)
  clear_output(wait=True)


def GJelim(matice, prava_strana):
  "Vrátí řešení soustavy rovnic A.x=b nalezené GJ eliminací. Parametry musí být typu ndarray."

  dim = len(prava_strana)

  A = np.zeros([dim,dim+1])

  # následující přiřazení provádí zároveň kontrolu velikosti vstupních polí
  A[0:dim, 0:dim] = matice
  A[0:dim,   dim] = prava_strana

  for r in range(dim):  # pro vsechny radky <r>
    print_A(A,sleep=2)

    # 1. částečná pivotace (řeší nulovost diagonálního prvku a zvyšuje přesnost)
    best_r = r
    for r1 in range(r+1,dim):
      if abs(A[r1,r]) > abs(A[best_r,r]):
        best_r = r1
    if best_r != r:
      # A[[best_r,r],r:] = A[[r,best_r],r:]  #alternativní varianta prohození
      print_A(A)
      A[best_r,r:], A[r,r:] = A[r,r:], A[best_r,r:].copy()
      print_A(A)

    # 2. normalizace radku
    # 2.1 najdeme prevracenou hodnotu prvku na diagonale
    faktor = 1 / A[r,r]
    # 2.2  "vsechny" sloupce radku <r> vynasobime timto faktorem
    A[r,r:] *= faktor

    print_A(A)

    # 3. eliminace
    for r1 in range(dim):  # pro vsechny radky <r1>
      if r != r1:
        faktor = A[r1,r];
        A[r1,r:] -= faktor * A[r,r:]
        print_A(A)


  print_A(A,sleep=2)
  return A[:,dim].copy()

matice = np.array( [
        [0, 2, 3, -1],
        [3, 2, 1, 1],
        [1, 9, 0, 4],
        [2, 3, 1, 1] ] )

vektor_ps = np.array( [  -4, 3, -5, 0 ] )

x = GJelim(matice, vektor_ps  )

print( '    řešení =', x )
with np.printoptions(precision=1):
  print( '   chyba x =', x - [1,-2,1,3] )
  print( 'chyba Ax-b =', matice.dot(x) - vektor_ps )
    řešení = [ 1. -2.  1.  3.]
   chyba x = [-5.6e-16 -8.9e-16  8.9e-16  1.8e-15]
chyba Ax-b = [-8.9e-16 -8.9e-16 -8.9e-16 -8.9e-16]

Dva technické detaily se nám dostaly i do tak jednoduchého programu, jako je řešení soustavy rovnic. Oba souvisí s tím že řezy podle numpy nejsou kopie prvků pole a jen jakési chytré odkazy na původní prvky pole.

Právě proto se v příkazu return A[:,dim].copy() vytváří kopie. Jde o to, že A[:,dim] je poslední sloupec obdélníkové matice A. Ta se zdá být lokální proměnnou funkce GJelim a tedy bychom očekávali, že po skončení běhu funkce paměť zabraná A je uvolněna. Pokud bychom ale funkci ukončili příkazem return A[:,dim] pak by matice A nemohla být zapomenuta, protože její poslední sloupeček představuje řešení soustavy rovnic. Teprve až to by nebylo potřeba, mohlo by se uklidit i místo pro matici A.

Při prohazování řádek se opět objevuje použití metody copy(). Přitom u prohazování např. celých čísel jsme mohli psát i,j = j,i a nedošlo k potížím. Řezy numpy jsou ale odkazy na pole a jak je vidět níže, zápis matice[2], matice[3] = matice[3], matice[2] způsobí, že řádek obsahující dvojky je přepsán dříve než se jeho hodnota může použít pro prohození. Proto se musí vytvořit "záloha" prostřednictvím copy().

Cvičení: Odstraňte z kódu prohazování řádků a vyzkoušejte, že pro matici

matice = np.array( [
        [1e-15, 2, 3, -1],
        [3, 2, 1, 1],
        [1, 9, 0, 4],
        [2, 3, 1, 1] ] )

Bude nalezené řešení soustavy rovnic naprosto špatně a že tedy ono prohazová řádků "za lepší" (částečná pivotace) je nezbytné.

Vzhledem k použití celých řezů můžeme eliminaci zapsat jediným příkazem A[r1,r:] -= A[r1,r] * A[r,r:], tedy pomocnou proměnnou faktor lze vynechat. Kdybychom ovšem prováděli tuto operaci cyklem, pak protože se hodnota prvku pole A[r1,r] během takového cyklu změní, je zapamatování původní hodnoty nezbytné.

Samozřejmě, řešení soustavy rovnic je k dispozici jako funkce numpy.linalg.solve.

Časová náročnost algoritmů

Je zřejmé, že některé výpočty mohou trvat dlouho. Uvažme následující obrázek úplného n-úhelníku pro \(n=7\).

n=6

tušíme, že věci se zesložití, pokud vezmeme větší počet vrcholů

n=6

Podobně tomu může být i počítačových programů - při růstu nějakého parametru může doba výpočtu překvapivě prudce. Jak něco takového odhadnout?

Pro každý problém zkusíme najít charakteristické \(N\) číslo udávající jeho velikost. Ne vždy je to možné ale pro představu pár příkladů:

  • Počet položek na skladu v programu pro skladové hospodářství.

  • Počet měst, které má navštívit obchodní cestující.

  • Počet neznámých v soustavě lineárních rovnic.

  • Počet hvězd tvořících hvězdokupu, jejíž vývoj zkoumáme.

  • Počet molekul v simulaci.

Nyní se můžeme ptát jak se bude program zabývající se výše uvedenými problémy chovat při růstu onoho charakteristického čísla \(N\). Samozřejmě, pokud neplánujeme růst skladu, zvyšování počtu neznámých atp., nemusíme se tím zabývat. I ve fyzice se ale projevuje tendence počítat složitější a složitější problémy (třeba ty jednodušší už někdo vyřešil) a tak na problém velkého \(N\) můžeme narazit.

Při porovnávání různých algoritmů máme možnost porovnat, jak se chovají pro různé velikosti vstupních dat přímo. Někdy závisí doba řešení problému na konkrétních vstupních hodnotách, často ale, jak jsme viděli, u soustav lineárních rovnic, závisí pouze na velikosti problému.

Vzpomeneme-li si na definici algoritmu, víme že jde o posloupnost elementárních operací. Předpokládejme nejdříve, že elementárními operacemi rozumíme základní operce jako jsou sčítání, násobení, větvení, přístup k jednoduché proměnné, přístup k prvku pole volání či návrat z podprogramu atp. Změříme-li dobu, kterou potřebuje počítač k provedení každé z těchto operací a budeme-li předpokládat, že celý program se vykoná za dobu odpovídající součtu těchto jednotlivých operací, víme jak spočíst dobu potřebnou k provedení programu. Uvažujme podprogram pro skalární součin dvou vektorů o délce \(N\)

def scalar_product(a,b):
   "Počítá skalární součin dvou vektorů"
   dim = len(a)

   s = 0.0
   for i in range(dim):
       s = s + a[i]*b[i]

   return s

Obsahuje různé operacem přičemž např. dim = len(a) se provádní jednou zatímco s = s + a[i]*b[i] se provádí \(N\times\). Mohli bychom tedy nějak zjistit čas provádení jednotlivých operací, ten vynásobit počtem provádění těchto operací, všechny časy sečíst a získat vzorec udávající závislost doby výpočtu na hodnotě \(N\) nap59klad v podobě

\[T[{\rm ns}] = 542 N + 1220\]

Jde o naivní přístup, na moderních počítačích je složité odhadnout dobu výpočtu i kdybychom znali prováděné instrukce procesoru, další komplikace přináší svojí povahou jazyk Python.

Můžeme také zkusit provést experiment:

[60]:
import time
import numpy as np
import matplotlib.pyplot as plt

def scalar_product(a,b):
    "Počítá skalární součin dvou vektorů"
    dim = len(a)

    s = 0.0
    for i in range(dim):
        s = s + a[i]*b[i]

    return s

def time_it(n):
    x = np.linspace(0,1,n)

    start = time.time()
    s = scalar_product(x,x)
    end = time.time()
    return end - start

n = [2**k for k in range(2, 24)]
t = [time_it(k) for k in n]

plt.plot(n,t)
plt.yscale('log')       # logaritmický průběh os
plt.xscale('log')
plt.xlabel('n')
plt.ylabel('čas')
plt.grid()
plt.show()
../_images/Prednaska_0090_Pole_a_linearni_algebra_29_0.png

Pohled na obrázek ukazuje, že jde o lineární závislot.

Podobně bychom pro násobení matice vektorem dostali třeba

\[T = 36 N^{2} + 44 N +26\]

Pokud nás ale zajímá chování pro větší \(N\), je rozhodující první člen (nejvyšší mocnina). Píšeme tedy

\[T \approx 36 N^{2}\]

Určit koeficient přesně je ale i tak velmi obtížné, nejsnazší možností je analýza měření závislosti času výpočtu na \(N\) (viz např. graf výše). Pokud chceme mluvit o výkonu algoritmu v okamžiku jeho návrhu ještě před jeho kódováním a nákupem počítače, ukazuje se, že nejjednoduuší je prostě psát

\[T = O(N^2)\]

Velké \(O(f)\) je označení pro libovolnou funkci \(g\), která splňuje vztah

\[0 < \lim_{N\rightarrow \infty} \left|\frac{g}{f}\right| <\infty\]

Následující tabulka má ilustrovat, že opravdu rozhodující je právě řád, nikoli konkrétní hodnota koeficientu u vedoucího členu.

hodnoty O(N)

Naše výpočetní možnosti podléhají principiálním omezením, od velkého třesku uplynulo asi \(4\times 10^{26}\) nanosekund a v pozorovaném vesmíru je \(\lessapprox 10^{83}\) atomů.

Vztah \(O(f) . O(g) = O(f . g)\) nám pak umožňuje místo rozkladu algoritmu na elementární kroky O(1), uvažovat strukturovaně.

Třeba Gauss-Jordanova eliminace (za předpokladu, že počet pravých stran je \(\le O(N)\) ) pro každý řádek ( tedy \({\it O(N) }\) krát) provádí

  • hledání největšího prvku na/pod diagonálou ${:nbsphinx-math:`it O(N) `} $

  • normalizace \({\it O(N) }\)

  • odečtení násobku řádku od všech ostatních ${:nbsphinx-math:`it O(N ^{2}) `} $

  • na začátku a na konci ještě nějaké kopírování \({\it O(N ^{2}) }\)

Odsud máme

\[O(N^{2}) + O(N)*[O(N)+ O(N)+ O(N^{2})] = O(N ^{3})\]

V tabulce si pak snadno najdeme, pro jaká \(N\) jde o zelený, oranžový či hnědý problém.

Podobně jako se s rostoucím \(N\) nějak chová čas potřebný k provedení výpočtu, nějak rostou i paměťové nároky. Protože nad velikostí datových struktur máme trochu lepší kontrolu, lze v praxi velmi dobře odhadnout, co se vejde a co ne. V rozvahách o schůdnosti různých algoritmů ale také můžeme použít Landauovu notaci velkého \(O\).