Wyniki testów numerycznych

X^(k)− G(k) X diag{(¯α(k) t )−1}^] +− X(k) ; // Kierunek poprawy

4 η¯^(k)= max{0, min{1, η(k)}} ; // gdzie η(k) dany jest w (4.271)

5 X^(k+1)= X^(k)+ P^(k)diag{¯η(k)};

6 α¯^(k+1) = max{αmin, min{αmax, α^(k+1)}} ; // gdzie α(k+1) wyra»ony jest przez (4.268)

Je±li αX = 0, z warunku stacjonarno±ci: ∂

∂η(k)Ψ (A, X^(k+1)), 0 oraz z (4.266) otrzymano reguª¦ aktualizacji dla η(k) w postaci analitycznej:

η^(k)=− ¹ T J [ P^(k)~ G(k) X ] 1^T_J [ P^(k)~ (ATAP^(k)) ]. (4.271)

Finalna posta¢ zmodykowanego algorytmu SPG wyra»ona jest przez algo-rytm 20.

Uwaga: Zªo»ono±¢ obliczeniow¡ pojedynczego kroku iteracyjnego algorytmu SPG dla aktualizacji macierzy X mo»na zgrubnie oszacowa¢ funkcj¡ O(IJT )+O(J2T ). Uwzgl¦dniaj¡c czªon regularyzuj¡cy oraz zakªadaj¡c, »e wyra»enie X(k)L_X wymaga kosztu obliczeniowego O(JT2), caªkowita zªo»ono±¢ zregularyzowanego algorytmu SPG wynosi O(IJT ) + O(J2T ) + O(J T²). Je±li macierz LX jest bardzo rzadka, zªo»ono±¢ t¦ mo»na jeszcze zmniejszy¢.

4.11. Wyniki testów numerycznych

Algorytmy numeryczne omawiane w niniejszym rozdziale byªy obszernie ba-dane i testowane w ró»nych pracach, z wykorzystaniem ró»nych zbiorów danych. W tym rozdziale przedstawiono jedynie wybrane wyniki bada« symulacyjnych,

które zostaªy wykonane na danych syntetycznych, pochodz¡cych gªównie z zada« ±lepej separacji sygnaªów nieujemnych. W badaniach tych macierz obserwacji Y generowana jest na podstawie znanych macierzy A i X, które nazwano macierzami oryginalnymi. Dla modelu dokªadnej faktoryzacji, macierz Y wyznaczona jest przez zale»no±¢ (1.1). Dane dla modelu przybli»onej faktoryzacji generowane s¡ wedªug modelu (1.3). Przyj¦to, »e N = [n1, . . . , nT], gdzie ∀t : nt∼ N (0, σ2

tII). Wariancje szumu σ2

t dobierano w taki sposób, aby uzyska¢ zamierzon¡ warto±¢ wspóªczynnika SNR.

Estymowane macierze ˆA i ˆX oceniano jako±ciowo wedªug miary SIR (ang. Signal-to-Interference Ratio) [89]: SIR^(A) = ¹ J J ∑ j=1

SIR^(A)_j , gdzie SIR^(A)_j =−20 log

( ||ˆaj− aj||2 ||aj||2 ) , (4.272) SIR^(X)= ¹ J J ∑ j=1

SIR^(X)_j , gdzie SIR^(X)_j =−20 log(||ˆxj− xj||2

||xj||2

)

, (4.273)

gdzie ˆaj  j-ty wektor kolumnowy macierzy ˆA, ˆxj  j-ty wektor wierszowy ma-cierzy ˆX.

4.11.1. Testy metod inicjalizacji faktorów

Metody inicjalizacji faktorów modelu NMF omówiono w rozdziale 4.1. W celu oszacowania ich efektywno±ci przeprowadzono eksperyment estymacji faktorów oryginalnych A i X w zadaniu ±lepej separacji sygnaªów nieujemnych. Do realiza-cji tego zadania wykorzystano dokªadny oraz zaburzony model NMF, odpowiednio (1.1) oraz (1.3). Macierze oryginalne A = [aij] ∈ RI×J

+ oraz X = [xjt] ∈ RJ×T

+

wygenerowano losowo w nast¦puj¡cy sposób:

∀i, j : aij ∼ U[0, 1], (4.274)

∀j, t : xjt = max{0, ξjt}, gdzie ξjt ∼ N (0, 1). (4.275) Przyj¦to: I = 100, J = 10 oraz T = 1000.

