Uczenie maszynowe i sztuczne sieci neuronowe/Wykład 6

Wstęp

W tym wykładzie zajmiemy się problemem konstrukcji Uogulnionych Modeli Liniowych (ang. generalized linear models - GML). Metodologia ta pozwala objąć w jednym formaliźmie zarówno problemy regresji jak i klasyfikacji. W pewnym sensie klasyfikacja jest podobna do regresji, z tą różnicą, że zmienne które chcemy przewidywać mogą przybierać tylko niewielką ilość dyskretnych wartości. Na początek skupimy się na dwóch nieco już oswojonych przykładach, pod koniec dołożymy jeden bardzo przydatny przykład klasyfikacji wielorakiej (z więcej niż dwoma klasami).

Zaczniemy od problemu klasyfikacji binarnej, czyli takiej w której wejściom mamy przypisywać jedną z dwóch klas (podobnie jak to było dla preceptronu Rosenblatta), np. oznaczonych 0 i 1. Na chwilę też odłożymy na bok sieci i zajmiemy się tym zagadnieniem bardziej z punktu widzenia statystyki.

Regresja logistyczna

Hipoteza

W tym podejściu ignorujemy fakt, że zbiór wartości jest dyskretny. Jako funkcję pełniącą rolę hipotezy wybierzemy sigmiodę uzyskaną z funkcji logistycznej, znaną nam już jako jedną z popularnych form nieliniowości neuronów:

[math]g(s) = \frac{1}{1+ \exp(-s)}[/math]

która wraz z parametrami [math]\theta[/math] i wejściami [math]x[/math] jest postaci:

[math]h_\theta(x) = g(\theta^T x) = \frac{1}{1+ \exp(-\theta^T x)}[/math]

Pod koniec wykładu okaże się dlaczego taki akurat wybór hipotezy jest bardzo naturalny.

Estymacja parametrów

Jak znaleźć parametry [math]\theta[/math]? W języku funkcji kosztu moglibyśmy oczywiście zapostulować odpowiednią funkcję kosztu i zastosować do niej minimalizację gradientową. Tu jednak spróbujemy pokazać, że analogiczny algorytm można też wyprowadzić z interpretacji probabilistycznej. Dzięki temu będziemy mogli nabrać nowego wglądu w proces doboru parametrów.

Hipotezę wybraliśmy tak, że:

[math] P(y=1|x;\theta) = h_\theta(x) [/math]

[math] P(y=0|x;\theta) = 1- h_\theta(x) [/math]

Zauważmy, że powyższe wyrażenia można zapisać w zwartej formie:

[math]P(y|x;\theta) = \left(h_\theta(x)\right)^y \left(1-h_\theta(x)\right)^{1-y}[/math]

Funkcja wiarygodności

Zakładając, że przykłady zbioru uczącego są niezależne od siebie prawdopodobieństwo zaobserwowania całego zbioru uczącego [math]{\left(X^{(j)},Y^{(j)}\right)}_{j=1,\dots,m}[/math] wynosi:

[math]P(Y|X;\theta) = \prod_{j=1}^m P(y^{(j)}|x^{(j)};\theta)[/math]

Możemy to prawdopodobieństwo potraktować jako funkcję parametrów [math]\theta[/math], nazywamy ją wówczas funkcją wiarygodności i oznaczamy [math]L(\theta)[/math].

[math]L(\theta)= \prod_{j=1}^m P(y^{(j)}|x^{(j)};\theta) = \prod_{j=1}^m \left(h_\theta(x^{(j)})\right) ^{y^{(j)}} \left(1-h_\theta(x^{(j)})\right)^{1-y^{(j)}}[/math]

Łatwiejsza w posługiwaniu się jest funkcja log-wiarygodności:

[math]l(\theta) = \log L(\theta) = \sum_{j=1}^m y^{(j)} \log h_{\theta}(x^{(j)}) + (1 - y^{(j)}) \log (1 - h_{\theta}(x^{(j)}))[/math]

Dobre parametry [math]\theta[/math] to te, dla których zaobserwowanie ciągu uczącego jest największe. Aby je znaleźć należy zmaksymalizować funkcję wiarygodności, czy też dowolną monotonicznie rosnącą funkcję funkcji wiarygodności np. log-wiarygodność. Robimy to modyfikując parametry zgodnie z jej pochodną:

