Podręcznik

Sieci samoorganizujące, dla których podstawą uczenia jest konkurencja między neuronami stanowią zwykle sieci jednowarstwowe, w których każdy neuron połączony jest ze wszystkimi składowymi $N$-wymiarowego wektora wejściowego $\mathbf{x}$, jak to przedstawiono schematycznie dla $N = 2$ na rys. 7.1 [46].

Wagi połączeń synaptycznych połączeń$i$-tego neuronu tworzą wektor $\mathbf{w}_i=[w_{i1}, w_{i2}, …, w_{iN}]^T$. Przy założeniu normalizacji wektorów wejściowych po pobudzeniu sieci wektorem $\mathbf{x}$ zwycięża we współzawodnictwie neuron, którego wagi najmniej różnią się od odpowiednich składowych tego wektora. Zwycięzca, neuron $w$-ty, spełnia relację

gdzie $d \left(\mathbf{x}, \mathbf{w}_i \right)$ oznacza odległość w sensie wybranej metryki między wektorem x i wektorem w_i, a n jest liczbą neuronów. W uczeniu (adaptacji wag) stosuje się bądź strategię typu Winner Takes All (WTA) lub Winner Takes Most (WTM).

W strategii WTM wokół neuronu zwycięzcy przyjmuje się topologiczne sąsiedztwo S_w(k) o określonym promieniu malejącym w czasie. Neuron zwycięzca i wszystkie neurony położone w obszarze sąsiedztwa podlegają adaptacji, zmieniając swoje wektory wag w kierunku wektora x, zgodnie z regułą Kohonena [46].

w której jest współczynnikiem uczenia i-tego neuronu z sąsiedztwa S_w(k) w k-tej chwili. W regule Kohonena wartość może być taka sama lub różna dla wszystkich neuronów z sąsiedztwa. Wagi neuronów spoza sąsiedztwa S_w(k) nie podlegają zmianom. Rozmiar sąsiedztwa oraz wartości współczynników uczenia poszczególnych neuronów są funkcjami malejącymi w czasie. Zostało wykazane, że przy takim sposobie uczenia funkcja gęstości rozkładu wektorów w_i poszczególnych neuronów jest zbliżona do zdyskretyzowanego rozkładu gęstości wektorów wymuszeń. Po przyłożeniu dwóch różnych wektorów x, np. x₁ i x₂, uaktywnią się dwa neurony sieci, każdy reprezentujący wagi najbliższe współrzędnym wektorów odpowiednio x₁ i x₂. Wagi te oznaczone w postaci wektorowej w₁ i w₂ mogą być zilustrowane w przestrzeni jako dwa punkty. Zbliżenie do siebie wektorów x₁ i x₂ powoduje podobne zmiany położeń wektorów w₁ i w₂. W granicy w₁=w₂ wtedy i tylko wtedy, gdy x₁ i x₂ są sobie równe lub zbliżone do siebie. Sieć spełniająca te warunki nazywa się mapą topograficzną (mapą Kohonena).

Inny sposób uczenia sieci reprezentuje strategia WTA. Zmiana wag zachodzi tu również według zależności (7.2), ale adaptacja dotyczy jedynie neuronu zwycięzcy. Neurony przegrywające konkurencję nie zmieniają swoich wag.

7.1.1 Miary odległości między wektorami

Proces samoorganizacji wymaga na każdym etapie wyłonienia zwycięzcy, czyli neuronu, którego wektor wagowy różni się najmniej od przyłożonego na wejściu sieci wektora x. Istotnym problemem staje się w tej sytuacji wybór metryki, w jakiej mierzona będzie odległość między wektorem x i wektorem w_i. Najczęściej używanymi miarami odległości są [10,16]:

• miara euklidesowa

• miara Mahalanobisa

gdzie S oznacza macierz kowariancji obu wektorów. W tej definicji wielkości nieskorelowane mają większy wpływ na końcowy wynik niż te, które są mocno lub średnio skorelowane.

• iloczyn skalarny

• miara według normy $L_1$ (Manhattan)

• miara według normy $L_ \infty$ (Czebyszewa)

• miara Minkowskiego $L_m$

Przy użyciu miary euklidesowej podział przestrzeni na regiony dominacji neuronów jest równoważny mozaice Voronoia, w której przestrzeń wokół punktów centralnych stanowi strefę dominacji danego neuronu. Zastosowanie innej miary przy samoorganizacji kształtuje podział stref wpływów poszczególnych neuronów nieco inaczej. Szczególnie zastosowanie iloczynu skalarnego bez normalizacji wektorów może prowadzić do niespójnego podziału przestrzeni, przy którym występuje kilka neuronów w jednym regionie, a w innym nie ma żadnego.

7.1.2 Normalizacja wektorów

Wykazano, że proces samoorganizacji prowadzi zawsze do spójnego podziału przestrzeni danych, jeśli wektory $\mathbf{x}$ podlegają normalizacji. Przy znormalizowanych wektorach uczących $\mathbf{x}$, wektory wag, nadążając za nimi, stają się automatycznie znormalizowane (norma wektora równa 1). Jednakże normalizacja wektora wagowego powoduje, że jeśli $\| \mathbf{x}_i\|=const$, wtedy dla wszystkich neuronów iloczyn $\| \mathbf{w}_i\| \| \mathbf{x}_i\|$ jest także stały przy określonej wartości $\mathbf{x}$. O zwycięstwie neuronu decyduje więc wartość kąta między wektorami $\mathbf{x}$ i $\mathbf{w}_i$ zgodnie z zależnością $\cos(\mathbf{x},  \mathbf{w}_i )$.

Badania eksperymentalne potwierdziły potrzebę stosowania normalizacji wektorów przy małych wymiarach przestrzeni, np. $N=2$, $N=3$, natomiast nie jest już tak istotna dla przestrzeni o bardzo dużych wymiarach. Normalizacja wektorów może być dokonana w dwojaki sposób [45]:

• redefinicja składowych wektora według wzoru

• zwiększenie wymiaru przestrzeni o jeden, przy takim wyborze dodatkowej składowej $(N+1)$-szej wektora $\mathbf{x}$, że

Przy wyborze tego sposobu normalizacji zachodzi zwykle konieczność wcześniejszego przeskalowania składowych wektora $\mathbf{x}$ w przestrzeni $R^N$ w taki sposób, aby możliwe było spełnienie równości (7.9). Przy zwiększaniu wymiaru wektora wejściowego efekt normalizacji staje się coraz mniej widoczny i przy dużych wymiarach sieci, $N > 100$, normalizacja nie odgrywa większej roli w procesie samoorganizacji.

7.1.3 Problem neuronów martwych

Przy losowej inicjalizacji wag sieci i ograniczeniu wielkości sąsiedztwa neuronów podlegających adaptacji może się zdarzyć, że część neuronów znajdzie się w strefie, w której nie ma danych lub ich liczba jest znikoma. Neurony takie mają niewielkie szanse na zwycięstwo i adaptację swoich wad, pozostając martwymi. W ten sposób dane wejściowe reprezentowane będą przez mniejszą liczbę neuronów (część neuronów martwa), a błąd reprezentacji danych, zwany również błędem kwantyzacji

(\ \mathbf{w}(\mathbf{x}_i) \) oznacza wektor wagowy neuronu zwycięzcy dla wektora wejściowego $\mathbf{x}_i$, natomiast $p$ jest liczbą danych) staje się większy. Istotnym problemem jest zatem uaktywnienie wszystkich neuronów sieci.

Można to osiągnąć, jeśli w algorytmie uczącym będzie uwzględniona liczba zwycięstw poszczególnych neuronów, a proces uczenia zostanie zorganizowany w taki sposób, aby dać szansę zwycięstwa neuronom mniej aktywnym. Sugestia takiej organizacji uczenia bierze się z obserwacji zachowania neuronów biologicznych, gdzie neuron wygrywający tuż po zwycięstwie pauzuje przez określony czas „odpoczywając” przed następnym współzawodnictwem [46]. Taki sposób uwzględniania aktywności neuronów nazywany jest mechanizmem zmęczenia. Oryginalna nazwa conscience mechanism jest w języku polskim określana również mianem mechanizmu świadomości.

Istnieje wiele mechanizmów uwzględniających aktywność neuronów w procesie uczenia. Jednym ze sposobów uaktywnienia wszystkich neuronów jest uwzględnienie liczby zwycięstw neuronu przy obliczaniu efektywnej odległości wag od wzorca uczącego $\mathbf{x}$, modyfikując ją proporcjonalnie do liczby zwycięstw danego neuronu w przeszłości (na początku wszystkim przypisuje się liczbę zwycięstw równą 1). Jeśli oznaczy się liczbę zwycięstw $i$-tego neuronu przez $N_i$, modyfikację taką można przyjąć w postaci

Neurony aktywne o dużej wartości $N_i$ są karane sztucznym zawyżeniem tej odległości. Należy zaznaczyć, że modyfikację odległości stosuje się jedynie przy wyłanianiu zwycięzcy. W momencie uaktualniania wagi bierze się pod uwagę odległość rzeczywistą. Modyfikacja odległości ma za zadanie uaktywnić wszystkie neurony przez wprowadzenie ich w obszar o dużej liczbie danych. Po wykonaniu zadania (zwykle po dwóch lub trzech cyklach uczących) wyłącza się ją, pozwalając na niezakłóconą konkurencję neuronów [46].

Celem uczenia sieci samoorganizujących się poprzez współzawodnictwo jest takie uporządkowanie neuronów (dobór wartości ich wag), które zminimalizuje wartość oczekiwaną zniekształcenia (błędu kwantyzacji), określaną jako błąd popełniany przy aproksymacji wektora wejściowego $\mathbf{x}$ wartościami wag neuronu zwyciężającego w konkurencji. Przy $p$ wektorach wejściowych $\mathbf{x}$ i zastosowaniu normy euklidesowej błąd kwantyzacji może być wyrażony w postaci zależności (7.10). Zadanie to zwane jest również kwantyzacją wektorową (ang. Vector Quantization - VQ). Numery neuronów zwyciężających przy kolejnych prezentacjach wektorów x tworzą tak zwaną książkę kodową.

7.2.1 Algorytmy WTA i WTM

W klasycznym rozwiązaniu kodowania stosuje sią algorytm K-uśrednień (ang. K-means) określany również jako uogólniony algorytm Lloyda [65]. W przypadku sieci neuronowych odpowiednikiem algorytmu Lloyda jest algorytm WTA (ang. Winner Takes All). W algorytmie tym po prezentacji wektora $\mathbf{x}$ określana jest jego odległość od wektorów wagowych wszystkich neuronów. Zwycięzcą staje się neuron o najmniejszej odległości. Zwycięzca ($w$-ty neuron) ma przywilej adaptacji swoich wag w kierunku wektora $\mathbf{x}$ według reguły [46]

Pozostałe neurony nie podlegają adaptacji, czyli

Algorytm umożliwia uwzględnienie zmęczenia neuronów poprzez uwzględnienie liczby zwycięstw neuronu i faworyzowanie jednostek o najmniejszej aktywności, dla wyrównania ich szans na zwycięstwo. Modyfikację taką, jak zaznaczono uprzednio, stosuje się głównie w początkowych fazach algorytmu, wyłączając ją po uaktywnieniu wszystkich neuronów. Ostatni sposób uczenia jest zaimplementowany w programie Kohon w postaci opcji CWTA. CWTA pochodzi od angielskiej nazwy algorytmu Conscience Winner Takes All i jest jednym z najlepszych i najszybciej działających algorytmów samoorganizacji.

Oprócz algorytmów WTA, w których tylko jeden neuron może podlegać adaptacji w każdej iteracji, w uczeniu sieci samoorganizujących stosuje się powszechnie algorytmy typu WTM (ang. Winner Takes Most), w których oprócz zwycięzcy uaktualniają swoje wagi również neurony z jego sąsiedztwa, przy czym im dalsza jest odległość neuronu od zwycięzcy, tym mniejsza jest zwykle zmiana wartości wag tego neuronu. Proces adaptacji wektora wagowego może być opisany uogólnioną zależnością [46] z funkcją sąsiedztwa $G(i, \mathbf{x} )$ $i$-tego neuronu względem wektora wejściowego $\mathbf{x}$

dla wszystkich neuronów $i$ należących do sąsiedztwa zwycięzcy $w$. We wzorze tym oddzielono współczynnik uczenia $\eta$ każdego neuronu od jego odległości względem neuronu zwycięzcy lub od prezentowanego wektora x, mającej wpływ na funkcję sąsiedztwa $G(i, \mathbf{x})$. Definiując tę funkcję w postaci

gdzie w oznacza numer zwycięzcy, otrzymuje się klasyczny algorytm WTA. Istnieje wiele różnych odmian algorytmów WTM, różniących się przede wszystkim postacią funkcji $G(i, \mathbf{x})$. Przedstawimy tu dwa rozwiązania: algorytm gazu neuronowego oraz algorytm klasyczny Kohonena dostosowany do tworzenia mapy.

7.2.2 Algorytm gazu neuronowego

Znaczącą poprawę samoorganizacji sieci w sensie mniejszego błędu kwantyzacji można uzyskać, stosując metodę zwaną przez autorów algorytmem gazu neuronowego ze względu na podobieństwo dynamiki zmian położenia neuronów do ruchu cząsteczek gazu.

W algorytmie tym wszystkie neurony podlegają sortowaniu w każdej iteracji w zależności od ich odległości ich wektorów wagowych względem wektora x. Neurony są ustawiane w kolejności odpowiadającej narastającej odległości $d_0$, gdzie $d_m=\left\|\mathbf{x}-\mathbf{w}_m\right\|$ oznacza odległość neuronu zajmującego w wyniku sortowania $m$-tą pozycję w szeregu za neuronem zwycięzcą, któremu przyporządkowano odległość $d_0$. Wartość funkcji sąsiedztwa dla $i$-tego neuronu jest określona zależnością [46]

w której zmienna $m(i)$ oznacza kolejność uzyskaną w wyniku sortowania $(m(i) = 0,1,\ldots, n-1)$, gdzie 0 oznacza neuron zwycięzcę, a $\lambda$ jest parametrem (analogicznie do promienia sąsiedztwa w algorytmie Kohonena), malejącym w czasie. Przy $\lambda=0$ tylko neuron zwycięzca podlega adaptacji i algorytm WTM przekształca się w zwykły algorytm WTA, natomiast przy $\lambda > 0$ adaptacji podlegają wagi wielu neuronów w stopniu uzależnionym od wartości funkcji $G(i, \mathbf{x})$.

Algorytm gazu neuronowego przypomina strategię zbiorów rozmytych, w której każdemu neuronowi przypisuje się wartość funkcji przynależności do sąsiedztwa, określoną zależnością (7.18). W celu uzyskania dobrych rezultatów samoorganizacji proces uczenia powinien rozpoczynać się z dużą wartością $\lambda$, po czym powinna ona maleć wraz z upływem czasu do wartości minimalnej (najczęściej zerowej). Zmiana $\lambda(k)$ może być liniowa bądź wykładnicza. Współczynnik uczenia i-tego neuronu $\lambda(k)$ może również zmieniać się bądź liniowo bądź wykładniczo. W praktyce lepsze wyniki organizacji uzyskuje się zwykle przy liniowej zmienności $\lambda(k)$.

Algorytm gazu neuronowego, obok CWTA uwzględniającego aktywność neuronów jest w praktyce najskuteczniejszym narzędziem samoorganizacji neuronów w sieci ze współzawodnictwem. Dobierając odpowiednio parametry sterujące procesem można uzyskać bardzo dobrą organizację sieci przy szybkości działania znacznie przewyższającej klasyczny algorytm Kohonena. Tym nie mniej należy zauważyć, że algorytm ten nie ma zastosowania przy tworzeniu mapy Kohonena, gdyż nie przenosi adaptacji wag na neurony położone w obszarze bliskiego sąsiedztwa (nie zachowuje spójności obszaru aktywności neuronów przy adaptacji).

7.2.3 Algorytmy uczące mapy Kohonena

Algorytm Kohonena należy do najstarszych metod uczenia WTM generujących mapę i w chwili obecnej istnieje wiele jego odmian. W klasycznym algorytmie Kohonena dostosowanym do tworzenia mapy inicjalizuje się sieć, przyporządkowując neuronom określone miejsce w przestrzeni i stowarzyszając je (łącząc) z sąsiadami. Stowarzyszenie to może być interpretowane na siatce prostokątnej (sąsiedztwo prostokątne) lub heksagonalnej (sąsiedztwo heksagonalne). Interpretacja obu rodzajów sąsiedztwa przedstawiona jest na rys. 7.2. Kontury wewnętrzne o kształcie prostokąta lub sześciokąta odpowiadają granicom sąsiedztwa o promieniu odpowiednio zerowym (jedynie zwycięzca – brak sąsiadów), pierwszym, drugim lub trzecim.

W momencie wyłonienia zwycięzcy uaktualnieniu podlegają nie tylko jego wagi, ale również wagi jego sąsiadów, pozostających w najbliższym sąsiedztwie. W ten sposób adaptacji podlega cała grupa neuronów wokół zwycięzcy zbliżając swoje wagi do prezentowanego wektora x, jak to pokazano w sposób graficzny na rys. 7.3.

Funkcja sąsiedztwa $G(i, \mathbf{x})$ w mapie Kohonena może zmieniać się w sposób skokowy (sąsiedztwo prostokątne) lub ciągły. W pierwszym sposobie zwanym przez Kohonena funkcją „bubble” jest ona definiowana w postaci

We wzorze tym $d(i,w)$ może oznaczać bądź odległość euklidesową między wektorami wag zwycięzcy (neuron $w$-ty) oraz neuronu $i$-tego, bądź odległość mierzoną w liczbie neuronów. Współczynnik $\lambda$ jest promieniem sąsiedztwa o wartości malejącej, poczynając od zadanej na wstępie wartości maksymalnej aż do zera, przy czym wielkość ta zmienia się w sposób skokowy, wyłączając z sąsiedztwa w kolejnych iteracjach neurony znajdujące się na krańcach obszaru sąsiedztwa. Zauważmy, że tak zdefiniowana funkcja sąsiedztwa $G(i, \mathbf{x})$ przyjmuje identyczne wartości dla wszystkich neuronów znajdujących się aktualnie w obszarze sąsiedztwa podlegającym adaptacji.

Drugim powszechnie stosowanym typem funkcji sąsiedztwa w mapach Kohonena jest sąsiedztwo typu gaussowskiego, w którym funkcja $G(i, \mathbf{x})$ określona jest wzorem [46]

O stopniu adaptacji neuronów z sąsiedztwa zwycięzcy decyduje nie tylko odległość euklidesowa $d(i,w)$ neuronu $i$-tego od zwycięzcy (neuronu w-tego), ale również promień sąsiedztwa $\lambda$. W odróżnieniu od sąsiedztwa typu prostokątnego, gdzie każdy neuron należący do sąsiedztwa zwycięzcy podlegał adaptacji w jednakowym stopniu, przy sąsiedztwie typu gaussowskiego stopień adaptacji jest zróżnicowany i zależy od wartości funkcji Gaussa. Dla zwycięzcy obowiązuje $G(w, w) = 1$, dla pozostałych neuronów $0 < G(i, w) < 1$. Sąsiedztwo gaussowskie prowadzi zwykle do lepszych rezultatów uczenia i lepszej organizacji sieci niż sąsiedztwo prostokątne.

Istotny wpływ na uczenie sieci Kohonena wywiera dobór wartości współczynnika uczenia $\eta$ w poszczególnych iteracjach. W ogólności wartość ta powinna maleć w miarę czasu uczenia. Na rys. 7.4 przedstawiono typowe sposoby zmian tego współczynnika.

W praktyce najlepsze wyniki samoorganizacji uzyskuje się przy zastosowaniu liniowej formy zmienności współczynnika uczenia.

W mapie Kohonena definiuje się z góry siatkę dwuwymiarową która będzie odzwierciedlała rzut położenia danych wielowymiarowych na płaszczyźnie. Należy z góry określić liczbę neuronów rozmieszczonych w osi poziomej i pionowej. Ogólna liczba neuronów w sieci będzie wówczas równa ich iloczynowi. Neurony te zajmują odpowiednie miejsce w przestrzeni dwuwymiarowej, zachowując powiązanie ze swoimi sąsiadami w procesie uczenia. W efekcie po zakończenia uczenia na etapie testowania zwycięstwo określonego neuronu jest przenoszone na pozycję tego neuronu w siatce dwuwymiarowej (mapie) odzwierciedlając w ten sposób jego rzut na mapę topograficzną. Mapa Kohonena stanowi więc sposób przedstawienia rozkładu danych wielowymiarowych na płaszczyźnie.

Rozkład ten może być przedstawiony na kilka sposobów. Rys. 7.5 przedstawia trzy najczęściej stosowane kształty: mapa płaska, cylindryczna oraz toroidalna. W mapie cylindrycznej dane na dwu przeciwległych krańcach zbliżone są do siebie. W mapie toroidalnej zbliżenie dotyczy dwu przeciwległych krańców w osi x i y.

Dane które w przestrzeni wielowymiarowej są odległe od siebie zajmują również odległe pozycje na mapie. Dane bliskie sobie w przestrzeni wielowymiarowej są również bliskie na mapie.

7.2.4 Przykład mapy Kohonena w analizie obciążeń elektroenergetycznych

Na rys. 7.6 przedstawiono mapę Kohonena o wymiarze 7×7 dla danych dotyczących obciążeń elektroenergetycznych w Polskim Systemie Elektroenergetycznym (PSE). Uwzględniono 12 rodzajów obciążeń odpowiadających 12 miesiącom roku (Jan – styczeń, Feb – luty, Mar – marzec, Apr – kwiecień, May – May, June – czerwiec, July – lipiec, Aug – sierpień, Sep – wrzesień, Oct – październik, Nov – listopad oraz Dec - grudzień). Dane rzeczywiste obciążeń dotyczą wektorów 24-wymiarowych (obciążenia 24 godzin doby). Dzięki mapie Kohonena zostały one zrzutowane na płaszczyznę reprezentując mapę obciążeń odpowiadających poszczególnym miesiącom roku [46].

Zauważmy, że dane dotyczące miesięcy letnich i zimowych są odległe od siebie, podczas gdy te związane z wiosną i jesienią są podobne i położone na mapie w bliskim sąsiedztwie. Podobnie można zrzutować obciążenia 24-godzinne na mapę odzwierciedlającą podział na dni tygodnia ( Mon – poniedziałki, Tue – wtorki, Wed – środy, Thu – czwartki, Fri – piątki, Sat – soboty, Sun – niedziele), jak to przedstawiono na rys. 7.7.

Wyraźnie widoczne są skupienia neuronów reprezentujących dni robocze (od poniedziałku do piątku) oraz oddzielnie świąteczne (Sun). Pomiędzy nimi znajdują się soboty (Sat).

Interesująca jest również mapa profili obciążenia reprezentowana przez wagi neuronów zwycięzców dla poszczególnych rejonów (rys. 7.8). Można zauważyć znaczne podobieństwo miedzy najbliższymi sąsiadami, zarówno w pionie jak i poziomie. Ponadto jest ewidentne, że wagi neuronów odległych od siebie na mapie (np. neurony górne i dolne lub neurony z prawej i lewej strony w tym samym wierszu mapy) są znacznie różniące się.

Rozkład danych wielowymiarowych uzyskanych w sieci samoorganizującej może być przedstawiony na płaszczyźnie lub w przestrzeni trójwymiarowej przy zastosowaniu rzutowań innych niż mapa Kohonena. Jednym ze znanych jest nieliniowe odwzorowanie Sammona [57]. Odwzorowanie to pozwala na rzutowanie danych z dowolnej przestrzeni N-wymiarowej w przestrzeń M-wymiarową (np. M=2 lub M=3) zachowując podstawowe cechy rozkładu danych z oryginalnej przestrzeni wielowymiarowej.

Niech będzie danych n wektorów N-wymiarowych x_i (i=1, 2, …n). Odpowiednio do nich definiuje się n wektorów w przestrzeni M-wymiarowej oznaczonych przez y_i. Odległości między poszczególnymi wektorami w przestrzeni N-wymiarowej oznaczane będą przez $d_{i j}^*=d\left(\mathbf{x}_i, \mathbf{x}_j\right)$ a w przestrzeni M-wymiarowej przez $d_{i j}=d\left(\mathbf{y}_i, \mathbf{y}_j\right)$. W określeniu odległości między wektorami można zastosować dowolną metrykę, w szczególności euklidesową. Zadanie odwzorowania nieliniowego Sammona polega na takim doborze wektorów y, aby zminimalizować funkcję błędu E zdefiniowaną wzorem [57]

gdzie

W zależnościach tych y_ij oznacza j-tą składową wektora y_i. W minimalizacji funkcji błędu (7.19) Sammon zastosował uproszczoną metodę optymalizacyjną Newtona, która pozwala wyrazić rozwiązanie z kroku na krok w sposób rekurencyjny w postaci

Wzór wyrażający poprawkę  reprezentuje iloraz odpowiedniej składowej gradientu przez diagonalny składnik hesjanu, określony w k-tej iteracji. Współczynnik $\eta$ jest odpowiednikiem stałej uczenia i przyjmowany jest z zakresu [0,3, 0,4]. Przy definicji funkcji błędu w postaci (7.19) odpowiednie składowe gradientu i hesjanu opisane są wzorami [57]

W wyniku wielu iteracji składowe wektorów y_i przyjmują wartości ostateczne minimalizujące wartość zdefiniowanej na wstępie funkcji błędu.

Na rys. 7.9 przedstawiono rzutowanie Sammona dla tych samych danych dotyczących obciążeń elektroenergetycznych w Polskim Systemie Elektroenergetycznym (PSE) przedstawionych na mapie Kohonena. Uwzględniono rodzaje obciążeń odpowiadających czterem porom roku (kolor zielony – wiosna, czerwony – lato, magenta – jesień oraz niebieski – zima). Dane rzeczywiste obciążeń dotyczą wektorów 24-wymiarowych (obciążenia 24 godzin doby).

Zauważmy, że również przy tym rzutowaniu dane dotyczące lata i zimy są odległe od siebie, podczas gdy te związane z wiosną i jesienią są bliskie sobie i położone na mapie w bliskim sąsiedztwie. Ponadto charakterystyczne jest, że dane dotyczące zimy są stosunkowo mało rozproszone w stosunku do innych pór roku, natomiast dane wiosenne charakteryzują się rozproszeniem największym (duże różnice między najwyższą i najniższą temperatura powodują duże zróżnicowanie w poborze energii elektrycznej).

Pakiet programów KOHON, napisany w języku Matlab wywoływany jest poleceniem kohon [46]. Umożliwia zarówno uczenie jak i testowanie sieci Kohonena. Może służyć do rozwiązywania różnorodnych problemów praktycznych. Oprócz programów uczących i testujących, pakiet zawiera szereg programów niezbędnych do wygenerowania danych uczących, konwersji obrazów graficznych i plików dźwiękowych na dane numeryczne, jak również programy do graficznego odzwierciedlenia danych uczących, mapy oraz wyników uczenia. Obsługa tego programu jest intuicyjnie prosta.

Na rys. 7.10. przedstawiono wygląd okna menu głównego programu. Pola Learning data, Testing data oraz Network z lewej strony menu służą do zadawania odpowiednio danych uczących, testujących oraz zapamiętywania i wczytywania zapamiętanej sieci Kohonena. Format danych powinien zawierać w pierwszej linii wymiar wektora x. Następne linie zawierają kolejne wektory x stanowiące zasadniczą treść danych.

Środkowa cześć menu (Network parameters) zawiera podstawowe informacje dotyczące definicji struktury sieci i sposobu tworzenia mapy odwzorowań. Są to:

• Network dimension – określa ilość neuronów umieszczonych w osi poziomej i pionowej. Należy podać dwie liczby, np. 5 7, z których pierwsza określa rozmiar mapy w osi x (równy 5) a druga w osi y (rozmiar 7). Liczba neuronów jest równa wówczas iloczynowi obu liczb (w przykładzie n=35).

• Neighborhood function – określa rodzaj zastosowanej funkcji sąsiedztwa. Dostępne są następujące typy: bubble –prostokątne, gaussian – gaussowskie, cutgaussian – zmodyfikowane gaussowskie z obciętą i podwyższoną podstawą oraz ep=max(0, 1-x²).

• Topology – rodzaj siatki określającej topologię sąsiedztwa. Dostępna jest siatka prostokątna (rect) oraz heksagonalna (hexa).

• Shape – określa sposób rysowania mapy Kohonena. Możliwe są 3 typy przedstawień mapy: sheet (struktura płaska mapy), cylinder (struktura cylindryczna) i toroid (struktura toroidalna).

• Initialization type - sposób inicjalizacji wartości wag początkowych sieci. Są możliwe dwie metody: randinit (przypisanie losowych wartości wag) oraz lininit – przypisanie wag według specjalnej procedury korzystającej z zależności liniowych odpowiadających danym uczącym.

• Initialize – przycisk do faktycznej inicjalizacji wag sieci według wybranej wcześniej procedury inicjalizacji.

Z prawej strony menu występują parametry definiujące sposób uczenia sieci (pole Learning). Na pole to składa się:

• Wybór algorytmu uczenia (Learning algorithm) w ramach którego można wybrać spośród następujących metod uczenia: WTA (metoda Winner Takes All), CWTA (metoda WTA z mechanizmem zmęczenia), WTM batch (algorytm Kohonena aktualizujący wagi po prezentacji wszystkich danych z pliku uczącego), WTM seq (algorytm Kohonena aktualizujący wagi sieci w sposób on-line) oraz Neural gas (algorytm gazu neuronowego).

• Training length – ustala liczbę cykli (epok) uczących

• Radius initial – podaje początkową wielkość sąsiedztwa (w liczbie neuronów) przy zastosowaniu algorytmu Kohonena)

• Alpha initial – początkowa wartość stałej uczenia $\alpha$

Po ustawieniu tych wielkości należy przycisnąć przycisk Learn, uruchamiający proces uczenia. Testowanie wytrenowanej sieci odbywa się na danych testujących po przyciśnięciu przycisku Test w polu Testing. Program umożliwia różny sposób wizualizacji wyników. Dostępne są następujące opcje (po wybraniu odpowiedniego przycisku):

• Graph – ilustracja położeń danych 2D na tle danych uczących.

• Map – przedstawienie położenia danych wielowymiarowych jako mapy Kohonena. Aby to uzyskać należy przeprowadzić uczenie sieci przy użyciu algorytmu Kohonena (możliwe jest krótkie douczenie przy zastosowaniu innych algorytmów).

• Sammon 2D i Sammon 3D – ilustracja położenia danych uczących odpowiednio w przestrzeni 2-wymiarowej (2D) lub trójwymiarowej (3D) przy zastosowaniu algorytmu Sammona.

• Voronoi – ilustracja położenia danych 2D i neuronów jak mozaiki Voronia.

W skład pakietu Kohon wchodzą również programy konwertujące dane obrazów graficznych zapisane w formacie PCX (pcxkoh.exe) na format plików pakietu Kohon i koh2pcx.exe dokonujące konwersji w drugą stronę. Program koherror.bat umożliwia natomiast obliczenie błędu PSNR (w dB) wyrażającego różnicę pomiędzy obrazem oryginalnym, a odtworzonym.

Przy ocenie jakości sieci neuronowej samoorganizującej ważną rolę odgrywa odwzorowanie danych jedno- i dwuwymiarowych, ze względu na czytelny i przejrzysty sposób interpretacji wyników na płaszczyźnie $x, y$. Biorąc pod uwagę, że wagi neuronów są odpowiednikiem współrzędnych punktów centralnych klastrów, na jakie dzielony jest zbiór danych, można każdemu wektorowi wagowemu przypisać odpowiedni punkt na płaszczyźnie. Łącząc te punkty z najbliższymi sąsiadami otrzymuje się regularną siatkę, odwzorowującą topograficzny rozkład danych (siatka prostokątna lub heksagonalna). Przy równomiernym rozkładzie wektorów uczących $x$ na płaszczyźnie spodziewane odwzorowanie wagowe poszczególnych neuronów przedstawione na płaszczyźnie powinno być równomierne. Jeśli rozkład danych jest nierównomierny, zagęszczenie wystąpi tam, gdzie prawdopodobieństwo wystąpienia wektorów uczących jest większe (patrz rys. 7.14).

Dobrym testem algorytmów uczących sieci jest odwzorowanie kształtu danych jednowymiarowych za pomocą sieci samoorganizującej. Przykładowo na rys. 7.15 zilustrowano proces samoorganizacji 40 neuronów odwzorowujących układ danych tworzących kształt eliptyczny (pozycje neuronów zaznaczone kółkami). Rys. 7.15a przedstawia stan wyjściowy wag neuronów (rozkład losowy), rys. 7.15b – stan po dwóch cyklach uczących, rys. 7.15c – stan po pięciu cyklach uczących, a rys. 7.15d – stan końcowy (po dziesięciu cyklach) na tle danych uczących tworzących kształt eliptyczny. Uczenie było przeprowadzone przy użyciu programu Kohon i algorytmu gazu neuronowego.

Na rys. 7.16 podano wyniki końcowe uporządkowania neuronów przy odwzorowaniu różnego rodzaju kształtów utworzonych przez dane na płaszczyźnie (x,y). W każdym przypadku wagi neuronów plasują się w strefach występowania danych w taki sposób, aby zminimalizować błąd kwantyzacji. Należy nadmienić, że podstawowa zaleta sieci samoorganizujących się jest widoczna w całej pełni dopiero przy danych wielowymiarowych, gdzie zawodzi zdolność człowieka do wyobrażenia sobie ich rozkładu przestrzennego.

Mechanizmy samoorganizacji wbudowane w algorytmy uczące takich sieci działają niezależnie od wymiarowości problemu. Istotną funkcją pełnioną przez sieć jest tutaj kwantowanie wektorowe, polegające na tym, że ogromna liczba danych tworzących klaster jest odwzorowana przez wektor wagowy jednego neuronu reprezentującego każdy punkt w przestrzeni. Rozkład przestrzenny neuronów pozwala więc określić rozkład skupień danych w przestrzeni wielowymiarowej oraz zasadnicze cechy statystyczne rozkładu danych, przydatne z punktu widzenia użytkownika (centra grup, rozproszenie mierzone wariancją danych w grupie, odległości między centrami różnych grup, itp.)

Podstawową cechą sieci samoorganizujących się jest kompresja danych, polegająca na tym, że duże grupy danych tworzących klaster reprezentowane są przez pojedynczy wektor wagowy neuronu zwycięzcy. Przy p danych podzielonych na P klastrów i reprezentacji każdego klastra przez jeden z n neuronów uzyskuje się znaczne zmniejszenie ilości informacji, zwane kompresją. Jest to kompresja typu stratnego, której towarzyszy pewien błąd kwantowania określony zależnością (7.10).

Przykładem wykorzystania kompresyjnych własności sieci Kohonena jest kompresja stratna obrazów, mająca za zadanie zmniejszenie ilości informacji reprezentującej dany obraz, przy zachowaniu błędu odwzorowania na określonym poziomie (zapewnienie odpowiednio dużej wartości współczynnika PSNR mierzącego stosunek sygnału do szumu). Zakłada się, że obraz o wymiarach $N_x \times N_y$ pikseli podzielony jest na równomierne ramki zawierające $n_x \times n_y$ pikseli. Piksele każdej ramki stanowią składowe wektorów wejściowych $\mathbf{x}$. Każdy zawiera $n_x n_y$ składników, reprezentujących stopień szarości poszczególnych pikseli w ramce. Przyporządkowanie pikselom wektora może się odbywać przez złączenie poszczególnych wierszy ramki w jeden ciąg bądź przez zastosowanie przyporządkowania typu zygzakowego, stosowanego między innymi w standardzie JPEG.

Sieć samoorganizująca zawiera n neuronów, każdy połączony wagami synaptycznymi ze wszystkimi składnikami wektora wejściowego $\mathbf{x}$. Uczenie sieci przy zastosowaniu jednego z algorytmów samoorganizacji polega na takim doborze wag poszczególnych neuronów, aby zminimalizować błąd kwantyzycji (7.10). W wyniku procesu uczenia następuje taka organizacja sieci, przy której wektorowi $\mathbf{x}$ każdej ramki odpowiada wektor wagowy neuronu zwycięzcy minimalizujący błąd kwantyzacji. Przy podobnym ukształtowaniu składników wektora $\mathbf{x}$ różnych ramek zwyciężać będzie ten sam neuron albo grupa neuronów o podobnych wektorach wagowych. Podczas kolejnej prezentacji ramek ustala się numery neuronu zwycięzcy dla poszczególnych ramek np. 1, 1, 3, 80 itd. Numery neuronów zwycięzców tworzą książkę kodową, a wagi tych neuronów reprezentują uśrednioną wartość odpowiadającą poszczególnym składowym wektora $\mathbf{x}$ (stopień szarości pikseli tworzących ramkę). Wektory te będą przy odtwarzaniu obrazu reprezentować poszczególne ramki. Biorąc pod uwagę, że liczba neuronów $n$ jest zwykle dużo mniejsza od liczby ramek $N_r$, otrzymuje się zmniejszenie ilości informacji przypisanej danemu obrazowi (kompresja).

Przy określaniu stopnia kompresji należy uwzględnić również określoną liczbę bitów użytych do zakodowania numerów neuronów zwycięzców dla poszczególnych ramek. Ostatecznie współczynnik kompresji obrazu definiuje się w postaci [46].

gdzie $n_x$ i $n_y$ oznaczają wymiary ramki w osiach $x , y, N_r$ - liczbę ramek, $n$ - liczbę neuronów, a $T$ i $t$ - liczbę bitów przyjętych w reprezentacji odpowiednio stopnia szarości pikseli i wartości wag. Jakość obrazu odtworzonego charakteryzuje się najczęściej przy użyciu miary PSNR, określonej wzorem

gdzie $MSE$ oznacza wartość błędu średniokwadratowego obrazu odtworzonego względem obrazu oryginalnego. Zastosowanie sieci samoorganizującej do kompresji umożliwia uzyskanie współczynnika kompresji obrazów rzędu nawet 16 przy współczynniku $PSNR$ około 25-28 dB.

Na rys. 7.17a przedstawiono wyniki uczenia sieci Kohonena dla obrazu Barbara o wymiarach $512 \times 512$ pikseli, podzielonego na ramki 16-elementowe (podobrazy o wymiarach $4 \times 4$) [43]. Sieć Kohonena zawierała 512 neuronów. Przy 8-bitowej reprezentacji danych uzyskano stopień kompresji równy: $K_r = 9.9$. Na rys. 7.17b pokazano obraz odtworzony na podstawie wag neuronów zwyciężających przy prezentacji kolejnych ramek. Współczynnik $PSNR$ dla obrazu odtworzonego wynosił 26,2dB. Różnice w jakości obrazu oryginalnego (rys. 7.17a) i odtworzonego (rys. 7.17b) są stosunkowo małe. Widoczne jest to na rys. 7.17c, oddającym obraz błędu, który jest różnicą między obrazem oryginalnym a odtworzonym po kompresji.

  

Ważną zaletą sieci neuronowych, z całą wyrazistością uwidaczniającą się przy kompresji obrazów, jest zdolność generalizacji, a więc możliwość skompresowania (przyporządkowania poszczególnym ramkom nowego obrazu numerów neuronów zwycięzców wytrenowanej wcześniej sieci) i zdekompresowania (przypisania odpowiednich wektorów wagowych zwycięzców poszczególnym ramkom). Jakość odtworzonego obrazu, który nie podlegał wcześniej procesowi uczenia, nie odbiega daleko od jakości obrazu poddanego uczeniu, pod warunkiem, że stopnie zróżnicowania obu obrazów są podobne.

Na rys. 7.18 przedstawiono przykładowy obraz linii papilarnych, poddany procesowi kompresji i dekompresji za pomocą sieci Kohonena, wytrenowanej przy zastosowaniu obrazu z rys. 7.17 [43]. Rysunek 9.18a jest obrazem oryginalnym, rys. 7.18b - obrazem zrekonstruowanym po kompresji, a rys. 7.18c obrazem różnicowym ilustrującym błąd kompresji. Stopień zniekształcenia obrazu odtworzonego jest porównywalny z obrazem poddanym uczeniu, a współczynnik zniekształcenia $PSNR = 26,4 dB$.

  

1. Określić i porównać ze sobą odległości między dwoma wektorami

$\mathbf{x}_1 =\begin{bmatrix} 0,2 \\ 0,6 \\ -0,5 \\ 0.4 \\ 0.9 \\ 0.7  \end{bmatrix},\;   \mathbf{x}_2 =\begin{bmatrix} -0,2 \\ 0,1 \\ 0,6 \\ -0.3 \\ 0.4 \\ -0.2  \end{bmatrix}.\;$

stosując różne normy: euklidesową, $L_1, L_\infty, L_M$ oraz w postaci iloczynu skalarnego.

2. Dwa neurony o wagach $\mathbf{w}_1 = [ 0.8 \; 0.9 \; 0.3]$,  $\mathbf{w}_2=[-0.1 \; 0.7 \; 0.5]$ otrzymały pobudzenia $\mathbf{x}_1 = [0.5 \; 0.7 \; 0]$ ,  $\mathbf{x}_2 = [0.1 \; 0.7 \; 0.4], \mathbf{x}_3 =[0.2 \; 0.5 \; 0.3]$. Wyłonić zwycięzców dla każdego pobudzenia.

3. Trzy neurony o wagach $\mathbf{w}_1 =[0.3 \; 0.8 \; 0.9]$, $\mathbf{w}_2 =[-0.2 \; 0.5 \; 0.1]$, $\mathbf{w}_3 =[-0.7 \; -0.4 \; 0.6]$ otrzymały pobudzenie w postaci wektora $\mathbf{x} =[0.2 \; 0.6 \; 0.5]$. Określić kolejność neuronów konkurujących ze sobą w algorytmie gazu neuronowego. Wyznaczyć wartość funkcji sąsiedztwa w metodzie gazu neuronowego przy założeniu $\lambda = 0.5$.

4. Dokonać normalizacji wektorów $\mathbf{x}$ danych w postaci $\mathbf{x}_1 = [1 \; 3 \; 8]$, $\mathbf{x}_2 = [5 \;9 \; 20]$ stosując rozszerzenie wektora i bez rozszerzania.

5. Korzystając z programu KOHON dokonać odwzorowania liniowego (jeden z wymiarów w osi $x$ lub $y$ równy $1$) danych jednowymiarowych tworzących różne kształty (spirala, koło, wielobok) stosując odwzorowanie Kohonena z uwzględnieniem dwu sąsiadów.

6. Korzystając z programu KOHON określić położenia wektorów wag neuronów mapy Kohonena odwzorowujących różne rozkłady danych na płaszczyźnie: rozkład równomierny w założonym kształcie maski, rozkład gaussowski, maski nieregularne. Zastosować odwzorowanie Kohonena z uwzględnieniem czterech sąsiadów.

7. Korzystając z programu KOHON określić położenia wektorów wag neuronów mapy Kohonena odwzorowujących rozkład danych na płaszczyźnie przy zadanym kształcie maski (koło, prostokąt trójkąt). Zastosować odwzorowanie jednowymiarowe Kohonena z uwzględnieniem tylko dwu sąsiadów (jeden z wymiarów $n_x$ lub $n_y$ równy $1$).

Analiza składników głównych PCA jest metodą statystyczną określającą przekształcenie liniowe $\mathbf{y} = \mathbf{Wx}$ transformujące opis stacjonarnego procesu stochastycznego dany w postaci zbioru N-wymiarowych wektorów $\mathbf{x}$ w zbiór wektorów $\mathbf{y}$ o zredukowanej wymiarowości $K \le N$ [13,46]. Przekształcenie to odbywa się za pośrednictwem macierzy $\mathbf{W}$ o wymiarach $K \times N$ w taki sposób, że przestrzeń wyjściowa o zredukowanym wymiarze zachowuje najważniejsze informacje dotyczące procesu. Innymi słowy, transformacja PCA zamienia dużą ilość informacji zawartej we wzajemnie skorelowanych danych wejściowych w zbiór statystycznie niezależnych składników, według ich ważności. Stanowi zatem formę kompresji stratnej, znanej w teorii komunikacji jako transformacja Karhunena--Loeve.

Dla zachowania maksimum informacji oryginalnej w zbiorze wektorów $\mathbf{y}$ o zredukowanym wymiarze macierz transformacji $\mathbf{W}$ powinna być dobrana w taki sposób, aby zmaksymalizować wartość wyznacznika $J$ [49]

W wyrażeniu tym $\mathbf{R}_{\mathbf{xx}}$ oznacza macierz kowariancji wektorów $\mathbf{x}$ (przy zerowych wartościach średnich zbioru $x$ macierz kowariancji jest równa macierzy korelacji). W praktyce centrowanie wektorów nie jest konieczne i można posługiwać się macierzą korelacji, niezależnie od zerowania się wartości średnich. Rozwiązanie powyższego problemu optymalizacyjnego uzyskuje się na podstawie rozkładu macierzy kowariancji zbioru wektorów $\mathbf{x}$ według wartości własnych.

Przyjmijmy, że $\mathbf{x} = [x_1, x_2, \ldots, x_N]^T$ oznacza wektor losowy o zerowej wartości średniej, a $\mathbf{R}_{\mathbf{xx}} = E[\mathbf{xx}]^T$ oznacza wartość oczekiwaną (średnią) macierzy autokorelacji (autokowariancji) po wszystkich wektorach $\mathbf{x}$. Macierz tę, przy skończonej liczbie $p$ wektorów $\mathbf{x}$, można estymować przy pomocy zależności

gdzie macierz danych $\mathbf{X}$ tworzą kolejne wektory uczące $\mathbf{X} = [\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_p]^T$. Oznaczmy przez $\lambda_i$ wartości własne macierzy autokorelacji $\mathbf{R_{xx}}$, a przez $\mathbf{w}_i$ ortogonalne wektory wartości własnych, skojarzone z nimi, przy czym $\mathbf{w} = [w_{i1}, w_{i2}, , \ldots, w_{iN}]^T$. Wartości własne oraz wektory własne macierzy $\mathbf{R_{xx}}$ powiązane są zależnością

dla $i=1, 2, …, N$. Wartości własne symetrycznej, nieujemnie określonej macierzy korelacji $\mathbf{R_{xx}}$ są rzeczywiste i nieujemne. Uporządkujmy je w kolejności malejącej poczynając od wartości największej $\lambda_1 > \lambda_2 > \ldots \lambda_N \ge 0$. W identycznej kolejności ustawimy wektory własne $\mathbf{w}_i$ skojarzone z odpowiednimi wartościami własnymi $\lambda_i$. Przy ograniczeniu się do $K$ największych wartości własnych macierz $\mathbf{W}$ przekształcenia PCA definiuje się w postaci

Macierz ta określa transformację PCA jako przekształcenie liniowe

Wektor y =[y₁, y₂,…, y_K]^T stanowi wektor składników głównych PCA, mających największy wpływ na rekonstrukcję oryginalnego wektora danych x =[x₁, x₂,…, x_N]^T. Transformacja PCA jest więc ściśle związana z rozkładem macierzy korelacji według wartości własnych. Uwzględnia zatem jedynie powiązania liniowe między wektorami danych.

Oznaczmy przez $\mathbf{L}$ macierz diagonalną utworzoną z wartości własnych $\lambda_i$ uwzględnionych w odwzorowaniu, to jest $\mathbf{L}  = diag\{ {\lambda_1, \lambda_2, \ldots, \lambda_K}\}$. Przy takich oznaczeniach macierz korelacji można przedstawić w postaci zdekomponowanej

gdzie znak przybliżenia wynika z uwzględnienia skończonej liczby (mniejszej niż wymiar macierzy korelacji) składników głównych. Przy $K=N$ zależność (8.6) będzie spełniona ze znakiem równości. Uwzględniając symetrię macierzy korelacji zależność tę można jednocześnie przedstawić w postaci

Z punktu widzenia statystycznego transformacja PCA określa zbiór $K$ wektorów ortogonalnych (kolejne wiersze macierzy $\mathbf{W}$), mających największy wkład w wariancję danych wejściowych. Pierwszy składnik główny odpowiadający wektorowi $\mathbf{x}$ jest równy iloczynowi skalarnemu wektora własnego $\mathbf{w}_1$ i wektora wejściowego $\mathbf{x}$ 

Wektor własny $\mathbf{w}_1$ określa zatem kierunek w przestrzeni wielowymiarowej w którym występuje największa wariancja (zmienność wartości) danych wejściowych zawartych w wektorach $\mathbf{x}$. Wariancja ta jest równa wartości własnej $\lambda_1$

Celem transformacji PCA jest określenie kierunków $\mathbf{w}_1, \mathbf{w}_2, \ldots, \mathbf{w}_K$ (zwanych kierunkami głównymi) w taki sposób, aby zmaksymalizować wielkość iloczynu skalarnego $\mathbf{E} (\mathbf{w}_i^T\mathbf{x})$ dla kolejnych wartości $i=1, 2, …, K$ przy spełnieniu warunku ortogonalności kolejnych wektorów $\mathbf{w}$ ze sobą, to jest $\mathbf{w}_i^T \mathbf{w}_j = 0$ oraz $\mathbf{w}_i^T \mathbf{w}_i = 1$ .

Rekonstrukcja wektora $\mathbf{x}$ na podstawie znajomości wektora składników głównych $\mathbf{y}$ oraz macierzy ortogonalnej $\mathbf{W}$ przekształcenia PCA odbywa się zgodnie z zależnością [46]

Macierz $\mathbf{W}$ dekompozycji PCA i macierz rekonstrukcji ($\mathbf{W}^T$) stanowią wzajemne transpozycje. PCA minimalizuje wartość oczekiwaną błędu rekonstrukcji danych, przy czym błąd ten określony jest wzorem ogólnym

Przy ograniczeniu się do $K$ największych wartości własnych ($K$ składników głównych), błąd ten jest proporcjonalny do sumy odrzuconych wartości własnych $\sum_{i=K+1}^N \lambda_i$.

Ze wzoru tego wynika, że minimalizacja błędu rekonstrukcji danych przy uwzględnieniu $K$ składników jest równoważna maksymalizacji wariancji rzutowania $\varepsilon_v$ na etapie rozkładu PCA

Zarówno $\varepsilon_r$, jak i $\varepsilon_v$ są nieujemne, gdyż wszystkie wartości własne macierzy korelacji, jako macierzy symetrycznej i nieujemnie określonej, są dodatnie bądź zerowe. Wnosimy stąd, że reprezentacja wektora danych $\mathbf{x}$ przez największe składniki główne $y_1, y_2, \ldots, y_K$ tworzące wektor $\mathbf{y}$ jest równoważna zachowaniu informacji o największej porcji energii zawartej w zbiorze danych.

Pierwszy (największy) składnik główny powiązany z $\lambda_1$ przez swój wektor własny $\mathbf{w}_1$ określa kierunek w przestrzeni wielowymiarowej, w którym wariancja danych jest maksymalna. Ostatni najmniejszy składnik główny (ang. Minor Principal Component) wskazuje kierunek, w którym wariancja jest najmniejsza. Na rys. 8.1 przedstawiono interpretację geometryczną najbardziej znaczącego i najmniej znaczącego składnika głównego transformacji PCA dla danych 2-wymiarowych. Pierwszy składnik główny odpowiada kierunkowi największej zmienności mierzonej poprzez wariancję (energię) sygnałów. Dokonując reprezentacji danych tylko za pomocą jednego składnika głównego oraz skojarzonego z nim wektora własnego i wybierając jako reprezentanta największy ze składników głównych ($y_1$), popełnia się najmniejszy błąd rekonstrukcji, maksymalizując jednocześnie wariancję transformacji. Najmniej znaczący składnik główny ma najmniejszy wpływ na dokładność odtworzenia danych. Stąd kompresja danych (zmniejszenie ilości informacji z najmniejszą stratą dla rekonstrukcji) wymaga reprezentowania tych danych przez zbiór największych składników głównych. Pominięcie składników najmniejszych ma najmniej znaczący wpływ na dokładność rekonstrukcji danych.

Transformacja PCA jest ściśle związana z korelacją zachodzącą między wieloma zmiennymi w zbiorze danych. Jeśli te zmienne są skorelowane ze sobą, to znajomość jedynie części z nich wystarczy do określenia pozostałych. Stąd taki zbiór danych może być reprezentowany przez mniejszą liczbę zmiennych. W przypadku gdy nie występuje korelacja między zmiennymi tworzącymi wektor x, predykcja części z nich na podstawie pozostałych jest niemożliwa.

8.1.2 Przykład zastosowania PCA w ekonomii

Jako przykład ilustrujący właściwości rozkładu danych na składniki główne rozpatrzone będą dane GUS-u dotyczące wielkości związanych z miesięcznymi wartościami odpowiadającymi wskaźnikowi cen i usług konsumpcyjnych (wcu), stopie bezrobocia (sb), wartości produkcji sprzedanej w przemyśle (wps) oraz średniej płacy miesięcznej (spm) w Polsce. Przykład ten wykorzystuje dane statystyczne GUS z 10.5 lat (126 wartości).

Wektory pomiarowe x tworzą w tym przypadku cztery składowe, kolejno: wcu, sb, wps oraz spm. Stosując definicję macierzy korelacji (wzór 10.2) w wyniku obliczeń otrzymano macierz korelacji R_xx w postaci

$\mathbf{R}_{\mathbf{xx}}=10^6\left[\begin{array}{llll} 0.0108 & 0.0017 & 0.0175 & 0.2564 \\ 0.0017 & 0.0003 & 0.0027 & 0.0390 \\ 0.0175 & 0.0027 & 0.0286 & 0.4226 \\ 0.2564 & 0.0390 & 0.4226 & 6.2953 \end{array}\right]$

Dokonując dekompozycji tej macierzy według wartości własnych uzyskuje się następujące wartości własne (w kolejności malejącej):

$\lambda_1  = 6334324.26$

$\lambda_2 = 590.69$

$\lambda_3 = 24.10$

$\lambda_4 = 7.48$

oraz skojarzone z nimi wektory własne

$\mathbf{w}_1=\left[\begin{array}{r} 0.0406 \\ 0.0062 \\ 0.0669 \\ 0.9969 \end{array}\right]$, $\mathbf{w}_2=\left[\begin{array}{r} -0.7609 \\ -0.1859 \\ -0.6173 \\ 0.0736 \end{array}\right]$, $\mathbf{w}_3=\left[\begin{array}{r} -0.6464 \\ 0.2778 \\ 0.7102 \\ -0.0231 \end{array}\right]$, $\mathbf{w}_4=\left[\begin{array}{r} -0.0402 \\ -0.9424 \\ 0.3316 \\ -0.0148 \end{array}\right]$

Na tej podstawie określa się pełną macierz transformacji PCA, zawierającą wszystkie wektory własne ułożone według największego znaczenia (w zależności od wielkości wartości własnych $\lambda_i$): $\mathbf{W}=\left[\mathbf{w}_1, \mathbf{w}_2, \mathbf{w}_3, \mathbf{w}_4\right]^T$ w postaci

$\mathbf{W}=\left[\begin{array}{rrrr}0.0406 & 0.0062 & 0.0669 & 0.9969 \\ -0.7609 & -0.1859 & -0.6173 & 0.0736 \\ -0.6464 & 0.2778 & 0.7102 & -0.0231 \\ -0.0402 & -0.9424 & 0.3316 & -0.0148\end{array}\right]$

oraz macierz diagonalną $\mathbf{L}$ złożoną z wartości własnych macierzy $\mathbf{R_{xx}}$ ułożonych według malejących wielkości $\mathbf{L} = diag \{6334324.26, 590.69, 24.10, 7.48 \}$. Największa wartość własna $\lambda_1 = 6334324.26$ skojarzona jest z pierwszym składnikiem głównym odpowiadającym wektorowi własnemu $\mathbf{w}_1$, stanowiącemu pierwszy wiersz macierzy $\mathbf{W}$. Składnik ten przy wektorze wejściowym $\mathbf{x}$ złożonym z czterech elementów ($wcu, sb., wps, spm$) opisany jest relacją $y_1 = \mathbf{w}_1^T\mathbf{x}$, która w tym przypadku przybiera konkretną postać: $y_1 = 0,0406*wcu + 0,0062*sb + 0,0.0669*wps + 0.9969*spm$. Jak widać największy wpływ na składnik główny $y_1$ ma zmienna $spm$. Każda z wartości własnych $\lambda_i$ odpowiada wariancji jaką reprezentuje dany składnik główny. Względny wkład poszczególnych składników głównych w łączną wariancję danych (energię) można określić wzorem:

$r_i=\frac{\lambda_i}{\sum_{j=1}^4 \lambda_j}$

Wartości te są następujące: $r_1 = 0.9999, r_2 = 9,32E-5, r_3 = 3.8E-6, r_4=1.18E-6$. Jak wynika z rozkładu tych wartości, największy składnik główny ma $99,99\%$ udziału w łącznej wariancji danych. Przy odtwarzaniu wszystkich składników ($wcu, sb, wps, spm$) na podstawie wektora $\mathbf{y}$ można ograniczyć się jedynie do jego największej składowej $y_1$, pomijając pozostałe jako nie wnoszące istotnego wkładu informacyjnego. Oznacza to 3-krotną redukcję ilości przetworzonej informacji. Po odtworzeniu w ten sposób danych uzyskano odtworzony zbiór danych, ze średnim błędem względnym równym $0.97\%$ (zdefiniowanym jako stosunek normy euklidesowej macierzy błędu do normy euklidesowej danych przyjętych w postaci macierzy $4 \times 126$).

Klasyczna transformacja PCA wymaga wyznaczenia na podstawie danego zbioru wektorów $\mathbf{x}$ najpierw macierzy korelacji $\mathbf{R_{xx}}$, a następnie rozkładu tej macierzy według wartości własnych. Jest to prosta i łatwa w zastosowaniu metoda (ciąg dekompozycji QR), jeśli problem jest dobrze uwarunkowany (współczynnik uwarunkowania $cond (\mathbf{R_{xx}})$ przyjmuje małe wartości). Wystarczy w tym celu zastosować funkcję eig.m w Matlabie. Przy bardzo dużych wymiarach wektorów x (powyżej kilku tysięcy) problem obliczeniowy wyznaczenia wartości własnych macierzy $\mathbf{R_{xx}}$ staje się zwykle źle uwarunkowany i trudny do przeprowadzenia. W takich przypadkach konkurencyjne może być zastosowanie metod neuronowych adaptacji pozwalających na wyznaczenie macierzy przekształcenia $\mathbf{W}$ bez tworzenia macierzy korelacyjnej, a jedynie poprzez bezpośrednie przetwarzanie wektorów wejściowych $\mathbf{x}$. Przetwarzanie sygnałów może być wówczas przedstawione w typowo sieciowy sposób jako liniowa sieć PCA zawierająca jedną warstwę.

8.2.1 Estymacja pierwszego składnika głównego

Wyznaczenie wektorów własnych $\mathbf{w_{i}}$ może przebiegać przy zastosowaniu metod adaptacyjnych, typowych dla sieci neuronowych. Są one oparte na uogólnionej regule Hebba i pozwalają na bezpośrednie przetworzenie wektorów wejściowych $\mathbf{x}$, bez potrzeby jawnego definiowania macierzy $\mathbf{R_{xx}}$. Metody te są niezastąpione w przypadkach akwizycji danych w trybie on-line, w którym utworzenie jawnej postaci macierzy korelacji byłoby niemożliwe. Powstały różne odmiany algorytmów, a w każdym z nich wykorzystuje się korelację zachodzącą między wektorami reprezentującymi dane wejściowe.

Bezpośrednie zastosowanie uogólnionej reguły Hebba dotyczy wyznaczenia jednego (największego) składnika głównego. Algorytm ten odpowiada zastosowaniu sieci liniowej jednowarstwowej o jednym neuronie wyjściowym przedstawionej na rys. 8.2.

Sygnał wyjściowy takiej sieci określony jest wzorem

Dobór wag wektora  $\mathbf{w_{1}}$ odbywa się według uogólnionej (znormalizowanej) reguły Hebba, zwanej regułą Oji, którą zapisać można w formie skalarnej [24,46]

We wzorach tych $\eta$ oznacza współczynnik uczenia malejący w czasie. W procesie uczenia sieci powtarza się wielokrotnie te same wzorce uczące, aż do ustabilizowania się wag sieci. Pierwszy składnik obu wzorów odpowiada zwykłej regule Hebba, a drugi zapewnia samo-normalizację wektorów wagowych, to jest $\| \mathbf{w}_1\|^2=1$. Dobór wartości $\eta$ ma istotny wpływ na zbieżność algorytmu. Dobre rezultaty uzyskuje się, przyjmując wartości $\eta(k)$ malejące wraz z upływem czasu uczenia.

Jakkolwiek metoda Oji pozwala wyznaczyć jedynie pierwszy składnik główny może być ona w prosty sposób przystosowana do wyznaczenia pozostałych składników. Przy uwzględnieniu wielu składników i wektorów własnych można sukcesywnie redukować wpływ wyselekcjonowanego składnika pierwszego odkrywając w ten sposób następny, traktowany w przekształconym wektorze $\mathbf{x}$ jako pierwszy składnik główny, dla którego można powtórzyć poprzednia procedurę. Przyjmując oznaczenie $i$-tego wektora własnego w postaci $\mathbf{w_{i1}} = [w_{i1}, w_{i2}, \ldots, w_{iN}]^T$ można zauważyć, że składowe wektora  $\mathbf{x}$ można wyrazić w postaci

dla $i=1, 2, \ldots, N$. Składnik pierwszy tego wzoru po wyznaczeniu pierwszego składnika głównego jest już znany i może zostać wyeliminowany poprzez utworzenie nowego wektora $\mathbf{x}^\prime$ o elementach zdefiniowanych następująco

Dla nowego wektora $x_i^{\prime}$ największą wartością własną jest teraz $\lambda_2$, gdyż wpływ $\lambda_1$ został wyeliminowany dzięki działaniu wyrażonemu wzorem (8.17). Procedura wyznaczenia tej wartości jest identyczna do omówionej wcześniej. Powtarzając te operacje $N$ razy można wyznaczyć wszystkie wektory własne stowarzyszone z kolejnymi wartościami własnymi.

8.2.2 Algorytm estymacji wielu składników głównych jednocześnie

Wyznaczanie wielu kolejnych składników PCA wymaga zastosowania wielu neuronów w warstwie wyjściowej sieci. Sieć neuronowa zawiera zatem tyle neuronów, ile składników głównych rozkładu ma być uwzględnionych. Są one ułożone w jednej warstwie - stąd sieć PCA jest siecią jednowarstwową o liniowych funkcjach aktywacji neuronów (rys. 8.3).

Istnieje wiele algorytmów uczenia takiej sieci. Tutaj ograniczymy się do reguły Sangera [24,46], będącej regułą lokalną, nie wymagającą w procesie uczenia rozwiązania układu równań. Przy $K$ neuronach liniowych w warstwie sygnały wyjściowe $y_i$ w kolejnych iteracjach $k$ określane są według wzoru

Adaptacja wag sieci w kolejnych iteracjach wykorzystuje aktualnie znane (zaadaptowane) wartości wag i przebiega według następującego wzoru

dla $j=1, 2, …, N$ oraz $i=1, 2, …, K$. Przyjmując oznaczenie

Jest oczywiste, że nawet przy istnieniu $K$ neuronów w warstwie wyjściowej reguła uczenia pozostaje nadal lokalną, pod warunkiem zmodyfikowania wartości sygnału wejściowego $x_j^\prime$. Zauważmy, że modyfikacja tego sygnału odbywa się przy wykorzystaniu wcześniej zaadaptowanych już wartości wag, a więc nie wymaga rozwiązania układu równań. Zależności skalarne (8.20) i (8.21) można zapisać w prostszej formie wektorowej wygodniejszej w działaniach

dla $i=1, 2,…, K$. Zauważmy, że dla neuronu pierwszego (pierwszy składnik główny PCA) mamy $\mathbf{x}^{\prime}(k)=\mathbf{x}(k)$. Dla neuronu drugiego otrzymuje się $\mathbf{x}^{\prime}(k)=\mathbf{x}(k)-\mathbf{w}_1(k) y_1(k)$ - wzór uzależniony od znanych już wag neuronu pierwszego. Podobnie dla trzeciego neuronu $\mathbf{x}^{\prime}(k)=\mathbf{x}(k)-\mathbf{w}_1(k) y_1(k)-\mathbf{w}_2(k) y_2(k)$ i wszystkich pozostałych, modyfikacja wektora $\mathbf{x}$ wyrażona jest przez wielkości wag wcześniej określone, w efekcie czego proces uczenia przebiega identycznie jak w przypadku algorytmu Oji, z samo-normalizującymi się wektorami $\mathbf{x}$, czyli $\left\|\mathbf{w}_i\right\|=1$.

Obecnie istnieje wiele różnych algorytmów neuronowych, pozwalających adaptacyjnie określić parametry transformacji PCA. Do najważniejszych, oprócz algorytmu Sangera, zalicza się algorytm Foldiaka, Rubnera oraz APEX (ang. Adaptive Principal component EXtraction). Szczegóły rozwiązań znaleźć można w opracowaniu książkowym Diamantarosa i Kunga.

Główne zastosowania transformacji PCA związane są z kompresją danych, która jest nieodłącznym składnikiem każdego przekształcenia PCA. Własność ta może być bezpośrednio wykorzystana do kompresji stratnej informacji – mówimy wtedy o kompresji sygnałów (dane jednowymiarowe) bądź obrazów (dane 2-D). Niezastąpionym zastosowaniem PCA jest ilustracja rozkładu danych wielowymiarowych na płaszczyźnie (układ 2-współrzędnych w postaci 2 najważniejszych składników głównych) lub w przestrzeni 3-D (układ 3-współrzędnych w postaci 3 najważniejszych składników głównych). Ponadto transformacja PCA reprezentująca $N$-wymiarowy wektor $\mathbf{x}$ przez $K$-wymiarowy ($K$) wektor $\mathbf{y}$ składników głównych pozwala traktować wektor $\mathbf{y}$ jako wektor cech diagnostycznych procesu reprezentowanego przez zbiór wektorów $\mathbf{x}$.

8.3.1 PCA w zastosowaniu do kompresji stratnej danych

Kompresja danych z zastosowaniem PCA polega na przekształceniu $N$-elementowego wektora wejściowego $\mathbf{x}$ w wektor $\mathbf{y}$ o zmniejszonym wymiarze $K$ ($K$). Redukcja wymiaru wektora $\mathbf{x}$ poprzez PCA zapewnia optymalność przekształcenia poprzez zachowanie w wektorze zredukowanym największej możliwie dawki informacji oryginalnej (przy założonej wartości $K$).

Na rys. 8.4 przedstawiono sieć PCA do kompresji (rys. 8.4a) oraz do rekonstrukcji (dekompresji) danych (rys. 8.4b). W sieci kompresyjnej wektor oryginalny $\mathbf{x}$ jest transformowany w wektor $\mathbf{y}$ o zredukowanym wymiarze $K$, przy czym $\mathbf{y} = \mathbf{Wx}$. Wektor $\mathbf{y}$ może podlegać bądź to transmisji na odległość bądź zapisaniu na dysku. W każdym przypadku możliwe jest odtworzenie wektora oryginalnego na podstawie jego zredukowanej formy $\mathbf{y}$, korzystając z sieci rekonstrukcyjnej z rys. 8.4b, wykonującej operację odwrotną $\mathbf{\hat{y}} = \mathbf{W^Tx}$. Biorąc pod uwagę pewną utratę informacji spowodowaną obcięciem wymiaru wektora odtworzenie to jest z pewnym przybliżeniem $\hat{\mathbf{x}} \simeq \mathbf{x}$.

O współczynniku kompresji decyduje liczba składników głównych $K$ uwzględnionych w przekształceniu PCA. Przy dużej liczbie wektorów $\mathbf{x}$ podlegających przekształceniu można pominąć liczbę bitów do kodowania wag sieci i współczynnik kompresji wyrazić wzorem przybliżonym

Im wyższy współczynnik kompresji tym większa oszczędność pamięci, ale gorsza jakość odtworzonego obrazu (większa porcja informacji utracona bezpowrotnie w wyniku redukcji wymiaru wektora oryginalnego).

Na rys. 8.5 przedstawiono obraz oryginalny (rys. 8.5a) oraz trzy obrazy zrekonstruowane na podstawie odpowiednio 1, 3 i 5 składników głównych PCA (rys. 8.5 b,c,d). Obraz poddany kompresji miał wymiar $512 \times 512$ pikseli i został podzielony na ramki o wymiarach $8 \times 8$ (wymiar wektora $\mathbf{x}$ równy $64$). Jakość odtworzonego obrazu jest ściśle uzależniona od liczby $K$ składników głównych uwzględnionych w odtwarzaniu. Im więcej jest tych składników, tym lepsza jakość obrazu, ale mniejszy współczynnik kompresji. Przy największym współczynniku kompresji (jeden składnik główny) wyraźnie widoczne są poszczególne ramki w obrazie. Obraz odtworzony na podstawie pięciu składników głównych nie różni się wzrokowo od obrazu oryginalnego. Współczynniki PSNR otrzymane dla poszczególnych obrazów są odpowiednio równe: 18,80 dB, 25,43 dB oraz 27,58 dB, przy czym $\mathbf{PSNR}$ określany jest wzorem

gdzie $MSE$ oznacza wartość błędu średniokwadratowego odtworzonego obrazu względem obrazu oryginalnego.

8.3.2 Przykład zastosowania PCA do ilustracji rozkładu danych wielowymiarowych

Oryginalne zastosowanie znalazło PCA w ilustracji graficznej rozkładu danych wielowymiarowych poprzez zrzutowanie ich w przestrzeń 2-D lub 3-D. Wybierając $K=2$ rzutujemy każdy $N$-wymiarowy wektor $\mathbf{x}$ w przestrzeń dwuwymiarową. W ten sposób każdy wektor $\mathbf{x}$ jest reprezentowany przez wektor $\mathbf{y} = [y_1 y_2]$, którego położenie może być bez problemu zilustrowane na płaszczyźnie, której oś poziomą stanowi teraz składnik główny $y_1$, a oś pionową składnik główny $y_2$. Wektory $\mathbf{x}$ „podobne” do siebie zajmą wówczas bliskie sobie położenia na płaszczyźnie, a wektory „dalekie” – położenia odległe.

Ta unikalna własność znajduje szerokie zastosowania w problemach klasyfikacyjnych, gdzie służy do badania jednorodności rozkładów danych w ramach poszczególnych klas oraz do określania średnich odległości między klasowych

Przykład takiego zastosowania pokażemy w przestrzeni 2-wymiarowej ilustrującej graficznie położenie poszczególnych województw Polski (na płaszczyźnie reprezentowanej przez 2 najważniejsze składniki główne). Wyniki dotyczą przykładowych danych GUS z jednego roku. Rzutowanie dotyczyło 13-wymiarowych elementów informacji GUS dla każdego województwa. Informacje dotyczą następujących elementów:

1. Procent ludności mieszkających w miastach

2. Odsetek zgonu niemowląt

3. Przyrost naturalny ludności

4. Stopa bezrobocia

5. Wynagrodzenie miesięczne brutto

6. Zasoby mieszkaniowe przypadające na 10000 ludności

7. Liczba osób hospitalizowanych w roku przypadających na 10000 ludności

8. Liczba ciągników rolniczych przypadająca na 100 ha gruntów

9. Produkcja sprzedana przemysłu przypadająca na głowę ludności

10. Produkcja sprzedana budownictwa przypadająca na głowę ludności

11. Liczba kilometrów dróg przypadających na 10km2 w przeliczeniu na głowę ludnosci

12. Wartość PKB per capita

13. Liczba firm zarejestrowanych w bazie REGON przypadająca na głowę ludności

Dla uproszczenia opisów przyjęto numeryczne oznaczenia poszczególnych województw w kolejności jak niżej:

1. Dolnośląskie

2. Kujawsko-pomorskie

3. Lubelskie

4. Lubuskie

5. Łódzkie

6. Małopolskie

7. Mazowieckie

8. Opolskie

9. Podkarpackie

10. Podlaskie

11. Pomorskie

12. Śląskie

13. Świętokrzyskie

14. Warmińsko-mazurskie

15. Wielkopolskie

16. Zachodniopomorskie

Tabela 8.1 przedstawia dane oryginalne dotyczące tych zagadnień. Wektor x dla każdego województwa tworzy wybranych 13 elementów informacji dotyczącej ekonomii, dostępności edukacji i opieki medycznej. Użyto następujących skrótów:

Tabela 8.1 Dane liczbowe GUS dotyczące 13 elementów informacji województw w Polsce. Wiersze reprezentują województwa (od 1 do 16), kolumny elementy uwzględnionej informacji (od 1 do 13)

Na rys. 8.6 przedstawiono ilustrację dwuwymiarową położenia względnego 16 województw Polski na podstawie danych statystycznych dotyczących rozwoju ekonomicznego, stanu edukacyjnego i opieki medycznej.

Bliskie położenie na płaszczyźnie oznacza podobieństwo cech charakteryzujących te województwa. Większe odległości miedzy poszczególnymi punktami oznaczają większe różnice w rozwoju tych województw (pod względem uwzględnionych 13 wskaźników ekonomicznych i społecznych). Najbardziej odstającym województwem, wyróżniającym się w grupie analizowanych jednostek okazało się województwo Mazowieckie. Najbardziej odległe od niego są województwa: Lubelskie, Podkarpackie, Podlaskie, Świętokrzyskie i Warmińsko-mazurskie. Najbliższe województwu Mazowieckiemu okazało się województwo Wielkopolskie i Śląskie. Zauważmy, że w tej analizie nie uwzględnia się wartości składników głównych, a jedynie względne odległości między położeniami na płaszczyźnie utworzonej przez dwa najważniejsze składniki główne. Należy dodatkowo zaznaczyć, że powyższy rozkład został uzyskany przy uwzględnieniu jedynie wybranych 13 parametrów, spośród wielu dostępnych w publikacjach GUS i dotyczy danych jednego roku.

1. Korzystając z oprogramowania Matlaba określić macierz korelacji $\mathbf{R}_{\mathbf{xx}}$ dla ciągu 20 wektorów 4-wymiarowych $\mathbf{x}$ utworzonych dla 20 kolejnych chwil czasowych w przedziale czasu $0 < t < 2$. Wektor ten opisują zależności

$\mathbf{x}=\left[\begin{array}{llll}\sin (3 t) & \exp (-0.1 * t) & \sin (3 * t) * \cos (0.5 * t) & \sin (3 * t-\cos (t))\end{array}\right]^T$

Wyznaczyć wartości i wektory własne tej macierzy.

2. Dany jest ciąg 20 wektorów losowych o wymiarze 10 i rozkładzie a) równomiernym, b) gaussowskim. Określić wartości własne i stowarzyszone z nimi wektory własne dla odpowiadających im macierzy korelacji $\mathbf{R}_{\mathbf{xx}}$. Porównać odpowiadające sobie wartości własne i ich rozkład w obu przypadkach przedstawiając je na wspólnym wykresie.

3. Dany jest zbiór wektorów $\mathbf{x}$ o postaci

$\mathbf{x}_1=\left[\begin{array}{r}1 \\ 2 \\ 0\end{array}\right], \; \mathbf{x}_2=\left[\begin{array}{r}0,8 \\ 1,8 \\ 0,5\end{array}\right], \; \mathbf{x}_3=\left[\begin{array}{r}1,3 \\ 2,5 \\ -0,5\end{array}\right], \; \mathbf{x}_4=\left[\begin{array}{r}0 \\ 1 \\ 2\end{array}\right], \; \mathbf{x}_5=\left[\begin{array}{r}0,5 \\ 1,5 \\ 3\end{array}\right], \; \mathbf{x}_6=\left[\begin{array}{r}0,8 \\ 2 \\ 4\end{array}\right]$

Określić macierz korelacji tych wektorów. Wyznaczyć wartości i wektory własne. Określić macierz PCA odwzorowującą te wektory w przestrzeń 2D. Zrekonstruować macierz korelacji na podstawie 2 najważniejszych składników głównych. Zrzutować wektory x_i w przestrzeń 2D. Narysować położenia tych rzutów w przestrzeni 2D.

4. Sygnał pomiarowy x tworzą funkcje czasu: $\sin(\omega_1 t)$, $\cos(\omega_2 t)$, $\exp(-\alpha t ) \cdot \sin(\omega_3 t)$, $\operatorname{sinc}(\omega_4 t)$ oraz dwie składowe szumowe o rozkładzie równomiernym i gaussowskim. Utworzyć zbiór 20 takich wektorów dla 20 różnych chwil czasowych. Każdy wektor x ma więc wymiar 6. Określić wartości własne macierzy korelacji tego zbioru wektorów. Dokonać przekształcenia PCA odpowiadającego 4 największym wartościom własnym macierzy korelacji. Odtworzyć te wektory i wykreślić je (zależność czasowa) na tle wartości oryginalnych.

5. Dla powyższego zbioru 20 wektorów zilustrować ich położenie na płaszczyźnie utworzonej przez dwa składniki główne: $y_1$ oraz $y_2$ ($K=2$).