Do oceny przybli»e« pocz¡tkowych A(0) oraz X(0) wybrano nast¦puj¡ce me-tody: inicjalizacja losowa (rozdz. 4.1.1), wielostartowa inicjalizacja losowa (algo-rytm 3) (K = N = 10), inicjalizacja sferyczn¡ metod¡ k-±rednich (algo(algo-rytm 4),

(a) (b) 1 2 3 4 5 6 12 14 16 18 20 22 24 SIR[dB] Algorytmy 1 2 3 4 5 6 10 15 20 25 30 SIR[dB] Algorytmy

Rys. 4.4. Statystyki rozkªadu wspóªczynników SIR macierzy ˆAestymowanych algorytmem HALS. Przybli»enia pocz¡tkowe A(0)i X(0)uzyskano za pomoc¡ nast¦puj¡cych metod inicjalizacji faktorów: 1  zwykªa losowa, 2  wielostartowa inicjalizacja losowa, 3  inicjalizacja

sferyczn¡ metod¡ k-±rednich, 4  NNDSVD, 5  inicjalizacja bidiagonalizacj¡ Lanczosa, 6  metoda CHV. Statystyki pokazano dla: (a) SNR = 20 dB, (b) SNR = 40 dB. Dla danych bez zaburze« szumowych metoda CHV estymuje macierz A ze wspóªczynnikiem SIR > 200 dB metoda NNDSVD [34] (algorytm 6), inicjalizacja bidiagonalizacj¡ Lanczosa (al-gorytm 8) oraz metoda CHV (al(al-gorytm 9) z funkcj¡ D wyra»on¡ przez kwadrat odlegªo±ci euklidesowej.

Ocen¦ efektywno±ci badanych metod przeprowadzono na podstawie wniosko-wania statystycznego z danej liczby prób. Wygenerowano 100 ró»nych macierzy

Y, gdzie macierze A i X otrzymywane byªy wedªug odpowiednich zale»no±ci (4.274) i (4.275). Macierze A(0) i X(0) dla ka»dej próby estymowane byªy za pomoc¡ badanych metod inicjalizacji faktorów. Nast¦pnie przeprowadzono fak-toryzacj¦ ka»dej macierzy Y algorytmem HALS. Dla ka»dej próby wykonywano 50 iteracji naprzemiennych. Estymowane macierze ˆA i ˆX oceniano wedªug miar SIR w (4.272) i (4.273). Statystyki rozkªadów warto±ci wspóªczynników SIR uzy-skanych dla próbek ˆA przedstawiono na rysunku 4.4. W eksperymencie tym sto-sowano dane zaburzone zakªóceniami szumowymi o ró»nej mocy. W tabeli 4.1 zamieszczono wyniki estymacji macierzy A(0) uzyskane tylko metod¡ CHV (bez ª¡czenia z jakimkolwiek algorytmem NMF).

Wniosek 4.2. Wyniki przedstawione w tabeli 4.1 pokazuj¡, »e dla danych fakto-ryzowalnych i niezakªóconych, metoda CHV dla p = 1 znajduje macierz A(0), dla której SIR > 200 dB. W takim przypadku nie ma potrzeby stosowania jakichkol-wiek algorytmów NMF. Niestety, nawet dla niezbyt silnych zakªóce« szumowych

Tabela 4.1. ‘rednie warto±ci wspóªczynnika SIR dla estymacji macierzy A(0), uzyskane algorytmem CHV dla 100 próbek (bez algorytmu NMF). Parametr p deniuje liczb¦

najbli»szych s¡siadów. W nawiasach podano odchylenia standardowe Dane dokªadne SIR = 40 dB SIR = 30 dB SIR = 20 dB p = 1 204,5 (49,95) 21,66 (2,61) 11,98 (2,27) 3,37 (0,92) p = 3 24,98 (18,07) 10,13 (2,74) 7,45 (1,92) 5,02 (0,92) p = 10 8,66 (0,51) 8,17 (0,5) 7,4 (0,44) 6,2 (0,43)

(SNR = 30 dB), wspóªczynnik SIR dla macierzy A(0) znacz¡co maleje i dlatego aproksymacj¦ A(0) mo»na jedynie traktowa¢ jako przybli»enie pocz¡tkowe dla metody NMF.

Na rysunku 4.4 pokazano, »e dla faktoryzacji przybli»onej, dla której SNR ≥ 40 dB oraz wystarczaj¡co rzadkiej macierzy X (wedªug def. 3.8), metoda CHV prowadzi do najlepszej jako±ci estymacji macierzy A(0). Jednak jako±¢ ta znacz¡co zale»y od mocy zakªóce« modelu faktoryzacji. Je±li SNR ∼= 40 dB oraz p = 1, macierz A(0) estymowana jest ze wspóªczynnikiem SIR ≥ 21 dB (tabela 4.1), a po zastosowaniu algorytmu HALS, wspóªczynnik SIR dla estymacji faktora A zwi¦ksza sie do warto±ci ok. 30 dB. Gdy dane obserwowane s¡ silnie zaburzone (SNR ≤ 20 dB), konieczna jest adaptacja parametru p w algorytmie 9 do okre-±lonego poziomu mocy zakªóce«. Ponadto, z rysunku 4.4 wynika, »e dla silnych zakªóce« metoda k-±rednich b¦dzie najlepszym wyborem.

4.11.2. Liniowe zadanie najmniejszych kwadratów

W podanym eksperymencie numerycznym oceniano przebieg zbie»no±ci wy-branych algorytmów iteracyjnych w liniowym zadaniu najmniejszych kwadratów:

min

X ||Y − AX||F, p.o. X ≥ 0. (4.276)

Testy wykonano, stosuj¡c dane wygenerowane z losowych macierzy A ∈ RI×J

+

i X ∈ R^J×T₊ . Zaªo»ono: I = 100, J = 10, T = 1000 oraz ∀i, j : aij ∼ U(0, 1), czyli a_ij > 0. Macierz X ∈ RJ×T

+ byªa równie» generowana z rozkªadu równomiernego. Nast¦pnie ok. 50% jej elementów losowo wybranych zast¡piono zerami, tak aby w ka»dym jej wektorze kolumnowym co najmniej dwa elementy byªy dodatnie.

Do bada« wybrano nast¦puj¡ce algorytmy:

• MUE  algorytm multiplikatywny, wyra»ony reguª¡ (4.68) do aktualizacji ma-cierzy X (rozdz.4.4.1),

• HALS  algorytm blokowo-sekwencyjny (rozdz. 4.6.1) do aktualizacji macie-rzy X, przedstawiony za pomoc¡ algorytmu 10,

• OPL  skalowany i projekcyjny algorytm Landwebera (rozdz. 4.7.1), wyra-»ony reguª¡ iteracyjn¡ w (4.174),

• LPG  algorytm rzutowania gradientu (rozdz. 4.7.2), zaproponowany przez Lina w [280],

• IPG  algorytm gradientów skalowanych, przedstawiony w rozdziale. 4.7.3, • OG  algorytm gradientów optymalnych, wyra»ony przez iteracje Nesterowa

w algorytmie 11 (rozdz. 4.7.4),

• SPG  algorytm spektralnego rzutowania gradientu, wyra»ony przez algorytm 20 w rozdziale 4.10.5,

• IIP  algorytm punktów wewn¦trznych, opisany w rozdziale 4.8.

Proces iteracyjny ka»dego algorytmu inicjalizowano macierz¡ X(0), któr¡ wy-generowano z rozkªadu równomiernego. Macierz ta byªa taka sama dla ka»dego algorytmu. W celu oszacowania ±rednich szybko±ci zbie»no±ci, eksperyment ten powtarzano 100 krotnie, gdzie za ka»dym razem losowo dobierano macierze A, X oraz X(0). ‘rednie warto±ci oraz odchylenia standardowe wska¹ników uwarunko-wania tych macierzy wynosz¡ odpowiednio: cond(A) = 7, 4, σ(cond(A)) = 0, 4, cond(X) = 2, 88, σ(cond(X)) = 0, 053.

Szybko±ci zbie»no±ci procesu iteracyjnego oceniano na podstawie u±rednio-nych przebiegów znormalizowanej odlegªo±ci euklidesowej (bª¦du residualnego):

Ψ (A, X^(k)) = ||Y−AX^(k)||F

||Y||F w funkcji liczby iteracji k oraz u±rednionego czasu wykonywania si¦ algorytmów. Przebiegi te pokazano na rysunku 4.5.

Wniosek 4.3. Na rysunku 4.5(a) pokazano, »e algorytm IIP charakteryzuje si¦ najwi¦ksz¡ szybko±ci¡ zbie»no±ci, znacznie ró»ni¡c¡ si¦ od pozostaªych algoryt-mów. Wniosek ten nie jest zaskakuj¡cy, poniewa» jest to algorytm oparty na ak-tualizacjach wzdªu» kierunku Newtona. W odró»nieniu od prostej, projekcyjnej metody Newtona (rozdz. 4.10), algorytm ten zapewnia monotoniczn¡ zbie»no±¢. Nale»y jednak zauwa»y¢, »e regularyzacja hesjanu i wprowadzenie parametru prze-suni¦cia τ(k)

t > 0w (4.196) ograniczaj¡ bª¡d residualny od doªu, co jest widoczne w postaci stagnacji bª¦du poni»ej warto±ci ok. 10−10. Z drugiej strony, taka stagna-cja jest ªatwiejsza do wykrycia przez kryteria zatrzymania procesów iteracyjnych i nie ma znacz¡cego wpªywu na nalny wynik aproksymacji rozwi¡zania. W tym rankingu na kolejnych miejscach znalazªy si¦ algorytmy: HALS, SPG, OG, LPG, OPL, IPG i MUE. Algorytm MUE okazaª si¦ najwolniejszy, co potwierdza analiz¦ teoretyczn¡ zamieszczon¡ w rozdziale 4.4.1.