[math] \begin{array}{lcl} \frac{\partial}{\partial \theta_i} l(\theta) &=& \sum_{j=1}^m\left(y^{(j)}\frac{1}{g(\theta^T x^{(j)})} - (1-y^{(j)})\frac{1}{1-g(\theta^Tx^{(j)})} \right) \frac{\partial}{\partial \theta_i} g(\theta^T x^{(j)})\\ &=& \sum_{j=1}^m \left(y^{(j)}\frac{1}{g(\theta^T x^{(j)})} - (1-y^{(j)})\frac{1}{1-g(\theta^Tx^{(j)})} \right) g(\theta^T x^{(j)})(1-g(\theta^T x^{(j)})) \frac{\partial}{\partial \theta_i} (\theta^T x^{(j)})\\ &=& \sum_{j=1}^m \left( y^{(j)} (1-g(\theta^T x^{(j)})) - (1-y^{(j)})g(\theta^T x^{(j)})\right)x_i^{(j)}\\ &=& \sum_{j=1}^m (y^{(j)}-h_\theta(x^{(j)}))x_i^{(j)} \end{array} [/math]

skorzystaliśmy po drodze z postaci pochodnej funkcji logistycznej. Zatem aby zwiększać funkcję wiarygodności powinniśmy parametry zmieniać zgodnie z obliczoną pochodną:

[math]\theta_i^{(j+1)} =\theta_i^{(j)} + \alpha \sum_{j=1}^m (y^{(j)} - h_\theta( x^{(j)}) )x_i^{(j)} [/math]

Może to się wydać dziwne, ale startując z zupełnie innych założeń i stosując optymalizację innej funkcji dostaliśmy taką samą regułę zmiany parametrów jak przy gradientowej minimalizacji funkcji (średniokwadratowej) kosztu

(proszę porównać)!

Uogólnione modele liniowe

Dotychczas rozważaliśmy przykłady regresji gdzie zwykle [math]y \in \mathcal{N}(\mu,\sigma^2)[/math] oraz klasyfikacji gdzie [math]y \in \text{Bernoulli}(\phi)[/math]. Wkrótce przekonamy się, że oba te problemy sa szczególnymi przypadkami większej rodziny modleli, tzw. uogólnionych modeli liniowych.

Rodzina wykładnicza

Mówimy, że dany rozkład należy do rodziny wykładniczej jeśli da się go zapisać w postaci:

[math] p(y;\eta) = b(y) \exp(\eta^T T(y) - a(\eta))[/math] (*)

tutaj:

• [math]\eta[/math] nazywana jest parametrem naturalnym lub kanonicznym dystrybucji;
• [math]T(y)[/math] jest tzw. statystyką wystarczającą (często [math]T(y) = y[/math]);
• wielkość [math]\exp(-a(\eta))[/math] jest czynnikiem normalizującym, takim aby rozkład [math]p(y;\eta)[/math] sumował/całkował się do 1.

Tak więc członka rodziny wykładniczej określamy podając konkretne postaci [math]\eta,\, T(y),\, a(\eta)\, \text{oraz} \, b(y)[/math].

Rozkład Bernouliego

Pokażemy teraz, że rozkład Bernouliego należy do rodziny wykładniczej. [math]\text{Bernouli}(\phi)[/math] to taki rozkład wartości [math]y \in \{0,1\}[/math], że [math]p(y=1;\phi) = \phi[/math] zaś [math]p(y=0;\phi) = 1-\phi[/math]. w sposób zwarty możemy napisać ten rozkład tak:

[math] \begin{array}{lcl} p(y;\phi) &=& \phi^y(1-\phi)^{1-y} \\ &=& \exp(y \log \phi + (1-y) \log (1-\phi))\\ &=& \exp\left( y \log \frac{\phi}{1-\phi} + \log(1-\phi)\right) \end{array} [/math]

Zatem w rozkładzie Bernouliego parametrem naturalnym jest [math] \eta = \log \frac{1}{1-\phi}[/math].

Warto zauważyć, że jeśli przekształcić to wyrażenie ze względu na [math]\phi[/math] to dostaniemy dobrzez znaną funkcję logistyczną: [math]\phi = \frac{1}{1+\exp(-\eta)} [/math].

Widzimy też, że:

[math]T(y) = y[/math]

[math]a(\eta) = -\log(1-\phi) = \log (1+exp(-\eta))[/math]

[math]b(y) = 1[/math]

Czyli możemy przedstawić ten rozkład w postaci (*).

Rozkład Gaussa

Teraz rozważymy rozkład Gaussa. Do estymacji parametrów regresji liniowej nie musieliśmy używać jego wariancji, więc dla uproszczenia obliczeń przyjmiemy, że nasz rozkład Gaussa ma wariancję [math]\sigma^2 = 1[/math]. Mamy:

[math] \begin{array}{lcl} p(y;\mu) &=& \frac{1}{\sqrt{2 \pi}} \exp\left( -\frac{1}{2} (y-\mu)^2\right)\\ &=& \frac{1}{\sqrt{2 \pi}} \exp\left(-\frac{1}{2}y^2 \right) \exp\left(\mu y - \frac{1}{2}\mu^2 \right) \end{array} [/math]

Widzimy więc, że rozkład Gaussa należy do rodziny wykładniczej z następującymi parametrami:

[math] \begin{array}{lcl} \eta &=& \mu \\ T(y) &=& y \\ a(\eta) &=& \mu^2/2 = \eta^2/2 \\ b(y) &=& \frac{1}{\sqrt{2 \pi}} \exp\left( -\frac{1}{2} y^2\right) \end{array} [/math]

Rodzina wykładnicza jest znacznie bogatsza. Zawiera w sobie rozkłady wielorakie, Poissona (do modelowania zliczeń), gamma i wykładnicze (np. interwałów czasowych) i wiele innych. W kolejnej sekcji podamy ogólny sposób na konstruowanie modeli, w których [math]y[/math] pochodzi z rozkładów wykładniczych.

Konstrukcja uogólnionego modelu liniowego

Załóżmy, że chcemy zbudować model służący do szacowania liczby (y) klientów odwiedzających sklep (lub witrynę) w dowolnej godzinie, na podstawie pewnych cech x, takich jak promocje, ostatnie reklamy, prognoza pogody, dzień tygodnia, itd. Wiemy, że rozkład Poissona zwykle daje dobry model zliczeń np. liczby odwiedzających. Wiedząc o tym, jak możemy wymyślić model dla naszego problemu? Na szczęście, rozkład Poissona należy do rodziny rozkładów wykładniczych, więc możemy zastosować uogólniony model liniowy (GLM). W tej sekcji opiszemy metodę konstruowania modeli GLM.

W ogólności chcielibyśmy przewidywać wartość zmiennej losowej (zależnej) y traktując ją jako funkcję zmiennej (niezależnej) x. Aby móc zastosować tu metodologię GLM musimy poczynić następujące założenia:

1. Zmienna y przy ustalonych x i [math]\theta[/math] podlega pewnemu rozkładowi wykładniczemu z parametrem [math]\eta[/math] tzn.: [math](y|x;\theta) \sim[/math] RodzinaWykładnicza([math]\eta[/math])
2. Naszym celem jest estymacja wartości oczekiwanej [math]T(y)[/math] mając dany [math]x[/math]. W wielu przypadkach mamy [math]T(y) = y[/math], co oznacza, że chcielibyśmy aby hipoteza otrzymana w wyniku uczenia spełniałą warunek [math]h(x)= E[y|x][/math]. (Zauważmy, że własność ta jest spełniona dla regresji liniowej i logistycznej. Np. dla regresji logistycznej mamy [math]h_\theta(x) = p(y=1|x;\theta) = 0 \cdot p(y=0|x;\theta) + 1\cdot p(y=1|x;\theta) = E[y|x;\theta][/math]).
3. Parametr naturalny [math]\eta[/math] jest liniowo związany z wejściem [math]x[/math]: [math]\eta = \theta^T x[/math].

Te trzy założenia pozwalają wyprowadzić klasę algorytmów uczących GLM. Poniżej przedstawimy trzy przykłady.

Regresja liniowa jako GLM

Przekonajmy się, że regresja liniowa jest szczególnym przykładem GLM. Zmienna zależna [math]y[/math] jest ciągła i jej prawdopodobieństwo warunkowe dla danego [math]x[/math] jest modelowane przez rozkład Gaussa [math]N(\mu,\sigma)[/math] ([math]\mu[/math] może zależeć od [math]x[/math]). Tak więc w tym wypadku wspomnianą w założeniu 1 RodzinąWykładniczą[math](\eta)[/math] jest rozkład Gaussa. Tak jak widzieliśmy wcześniej w tym przypadku [math]\eta = \mu[/math]. Dalej mamy:

[math] \begin{array}{lcl} h_\theta(x) &=& E[y|x;\theta]\\ &=& \mu \\ &=& \eta \\ &=& \theta^T x \end{array} [/math]

• pierwsza równość wynika z założenia 2,
• druga równość wynika z tego, że [math](y|x;\theta) \sim N(\mu,\sigma^2)[/math], tak więc wartość oczekiwana wynosi [math]\mu[/math]
• trzecia równość wynika z założenia 1
• ostatnia równość wynika z założenia 3

Zauważmy jak przyjęcie założeń co do postaci rozkładu zmiennej zależnej i metodologi GLM (trzy założenia) prowadzą do konkretnej postaci hipotezy.

Regresja logistyczna jako GLM

Rozważmy teraz regresję logistyczną. Interesuje nas tutaj klasyfikacja binarna, więc [math]y \in \{0,1\}[/math]. Ponieważ [math]y[/math] przyjmuje wartości binarne to naturalnym rozkładem prawdopodobieństwa do modelowania warunkowego rozkładu [math](y|x)[/math] jest rozkład Bernoulliego z parametrem [math]\phi[/math] ([math]\phi[/math] jest prawdopodbieństwem tego, że y=1). Pokazaliśmy już wcześniej, że rozkład ten należy do rodziny RozkładówWykładniczych i, że w tym przypadku [math]\phi = \frac{1}{1+\exp(-\eta)}[/math]. Dalej zauważmy, że jeśli [math](y|x;\theta) \sim \text{Bernoulli}(\phi)[/math], to [math]E[y|x;\theta] = \phi[/math], więc analogicznie jak dla regresji liniowej mamy:

[math] \begin{array}{lcl} h_\theta(x) &=& E[y|x;\theta]\\ &=& \phi \\ &=& \frac{1}{1+\exp(-\eta)} \\ &=& \frac{1}{1+\exp(-\theta^T x)} \end{array} [/math]

Otrzymaliśmy więc funkcję hipotezy w postaci [math]h_\theta(x) = \frac{1}{1+\exp(-\theta^T x)}[/math]. To jest wytłumaczenie dlaczego stosuje się funkcje logistyczne w problemach klasyfikacji: jak tylko założymy, że zmienna zależna podlega rozkładowi Bernoulliego to funkcja logistyczna jest konsekwencją definicji uogólnionych modeli liniowych i rodziny rozkładów wykładniczych.

Wprowadzimy tu jeszcze dwa pojęcia:

• funkcja [math]g[/math] dająca średnią rozkładu jako funkcję parametru naturalnego, [math]g(\eta) = E[T(y);\eta][/math] nazywamy kanoniczną funkcją odpowiedzi (canonical response function)
• jej odwrotność, [math]g^{-1}[/math] nazywamy kanoniczną funkcją łączącą (canonical link function.)

Regresja wieloraka (softmax) jako GLM

Rozważmy jeszcze jeden przypadek modelu GLM. Tym razem nasz problem polega na przydzieleniu zmiennych niezależnych do jednej z k klas, czyli zmienna zależna nadal jest dyskretna, ale może przyjmować jedną z k wartości: [math]y \in \{1,2,\dots,k\}[/math]. Mówimy, że zmienne y podlegają rozkładowi wielorakiemu (ang. multinomial).

Wyprowadzimy teraz GLM do modelowania takich wielorakich danych. Zaczniemy od wyrażenia rozkładu wielorakiego jako rozkładu należącego do rodziny rozkładów wykładniczych.

Aby sparametryzować wielorakość z k możliwymi wynikami, można by zacząć od k parametrów [math]\phi_1,\dots,\phi_k[/math] określających prawdopodobieństwo każdego z wyników. Taka parametryzacja jest jednak redundantna, tzn. parametry te nie są niezależne (znając [math]k-1[/math] spośród [math]\phi_i[/math] parametrów ostatni, k-ty, parametr jest jednoznacznie określony bo musi być spełniona równość [math]\sum_{i=1}^k \phi_i =1[/math].) Tak więc sparametryzujemy rozkład przez k-1 parametrów:

[math]\phi_1,\dots,\phi_{k-1}[/math],

gdzie:

• [math]\phi_i = p(y=i;\phi)[/math]
• [math]p(y=k;\phi) = 1 - \sum_{i=1}^{k-1} \phi_i[/math].

