In [95]:
# Import potrebných balíčkov
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy.utilities.lambdify import lambdify
import sympy as sp

Metóda klesajúceho gradientu na jednoduchom paraboloide

Ako prvý krok si vyskúšajme metódu klesajúceho gradientu aplikovať na minimalizáciu jednoduchej funkcie – paraboloidu podľa:

\begin{equation} z = f(x, y) = x^2 + y^2 \end{equation}

Vizualizácia funkcie

Na začiatok si funkciu zadefinujme a vizualizujme:

In [96]:
def f(x, y):
    return x**2 + y**2

Vygenerujeme si všetky kombinácie bodov $x, y$ z určitého rozsahu a vypočítame pre ne $z = f(x, y)$:

In [97]:
xx, yy = np.mgrid[-10:10.2:0.2, -10:10.2:0.2]
zz = f(xx, yy)

Výsledok vizualizujeme ako 3D graf:

In [98]:
plt.figure()
ax = plt.gca(projection='3d')
ax.plot_surface(xx, yy, zz, cmap='Spectral',
                linewidth=0, antialiased=True)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
plt.savefig("gradient_3d_plot.pdf")

V 3D grafoch sa budú niektoré veci navzájom prekrývať – aby boli naše grafy ľahšie čitateľné, budeme nižšie namiesto 3D grafov používať 2D vrstevnicové grafy:

In [99]:
plt.figure()
plt.contour(xx, yy, zz, cmap='Spectral')
# obe osi v rovnakej škále + legenda
plt.gca().set_aspect('equal')
plt.xlabel('x'); plt.ylabel('y')
plt.colorbar(label='z')
plt.savefig("gradient_func_contour.pdf",
            bbox_inches="tight",
            pad_inches=0)

Obaľme si tento kód do pomocnej funkcie, aby sme ho nižšie nemuseli zakaždým opakovať:

In [100]:
def plot_func(xx, yy, f, X=None):
    if not X is None:
        Xmin, Xmax = X[:, 0].min(), X[:, 0].max() 
        Ymin, Ymax = X[:, 1].min(), X[:, 1].max()
        
        if (Xmin < xx.min() or Xmax > xx.max() or
                Ymin < yy.min() or Ymax > yy.max()):            
            xx = np.linspace(Xmin, Xmax, 100)
            yy = np.linspace(Ymin, Ymax, 100)
            xx, yy = np.meshgrid(xx, yy)
            
        plt.scatter(X[:, 0], X[:, 1], zorder=10)
        
    zz = f(xx, yy)
    plt.contour(xx, yy, zz, cmap='Spectral')
    # obe osi v rovnakej škále + legenda
    plt.gca().set_aspect('equal')
    plt.xlabel('x'); plt.ylabel('y')
    plt.colorbar(label='z')

Gradient funkcie

Aby sme mohli funkciu $f(x, y)$ minimalizovať pomocou metódy klesajúceho gradientu, musíme určite jej gradient. Pripomeňme si, že gradient $\nabla f(x, y)$ funkcie $f(x, y)$ je vektor jej prvých parciálnych derivácií. V našom prípade teda:

\begin{equation} \nabla f(x, y) = \left( \frac{\partial f}{\partial x}, \frac{\partial f}{\partial y} \right) \end{equation}

Pripomeňme, že naša funkcia má tvar $f(x, y) = x^2 + y^2$. Určiť parciálne derivácie je teda ľahké. Pri parciálnej derivácii podľa $x$ sa bude člen $y^2$ brať ako konštanta a derivovať budeme len $x^2$. Pri derivácii podľa $y$ to bude naopak. Máme teda:

\begin{align} \frac{\partial f}{\partial x} &= 2x \\[0.75em] \frac{\partial f}{\partial y} &= 2y \end{align}

Úloha 1: Výpočet gradientu

Doplňte teraz do nasledujúcej bunky kód tak, aby funkcia grad_f navrátila gradient funkcie $f(x, y)$:

In [101]:
def grad_f(x, y):
    # sem treba doplniť kód
    return np.array([      ,      ])

Vizualizácia gradientu

Ako vieme, gradient funkcie indikuje smer jej najväčšieho rastu. Aby sme si vedeli predstaviť, čo to znamená, môžeme si vizualizovať gradient funkcie, ktorý sme si práve zadefinovali:

In [102]:
xxg, yyg = np.mgrid[-10:11:1.5, -10:11:1.5]
gg = np.array(
    [[grad_f(x, y) for x, y in zip(rx, ry)] 
          for rx, ry in zip(xxg, yyg)]
)
In [103]:
plt.figure(figsize=[8,8])
plot_func(xx, yy, f)
plt.quiver(xxg, yyg, gg[..., 0], gg[..., 1],
   #np.sqrt(gg[..., 0]**2 + gg[..., 1]**2), cmap='jet'
)
plt.savefig("gradient_func_quiver.pdf",
            bbox_inches="tight",
            pad_inches=0)

Šípky ukazujú smer gradientu. Ako vidno, všetky ukazujú von zo stredu – kde paraboloid stúpa. Veľkosť šípok indikuje veľkosť gradientu. Šípky blízko stredu sú maličké (derivácia funkcie v minime je nulová) a smerom ku okrajom sa zväčšujú: pretože hodnota funkcie rastie čoraz strmšie.

Metóda klesajúceho gradientu

Ďalej budeme pokračovať aplikáciou metódy klesajúceho gradientu, aby sme funkciu minimalizovali. Na začiatok si definujeme zopár parametrov: počet krokov a rýchlosť učenia:

In [104]:
num_steps = 20
learning_rate = 0.1