6. Dana jest transmitancja operatorowa drugiego rzędu $H(s)=\frac{b_0}{s^2+a_1 s+a_0}$. Założyć określone wartości współczynników $a_i$, $b_i$ dla $i = 0, 1$ (na przykład $a_0 = 1$, $a_1 = 0,7$, $b_0 = 1$. Wykreślić charakterystykę amplitudową i fazową układu dla 10 wybranych pulsacji (np. wartości zmieniające się od $\omega=0.1$ do $\omega=2$). Utworzyć 20-wymiarowy wektor $\mathbf{x}$ zawierający wartość modułu $\left|H\left(j \omega_i\right)\right|$ fazy $arg(H(j\omega_1))$ dla przyjętych wartości częstotliwości (każdy punkt częstotliwości odpowiada dwu wartościom: amplitudy i fazy). Utworzyć ciąg takich wektorów dla wartości współczynników $a_0$, $a_1$, $b_0$, $b_1$ zmieniających się od zera do wartości nominalnej (w każdym przypadku zmiana wartości jednego współczynnika). Dokonać przekształcenia PCA dla takiego zbioru wektorów. Narysować wykres zmian położeń tych wektorów na płaszczyźnie utworzonej przez dwa pierwsze składniki główne: $y_1$ i $y_2$.

7. Przedstawić na płaszczyźnie utworzonej przez dwa najważniejsze składniki główne rozkład danych odpowiadających rozwojowi ekonomicznemu województwa wielkopolskiego i podlaskiego w latach 2004-2008 według danych zgromadzonych w pliku pca_zad5.mat.

Słownik opanowanych pojęć

Wykład 8

PCA – transformacja liniowa stosowana do redukcji wymiaru danych (bazuje na transformacji Karhunena-Loewe).

Macierz kowariancji – macierz stanowiąca uogólnienie pojęcia wariancji na przypadek wielowymiarowy

Wartości własne – pierwiastki $\lambda$ równania $det(\lambda \mathbf{1} - \mathbf{R}_{\mathbf{xx}} ) = 0$, gdzie $\mathbf{R}_{\mathbf{xx}}$ jest macierzą kowariancji.

Wektory własne – wektory $\mathbf{w}_i$ stowarzyszone z wartościami własnymi i macierzą $\mathbf{R}_{\mathbf{xx}}$ poprzez relację $\mathbf{R}_{\mathbf{xx}} \mathbf{w}_i = \lambda_i \mathbf{w}_i$.

Składniki główne – elementy wektora wyjściowego $\mathbf{y}$ po transformacji PCA, $\mathbf{y} = \mathbf{Wx}$.

Reguła Oji – uogólnienie reguły Hebba w implementacji on-line transformacji PCA.

Kompresja stratna danych – reprezentacja przybliżona danych poprzez zredukowaną (w stosunku do oryginału) liczbę elementów.

Oryginalne rozwiązanie Heraulta-Juttena dotyczyło problemu separacji sygnałów $s_j (t)$ zmiennych w czasie na podstawie informacji zawartej w ich liniowej superpozycji. Przyjmiemy za ich twórcami, że danych jest n niezależnych (nieznanych) sygnałów $s_j (t)$ oraz macierz mieszająca $\mathbf{A}$ (również nieznana) o wymiarze $n \times n$. Dla pomiarów dostępne są jedynie sygnały $x_i (t)$ będące liniową superpozycją $s_j (t)$ przy czym [8]

dla $i = 1, 2, \ldots, n$. Główna trudność polega na tym, że zarówno $a_{ij} (t)$ jak i $s_j (t)$ nie są znane. Przy założeniu statystycznej niezależności sygnałów J. Herault i C. Jutten zaproponowali rozwiązanie problemu z wykorzystaniem sieci neuronowej [8].

Schemat ogólny włączenia tej sieci w system pomiarowy przedstawiono na rys. 9.1. Istotnym założeniem w ich rozwiązaniu jest statystyczna niezależność sygnałów źródłowych.

Niezależność statystyczna sygnałów losowych jest pojęciem ogólniejszym niż dekorelacja. Pojęcie korelacji dotyczy jedynie pojęcia liniowej zależności występującej między sygnałami, natomiast niezależność statystyczna dowolnej zależności, w tym nieliniowej, istniejącej między sygnałami. W ogólnym przypadku dwie zmienne losowe $y_i$ oraz $y_j$ są statystycznie niezależne, jeśli informacja o jednej zmiennej nie wnosi żadnej wiedzy o zachowaniu drugiej zmiennej. Z matematycznego punktu widzenia niezależność statystyczna oznacza, że dwuwymiarowa łączna gęstość prawdopodobieństwa $p\left(y_i, y_j\right)$ jest równa iloczynowi jednowymiarowych funkcji gęstości zmiennej $y_i$ oraz $y_j$

Dla sygnałów statystycznie niezależnych uogólniona macierz kowariancji funkcji $f(y_i)$ oraz $f(y_j)$ (obie funkcje muszą być nieparzyste) tworzy nieosobliwą macierz diagonalną, mającą postać

W równaniu tym symbol E oznacza wartość oczekiwaną. Z warunku niezależności statystycznej wynika, że wszystkie uogólnione kowariancje wzajemne są zerowe, a zatem $\mathrm{E}\left[\mathbf{f}(\mathbf{y}) \mathbf{g}^T(\mathbf{y})\right]-\mathrm{E}[\mathbf{f}(\mathbf{y})] \cdot \mathrm{E}\left[\mathbf{g}^T(\mathbf{y})\right]=0$, natomiast kowariancje własne są niezerowe  $\mathrm{E}\left[\mathbf{f}(\mathbf{y}) \mathbf{g}^T(\mathbf{y})\right]-\mathrm{E}[\mathbf{f}(\mathbf{y})] \cdot \mathrm{E}\left[\mathbf{g}^T(\mathbf{y})\right] \neq 0$. Warunek statystycznej niezależności sygnałów jest utożsamiony z zerowaniem się kumulantów wzajemnych wyższego rzędu [8].

W rozwiązaniu problemu separacji sygnałów statystycznie niezależnych J. Herault i C. Jutten zaproponowali sieć neuronową liniową ze sprzężeniem zwrotnym, przedstawioną na rys. 9.2.

Sieć zawiera $n$ liniowych neuronów, powiązanych ze sobą przez wzajemne sprzężenia zwrotne. Wagi synaptyczne $w_{ij}$ w rozwiązaniu oryginalnym Heraulta i Juttena są różne od zera tylko przy sprzężeniach wzajemnych (zasada ta została odrzucona przez prof. Cichockiego pozwalając na lepsze działanie systemu). Każdy neuron w sieci generuje sygnał wyjściowy

Przy oznaczeniu przez $\mathbf{A}$ macierzy mieszającej, przez $\mathbf{W}$ macierzy wagowej (obie kwadratowe o tych samych wymiarach)

a przez $\mathbf{x}(t)$, $\mathbf{s}(t)$ oraz $\mathbf{y}(t)$  wektorów odpowiednio obserwowanych sygnałów $x_i(t)$ przetworzonych według zależności (9.1), wektora sygnałów źródłowych $s_i(t)$ oraz wektora sygnałów wyjściowych $y_i(t)$ sieci, gdzie

działanie sieci z rys. 9.2 można opisać równaniem macierzowym

Przy nieznanej macierzy mieszającej $\mathbf{A}$ i wektorze $\mathbf{s}(t)$ oraz założeniu statystycznej niezależności składników wektora $\mathbf{s}(t)$, zadanie sieci sprowadza się do takiego określenia wektora rozwiązania $\mathbf{y}(t)$, które umożliwi odtworzenie sygnałów pierwotnych $s_i(t)$ tworzących wektor $\mathbf{s}(t)$. Z równania (9.7) wynika, że rozwiązanie takie musi spełniać warunek

Wektor $\mathbf{y}(t)$ będzie odtwarzał wektor poszukiwanych sygnałów źródłowych $\mathbf{s}(t)$. Odtworzenie to jest możliwe z dokładnością do pewnej, bliżej nieokreślonej skali $d_i$, czyli $\mathbf{y}(t) = \mathbf{Ds}(t)$ gdzie $\mathbf{D}$ jest macierzą diagonalną, $\mathbf{D}=\operatorname{diag}\left\{d_1, d_2, \ldots, d_n\right\}$, przy praktycznie dowolnej kolejności występowania poszczególnych składników $\mathbf{y}(t)$ w wektorze $\mathbf{y}(t)$ (przy nieznanej postaci składników wektora sygnałów źródłowych nie ma to praktycznie żadnego znaczenia). Stąd podstawowym celem jest rekonstrukcja „kształtu” sygnałów źródłowych. Innymi słowy, poszukiwanym systemem separującym jest taka macierz $\mathbf{W}$, że zachodzi relacja

Można przyjąć, że macierz estymowana $\hat{\mathbf{A}}=\mathbf{W}^{+}$ jest w rzeczywistości określona jako $\hat{\mathbf{A}}=\mathbf{A P D}$ oraz $\mathbf{W A} \hat{\mathbf{A}}=\mathbf{W} \mathbf{W}^{+}=\mathbf{1}_n$, gdzie $\mathbf{P}$ jest macierzą permutacji odpowiedzialnej za przestawienie kolejności składników wynikowych $y_i$, $\mathbf{D}$ – macierz diagonalna skalująca wartości sygnałów wynikowych, $\mathbf{W}^{+}$ jest macierzą pseudo-odwrotną do macierzy $\mathbf{W}$.

Rozwiązanie określające wektor $\mathbf{y}(t)$ spełniający warunek (9.8) jest możliwe do osiągnięcia przy dowolnej liczbie $n$ sygnałów. Przy większej niż dwa liczbie źródeł i nieznanych z góry wartościach współczynników $a_{ij}$ macierzy mieszającej można osiągnąć separację sygnałów, stosując jedynie metody algorytmiczne adaptacyjnego doboru wag sieci neuronowej.

Rozwiązanie problemu separacji sygnałów przy zastosowaniu sieci rekurencyjnej zostało zaproponowane przez J. Heraulta i C. Juttena i ulepszone przez Cichockiego [8] sprowadza się do rozwiązania układu równań różniczkowych, opisujących zmiany wag sieci neuronowej dla każdej chwili czasowej $t$.

Zgodnie z modyfikacja Cichockiego wprowadza się sprzężenie zwrotne własne neuronu z wagą $w_{i i} \neq 0$. Sprzężenie to powoduje samo-normalizację sygnałów wyjściowych, sprowadzając je wszystkie do takiego samego poziomu liczbowego (wartości znormalizowanej) i ułatwiając w ten sposób proces separacji. Wówczas zależności adaptacyjne wag opisane są wzorami (przy założeniu, że $y_i(t)$ nie zawierają składowej stałej)

dla wag łączących różne neurony $i \neq j$, oraz

dla kolejnych wag sprzężenia własnego (wagi diagonalne macierzy $\mathbf{W}$). Obie zależności można zapisać w postaci macierzowej wspólnej dla obu rodzajów wag

w której $\mathbf{W}$ oznacza macierz wag sieci o wymiarze $n \times n$, a $\mathbf{f}(\mathbf{y})$ i $\mathbf{g}(\mathbf{y})$ oznaczają wektory sygnałowe, $\mathbf{f}(\mathbf{y})=\left[f\left(y_1(t)\right), f\left(y_2(t)\right), \ldots, f\left(y_n(t)\right)\right]^T$, $\mathbf{g}(\mathbf{y})=\left[g\left(y_1(t)\right), g\left(y_2(t)\right), \ldots, g\left(y_n(t)\right)\right]^T$.

W praktyce stosuje się różne rodzaje funkcji $f(x)$ i $g(x)$, najczęściej przyjmując jedną z nich typu wypukłego, a drugą typu wklęsłego. Przykładowy wybór funkcji to $f(x)=x^3$, $f(x)=x^5$. W przypadku funkcji $g(x)$ dobre rezultaty uzyskuje się przy $g(x)=tanh(x)$, $g(x)=arctg(x)$, $g(x)=x$, $g(x)=sgn(x)$, itp. Jak zostało udowodnione w pracy [8] obie funkcje $f(x)$ oraz $g(x)$ odpowiadają momentom statystycznym wyższych rzędów, które przy założeniu statystycznej niezależności sygnałów automatycznie zapewniają wartość średnią $\overline{f\left(y_i(t)\right) g\left(y_i(t)\right)}$ równą zeru, co jest warunkiem zbieżności algorytmu uczącego. Wybór współczynnika uczenia $\eta(t)$ odgrywa istotną rolę w procesie adaptacji. W przypadku sygnałów stacjonarnych na początku uczenia jest to wartość stała, a następnie maleje do wartości minimalnej w miarę upływu czasu.

Z rozwiązania układu równań różniczkowych opisanych zależnością (9.12) otrzymuje się aktualne wartości wag, służące do określenia sygnałów wyjściowych opisanych wektorem $\mathbf{y}(t)$, przy czym przy $\mathbf{y}(t)=(\mathbf{1}+\mathbf{W})^{-1} \mathbf{x}(t)$.

Głównym źródłem zwiększonej efektywności algorytmu jest samo-normalizacja (do wartości jednostkowej) sygnałów wyjściowych $y_i(t)$. Mianowicie w stanie ustalonym $\frac{d \mathbf{W}}{d t}=\mathbf{0}$, skąd wynika, że $\overline{f\left(y_i(t)\right) g\left(y_i(t)\right)} = 1$, co oznacza, że niezależnie od aktualnego poziomu sygnałów $s_i(t)$ następuje automatyczne skalowanie wszystkich sygnałów w sieci do poziomu jednostkowego. Badania symulacyjne sieci o zmodyfikowanym algorytmie uczenia wykazały, że możliwa jest separacja sygnałów o amplitudach różniących się nawet w stosunku  $1 : 10^9$.

Program został przetestowany przy separacji wielu różnorodnych sygnałów, wykazując się bardzo dobrą skutecznością.

Na rys. 9.3 przedstawiono ilustracje procesu separacji czterech sygnałów s_i(t) o znacznie różniących się amplitudach

$\begin{aligned} & s_1(t)=0,001 \sin (300 t+6 \cos (60 t)) \\ & s_2(t)=\operatorname{rand}(0,00001, t) \\ & s_3(t)=0,001 \operatorname{sgn}(\cos (155 t)) \\ & s_4(t)=0,00001 \sin (1200 t) \sin (50 t)\end{aligned}$

Rys. 9.3a przedstawia (u góry) sygnały oryginalne $s_i(t)$, rys. 9.3b - sygnały zmieszane $x_i(t)$, a rys. 9.3c (na dole) - sygnały estymowane przez sieć neuronową w procesie separacji (sygnały niezależne będące repliką sygnałów oryginalnych). Sygnałami wejściowymi sieci neuronowej są sygnały zmieszane (trzy środkowe wykresy czasowe na rys. 9.3b). Ze względu na ogromną różnicę amplitud sygnałów oryginalnych w sygnałach zmieszanych podawanych na sieć widoczne są jedynie największe sygnały szumu, natomiast sygnały o małej amplitudzie są niezauważalne.

Proces separacji przeprowadzony został przy zastosowaniu funkcji nieliniowych $f(x) = x^3$ i $f(x) = tgh(10x)$ oraz współczynniku uczenia $\eta(t)$ zmienianym adaptacyjnie przy wartości startowej równej 2000. Taki dobór parametrów pozwolił na separację wszystkich sygnałów, niezależnie od poziomu amplitud (trzy dolne wykresy czasowe na rys. 9.3c). Odseparowane przez sieć sygnały charakteryzowały się jednakowym poziomem amplitud, uzyskanym dzięki wprowadzeniu sprzężenia własnego neuronów. Ich kolejność występowania jest inna niż kolejność sygnałów oryginalnych.

Wadą sieci rekurencyjnej jest potrzeba zapewnienia bezwzględnej stabilności przy separacji sygnałów, szczególnie wtedy, gdy macierz mieszająca $\mathbf{A}$ jest bardzo źle uwarunkowana, a sygnały źródłowe różnią się znacznie pod względem amplitud. Należy ponadto zauważyć, że w sieci rekurencyjnej w każdym kroku występuje konieczność inwersji macierzy wagowej (wzór 9.8), co wydatnie zwiększa złożoność obliczeniową algorytmu. Zastosowanie sieci jednokierunkowej bez sprzężeń zwrotnych pozwala wyeliminować te problemy. Większość współczesnych rozwiązań ślepej separacji unika tych problemów stosując przekształcenie liniowe

odpowiadające sieci liniowej o jednym kierunku przepływu sygnałów. Macierz $\mathbf{W}$ występująca w tej zależności jest macierzą pełną (nie występują założenia co do wartości zerowych pewnych wag). Sieć jednokierunkowa odpowiadające tej zależności pełni identyczną funkcję jak sieć sprzężeniem zwrotnym. Przyjmując, że wektor wyjściowy $\mathbf{y}$ w obu rodzajach sieci jest równy sobie, można zapisać 

W równaniu tym $\mathbf{W}_r$ oznacza macierz wag sieci rekurencyjnej, a $\mathbf{W}$ macierz sieci jednokierunkowej. Oznacza to, że obie sieci będą sobie równoważne, jeśli

W praktyce odchodzi się od równań różniczkowych na rzecz równania różnicowego, jako powszechnie stosowanej metody rozwiązania równań różniczkowych. Powstało wiele odmian algorytmów uczących implementujących relacje (9.13), charakteryzujących się szczególnie dobrymi własnościami separacji przy złym uwarunkowaniu macierzy mieszającej $\mathbf{A}$ lub przy dużym zróżnicowaniu amplitud sygnałów źródłowych s_i(t). Można je w ogólności zaliczyć do algorytmów bazujących na statystykach drugiego rzędu (ang. Second Order Statistics - SOS) oraz na statystykach wyższych rzędów (ang. Higher Order statistics - HOS) [8]. Większość z tych algorytmów rozwinęła się na bazie teorii statystycznego przetwarzania sygnałów i nie wykorzystuje bezpośrednio zależności odnoszących się do sieci neuronowych.

Jednym z najprostszych algorytmów klasy SOS jest AMUSE (ang. Algorithm for Multiple Unknown Source Extraction). Algorytm ten bazuje na diagonalizacji macierzy kowariancji (9.3). Wykorzystuje dwie podstawowe operacje: wybielanie (ang. prewhitening) oraz diagonalizację macierzy kowariancji dla jednego opóźnienia czasowego. Można go traktować jako dwukrotne zastosowanie metody składników głównych PCA [8].

• W pierwszym kroku estymowana jest macierz kowariancji oryginalnych wektorów $\mathbf{x}$ bez opóźnień

• Następnie obliczana jej dekompozycja macierzy kowariancji według wartości własnych (ang. Eigen Value Decomposition - EVD),

w której $\mathbf{V}$ jest macierzą ortogonalną, a $\mathbf{L}$ macierzą diagonalną.

• Następnie obliczana jest macierz $\mathbf{Q}$ liniowej transformacji wybielającej

Wektor sygnałów obserwowanych $\mathbf{x}(k)$ w wyniku wybielenia jest transformowany w wektor sygnałów wybielonych $\mathbf{z}(k)$

• Dla sygnałów wybielonych (zdekorelowanych) $\mathbf{z}(k)$ obliczana jest nowa macierz kowariancji, tym razem z opóźnieniem czasowym równym jeden, z jednoczesną dekompozycją SVD

gdzie $\mathbf{S}$  jest macierzą diagonalną wartości osobliwych macierzy $\mathbf{R}_\mathbf{zz}(1)$ a macierze $\mathbf{U}$, $\mathbf{V}$ są ortogonalne.

• Macierz separującą $\mathbf{W}$ stanowiącą inwersję macierzy mieszającej $\mathbf{A}$ ($\mathbf{W} = \mathbf{A}^{-1}$) określa wówczas wzór

Nieznana macierz mieszająca $\mathbf{A}$ może więc być estymowana w postaci odwrotnej do $\mathbf{W}$, czyli ($\mathbf{A} = \mathbf{Q}^{T}\mathbf{U}$). Algorytm AMUSE jest zwykle dość wrażliwy na szum zakłócający pomiary. Dla zwiększenia jego odporności można w drugim kroku zamiast macierzy kowariancji dla jednego opóźnienia czasowego użyć liniowej kombinację wielu macierzy kowariancji dla różnych opóźnień czasowych. Można zastosować równoczesną diagonalizację kilkunastu lub nawet kilkudziesięciu macierzy kowariancji. W takim przypadku macierz kowariancji $\mathbf{R}_{zz}$ jest wartością oczekiwaną (średnią) z macierzy kowariancji obliczonych dla wielu ($K$) realizacji.

dla $p=1,2,…,K$. Taki sposób postępowania został zaimplementowany miedzy innymi w algorytmie SOBI (ang. Second Order Blind Identification).

W przypadku algorytmów bazujących na statystykach wyższych rzędów wykorzystuje się niezależność statystyczną sygnałów (ICA). We współczesnych rozwiązaniach algorytmów uczących macierz $\mathbf{W}$ otrzymuje się z iteracji na iterację w postaci dyskretnej

unikając w ten sposób potrzeby rozwiązania układu równań różniczkowych. Istnieje wiele różnych rozwiązań ślepej separacji należących do rodziny ICA. Spośród nich można wymienić [8]:

• algorytm Cichockiego-Amari-Younga

gdzie  $\mathbf{L}$ jest macierzą diagonalną o elementach $\lambda_{ii} \geq  0$ (najczęściej $\lambda_{ii} = 1$).

• algorytm opierający się na gradiencie naturalnym prof. Amari

Macierz $\mathbf{L}$ dobiera się podobnie jak w algorytmie Cichockiego-Amari-Younga.

• algorytm Cardossa

gdzie $\alpha$ i $\beta$ są współczynnikami liczbowymi z zakresu $[0, 1]$.

Podobnie jak w sieciach rekurencyjnych, współczynnik uczenia $\eta(k)$ stanowi funkcję malejącą w czasie do zera. Zwykle jest to funkcja wykładnicza postaci $\eta(t)=A e^{-\alpha k / \tau}$, o wartości amplitudy $A$ i stałej czasowej $\tau$ dobieranej indywidualnie dla poszczególnych przypadków.

Przykład zastosowania programu ICALAB dotyczy separacji wypowiedzi czterech wyrażeń: Politechnika Warszawska, Wojskowa Akademia Techniczna, Sieci neuronowe oraz Ala ma kota. Sygnały te zostały uzupełnione przez sygnał losowy. Wszystkie 5 sygnałów zostały zmieszane przy pomocy macierzy losowej A. W następnym etapie tak zmieszane wektory podlegały ślepej separacji przy użyciu jednego z wybranych algorytmów typu SOS (AMUSE) [74]. Na rys. 9.5 przedstawiono wyniki wygenerowane przez program przedstawiające: sygnały oryginalne przed zmieszaniem, sygnały po zmieszaniu, sygnały odtworzone oraz macierz przedstawiającą iloczyn A*W (przy idealnej separacji powinna to być postać macierzowa o wartości 1 na jednej z 5 pozycji w każdym kanale odtworzonym). Wynik działania jest zbliżony do idealnego (przesłuchy z innych kanałów są minimalne, co odzwierciedlone jest na niediagonalnych pozycjach macierzy $\mathbf{AW}$ w postaci bardzo małych wartości, największa wartość niediagonalna to 0.1130). Odtworzenie słuchowe sygnałów wynikowych potwierdza bardzo dobre działanie systemu separacji.

1. Wykorzystując funkcję xcorr Matlaba wyznaczyć funkcje autokorelacji i korelacji wzajemnej dwu losowych wektorów 100 elementowych o rozkładzie gaussowskim (wartości znormalizowane). Narysować przebieg tej funkcji dla różnych opóźnień.

2. Wykorzystując funkcję xcorr Matlaba wyznaczyć funkcje autokorelacji i korelacji wzajemnej dwu wektorów 100 elementowych z których jeden składa się z elementów funkcji sinusoidalnej $x_1(t)=\sin(2t)$ a drugi funkcji wykładniczej $x_2(t)=\exp(-0.1t)$. Narysować przebieg tej funkcji dla różnych opóźnień.

3. Wykorzystując program ICALAB dokonać separacji dwu sygnałów: prostokątnego (np. $A_1 \cdot \textrm{sign}( \sin( \omega_1 t))$ oraz sinusoidalnego (np. $A_2 \cdot  \sin( \omega_2 t)$ dla wybranych wartości $\omega_1$ oraz $\omega_2$ i zróżnicowanych wartości  $A_1$ i $A_2$, np. $A_1 = 1$, $A_2 = 100$ oraz $A_1 = 1$, $A_2 = 10^8$.

4. Sprawdzić działanie sieci separującej (program ICALAB) dla 5 sygnałów wygenerowanych w sposób dowolny. Rozważyć przypadek sygnałów zależnych i niezależnych oraz sygnałów będących superpozycją sygnału deterministycznego i szumu losowego.

5. Pokazać związek występujący między macierzą wag $\mathbf{W}$ sieci separującej rekurencyjnej i jednokierunkowej.

6. Korzystając z funkcji audiorecorder Matlaba dokonać nagrania 3 różnych wypowiedzi tworząc 3 wektory o tej samej długości. Utworzyć z nich macierz X dodając czwarty sygnał szumu losowego (funkcja randn Matlaba) o dużej wariancji (amplituda szumu wielokrotnie większa od amplitudy sygnałów mowy). Wykorzystując program ICALAB sygnały te poddać zmieszaniu przy założeniu różnej macierzy mieszającej $\mathbf{A}$, a następnie separacji. Odtworzyć dźwiękowo sygnały zmieszane oraz sygnały odseparowane używając funkcji soundsc.

7. Sprawdzić i porównać działanie różnych algorytmów ślepej separacji sygnałów zaimplementowanych w programie ICALAB wykorzystując sygnały wbudowane w bazie danych tego pakietu.

Słownik opanowanych pojęć

Wykład 9

ICA – ślepa separacja sygnałów niezależnych (ang. Independent Component Analysis).

BSS – uogólnienie ICA na inne rodzaje dekompozycji sygnałów (BSS - ang. Blind Signal Separation) .

Niezależność statystyczna sygnałów – własność sygnałów oznaczająca, że dwuwymiarowa łączna gęstość prawdopodobieństwa $p(y_i,y_j)$ jest równa iloczynowi jednowymiarowych funkcji gęstości zmiennej $y_i$ oraz $y_j$.

SOS – statystyki drugiego rzędu (ang. Second Order Statistics).

HOS – statystyki wyższych rzędów (ang. Higher Order Statistics).

EVD – dekompozycja macierzy według wartości własnych (ang. Eigen Value Decomposition).

AMUSE – jeden z algorytmów dekompozycji BSS wykorzystujący SOS.

Toolbox ICALAB – program komputerowy z interfejsem graficznym w Matlabie do ślepej separacji sygnałów.

Preprocessing – wstępne przetwarzanie sygnałów.

Postprocessing – końcowe przetwarzanie sygnałów po właściwej operacji.