Dla wygody notacji zapiszemy, że [math]\phi_k = 1-\sum_{i=1}^{k-1} \phi_i[/math], ale będziemy pamiętać, że to nie jest parametr, i że nasz rozkład wieloraki jest w pełni określony przez parametry: [math]\phi_1,\dots,\phi_{k-1}[/math].

Aby wyrazić rozkład wieloraki w języku rodziny rozkładów wykładniczych zdefiniujmy [math]T(y) \in \mathcal{R}^{k-1}[/math] w następujący sposób:

[math] T(1) = \left[ \begin{array}{c} 1\\ 0\\ 0\\ \vdots\\ 0 \end{array} \right], T(2) = \left[ \begin{array}{c} 0\\ 1\\ 0\\ \vdots\\ 0 \end{array} \right], T(3) = \left[ \begin{array}{c} 0\\ 0\\ 1\\ \vdots\\ 0 \end{array} \right],\dots, T(k-1) = \left[ \begin{array}{c} 0\\ 0\\ 0\\ \vdots\\ 1 \end{array} \right], T(k) = \left[ \begin{array}{c} 0\\ 0\\ 0\\ \vdots\\ 0 \end{array} \right] [/math]

W odróżnieniu od poprzednich przykładów nie mamy tu [math]T(y) = y[/math], ale [math]T(y)[/math] jest k-1 wymiarowym wektorem a nie skalarem. Aby oznaczyć i-ty element tego wektora będziemy pisać [math](T(y))_i[/math].

I jeszcze jedna użyteczna konwencja. Wprowadźmy funkcję [math]1\{\cdot\}[/math], przyjmuje ona wartość 1 gdy jej argument jest prawdziwy i 0 gdy jest fałszywy, np.: [math]1\{2==3\} = 0[/math] zaś [math]1\{3==7-4\} = 1[/math]. Tak więc [math](T(y))_i = 1\{y==i\}[/math]. Dalej mamy [math]E[(T(y))_i] = P(y=i) = \phi_i[/math].

Teraz możemy pokazać jak rozkład wieloraki wyrazić w języku rozkładu wykładniczego i sformułować model klasy GLM.

[math] \begin{array}{lcl} p(y;\phi) &=& \phi_1^{1\{y==1\}}\phi_2^{1\{y==2\}} \dots \phi_k^{1\{y==k\}}\\ &=& \phi_1^{1\{y==1\}}\phi_2^{1\{y==2\}} \dots \phi_k^{1-\sum_{i=1}^{k-1}1\{y==i\}}\\ &=& \phi_1^{(T(y))_1}\phi_2^{(T(y))_2}\dots\phi_k^{1-\sum_{i=1}^{k-1}(T(y))_i}\\ &=& \exp\left[ (T(y)_1\log(\phi_1)) + (T(y)_2\log(\phi_2)) + \dots + (1-\sum_{i=1}^{k-1} (T(y))_i)\log(\phi_k)\right]\\ &=& \exp \left[ (T(y))_1 \log \frac{\phi_1}{\phi_k} + (T(y))_2 \log \frac{\phi_2}{\phi_k} + \dots + (T(y))_{k-1} \log \frac{\phi_{k-1}}{\phi_k} + \log(\phi_k)\right]\\ &=& b(y) \exp(\eta^T(y) - a(\eta)) \end{array} [/math]

gdzie:

[math] \begin{array}{lcl} \eta &=& \left[ \begin{array}{c} \log \frac{\phi_1}{\phi_k}\\ \log \frac{\phi_2}{\phi_k}\\ \vdots \\ \log \frac{\phi_{k-1}}{\phi_k} \end{array} \right]\\ a(\eta)&=& -\log(\phi_k) \\ b(y) &=& 1 \end{array} [/math]

To kończy prezentację rozkładu wielorakiego jako członka rodziny rozkładów wykładniczych.

Funkcja łącząca (dla [math]i = 1,\dots,k[/math]) dana jest przez:

[math]\eta_i = \log \frac{\phi_i}{\phi_k}[/math]

dla wygody zdefiniowaliśmy także [math]\eta_k = \log \frac{\phi_k}{\phi_k}[/math] Stąd mamy:

[math]\exp(\eta_i) = \frac{\phi_i}{\phi_k}[/math] (**)

[math]\phi_k \exp(\eta_i) = \phi_i[/math]