Aby sme si mohli neskôr postup minimalizácie vizualizovať, uložíme si všetky vypočítané body do matice. Matica má 2 stĺpce – jeden pre $x$ a druhý pre $y$:

In [105]:
X = np.zeros((num_steps + 1, 2))

Počiatočný bod zvolíme buď náhodne:

In [106]:
X[0] = np.random.uniform(-10, 10, (2,))

alebo si môžeme vybrať nejaký fixný bod, aby sme zakaždým dostali ten istý výsledok:

In [107]:
X[0] = [-9, -8]

Metódu klesajúceho gradientu budeme aplikovať tak, že začneme z počiatočného bodu $\mathbf{x}_0$ a budeme iteratívne aplikovať známe pravidlo, pomocou ktorého vypočítame vždy nasledujúci bod: \begin{equation} \mathbf{x}_{i+1} = \mathbf{x}_i - \gamma \nabla f(\mathbf{x}_i) \end{equation} kde $\gamma$ je rýchlosť učenia a $\nabla f(\mathbf{x}_i)$ je gradient minimalizovanej funkcie.

Pravidlo sme zapísali vo vektorovom tvare, takže $\mathbf{x}$ označuje vektor – v našom prípade $\mathbf{x} = (x, y)$.

Úloha 2: Implementácia metódy klesajúceho gradientu

Doplňte do nasledujúcej bunky chýbajúcu časť kódu:

In [108]:
for i in range(num_steps):
    
    # sem treba doplniť kód:
    
    X[i+1] =

Vypočítané body si nakoniec vizualizujeme a skontrolujeme, či postupnosť naozaj smeruje ku minimu paraboloidu:

In [120]:
plot_func(xx, yy, f, X)
plt.savefig("gradient_mini_steps.pdf",
            bbox_inches="tight",
            pad_inches=0)

Úloha 3: Obalenie kódu do funkcie

Kód si teraz znovu obaľme do funkcie, aby sme ho mohli opätovne použiť:

In [110]:
def grad_desc(grad_func, init_point,
              num_steps, learning_rate):
    X = np.zeros((num_steps + 1, 2))
    X[0] = init_point
    
    for i in range(num_steps):
        
        # sem treba doplniť kód:
        
        X[i+1] = 
        
        
    
    return X

Novú funkciu môžeme otestovať takto:

In [111]:
X = grad_desc(grad_f, init_point=[-9, -8],
              num_steps=20, learning_rate=0.1)
plot_func(xx, yy, f, X)

Testovanie iných rýchlostí učenia

Pre ilustráciu si teraz metódu klesajúceho gradientu otestujeme aj pri iných rýchlostiach učenia.

Začnime s $\gamma = 0.45$:

In [112]:
X = grad_desc(grad_f, init_point=[-9, -8], num_steps=20,
      learning_rate=0.45
)
plot_func(xx, yy, f, X)
plt.savefig("gradient_lr_0_45.pdf",
            bbox_inches="tight",
            pad_inches=0)

Pre $\gamma = 1.0$ bude už algoritmus oscilovať bez toho, aby sa ku minimu priblížil:

In [113]:
X = grad_desc(grad_f, init_point=[-9, -8], num_steps=20,
      learning_rate=1.0
)
plot_func(xx, yy, f, X)
plt.savefig("gradient_lr_1_0.pdf",
            bbox_inches="tight",
            pad_inches=0)

Pre $\gamma > 1$ dôjde ku divergencii a algoritmus sa bude od minima naopak vzďaľovať:

In [114]:
X = grad_desc(grad_f, init_point=[-9, -8], num_steps=20,
      learning_rate=1.02
)
plot_func(xx, yy, f, X)
plt.savefig("gradient_lr_1_02.pdf",
            bbox_inches="tight",
            pad_inches=0)

-19.720108287300768 18.961642583943046 -17.528985144267384 16.854793407949405

Automatický výpočet symbolického gradientu

Vyššie sme gradient funkcie vypočítali analyticky a ručne prepísali do kódu. Symbolický výpočet gradientu je však v jazyku Python možné vykonať aj automaticky – pomocou balíčka sympy. V ďalšom príklade si teda ukážeme, ako sa to dá urobiť.

Najprv si definujeme symbolické premenné, ktoré budeme potrebovať – v našom prípade $x$ a $y$:

In [115]:
symx, symy = sp.symbols('x y')

Analyticky definujeme funkciu $f(x,y)$:

In [116]:
symf = symx**2 + symy ** 2
symf
Out[116]:
$\displaystyle x^{2} + y^{2}$

Pri výpočte symbolického gradientu použijeme drobný trik. Zo skalárnej funkcie si najprv spravíme maticu a potom vypočítame jej jakobián. Výsledkom bude riadkový vektor zodpovedajúci gradientu:

In [117]:
symg = sp.Matrix([symf]).jacobian([symx, symy])
symg
Out[117]:
$\displaystyle \left[\begin{matrix}2 x & 2 y\end{matrix}\right]$

Aby sme výslednú symbolickú reprezentáciu gradientu mohli použiť na výpočet pre konkrétne hodnoty, ešte ju konvertujeme na bežnú numerickú funkciu. To isté spravíme pre funkciu $f$:

In [118]:
f = lambdify((symx, symy), symf, "numpy")
grad_f = lambdify((symx, symy), symg, "numpy")

Metódu klesajúceho gradientu môžeme teraz aplikovať rovnako, ako v predchádzajúcom prípade – akurát sme nemuseli ručne vypočítať gradient:

In [119]:
xx, yy = np.mgrid[-10:10.2:0.2, -10:10.2:0.2]
zz = f(xx, yy)

X = grad_desc(grad_f, init_point=[-9, -8],
              num_steps=20, learning_rate=0.1)
plot_func(xx, yy, f, X)
In [ ]: