machine learning

Broyden–Fletcher–Goldfarb–Shanno – AI Tutorial 4

Poprzedni post z cyku.

Miało być o racjonalnych decyzjach ale podjąłem takową dzisiaj i będzie o czymś innym.

Wróćmy jeszcze do uczenia maszynowego i wprowadźmy kilka usprawnień.
Na początek zobaczmy dane na jakich będziemy pracować. (X – 0, gwiazdka – 1).

raw

Pierwszą zmianą będzie użycie innej metody minimalizowania. Gradient prosty jest prosty i tu w zasadzie kończą się jego zalety. Jednakże zamiast samemu implementować inny mechanizm skorzystajmy z tego co dostarcza biblioteka do obliczeń symbolicznych SymPy.
Wybór padł na algorytm Broyden–Fletcher–Goldfarb–Shanno. Nie będę go opisywać gdyż nie jest to blog o analizie numerycznej (jeszcze nie ( ͡° ͜ʖ ͡°) ).
Przegląd innych metod minimalizacji jest dostępny w dokumentacji.
Użycie jest bardzo proste, wystarczy zaimportować odpowiednią funkcję.

from scipy.optimize import fmin_bfgs

A następnie wywołać ją, przekazując jej funkcję do zoptymalizowania oraz miejsce od którego zacząć:

theta = fmin_bfgs(J, initialTheta)

Znajomość pochodnych nie jest konieczna, wystarczy nam sama funkcja J

#ostatecznie uznałem, że tak będzie czytelniej
def single_penalty(Theta, x, y):
    return -np.log(h(Theta, x)) if y == 1 else -np.log(1-h(Theta, x))                
        
def J(Theta):
    sum = 0.
    for i in range(m):
        sum += single_penalty(Theta, TrainingData[i,:], Y[i])
    return sum/m

Funkcja fmin_bfgs wypisze nam jak przebiegła optymalizacja oraz zwróci nam wynik (poszukiwane przez nas \theta )Np:

Optimization terminated successfully.
Current function value: 0.000002
Iterations: 30
Function evaluations: 487
Gradient evaluations: 37

Takie podejście do sprawy daje nam komfort swobodnych eksperymentów z funkcjami g_\theta(x) i J(\theta). Zacznijmy od najprostszej wersji hipotezy:
g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2

Możemy się spodziewać, że dla testowych danych trudno będzie poprowadzić sensownie prostą linię. I tak jest w istocie.

1

Wygląda źle, ale nie oczekiwaliśmy niczego innego. Dodajmy jeden wyraz do hipotezy.

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2

2

Eeee nie. Coś się zaczęło co prawda „wyginać” ale daleko nam do celu. Kolejne dwa wyrazy:

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2

3Znacznie lepiej. 3 kolejne okręgi (od środka) oznaczają kolejno P(y=1) = 1%, P(y=1) = 50%, P(y=1) = 99%. Im będą bliżej siebie, tym bardziej wyraźna granica decyzji i mniej obszarów niepewności. Czy dodanie kolejnego wyrazu może coś poprawić? Sprawdźmy

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2 + \theta_6 x_1^2x_2^2

4

Kształt zmienił się zdecydowanie, obszar niepewności również. Czy jest to lepsze dopasowanie niż poprzednio to już kwestia dyskusyjna. Faktem jest jednak że algorytm wytrenował się poprawnie do podanych danych. Dodajmy jeszcze coś.

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2 + \theta_6 x_1^2x_2^2 + \theta_7 x_1^2 x_2 + \theta_8 x_1 x_2^2

5

Sytuacja jest podobna. Ostatni krzyk rozpaczy, dodajemy trzecie potęgi.

g_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_1 x_2 + \theta_4 x_1^2 + \theta_5 x_2^2 + \theta_6 x_1^2x_2^2 + \theta_7 x_1^2 x_2 + \theta_8 x_1 x_2^2 + \theta_9 x_1^3 + \theta_{10} x_2^3

6
Wystarczy. Dodawanie kolejnych potęg nic nie wnosi, a wręcz może zaszkodzić. W przyszłości poznamy mechanizm zniechęcania do wykorzystywania zbyt skomplikowanych krzywych.

Podsumowanie

Teraz już na pewno odstawiamy na bok (tymczasowo) uczenie maszynowe i przechodzimy do czegoś innego. O tym jak znaleźć drogę w labiryncie już niedługo.

Całość kodu

Regresja logistyczna (logistic regression) – AI tutorial 3

Pierwszy post z cyklu
Drugi post z cyklu

Zostawmy na jakiś czas egzaminy, przejdźmy do innych problemów.
Idąc do lekarza mówimy co nas boli, robimy badania (jeśli NFZ przypadkiem jeszcze nie zbankrutował) a wyniku pracy inteligentnej osoby dostajemy odpowiedź czy jesteśmy chorzy. Załóżmy, że chcemy zastąpić lekarzy maszynami. Doktor na podstawie zdobytej wiedzy oraz doświadczenia (czyli de facto również wiedzy) podejmuje decyzję, czy jesteśmy zdrowi czy chorzy. Jak to się ma do poprzedniego problemu? Trochę podobnie, trochę inaczej. Na wejściu jest to samo – n zmiennych wejściowych (po ang. zwane features). Wyjście jest inne – zamiast dostać wartość z ciągłego zbioru (punkty) dostajemy binarną odpowiedź chory/zdrowy. Mamy tutaj do czynienia z problemem klasyfikacji – efekt pracy może należeć tylko do pewnego skończonego zbioru. Wymaga to nieco innego podejścia. Inne przykłady takich problemów:

X oznacza wektor zmiennych wejściowych

  • Email klasyfikowany jako spam X → {0, 1}
  • Email klasyfikowany do pewnej kategorii (praca, prywatne, newslettery) X → {0, 1, 2}
  • Kategoria aut (sedan, pickup etc) na podstawie zdjęcia X → {0, 1, 2, 3, 4…}
  • Klasyfikacja płci X → {0, 1}

Testowany obiekt musi być w dokładnie jednej kategorii.

Na razie ograniczmy się do odpowiedzi tak/nie: y \in \{0,1\}. Postawmy taki problem: mamy chorego na określoną chorobę i chcemy określić czy może zostać wyleczony czy też nie. Do dyspozycji mamy takie dane jak jego wiek, waga, wyniki badań, płeć etc. Na początku załóżmy że rozważamy tylko wiek. Przykładowe dane mogą wyglądać następująco:

pic1

Kółko oznacza 0, gwiazdka 1 – będę dalej używać tej notacji.

Widzimy, że tendencja jest bardzo prosta – im pacjent starszy tym gorsze rokowania. Gdybyśmy chcieli wyznaczyć przedziały dla y=0 (wyzdrowienie) oraz y=1 (śmierć) mogło by wyglądać to tak:

pic2

Widać tutaj pewien problem – w okolicach miejsca podziału pewnie obiekty zostaną źle sklasyfikowane, nie da się tego całkowicie uniknąć. Co możemy z tym zrobić? Zamiast mówić „Hej, umrzesz” lepiej będzie powiedzieć „Umrzesz z prawdopodobieństwem 87%”. Oczywistą zaletą jest, że w niepewnej okolicy będziemy dostawać wyniki oscylujące wokół 50% co wskaże fakt niepewności klasyfikacji. Rzućmy okiem na przykład z dwiema zmiennymi wejściowymi:

db1

Jedna zmienna jest na osi X, druga na Y. Zaznaczona linia wskazuje miejsce w którym tniemy płaszczyznę (czy to zawsze i wszędzie muszą być linie proste? Nie, pokażemy wkrótce inne warianty).

Hipoteza h(x)

Ostatnio szukaliśmy takiej funkcji h która dla pewnego zestawu danych zwróci nam przewidywany wynik egzaminu. Teraz szukamy takiej która zwróci nam prawdopodobieństwo tego, że y = 1. W związku z tym funkcja musi być ograniczona 0 \le h_\theta(x) \le 1 . Konstrukcja takiej funkcji może być bardzo problematyczna.
Prościej było by uznać że jeśli h(x) > 0 – wynik pozytywny (y = 1). W takim rozumieniu tracimy jednak informację o prawdopodobieństwie chociaż intuicyjnie czujemy że im większa wartość tym większa pewność wyniku pozytywnego, a im mniejsza negatywnego. Potrzebujemy więc pewnej funkcji (nazwijmy ją S) która zamieni wygodną dla nas reprezentację na tą bardziej użyteczną.

S: (-\infty; +\infty) \mapsto (0;1)

Taka funkcja istnieje, nazywa się ona sigmoid function (tłumaczenia się nie podejmę) i jest zdefiniowana następująco:

S(t) = \frac{1}{1+e^{-t}}

Spójrzmy na wykres

 

Nada się znakomicie – dla t = 0 zwraca dokładnie 0,5. im dalej w kierunku nieskończoności tym bliżej 0 lub 1.

Uporządkujmy nieco oznaczenia. Funkcja hipotezy h_\theta(x) będzie teraz korzystać z funkcji S aby „przetłumaczyć” wynik z formy wygodnej dla nas na tą użyteczną. Do kreślenia prostych używaliśmy ostatnio h_\theta(\boldsymbol{X}) = \boldsymbol{\theta}^T \boldsymbol{X} . Teraz też będzie odpowiednia. Oznaczmy ją tym razem jako g(X):
g(\boldsymbol {X}) = \boldsymbol{\theta}^T \boldsymbol{X}
pamiętając że mamy do dyspozycji funkcję S:
S(t) = \frac{1}{1+e^{-t}}

Otrzymujemy ostatecznie:

h_\theta(\boldsymbol{X}) = S(g(\boldsymbol{X}))
lub równoważnie:
h_\theta(\boldsymbol{X}) = \frac{1}{1+e^{-\boldsymbol{\theta}^T \boldsymbol{X}}}

Kiedy będziemy chcieli zmienić zachowanie funkcji zmodyfikujemy tylko g(X)

Właśnie, jak mamy interpretować wartość zwracane przez g? Przy regresji liniowej była to wartość oczekiwana dla danego wejścia (co zresztą było bardzo dobrze widać na wykresie). Do tej pory ustaliliśmy: funkcja h zwraca prawdopodobieństwo wyniku pozytywnego (y = 1). Jeżeli jest ono większe niż 50% klasyfikujemy wynik jako pozytywny, w przeciwnym razie jako negatywny. Progiem „przełączania” jest wartość 0.5 funkcji h czyli wartość 0 funkcji g. Wynika z tego, że równanie g(x) = 0 pokazuje linię (powierzchnię, hiperpłaszczyznę etc) po której jednej stronie wartości są klasyfikowane jako pozytywne, a po drugiej – negatywne. Spójrzmy na przykłady

db2

 

Dotychczas używana przez nas metoda (konkretnie funkcja g) pozwalała tylko na oddzielanie elementów linią prostą. Nic jednak nie stoi na przeszkodzie aby była nieco bardziej skomplikowana, np:

g(x) = \theta_0 x_0 +\theta_1 x_1+\theta_2 x_2+ \theta_3 x_1^2 +\theta_4 x_2^2 +\theta_5 x_1 x_2 + \dots

Spowoduje to dopasowanie pewną krzywą, im dłuższe wyrażenie tym potencjalnie bardziej skomplikowaną. Rozwiniemy temat przy okazji sieci neuronowych (przykład na coś innego niż prosta pojawi się jeszcze w tym poście).

Funkcja kosztu J

Funkcję hipotezy mamy za sobą, co z funkcją kosztu J? Poprzednio miała formę:

J(\boldsymbol{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m \left ( h_\theta(x^{(i)}) - y^{(i)} \right )^2

Wyizolujmy z niej część odpowiedzialną za obliczenie kosztu pojedynczego zestawu danych. Będzie to:
h_\theta(x^{(i)}) - y^{(i)}

Zastanówmy się jak powinna wyglądać obecnie. Z całą pewnością chcemy, aby poprawny wynik był obarczony zerowym kosztem. Z drugiej strony jeżeli nasza hipoteza stwierdzi z dużą pewnością nieprawidłowy wynik, penalizacja powinna być wysoka. Pamiętajmy że poruszamy się w przedziale (0;1).

Na początek załóżmy, że y = 0. Jak ma kształtować się kara dla różnych wyników h(x)? Na przykład tak:

j0

Dla prawidłowego rezultatu kara jest równa zero. Im bardziej algorytm był pewny złej odpowiedzi tym większą karę ponosi (y dąży do nieskończoności wraz z x dążącym do 1).

Dla y = 1 sytuacja jest odwrotna:

j1

Wyrażenie opisujące pojedynczy zestaw danych będzie więc wyglądać tak:

\text{single penalty}(\boldsymbol{X}) =  \begin{cases}  -\log(h_\theta(\boldsymbol{X})) & \text{if } y = 1 \\  -\log(1-h_\theta(\boldsymbol{X})) & \text{if } y = 0  \end{cases}

Możemy pozbyć się niepożądanego warunku zapisując całość jako: (ze względu na y \in \{0,1\})

\text{single penalty}(\boldsymbol{X}) =-y \log(h_\theta(\boldsymbol{X})) - (1-y) \log(1-h_\theta(\boldsymbol{X}))

Podstawiając to do poprzedniego wyrażenia na wartość funkcji J(\boldsymbol{\theta}) otrzymujemy ostatecznie:

J(\boldsymbol{\theta}) = -\frac{1}{m} \left( \sum\limits_{i=1}^m \left (  y \log(h_\theta(\boldsymbol{X})) - (1-y) \log(1-h_\theta(\boldsymbol{X})    \right ) \right )

Optymalizacja

Tradycyjnie, kiedy mamy już opisną funkcję J problem jest sprowadzony do jej minimalizacji. Wyznaczmy pochodne kierunkowe względem wektora \boldsymbol{\theta}.

\frac{\partial}{\partial \theta_j} J(\boldsymbol{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_j^{(i)}

???

Nie, to nie jest błąd. Czy dostaliśmy to samo (z dokładnością do mnożenia przez 2)? Prawie. Pamiętajmy że obecne h_\theta(x) wygląda inaczej niż w regresji liniowej.

Całość algorytmu gradientu prostego jest identyczna, więc nie będę go opisywał ponownie.

Implementacja

Zaimplementujmy wersję w której funkcja g(x) będzie prezentować się następująco:

g(x) = \theta_0 x_0 + \theta_1 x_1^2 + \theta_2 x_2^2

Jest to nic innego jak elipsa o środku w O=(0,0).

Funkcje S oraz h w kodzie Pythona:

def sigmoid(x):
    return 1/(1 + np.e**(-x))

def h(Theta, X):
    z = Theta[0] * X[0] + Theta[1] * X[1] **2 + Theta[2] * X[2]**2
    return sigmoid(z)

Obliczenia pochodnych cząstkowych są identyczne jak poprzednio podobnie jak i sam gradient. Spójrzmy na przykładowy wynik działania. Mając dostępne następujące dane:

raw

Algorytm wytyczył krzywą dla której prawdopodobieństwo obu wyników jest takie samo (P(y=1)=0,5 w sposób następujący:

line

Dodatkowo krzywe dla 75% prawdopodobieństwa y=1 oraz y=0.

imgNastępnym razem zostawimy uczenie się na pewien czas i zastanowimy się nad podejmowaniem rozsądnych decyzji.

Całość kodu

Regresja liniowa wielu zmiennych (AI tutorial 2)

Poprzedni post z cyklu

Założenie, że wynik egzaminu zależy tylko od obecności na wykładach jest bardzo naiwne. Bez problemu możemy wskazać inne czynniki – ocena z ćwiczeń, IQ, nastrój prowadzącego (im lepszy tym łagodniej ocenia), ilość godzin spędzonych na nauce, etc. Uogólnijmy więc metodę regresji liniowej na wiele zmiennych. Przykłady będę prowadził na dwóch, ale podane algorytmy działają dla dowolnej liczby (w granicach mocy obliczeniowej).

Powróćmy do przykładu z egzaminem. Tym razem do dyspozycji mamy również ocenę z laboratoriów. Dostępne dane treningowe prezentują się następująco:

training

Hipoteza h(x)

Trend wydaje się następujący: im więcej nieobecności tym gorszy wynik, im lepsza ocena tym więcej punktów.
Metoda postępowania jest bardzo podobna do poprzedniej, właściwie jedyną rzeczą którą będziemy robić będzie uogólnianie. Na początek przyjrzyjmy się funkcji h_\theta(x) . Poprzednio miała postać:

h_\theta(x) = \theta_0 + \theta_1 x

Oznaczmy opuszczone wykłady przez x_1 a ocenę z laborek przez x_2. W ogólnym przypadku x_k oznacza k-tą zmienną wejściową. Funkcja będzie wyglądać następująco:

h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2

W przypadku ogólnym (dla n zmiennych) będzie to:

h_\theta(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

Spróbójmy uprościć nieco zapis (nie mam tu na myśli wstawienia znaku \sum ). Zmiennych wejściowych jest wiele, zapiszmy je jako wektor.

\boldsymbol{X} = \left [ x_1, x_2 \dots x_n \right ] n-elementowy wektor wierszowy \boldsymbol{X} \in \mathbb{R}^n

Parametrów θ jest również wiele – zróbmy z nich wektor \boldsymbol{\theta} \in \mathbb{R}^{n+1}. n + 1 ponieważ mamy \theta_0 nie powiązane z żadnym x.

\boldsymbol{\theta} = \left [ \theta_0, \theta_1 \dots \theta_n \right ]

Na razie nie pomogło nam to zbyt mocno, ilosć elementów jest różna. Wystarczy jeden prosty trick aby ją unormować. Dodajmy „fałszywą” zmienną x_0 = 1 i przypiszmy ją do pary razem z \theta_0

h_\theta(x) = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n

Mnożenie przez jeden oczywiście nic nie zmieni, zostaje to co było ale w wygodniejszej formie. X wygląda teraz tak:

\boldsymbol{X} = \left [ x_0 ,x_1, x_2 \dots x_n \right ], \quad \boldsymbol{X} \in \mathbb{R}^{n+1}, x_0 = 1

Dochodzimy do konkluzji. Przypominając sobie coś algebry całą funkcje zapisujemy w postaci:

h_\theta(\boldsymbol{X}) = \boldsymbol{\theta}^T \boldsymbol{X}

Zapis matematyczny jest zwięzły, ale czy przekłada się to na kod? Tak (zakładając że korzystamy z sensownego środowiska). Używam Pythona w wersji 3.4 oraz NumPy, funkcja h prezentuje się następująco:

def h(Theta, X):
    return Theta.transpose().dot(X)

Gdzie Theta oraz X są wektorami wierszowymi. transpose() zwraca transpozycję macierzy (pamiętajmy, że wektor to szczególny przypadek macierzy gdzie jeden z wymiarów jest równy 1). dot() dokonuje mnożenia macierzowego.

Funkcja kosztu J

Nie wiem czy funkcja kosztu to poprawne tłumaczenie ale uważam je za całkiem zręczne – koszt podobnie jak funkcję J chcemy zminimalizować.

Co z funkcją J? Nic. Mamy jedynie kosmetyczną zmianę, zamiast dwóch parametrów \theta_0 oraz \theta_1 mamy teraz jeden – wektor \boldsymbol{\theta}.

J(\boldsymbol{\theta}) = \frac{1}{m} \sum\limits_{i=1}^m \left ( h_\theta(x^{(i)}) - y^{(i)} \right )^2

Również pochodne cząstkowe są praktycznie identyczne. Kontynuując nasz przykład dla dwóch zmiennych mamy trzy pochodne kierunkowe (pamiętajmy że \boldsymbol{\theta} to wektor (n+1)-elementowy gdzie n to ilość zmiennych wejściowych).

\frac{\partial}{\partial \theta_0} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right)
\frac{\partial}{\partial \theta_1} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_1^{(i)}
\frac{\partial}{\partial \theta_2} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_2^{(i)}

Pojawiła się nowa notacja x_j^{(i)} Oznacza ona wartość zmiennej x_j w i-tym zestawie danych. Spójrzmy na przykład. Weźmy i = 8, j = 2. Na podstawie tabelki z początku postu odczytujemy że jest to 4,5 (wiersz 8, druga zmienna, dla przypomnienia wartości w kolumnie „punkty na egzaminie” to nasze y).

Pochodne po \theta_1 i \theta_2 możemy swobodnie uogólnić.
Czy to samo możemy zrobić z pochodną po \theta_0? Tak, x_0 jest przecież zawsze równe jeden. Ostatecznie wszystkie pochodne wyznaczamy ze wzoru:
\frac{\partial}{\partial \theta_j} J(\boldsymbol{\theta}) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x_j^{(i)}

Kod funkcji obliczajacej wartość pochodnej funkcji J po j-tej zmiennej dla danego \boldsymbol{\theta}

def Jpartial(j,Theta):
    sum = 0
    for i in range(m):
        sum += (h(Theta, TrainingData[i]) - Y[i]) * TrainingData[i][j]
    return 2 * sum / m

Pętlą for przechodzimy po wszystkich zestawach danych. TrainingData jest macierzą z danymi treningowymi. Dla rozważanych danych (tabela) wygląda ona tak:

TrainingData = np.array([     
	[1.,1.	,4.5],[1.,1.	,3. ],[1.,1.	,4. ],
	[1.,2.	,3.5],[1.,2.	,5. ],[1.,2.	,3. ],
	[1.,3.	,3. ],[1.,3.	,4.5],[1.,3.	,4. ],
	[1.,4.	,3. ],[1.,4.	,3.5],[1.,4.	,5. ],
	[1.,5.	,4. ],[1.,5.	,3. ],[1.,5.	,3.5],
	[1.,6.	,4.5],[1.,6.	,4.5],[1.,6.	,4.5],
	[1.,7.	,4. ],[1.,7.	,3. ],[1.,7.	,5. ],
	[1.,8.	,3. ],[1.,8.	,4. ],[1.,8.	,3.5],
	[1.,9.	,3. ],[1.,9.	,4. ],[1.,10.	,4.5],
	[1.,10.	,5. ],[1.,10.	,3. ]])
m,n = TrainingData.shape

Jest to macierz (w terminologii NumPy tablica) o m wierszach i n + 1 kolumnach. Kolumna 0 to oczywiście sztuczne x_0.
Zapis

TrainingData[i][j]

Oznacza wybranie wartości z i-tego wiersza i j-tej kolumny.

TrainingData[i]

to wybór całego i-tego wiersza.
Y to wektor z wynikami (trzecia kolumna w tabelce).

Metoda gradientu prostego dla wielu zmiennych

Już poprzednio szukanie minimum wymagało modyfikowania dwóch zmiennych. Dwie pętle for i jesteśmy w domu. Zanim to zrobimy zastanówmy się nad warunkiem stopu. Poprzednio podałem bardzo ogólne „skończ kiedy uznasz za stosowne”. Teraz zastanówmy się kiedy stosownie będzie zakończyć.

 

Analiza matematyczna mówi nam, że w optimum lokalnym wartości pochodnych są równe zero. Szansa na to, że dostaniemy *dokładnie* zero jest niewielka (a właściwie to zerowa – im bliżej jesteśmy tym wolniej idziemy, do zera zbliżamy się asymptotycznie). Pamiętajmy że obliczenia numeryczne są mniej lub bardziej przybliżone. Wybierzmy zatem pewną liczbę ε poniżej której uznamy, że funkcja się zbiegła. Jeżeli wszystkie pochodne są mniejsze od niej stwierdzamy osiągnięcie minimum i kończymy.

converge

Zbieżność do zera. X – iteracje, Y = J(θ)

notcoverage

Rozbieżność. Wartość α jest zbyt duża. X – iteracje, Y = J(θ)

 

def gradient_descent(startTheta):
    theta = startTheta
    converge = False
    i = 0
    while not converge:
        differentials = np.zeros(n)
        for j in range(n):
            differentials[j] = Jpartial(j, theta)
        for j in range(n):
            theta[j] -= alpha * differentials[j]
			
        if i % 500 == 0:
            print (np.amax(differentials))
        converge = np.amax(differentials) < epsylon
    return theta

startTheta to punkt początkowy. Dopóki nie stwierdzimy zbieżności obliczamy wartości wszystkich pochodnych. np.zeros(n) tworzy tablicę w której będziemy zapisywać ich wartości. Po ich obliczeniu przesuwamy wszystkie θ. Jeśli wszystkie są mniejsze od ε kończymy.
Dodatkowo co 500 iteracji wypisujemy informację która pozwoli nam ocenić jak długo musimy jeszcze czekać.

Wyniki

Zobaczmy jakich wartości θ wyuczył się algorytm.
\theta_0 = -0,173
\theta_1 = -0,396
\theta_2 = 1,935

Opuszczone wykłady są skorelowane z modyfikatorem -0,4. Im więcej opuszczonych tym mniej punktów. Ocena z laborek ma mnożnik ~2. Widać że komputer doszedł do „rozsądnych” wniosków (na podstawie zadanych wartości treningowych). Zauważyć można również fakt, że ocena jest ~4x bardziej istotna przy przewidywaniu wyniku. Zapytajmy teraz o kilku studentów i sprawdźmy co na temat ich egzaminu powie algorytm.

results

Pamiętajmy, że jakość otrzymywanych wyników zależy w kluczowym stopniu od podanego zestawu treningowego!

Całość kodu

It’s been a long time (AI tutorial 1)

Sporo czasu upłynęło od ostatniego razu. Po raz kolejny wznawiam bloga, tym razem mając konkretne cele. Blog ma mi towarzyszyć w moich zmaganiach na drodze do zostania wymiataczem w AI. Temat bardzo ciekawy, bardzo trudny, bardzo nęcący ale zawsze odstawiany przeze mnie na bok ze względu na różne wymówki. Wystarczy, w końcu coś pęka i człowiek zaczyna robić to co powinien zacząć już dawno. Każda kolejna rzecz opisana na blogu w celu raz usystematyzowania wiedzy, dwa posiadania notatek.

Na początek odświeżenie tego co już wiadomo.

Wyobraźmy sobie taką sytuację. Jest egzamin, jedni zdają inni nie zdają, wszystko w normie. Każdy jednak woli być w tym zbiorze studentów który zdał. A najlepiej było by wiedzieć czy się zdało jeszcze przed. Można wykorzystać dane historyczne aby spróbować przewidzieć przyszłość. Czy można to nazwać AI? IMO tak. Naturalne jest że wiele decyzji podejmujemy na postawie poprzednich doświadczeń. Niech komputer na podstawie tego co zdarzyło się wcześniej przewidzi przyszłość. Możemy to nazwać uczeniem maszynowym. Pewne dane są dostępne i zrozumiane. Przekazujemy je do programu, który na ich podstawie uczy się pewnych wzorców. Kiedy skończy podajemy nowe dane a on przekazuje przewidywany wynik. Brzmi nieźle.

Drawing1

Zacznijmy od czegoś prostego.

Mamy dostępne dane na temat zdobytych punktów z egzaminu oraz opuszczonych wykładów. Na tej podstawie spróbujmy przewidzieć jaki wynik uzyska ktoś kto był dokładnie na połowie wykładów.
Zobaczmy jak to wygląda na wykresie.

dd

Jak widać dane układają się mniej więcej wzdłuż pewnej linii. Spróbujmy ją naszkicować na oko.

łan

Widać że dla 50% obecności możemy oczekiwać wyniku w okolicach 5 punktów. No dobrze, ale jak komputer ma do tego dojść?

Regresja liniowa jednej zmiennej.

Podsumujmy co mamy do tej pory: zbiór par wartości (obecność, punkty). Szukamy pewnej funkcji liniowej która pozwoli przewidywać wynik. Postać ogólną, którą wszyscy znają wygląda tak:
y = ax + b
Dalej będę używał formy:
h_\theta(x) = \theta_0 + \theta_1 x

Potrzebujemy wyznaczyć dwa współczynniki \theta_0 i \theta_1 . Tym zajmie się tytułowy algorytm.

Aby stwierdzić czy rozwiązanie jest dobre zdefiniujmy sobie coś co pozwoli ocenić jego jakość. Taka funkcja musi sprawdzać jak bardzo rzeczywiste dane różnią się od przewidywanych. Najprościej będzie użyć sumy kwadratów różnic, innymi słowy metody najmniejszych kwadratów.

Dla pojedynczej pary wyrażenie będzie wyglądać tak:
(h_\theta(x) - y)^2
Podnosimy do kwadratu różnicę między tymi wartościami. Gratisowo pozbywamy się ew. wartości ujemnych.
Jak będzie to wyglądać dla całego zbioru danych?

\sum\limits_{i=1}^m (h_\theta(x^{(i)}) - y^{(i)})^2
Gdzie x^{(i)} i y^{(i)} oznaczają po prostu i-tego x i y, a m liczność zbioru z danymi.

Ostatecznie zapiszmy to jako średnią powyższych odchyleń:
J(\theta_0, \theta_1) = \frac{1}{m} \sum\limits_{i=1}^m \left ( h_\theta(x^{(i)}) - y^{(i)} \right )^2

Warto zauważyć, że h_\theta(x^{(i)}) to to samo co \theta_0 + \theta_1 x

Będziemy starali się tak manipulować parametrami \theta_0 i \theta_1 aby zminimalizować funkcję J i w efekcie uzyskać lepsze dopasowanie.

Minimalizacja

Pamiętajmy, że rozwiązanie problemu sprowadza się teraz do znalezienia minimum funkcji J.

Funkcja J jest dwuargumentowa. Do wyznaczenia minimum (niekoniecznie globalnego!!!) posłużymy się algorytmem gradientu prostego. Nazwa sugeruje że jest to prosta metoda i tak jest w istocie.

Pomysł jest następujący:
1. Wybierz punkt początkowy (wartości dla \theta_0 i \theta_1 ).
2. Rozejrzyj się dookoła gdzie możesz zejść w dół (i zejdź).
3. Powtarzaj tak długo aż znajdziesz optimum/skończy się czas/etc.

Krok 1 jest oczywisty. Zastanówmy się nad drugim. Aby zdecydować w którym kierunku iść, będziemy potrzebowali pochodnych cząstkowych (jak zmienia się funkcja w obecnym położeniu). Wyglądają one następująco:

\frac{\partial}{\partial \theta_0} J(\theta_0, \theta_1) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right)
\frac{\partial}{\partial \theta_1} J(\theta_0, \theta_1) = \frac{2}{m} \sum\limits_{i=1}^m \left (h_\theta \left (x^{(i)} \right ) - y^{(i)}\right) \cdot x^{(i)}

W jaki sposób pochodne mają nam pomóc? Aby zobaczyć jak to działa weźmy funkcję kwadratową.
y=x^2-2x+3
Jej pochodna
y'=2x-2
Weźmy punkt startowy x = 10 i poszukajmy minimum. Będzie nam potrzebny jeszcze współczynnik \alpha . O nim za chwilę, przyjmijmy że ma wartość \alpha=0.1

Pochodna w tym miejscu jest równa y'(10) = 18
Zmodyfikujmy x następująco:
x_{i+1} = x_i - \alpha \cdot 18 1iter
Przesunęliśmy się we właściwym kierunku!
Wykonajmy więcej iteracji.

10iter

Jak widać osiągnęliśmy optimum. Czym jest \alpha? Jest to współczynnik mówiący jak ostro mamy iść w kierunku rozwiązania. Im mniejszy tym więcej iteracji trzeba wykonać, im większy tym szybciej ale musimy uważać na zbyt dużą wartość. Ustawmy \alpha=1.2, 20 iteracji.

overshoot

Coś nie zagrało. Nad przyczyną polecam zastanowić się samemu.

OK, sprawdziliśmy że metoda działa, przełóżmy ją na dwie zmienne.
Rozumowanie jest identyczne, jedyne o czym musimy pamiętać to jednoczesna zmiana obu parametrów (liczymy krok dla \theta_0 oraz \theta_1 a dopiero potem zmieniamy wartości tych zmiennych).

Co zobaczymy po implementacji?

lech

Seems legit

 

Zobaczmy jak będzie wyglądał kod.

X = [10,9,8,7,6,5,4,3,2,1,0]
Y = [1,2,4,3,5,4,6,7,4,8,9]

m = len(X)
alpha = 0.01

Wprowadzamy dane, ustawiamy m oraz α.

def h(x, theta):
    return theta[0] + theta[1] * x


def J_partial0(theta):
    sum = 0
    for i in range(m):
        sum += h(X[i], theta) - Y[i]

return 2 * sum / m;

def J_partial1(theta):
    sum = 0
    for i in range(m):
        sum += (h(X[i], theta) - Y[i]) * X[i]

    return 2 * sum / m;

Definiujemy fukcję h oraz pochodne funckcji J.

def gradient_descent(startTheta):
    theta = startTheta
    for i in range(5000):
        d0 = J_partial0(theta)
        d1 = J_partial1(theta)
        theta[0] -= alpha * d0
        theta[1] -= alpha * d1

    return theta

Pamiętamy aby obliczyć d0 oraz d1 przed ich dodaniem do aktualnego wyniku.

Podsumowanie

To wszystko, metoda jest bardzo prosta ale to dopiero pierwszy krok. Następnym razem: co jeśli mamy dodatkowe dane i chcemy je uwzględnić w uczeniu.

Cały kod dostępny tutaj: https://gist.github.com/Xevaquor/e128c3c31c48dafd5766

DISCLAIMER: Używana notacja jest luźna, niekoniecznie musi być w pełni formalnie poprawna.