Suma po wszystkich możliwych zdarzeniach daje zdarzenie pewne:

[math]\phi_k \sum_{i=1}^{k} \exp(\eta_i) = \sum_{i=1}^k \phi_i = 1[/math]

Stąd:

[math]\phi_k = \frac{1}{ \sum_{i=1}^{k} \exp(\eta_i)}[/math]

Podstawiając to wyrażenie do (**) otrzymujemy funkcję odpowiedzi postaci:

[math]\phi_i = \frac{\exp(\eta_i)}{\sum_{j=1}^k \exp(\eta_j)}[/math]

Ta funkcja mapująca [math]\eta[/math] na [math]\phi[/math] nazywa się funkcją softmax. Aby dokończyć formulację modelu użyjemy założenia 3, że [math]\eta[/math] jest liniowo związana ze zmienną niezależną [math]x[/math]. Tak więc mamy:

[math]\eta_i = \theta_i^T x[/math] dla [math]i = 1, \dots,k-1[/math] gdzie [math]\theta[/math] to parametry modelu. Dla wygody notacji definiujemy [math]\theta_k = 0[/math]. Wynika stąd, że [math]\eta_k = \theta_k^T x = 0[/math].

Zatem model nasz zakłada, że rozkład warunkowy [math](y|x)[/math] dany jest przez:

[math] \begin{array}{lcl} p(y=i|x;\theta) &=& \phi_i \\ &=&\frac{\exp(\eta_i)}{\sum_{j=1}^k \exp(\eta_j)} \\ &=& \frac{\exp(\theta_i^T x)}{\sum_{j=1}^k \exp(\theta_j^T x)} \end{array} [/math]

W wyprowadzonym powyżej modelu regresji softmax hipoteza ma postać:

[math] \begin{array}{lcl} h_\theta(x) &=& E[T(y)|x;\theta]\\ &=& E \left[ \begin{array}{lcl} 1\{y==1\} &|& \\ 1\{y==2\} &|& \\ \vdots &|& \\ 1\{y==k-1\} &|& \end{array} x;\theta \right]\\ &=& \left[ \begin{array}{l} \phi_1 \\ \phi_2 \\ \vdots \\ \phi_{k-1} \end{array} \right]\\ &=& \left[ \begin{array}{l} \frac {\exp(\theta_1^T x)}{\sum_{j=1}^k \exp(\theta_j^T x)}\\ \frac {\exp(\theta_2^T x)}{\sum_{j=1}^k \exp(\theta_j^T x)} \\ \vdots \\ \frac {\exp(\theta_{k-1}^T x)}{\sum_{j=1}^k \exp(\theta_j^T x)} \end{array} \right] \end{array} [/math]

Wyrażając to słowami: nasza hipoteza zwróci prawdopodobieństwo warunkowe przynależności danego x do każdej z klas i:

[math]p(y=i|x;\theta) [/math] dla [math]i = 1, \dots,k[/math],

przy czym prawdopodobieństwo przynależności do ostatniej klasy dane jest przez:

[math]p(y=k|x;\theta) = 1-\sum_{j=1}^{k-1} \phi_j[/math].

Na koniec zastanówmy się jak estymować parametry tego modelu. Podobnie jak w przypadku regresji liniowej i regresji logistycznej potrzebny jest nam zbiór uczący postaci [math]\left\{(x^{(j)},y^{(j)})\right\}_{ j =1,\dots,m }[/math]. Można ponownie skorzystać z zasady największej wiarygodności i wyznaczyć parametry [math]\theta[/math], które maksymalizują prawdopodobieństwo zaobserwowania całego zbioru uczącego. Funkcja log-wiarygodności ma postać:

[math] \begin{array}{lcl} l(\theta) &=& \sum_{j=1}^m \log p(y{(j)}|x^{(j)};\theta) \\ &=& \sum_{j=1}^m \log \prod_{i=1}^k \left( \frac{\exp(\theta_i^T x^{(j)} )}{\sum_{n=1}^k \exp (\theta_n^T x^{(j)})} \right)^{1\{y^{(j)}==i\}} \end{array} [/math]

Teraz maksymalizację [math]l(\theta)[/math] można przeprowadzić np. za pomocą algorytmu gradientowego (tzn. zmieniamy iteracyjnie parametry w kierunku zgodnym z gradientem funkcji log-wiarygodności).