Dowód 4.2. Z warunku stacjonarno±ci funkcji G(x t , x (n) t ) wynika:
4.7. Gradientowe algorytmy kierunków poprawy
Formuª¦ iteracyjn¡ gradientowych algorytmów kierunków poprawy dla aktu-alizacji wektora xt ∈ RJ w k-tym kroku iteracji wewn¦trznych mo»na wyrazi¢ nast¦puj¡co:
x(k+1)t = x(k)t + η(k)t p(k)t , (4.167) gdzie p(k)
t ∈ RJ jest kierunkiem poprawy, okre±laj¡cym kierunek poszukiwa« mini-mum funkcji celu, a η(k)
t > 0jest dªugo±ci¡ kroku poszukiwa« lub wspóªczynnikiem uczenia.
Analogiczn¡ formuª¦ iteracyjn¡ mo»na zapisa¢ dla aktualizacji i-tego wektora wierszowego macierzy A. Dla tego typu algorytmów, k-ty krok iteracyjny jest rozumiany jako krok iteracji wewn¦trznych, natomiast n-ty krok iteracyjny jako krok iteracji zewn¦trznych (naprzemiennych). W ka»dym n-tym kroku iteracji naprzemiennych zakªada si¦, »e aktualizacja zarówno faktora A, jak i faktora X, wymaga wykonania okre±lonej liczby (zwykle wi¦kszej od jeden) kroków iteracji wewn¦trznych.
Aby sekwencje iteracyjne generowane przez (4.167) byªy nieujemne, trzeba okre±li¢ strategi¦ wymuszania nieujemno±ci. Typowym i najprostszym rozwi¡za-niem jest zastosowanie rzutowania (4.137) do (4.167), co w efekcie prowadzi do formuªy iteracyjnej algorytmu rzutowania gradientu dla ogranicze« przedziaªo-wych [153]: x(k+1)t = [ x(k)t + ηt(k)p(k)t ] +. (4.168)
Niestety takie podej±cie nie gwarantuje zbie»no±ci procesu iteracyjnego do punktu stacjonarnego okre±lonego warunkami KKT . Niemniej jednak reguª¦ ak-tualizacji (4.168) mo»na zastosowa¢ do modelu NMF lub NTF, je±li poª¡czy si¦ j¡ z odpowiednim kryterium zatrzymania, np. z kryterium gradientów rzutowanych w (4.42).
Nieujemno±¢ iteracyjnych aktualizacji mo»na te» wymusza¢ poprzez odpo-wiedni dobór dªugo±ci kroku η(k)
t . Je±li aktualizowany faktor X jest rzadki, mo»na zaªo»y¢: ∀t : minjx(k+1)jt = 0. W rezultacie z reguªy (4.167) otrzymuje si¦:
ηt(k)= min j { −x (k) jt p(k)jt dla p (k) jt < 0 } , (4.169) gdzie p(k)
t = [p(k)jt ]∈ RJ. W przeciwnym razie wspóªczynnik η(k)
t mo»na otrzyma¢ z zaªo»enia ∃t : minjx(k+1)jt = 0.
Wektor kierunku poprawy p(k)
t mo»e by¢ de niowany na ró»ne sposoby. Wspóln¡ cech¡ gradientowych metod kierunków poprawy jest zaªo»enie, »e wek-tor ten jest wyznaczany na podstawie gradientu funkcji celu w pewnym punkcie niekoniecznie punkcie x(k)
t . W wielu metodach omawianych w niniejszej pracy wektor kierunku poprawy mo»na zapisa¢ nast¦puj¡co:
p(k)t =−B(k)
t g(k)t , (4.170)
gdzie g(k)
t = ∇xtΨ (A, ¯x(k)t ) ∈ RJ gradient funkcji Ψ(A, xt) w punkcie ¯x(k)
t ,
B(k)t ∈ RJ×J macierz skaluj¡ca, która mo»e przybiera¢ ró»ne formy, w zale»no±ci od przyj¦tej metody.
Punkt ¯x(k)
t mo»e by¢ wybierany jako punkt x(k)
t (metoda gradientów prostych) albo mo»e by¢ de niowany przez funkcj¦ okre±lon¡ na zbiorze punktów {x(k)
t }
(metoda gradientów optymalnych). Je±li B(k)
t jest macierz¡ dodatnio okre±lon¡, to (g(k)
t )Tp(k)t =−(g(k)
t )TB(k)t g(k)t < 0. Sugeruje to, »e przesuwaj¡c si¦ z punktu
x(k)t w kierunku wektora p(k)
t , mo»liwe jest zmniejszenie warto±ci funkcji celu, tzn. Ψ(A, x(k+1)
t ) < Ψ (A, x(k)t )dla dowolnego x(0)
t ∈ RJ. Je±li B(k)
t jest macierz¡ nieujemnie okre±lon¡, to Ψ(A, x(k+1)
t )≤ Ψ(A, x(k)
t ).
W niniejszym rozdziale przedstawiono wybrane algorytmy NMF, które nale»¡ do rodziny gradientowych algorytmów kierunków poprawy. Formuªy iteracyjne bazuj¡ wprost na reguªach (4.167) lub (4.168). Nale»y jednak zauwa»y¢, »e algo-rytmy multiplikatywne przedstawione w rozdziale 4.4, algoalgo-rytmy projekcyjne naj-mniejszych kwadratów (rozdz. 4.5) oraz algorytmy blokowo-sekwencyjne (rozdz. 4.6) mo»na równie» w jakim± sensie zakwali kowa¢ do rodziny gradientowych algorytmów kierunków poprawy. Zgodnie z reguª¡ aktualizacji (4.50) algorytmy multiplikatywne wynikaj¡ bezpo±rednio z formuªy iteracyjnej algorytmu gradien-towego. Mo»na ªatwo pokaza¢, »e algorytmy ALS, omówione w rozdziale 4.5, s¡ szczególnymi przypadkami projekcyjnej metody Newtona, która bazuje na sche-macie (4.168) [89]. Analogicznie, zgodnie z (4.165) metody HALS i SCWA, mo»na wyrazi¢ formuª¡ (4.168). Niemniej jednak wspomniane algorytmy zaliczono do od-dzielnych grup algorytmów, ze wzgl¦du na ich nalne reguªy aktualizacji faktorów.
4.7.1. Projekcyjna metoda Landwebera
Metoda Landwebera [20], nazywana te» metod¡ Richardsona pierwszego ro-dzaju [30], bazuje na regule iteracyjnej (4.167), w której kierunek poprawy wy-ra»ony jest przez (4.170), gdzie B(k)
t = IJ, zatem: p(k)
Ψ (A, xt) = 12D(Y||Axt), gdzie D(Y ||Axt) jest funkcj¡ odlegªo±ci euklidesowej w (2.4), to p(k)
t =−AT(Ax(k)t −yt). Sekwencje iteracyjne {x(k)
t } generowane tak¡
reguª¡ dla dowolnego x(0)
t ∈ RJ s¡ zbie»ne do rozwi¡zania w sensie kryterium naj-mniejszych kwadratów, tzn. limk→∞x(k)t = x∗
t, je±li η(k)
t ∈ (0, ηmax). Aby wyzna-czy¢ ηmax, reguª¦ iteracyjn¡ w (4.167) mo»na przeksztaªci¢ do postaci: x(k+1)
t = x(k)t − η(k) t AT(Ax(k)t − yt) = (IJ− η(k) t ATA)x(k)t + ηt(k)ATyt= G(k)t x(k)t + bt. Se-kwencje {x(k) t } s¡ zbie»ne, je±li ρ(G(k) t ) < 1, gdzie ρ(G(k) t ) = max1≤j≤J|λj(G(k)t )|
jest promieniem spektralnym macierzy G(k)
t , a λj(G(k)t ) jest jej j-t¡ warto±ci¡ wªasn¡. atwo zauwa»y¢, »e je±li ηmax = 2
λmax(ATA), gdzie λmax(ATA) jest maksymaln¡ warto±ci¡ wªasn¡ macierzy ATA, to ρ(G(k)
t ) < 1.
Niestety estymacja maksymalnej dªugo±ci kroku ηmax bezpo±rednio z de nicji wymaga zbyt du»ego kosztu obliczeniowego. Dlatego takie podej±cie nie mo»e by¢ zastosowane do omawianych metod faktoryzacji macierzy lub tensorów. Po-niewa» macierz A jest macierz¡ nieujemn¡, λmax(ATA) mo»e by¢ oszacowana na podstawie twierdzenia Perrona-Frobeniusa [17]. Je±li ATA ∈ RJ×J+ jest ma-cierz¡ nieredukowaln¡3, promie« spektralny macierzy ATA speªnia nierówno±ci: min1≤u≤J∑J
v=1[ATA]uv≤ ρ(ATA)≤ max1≤u≤J∑J
v=1[ATA]uv. Z tego wynika:
ηmax= 2
max1≤u≤J∑J
v=1[ATA]uv. (4.171)
Stosuj¡c ograniczenia nieujemno±ci wedªug reguªy w (4.168), projekcyjny al-gorytm Landwebera mo»na wyrazi¢ nast¦puj¡c¡ formuª¡ iteracyjn¡:
X(k+1)= [
X(k)− ηXAT(AX(k)− Y )]
+, ηX = ξηmax, (4.172) gdzie 0 < ξ . 1.
Johansson i inni wspóªautorzy pracy [213] pokazali, »e proces iteracyjny w (4.172) mo»na dodatkowo przyspieszy¢, je±li macierz B(k)
t w (4.170) zde nio-wana jest jako diagonalna macierz skaluj¡ca, tzn.
B(k)t = D = diag ( [ATA1J]−1 j ) . (4.173)
3 Dowolna macierz kwadratowa jest macierz¡ nieredukowaln¡, wtedy i tylko wtedy, gdy nie mo»na j¡ sprowadzi¢ do macierzy blokowej górnotrójk¡tnej przez jednoczesne permutacje jej wierszy i kolumn.
Dla tak zde niowanej macierzy skaluj¡cej, projekcyjny algorytm Landwebera w (4.172) sprowadza si¦ do algorytmu OPL (ang. Oblique Projected Landweber):
X(k+1)= [
X(k)− DAT(AX(k)− Y )]
+, (4.174)
gdzie D jest wyra»one przez (4.173).
atwo zauwa»y¢, »e wspóªczynnik ηmax zawarty jest w de nicji macierzy D. W pracy [213] udowodniono, »e sekwencje iteracyjne {X(k)} generowane przez
(4.174) dla A ∈ RI×J
+ i X(0) ≥ 0 s¡ zbie»ne do rozwi¡zania liniowego
zada-nia najmniejszych kwadratów przy ograniczezada-niach nieujemno±ci. Algorytm OPL zostaª zastosowany po raz pierwszy do modelu NMF w pracy [499].
4.7.2. Algorytm rzutowania gradientu
Typowym przedstawicielem algorytmów rzutowania gradientu w metodzie NMF jest algorytm zaproponowany przez Lina w [280]. Oznaczono go skrótowcem: LPG (ang. Lin's Projected Gradient). Algorytm ten bazuje na formule iteracyjnej, podanej w (4.168), gdzie dªugo±¢ kroku η(k)
t estymowano w sposób przybli»ony na podstawie reguªy Armijo [23]. Reguªa ta okre±la jeden z podstawowych warunków Wolfe'a [329], wymaganych w celu zapewnienia wystarczaj¡cej poprawy wzdªu» kierunku p(k)
t . Dªugo±¢ kroku η(k)
t dobierana jest z ci¡gu liczb: 1, β, β2, β3, . . ., gdzie β ∈ (0, 1). Optymalna dªugo±¢ kroku speªnia warunek: η(k)
t = βmk dla pewnej nieujemnej liczby caªkowitej mk = 0, 1, 2, . . ., dobranej w taki sposób, aby: Ψ (A, x(k+1)t )− Ψ(A, x(k) t )≤ σ(∇xtΨ (A, x(k)t ) )T (x(k+1)t − x(k) t ), (4.175) gdzie σ ∈ (0, 1). W praktyce najcz¦±ciej: σ = 0, 01. Bertsekas [22] udowodniª, »e istnieje takie η(k)
t > 0, dla którego punkt graniczny aproksymacji (4.168) jest punktem stacjonarnym.
Lin [280] zastosowaª do metody NMF równie» zmody kowan¡ wersj¦ reguªy Armijo. Lin oraz Moré [281] zauwa»yli, »e η(k)
t i η(k−1)
t mog¡ mie¢ podobne war-to±ci. Dlatego zamiast rozpoczyna¢ proces selekcji dªugo±ci kroku od warto±ci
ηt(k) = β0 = 1, mo»na zaczyna¢ od warto±ci η(k−1)
t , a nast¦pnie krok ten zwi¦k-sza¢ lub zmniejzwi¦k-sza¢ w zale»no±ci od warunku (4.175). Je±li w kroku pocz¡tkowym speªniony jest warunek (4.175), η(k)
t nale»y zwi¦ksza¢. W przeciwnym wypadku, parameter η(k)
t jest zmniejszany, zgodnie z pierwotn¡ zasad¡ Armijo. Takie podej-±cie znacz¡co obni»a koszt obliczeniowy procesu aktualizacji.
4.7.3. Algorytm gradientów skalowanych
Algorytm gradientów skalowanych bazuje na algorytmie gradientowych punk-tów wewn¦trznych (IPG ang. Interior-Point Gradient), który ukazaª si¦ w pracy Merritta i Zhanga [302]. Przeznaczony jest do rozwi¡zywania liniowych zada« najmniejszych kwadratów z ograniczeniami nieujemno±ci. Do metody NMF za-stosowano go w pracach [85, 89, 499]. Algorytm ten realizuje aktualizacje addy-tywne wedªug reguªy (4.167), gdzie kierunek poprawy p(k)
t okre±lony jest przez (4.170). Macierz B(k)
t zde niowana jest przez diagonaln¡ macierz skaluj¡c¡ st¡d algorytm skalowanych gradientów. Przyj¦to wi¦c: B(k)
t = diag { x(k)jt [AT AX(k) ]jt } . W rezultacie kierunek poprawy ma posta¢:
P(k)=− ˜D(k)~ ∇XΨ (A, X(k)), (4.176)
gdzie ˜D(k)= X(k)⊘ ATAX(k) macierz skaluj¡ca. Je±li ∀t : η(k)
t = 1w (4.167), algorytm ten sprowadza si¦ do algorytmu MUE, przedstawionego w rozdziale 4.4.1. W algorytmie IPG [302] faktor X aktualizo-wany jest wedªug reguªy:
X(k+1) = X(k)+ η(k)P(k). (4.177) Dªugo±¢ kroku η(k)dobierana jest w taki sposób, aby aktualizacje byªy nieujemne i je±li to mo»liwe, zgodne z reguª¡ najszybszego spadku (ang. steepest descent). Tak wi¦c η(k) = min{τ ˆη(k), ¯η(k)}, gdzie ˆη(k) = arg maxη
{
X(k)+ ηP(k)≥ 0},
τ . 1 oraz ¯η(k)= arg minηΨ (A, X(k)+ ηP(k)). Krok ˆη(k) zapewnia nieujemno±¢ aktualizacji i mo»e by¢ dobierany na podstawie reguªy (4.169), zatem: ˆη(k) = minj,t { −x(k)jt p(k)jt dla p(k)jt < 0 } . Je±li Ψ(A, X) = 1
2D(Y||AX), gdzie D(Y ||AX)
jest funkcj¡ odlegªo±ci euklidesowej, otrzymuje si¦:
Ψ (A, X(k)+ ηP(k)) = 1 2||Y − AX(k)+ ηAP(k)||2 F = 1 2η 2tr ( (P(k))TATAP(k) ) + η tr ( (P(k))TG(k) ) + const, (4.178)
gdzie G(k)=∇XΨ (A, X(k)) = AT(AX(k)−Y ). Z warunku ∂
∂ηΨ (A, X(k+1)), 0 oraz z (4.178) wynika:
¯
η(k)=− (vec(P(k)))Tvec(G(k))
Dªugo±ci kroków ˆη(k) i ¯η(k) mog¡ by¢ te» de niowane, oddzielnie dla ka»dego t i wówczas faktor X aktualizowany jest wedªug reguªy (4.167), gdzie η(k)
t =
min{τ ˆη(k)
t , ¯ηt(k)}. Parametr ˆη(k)
t estymowany jest na podstawie (4.169), natomiast:
¯ ηt(k)=− (p (k) t )Tg(k)t (Ap(k)t )TAp(k)t , (4.180) gdzie g(k)
t = ∇xtΨ (A, x(k)t ). Takie podej±cie zapewnia lepsz¡ jako±¢ estymacji, ale koszt obliczeniowy jest nieco wi¦kszy.
4.7.4. Algorytm gradientów optymalnych
Algorytm gradientów optymalnych opracowano na podstawie metody Neste-rowa [326], przeznaczonej do rozwi¡zywania zada« optymalizacji nieliniowej bez ogranicze« dziedziny. Jest to metoda gradientowa, której wspóªczynnik zbie»no-±ci wynosi O(1/k2), je±li minimalizowana funkcja f(·) ∈ C1,1
L : RM → R jest
wypukª¡ funkcj¡ klasy C1,1
L (z ci¡gªym gradientem i jednostajnie ci¡gª¡ z wa-runkiem Lipschitza). atwo zauwa»y¢, »e funkcja odlegªo±ci euklidesowej w (2.4) jest funkcj¡ klasy C1,1
L , zarówno ze wzgl¦du na argument A, jak i X. Niech
Ψ (A, X) = 12D(Y||AX), gdzie D(Y ||AX) wyra»one jest poprzez (2.4), zatem: ||∇AΨ (A, X)− ∇AΨ ( ¯A(k), X)||F ≤ LA||A − ¯A(k)||F, (4.181)
||∇XΨ (A, X)− ∇XΨ (A, ¯X(k))||F ≤ LX||X − ¯X(k)||F, (4.182) gdzie {A, ¯A(k)} ∈ RI×J, {X, ¯X(k)} ∈ RJ×T. Czynniki LA = ||XXT||2 i LX =
||ATA||2 s¡ staªymi Lipschitza.
Z nierówno±ci (4.181) i (4.182) wynika: Ψ (A, X) ≤ Ψ( ¯A(k), X) + ⟨ ∇AΨ ( ¯A(k), X), A− ¯A(k) ⟩ + L˜A 2 ||A − ¯A(k)||2 F = FA(A, ¯A(k)), (4.183) Ψ (A, X) ≤ Ψ(A, ¯X(k)) + ⟨ ∇XΨ (A, ¯X(k)), X− ¯X(k) ⟩ + L˜X 2 ||X − ¯X(k)||2 F = FX(X, ¯X(k)). (4.184) je±li LA≤ ˜LAi LX ≤ ˜LX.
W celu minimalizacji funkcji Ψ(A, X) ze wzgl¦du na X przyj¦to, »e
w de nicji 4.1. Jest to funkcja wypukªa i okre±lona na zbiorze RJ×T
+ dla LX ≤ ˜LX. Stosuj¡c metod¦ gradientów proksymalnych (ang. proximal gradient) [11, 92, 326], uzyskuje si¦:
X(k+1)= proxh( ¯X(k)) = arg min
X
(
h(X) + FX(X, ¯X(k)) )
, (4.185)
gdzie proxh( ¯X(k)) operator proksymalny pewnej wypukªej funkcji h(X). Dla
h(xjt) = {
0, je»eli xjt∈ R+,
∞, w przeciwnym razie (4.186)
reguªa aktualizacji faktora X jest nast¦puj¡ca:
X(k+1) = [ ¯ X(k)− ˜L−1X GX( ¯X(k)) ] +, (4.187) gdzie GX( ¯X(k)) =∇XΨ (A, ¯X(k)) = AT(A ¯X(k)− Y ).
Punkt ¯X(k) jest punktem linearyzacji funkcji pomocniczej FX(X, ¯X(k)). Mo»na go zde niowa¢ na ró»ne sposoby. Je±li ¯X(k)= X(k), to reguªa aktualizacji w (4.187) mo»e by¢ interpretowana jako reguªa iteracyjna projekcyjnej metody Landwebera. W metodzie Nesterowa ¯X(k) jest wyznaczany jako punkt ekstrapo-lacji kierunków {X(k−1), X(k)}. St¡d:
¯
X(k)= X(k)+ β(k)(X(k)− X(k−1)). (4.188) Stosuj¡c podej±cie Nesterowa [326], optymalny wspóªczynnik zbie»no±ci se-kwencji iteracyjnej w (4.187) uzyskuje si¦, gdy czynnik β(k) wyra»ony jest przez
β(k) = γ(k−1)−1
γ(k) . Parameter γ(k) obliczany jest z równania (γ(k))2 − γ(k) −
(γ(k−1))2 = 0.
Finalna wersja algorytmu gradientów optymalnych dla aktualizacji faktora
X wyra»ona jest algorytmem 11, gdzie kmax oznacza maksymaln¡ liczb¦ iteracji wewn¦trznych. Algorytm aktualizacji faktora A jest analogiczny do algorytmu 11. Twierdzenie 4.2. Je±li sekwencje iteracyjne {X(k)} i {Z(k)} generowane s¡
al-gorytmem 11 dla k ≥ 1, oraz X∗= limk→∞X(k), to
Ψ (A, X(k))− Ψ(A, X∗)≤ 2LX||X(k)− X∗||2
F
(k + 2)2 . (4.189)
Algorytm 11. Iteracje Nesterowa
Wej±cie: A ∈ RI×J, Y ∈ RI×T, X(0)∈ RJ×T
+ , kmax maksymalna liczba iteracji wewn¦trznych
Wyj±cie: X estymowany faktor
1 Inicjalizacja: Z(0) = X(0), LX =||ATA||2, γ(0)= 1;
2 for k = 1, 2, . . . , kmax do
3 G(k)X = AT(AZ(k−1)− Y ) ; // Gradient w punkcie Z(k)
4 X(k)= [ Z(k−1)− L−1 X G(k)X ] + ; // Aktualizacje projekcyjne 5 γ(k)= 1+ √ 4(γ(k−1))2−1 2 , β(k)= γ(k−1)−1 γ(k) ; Z(k)= X(k)+ β(k)(X(k)− X(k−1)) ; // Kierunek poprawy Maksymalna liczba iteracji kmax w algorytmie 11 mo»e by¢ okre±lana przez odpowiednie kryterium zatrzymania, np. kryterium gradientów rzutowanych w (4.42). Mo»liwe te» s¡ inne podej±cia. Przykªadowo, w pracy [504] liczb¦ t¦ dobierano wedªug zasady: kmax = min{10, n}, gdzie n jest krokiem iteracji
ze-wn¦trznych. To podej±cie motywowane jest reguªami (4.146) i (4.147). Silniejszy wpªyw parametrów regularyzacji w pocz¡tkowej fazie naprzemiennej aktualizacji faktorów jest równowa»ny w kontek±cie regularyzacji wcze±niejszemu przerywaniu iteracji wewn¦trznych. Jest to tzw. regularyzacja przez wcze±niejsze przerywanie procesu iteracyjnego (ang. truncated iterations).
Metoda Nesterowa wykorzystywana jest w ró»nych dziedzinach przetwarzania sygnaªów i obrazów [11, 92, 529], a tak»e w grupowaniu dokumentów teksto-wych [162]. W pracy [504] pokazano, »e metoda ta jest równie» bardzo efektywna w klasy kacji obrazów tekstur. Badano j¡ równie» w kontek±cie ekstrakcji wekto-rów cech z obrazów bazy ORL i z obrazów hiperspektralnych [462]. Nowe algo-rytmy NMF bazuj¡ce na metodzie Nesterowa mo»na te» odnale¹¢ w pracy [159].