Algorytmy projekcyjne najmniejszych kwadratów

Algorytmy projekcyjne najmniejszych kwadratów wykorzystywane do estyma-cji nieujemnych faktorów zarówno w modelu NMF, jak i NTF s¡ projekcyjnymi wersjami algorytmu ALS (ang. Alternating Least Squares). Jest to najstarszy i za-razem najprostszy algorytm naprzemiennej estymacji zmiennych wedªug kryte-rium najmniejszych kwadratów. Koncepcja zastosowania algorytmu ALS do esty-macji faktorów w dekompozycji PARAFAC pojawiªa si¦ w pracy Harshmana [173], a tak»e niezale»nie w pracy Carrolla i Changa [53] jako model zwany CANDE-COMP. Nast¦pnie, idea ta byªa rozwijana i stosowana w wielu obszarach bada« naukowych [104, 235, 240, 241, 340, 408].

Podstawowa wersja algorytmu ALS nie gwarantuje znalezienia faktorów o ujemnych elementach. Aby jednak algorytm ten byª u»yteczny do estymacji nie-ujemnych faktorów, konieczne s¡ usprawnienia, wymuszaj¡ce nieujemne rozwi¡-zania. Najprostsz¡ i do±¢ popularn¡ technik¡ wymuszania nieujemno±ci jest rzu-towanie estymowanych faktorów na nieujemny ortant, czyli zast¡pienie ujemnych elementów w estymowanym faktorze warto±ciami zerowymi. Rzutowanie to dla elementu ξ ∈ R realizowane jest przeksztaªceniem:

[ξ]₊= max{0, ξ}. (4.137)

W praktycznych zastosowaniach najcz¦±ciej wykonuje si¦ rzutowanie na dodatni ortant, czyli [ξ]_ϵ= max{ϵ, ξ}, gdzie ϵ > 0 i ξ ∈ R. Takie podej±cie stosowane jest

we wszystkich algorytmach projekcyjnych najmniejszych kwadratów, omawianych w rozdziale 4.5, a tak»e w wielu pracach na temat nieujemnej faktoryzacji macie-rzy, np. w [18, 19, 89, 194, 395].

Podstawowe algorytmy projekcyjne ALS pomimo swej prostoty w implemen-tacji i niskiego kosztu obliczeniowego maj¡ bardzo istotn¡ wad¦, a mianowicie nie

gwarantuj¡ monotoniczno±ci procesu iteracyjnego ani te» zbie»no±ci do punktu stacjonarnego, okre±lonego warunkami optymalno±ci KKT. Pomimo tych wad, ba-dania eksperymentalne [81, 489] pokazuj¡, »e dla specycznych zada« algorytmy te zachowuj¡ monotoniczn¡ zbie»no±¢. Nie mo»na wi¦c ich wykluczy¢ z obszaru zastosowania modelu NMF lub NTF.

4.5.1. Algorytm ALS

Algorytm ALS w wersji podstawowej minimalizuje funkcj¦ odlegªo±ci eukli-desowej (2.4), naprzemiennie ze wzgl¦du na argumenty A i X. Z warunków stacjonarno±ci funkcji celu w (3.38), gdzie D(Y ||AX) jest funkcj¡ w (2.4) oraz

α_A= α_X = 0, wynika:

∇AΨ (A, X) = 2(AX− Y )XT , 0, (4.138)

∇XΨ (A, X) = 2A^T(AX− Y ) , 0. (4.139)

Zakªadaj¡c rank(A) = rank(X) = J dla J ≤ min{I, T }, z równa« (4.138) i (4.139) otrzymano algorytm ALS:

A = Y X^T(XX^T)−1_{= Y X}†_{, X = (A}TA)−1_ATY = A†_{Y ,} (4.140) gdzie X†i A† s¡ pseudoodwrotno±ciami Moore'aPenrose'a odpowiednich macie-rzy X i A. Po zastosowaniu przeksztaªcenia (4.137), uzyskuje si¦ naprzemienny algorytm projekcyjny najmniejszych kwadratów2:

A = [ Y X†] +, X = [ A†_Y] +. (4.141)

Zªo»ono±¢ obliczeniow¡ algorytmu (4.140) dla aktualizacji faktora A lub X mo»na zgrubnie oszacowa¢ jako O(J3+ J²(I + T ) + IJ T ), przyjmuj¡c »e koszt obliczania macierzy odwrotnych w (4.140) wynosi O(J3).

Algorytm postaci (4.141) jednak rzadko stosowany jest w praktyce. Pomimo zaªo»e« rank(A∗_{) = rank(X}∗_{) = J} oraz rank(A(0)) = rank(X⁽⁰⁾) = J dla

J ≤ min{I, T }, naprzemienna estymacja faktorów mo»e spowodowa¢, »e w n-tym

kroku iteracyjnym rankϵ(A⁽ⁿ⁾) < J lub rankϵ(X⁽ⁿ⁾) < J. Funkcja rankϵ(·)

ozna-cza numeryczny rz¡d macierzy, szacowany jako liczba warto±ci osobliwych, które s¡ wi¦ksze od przyj¦tej tolerancji ϵ. Niech rankϵ(A⁽ⁿ⁾) < J w n-tym kroku itera-cyjnym. W takim przypadku estymacja faktora X(n) odbywa si¦ z utrat¡ cz¦±ci 2 W dalszej cz¦±ci pracy naprzemienny algorytm projekcyjny najmniejszych kwadratów jest oznaczany skrótowcem ALS, zgodnie z powszechnie przyj¦t¡ nomenklatur¡.

informacji widmowej, co mo»e przekªada¢ si¦ na gorszy wynik estymacji faktora

A⁽ⁿ⁺¹⁾. W praktyce lepsze efekty uzyskuje si¦, je±li proces estymacji faktorów podlega kontroli, np. przez odpowiedni¡ regularyzacj¦.

4.5.2. Algorytm RALS

Zregularyzowane algorytmy ALS (RALS  ang. Regularized ALS) minimalizuj¡ funkcj¦ celu w (3.38), gdzie D(Y ||AX) jest funkcj¡ w (2.4), a czªony regularyzu-j¡ce s¡ lub mog¡ by¢ modelowane funkcjami kwadratowymi. Niech Ψ(A, X) jest postaci (3.62) dla αA> 0 i αX > 0, wówczas

∇AΨ (A, X) = 2(AX− Y )XT + 2αAA, 0, (4.142)

∇XΨ (A, X) = 2A^T(AX− Y ) + 2αXX , 0. (4.143)

Po przeksztaªceniu równa« (4.142) i (4.143) oraz uwzgl¦dnieniu rzutowania (4.137) otrzymuje si¦ projekcyjny algorytm RALS:

A =[

Y X^T(XX^T + αAIJ)−1]

+, X =[

(A^TA + αXIJ)−1_ATY]

+. (4.144) Je±li αA> 0 i αX > 0, macierze odwrotne w (4.144) istniej¡ niezale»nie od rz¦du macierzy A i X.

Gdyby zaªo»y¢, »e rank(A) = J oraz faktor X jest rzadki, a informacja a priori o rzadko±ci tego faktora modelowana jest czªonem UX(X) w (3.46), wówczas reguªa aktualizacji faktora X miaªaby posta¢:

X =[

(A^TA + αX1J×J)−1_ATY]

+, (4.145)

gdzie αX > 0 oraz 1J×J ∈ RJ×J

+ jest macierz¡ samych jedynek. Je±li rzad-ko±¢ wymuszana jest w faktorze A, rank(X) = J i αA > 0, to po uwzgl¦d-nieniu funkcji (3.47) reguªa aktualizacji dla macierzy A ma posta¢: A = [

Y X^T(XX^T + αA1J)−1]

+.

Zgodnie z modelem (3.50) wymuszanie rzadko±ci powinno odbywa¢ si¦ tylko w jednym z estymowanych faktorów. Je±li faktor X aktualizowany jest wedªug (4.145), to faktor A powinien by¢ estymowany wedªug reguªy (4.144) lub innej wymuszaj¡cej gªadko±¢ estymowanego faktora.

Stosuj¡c algorytm RALS nale»y tak»e mie¢ na uwadze wªa±ciwy dobór parame-trów αAi αX. Korzystaj¡c z reguªy aktualizacji w (4.144) oraz je±li rank(A) < J lub rank(X) < J, zadania estymowane s¡ ¹le uwarunkowane, a parametry αA

i αX maj¡ charakter parametrów regularyzacji. W takim przypadku mo»na zasto-sowa¢ typowe narz¦dzia do estymacji parametrów regularyzacji Tichonowa. Mog¡ to by¢ metody krzywej L (ang. L-curve) lub metody uogólnionej walidacji krzy-»owej (ang. Generalized Cross-Validation  GCV) [30, 172]. Je±li jednak faktory

A i X s¡ peªnego rz¦du i dobrze uwarunkowane, to wspomniane metody esty-macji parametrów regularyzacji nie s¡ u»yteczne. W takim przypadku zadaniem czªonów kary w funkcji celu (3.38) nie jest regularyzacja ¹le postawionych lub ¹le uwarunkowanych zada«, a wymuszanie okre±lonych cech estymowanych faktorów. Parametry te powinny by¢ tak dobierane, aby uzyska¢ zadawalaj¡cy charakter estymowanych faktorów. Niestety, tak sformuªowane zadanie jest bardzo trudne do realizacji w praktycznych zastosowaniach. Zwykle nie mo»na z góry oszaco-wa¢ jaki jest stopie« rzadko±ci lub gªadko±ci estymowanego faktora. W praktyce, najcz¦±ciej parametry te dobierane s¡ na podstawie symulacji.

Mo»liwe jest te» alternatywne rozwi¡zanie, polegaj¡ce na adaptacyjnym lub zaplanowanym wspomaganiu procesu aktualizacji faktorów dan¡ informacj¡ a priori. Tak wi¦c, parametry kary αAlub αX mog¡ zmienia¢ si¦ z krokami itera-cji naprzemiennych. Sposób ich zmian mo»e by¢ okre±lony przez narzucon¡ reguª¦ deterministyczn¡ lub poprzez reguª¦ adaptacyjn¡, zale»n¡ od oceny pewnych cech. Wyniki bada« eksperymentalnych przedstawione w wielu pracach [81, 89, 484 486, 497, 502, 511] pokazuj¡, »e parametry kar powinny przyjmowa¢ du»e warto±ci na pocz¡tku iteracyjnego procesu aktualizacji faktorów, a nast¦pnie warto±ci tych wspóªczynników powinny asymptotycznie male¢ w funkcji liczby iteracji do pewnej maªej warto±ci progowej. Podej±cie to b¦dzie analizowane i wyja±niane w dalszej cz¦±ci pracy. Przyj¦to wi¦c nast¦puj¡ce reguªy:

• eksponencjalna:

α⁽ⁿ⁾= ¯α + α⁽⁰⁾exp{−τn}, (4.146)

• bisekcyjna:

α⁽ⁿ⁾= ¯α + α⁽⁰⁾2−n_, (4.147)

gdzie α(n)  parametr αAlub αX w n-tym kroku iteracyjnym, α(0) > 0  pocz¡t-kowa warto±¢ parametru, 0 < τ ≤ 1  szybko±¢ zmian parametru, ¯α > 0  warto±¢ progowa parametru.

Omówione algorytmy RALS nie wyczerpuj¡ w peªni tej tematyki. Obszerny przegl¡d ró»nych wersji algorytmu RALS mo»na znale¹¢ w pracach [81, 89].