Shlukování a snížení dimenze: k-means, PCA a t-SNE

Nechť d ≥ 1 a n ≥ 1 jsou přirozená čísla, x₁, …, x_n ∈ ℝ^d body a k ∈ ℕ, 1 ≤ k ≤ n, počet shluků. Přiřazení je zobrazení c : {1, …, n} → {1, …, k}, které bodu x_i přiřadí číslo jeho shluku. Pro středy μ₁, …, μ_k ∈ ℝ^d a přiřazení c je inercie

J(μ₁, …, μ_k, c) = ∑_i=1ⁿ ‖x_i − μ_c(i)‖²,

kde ‖x − y‖ = √∑_r=1^d (x_r − y_r)² je euklidovská vzdálenost bodů x, y ∈ ℝ^d. Úloha k-means je najít středy a přiřazení, pro které je J nejmenší možné. (Díl Fuzzy shlukování zapisuje tutéž úlohu s počtem shluků c, středy v_j a přiřazením „shluk“; tady držíme písmeno k z názvu metody.)

Lloydův algoritmus (přibližné řešení úlohy k-means). Vstup: body x₁, …, x_n ∈ ℝ^d a počet shluků k. Výstup: středy μ₁, …, μ_k a přiřazení c. Kroky:

1. Za počáteční středy zvol k navzájem různých bodů z x₁, …, x_n (náhodně).
2. Přiřazení: pro každé i = 1, …, n zvol c(i) ∈ argmin_{j ∈ {1, …, k}} ‖x_i − μ_j‖²; je-li nejbližších středů víc a je mezi nimi dosavadní střed bodu, ponech bod v jeho shluku.
3. Zastavení: pokud se přiřazení oproti předchozímu průchodu nezměnilo, skonči a vrať středy a přiřazení.
4. Středy: pro každé j = 1, …, k označ S_j = {i : c(i) = j} a je-li S_j ≠ ∅, polož μ_j = 1/|S_j| ∑_{i ∈ Sj} x_i (průměr bodů shluku; prázdný shluk střed nemění).
5. Pokračuj krokem 2.

Algoritmus tedy střídá dva kroky — každý bod přiřadí nejbližšímu středu a každý střed přesune do průměru svých bodů —, dokud se přiřazení mění.

Nechť x₁, …, x_n ∈ ℝ^d, 1 ≤ k ≤ n a Lloydův algoritmus běží podle definice výše (včetně pravidla, že při shodné vzdálenosti zůstane bod ve svém shluku). Pak

1. žádný krok přiřazení ani krok středů nezvětší inercii J;
2. algoritmus skončí po konečně mnoha průchodech.

Důkaz. Krok přiřazení: člen ‖x_i − μ_c(i)‖² se nahradí nejmenším ze všech ‖x_i − μ_j‖², takže se nezvětší.

Krok středů: pro neprázdnou množinu indexů S jednoho shluku s průměrem x̄ = 1/|S| ∑_i∈S x_i a libovolné μ ∈ ℝ^d platí

∑_i∈S ‖x_i − μ‖² = ∑_i∈S ‖x_i − x̄‖² + |S| · ‖x̄ − μ‖²

(smíšený člen 2 ∑ (x_i − x̄) · (x̄ − μ) je nulový, protože ∑ (x_i − x̄) = 0). Druhý sčítanec je ≥ 0 a nulový právě pro μ = x̄, takže průměr je nejlepší možný střed. Prázdný shluk do J nepřispívá a jeho střed se nemění.

Přiřazení c : {1, …, n} → {1, …, k} je konečně mnoho (kⁿ) a při každé změně přiřazení inercie ostře klesne (při shodné vzdálenosti bod zůstává ve svém shluku), takže se žádné nezopakuje a algoritmus skončí. ∎ Věta neříká, že konec je globální minimum — jen že algoritmus nebude cyklit. Podrobnější verzi s fuzzy příslušnostmi má díl Fuzzy shlukování.

Body na přímce 1, 2, 4, 7, 9, 10, k = 2, start μ₁ = 1, μ₂ = 2 (záměrně špatný — oba středy vlevo). Řádek tabulky = středy, přiřazení k nim a inercie (barvy podle vzorce výše: středy, sčítance ‖x_i − μ_c(i)‖² a výsledné J):

Středy jsou průměry předchozího řádku: 2 + 4 + 7 + 9 + 10/5 = 6,4; pak 1 + 2/2 = 1,5 a 4 + 7 + 9 + 10/4 = 7,5; nakonec 1 + 2 + 4/3 = 7/3 ≈ 2,33 a 7 + 9 + 10/3 = 26/3 ≈ 8,67. Ve 3. kroku je součet čtverců v každém shluku 42/9. Bod 4 přeběhl ve 2. kroku, protože |4 − 1,5| = 2,5 < |4 − 7,5| = 3,5. Ve 3. kroku se přiřazení nezměnilo — konec. Inercie klesala 142 → 26,84 → 15,5 → 9,33, jak věta slibuje.

Nechť x₁, …, x_n ∈ ℝ^d jsou body, mezi kterými je aspoň k navzájem různých, a k ∈ ℕ, 1 ≤ k ≤ n. Jsou-li už vybrány středy μ₁, …, μ_j−1, označme pro bod x ∈ ℝ^d

D(x) = min_{l ∈ {1, …, j−1}} ‖x − μ_l‖

jeho vzdálenost k nejbližšímu už vybranému středu. k-means++ (Arthur a Vassilvitskii, 2007) je tento způsob volby počátečních středů. Vstup: body x₁, …, x_n a k. Výstup: středy μ₁, …, μ_k. Kroky:

1. μ₁ = x_i, kde index i ∈ {1, …, n} je vybrán náhodně, každý s pravděpodobností 1/n.
2. Pro j = 2, …, k: μ_j = x_i, kde index i je vybrán náhodně s pravděpodobností P(x_i) = D(x_i)²/∑_l=1ⁿ D(x_l)² (jmenovatel je kladný, protože dokud je vybráno méně než k středů, leží aspoň jeden z různých bodů mimo ně).
3. Středy μ₁, …, μ_k použij místo kroku 1 Lloydova algoritmu a pak ho spusť obvyklým způsobem.

Každý další střed se tedy volí náhodně, ale body daleko od všech dosavadních středů mají větší šanci.

Ať první střed padl na bod 1. Pak D(x)² pro ostatní body je 1, 9, 36, 64, 81 a součet 191 (barvy podle vzorce výše). Bod 10 bude druhým středem s pravděpodobností P(10) = 81/191 = 42,4 %, bod 2 jen s P(2) = 1/191 = 0,5 %. Že druhý střed padne do pravé skupiny {7, 9, 10}, má pravděpodobnost 36 + 64 + 81/191 = 94,8 %. Náhodný start (každý ze zbylých pěti bodů se stejnou šancí) by tam trefil jen 3/5 = 60 %.

Metoda lokte. Nechť pro každé k = 1, 2, …, k_max je J^*(k) nejmenší inercie, kterou jsme pro k shluků našli (třeba nejlepší z několika běhů k-means++). Metoda lokte vynese body (k, J^*(k)) do grafu a zvolí to k, za kterým se pokles J^*(k) − J^*(k + 1) výrazně zmenší (graf tam má „loket“). Je to vizuální pravidlo, ne přesné kritérium — co je „výrazně“, matematicky určeno není.

Silueta. Nechť x₁, …, x_n ∈ ℝ^d jsou navzájem různé body, k ≥ 2 a C₁, …, C_k ⊆ {1, …, n} neprázdné, navzájem disjunktní množiny indexů s C₁ ∪ … ∪ C_k = {1, …, n} (shluky). Pro bod i ∈ C_j definujeme:

1. soudržnost (pro |C_j| ≥ 2) a(i) = 1/|C_j| − 1 ∑_{l ∈ Cj, l ≠ i} ‖x_i − x_l‖, průměrnou vzdálenost k ostatním bodům vlastního shluku;
2. oddělení b(i) = min_{h ≠ j} 1/|C_h| ∑_{l ∈ Ch} ‖x_i − x_l‖, průměrnou vzdálenost k bodům nejbližšího cizího shluku;
3. siluetu s(i) = b(i) − a(i)/max(a(i), b(i)) pro |C_j| ≥ 2,   s(i) = 0 pro |C_j| = 1 (jmenovatel je kladný, protože body jsou různé, a tedy b(i) > 0);
4. průměrnou siluetu rozkladu s̄ = 1/n ∑_i=1ⁿ s(i).

Vždy −1 ≤ s(i) ≤ 1. Silueta bodu je tedy kladná, právě když má bod k vlastnímu shluku v průměru blíž než k nejbližšímu cizímu. Ze všech zkoušených k se volí to s největší průměrnou siluetou s̄.

Bod 4 (barvy podle vzorce výše): a = |4 − 1| + |4 − 2|/2 = 2,5; b = 3 + 5 + 6/3 ≈ 4,67; s = 4,67 − 2,5/4,67 ≈ 0,46. Leží na kraji svého shluku, a tak má siluetu nejmenší ze všech šesti bodů.

Bod 1: a = 1 + 3/2 = 2; b = 6 + 8 + 9/3 ≈ 7,67; s = 7,67 − 2/7,67 ≈ 0,74. Průměrná silueta všech šesti bodů vyjde 0,66.

Nechť X ⊆ ℝ^d je konečná množina bodů, ε > 0 reálné číslo (poloměr) a m ∈ ℕ, m ≥ 1 (minimální počet bodů). Pro x, y ∈ X:

1. ε-okolí bodu x je N_ε(x) = {y ∈ X : ‖x − y‖ ≤ ε} (obsahuje i x samotný);
2. x je jádrový, právě když |N_ε(x)| ≥ m;
3. y je z x přímo dosažitelný, právě když x je jádrový a y ∈ N_ε(x);
4. y je z x dosažitelný, právě když existují r ≥ 2 a body z₁, …, z_r ∈ X takové, že z₁ = x, z_r = y a z_l+1 je přímo dosažitelný ze z_l pro každé l = 1, …, r − 1;
5. x a y jsou propojené, právě když existuje o ∈ X, ze kterého je dosažitelný x i y.

Shluk (Ester a kol., 1996) je neprázdná množina C ⊆ X, která splňuje

1. maximalitu: ∀ x, y ∈ X : (x ∈ C ∧ y je dosažitelný z x) ⇒ y ∈ C;
2. souvislost: ∀ x, y ∈ C : x a y jsou propojené.

Body z X, které nepatří do žádného shluku, jsou šum; nejádrový bod, který do nějakého shluku patří (leží v ε-okolí jádrového bodu), je okrajový. DBSCAN je algoritmus níže, který tyto shluky najde. Shluk tedy tvoří jádrové body, mezi kterými se dá přejít po krocích délky nejvýš ε přes další jádrové body, a k tomu všechny body v jejich okolí. Okrajový bod dosažitelný ze dvou různých shluků patří podle definice do obou; algoritmus ho přiřadí tomu, který ho najde dřív.

Nechť n ≥ 1 je počet bodů, G₁, …, G_q ⊆ {1, …, n} skutečné skupiny (navzájem disjunktní, se sjednocením {1, …, n}) a C₁, …, C_p ⊆ {1, …, n} nalezené shluky (navzájem disjunktní; indexy mimo všechny shluky jsou šum). Shoda nalezených shluků se skupinami je číslo

shoda = 1/n ∑_l=1^p max_{h ∈ {1, …, q}} |C_l ∩ G_h| ∈ [0, 1].

Je to podíl bodů, které leží ve shluku spolu s většinou své skupiny — každý shluk „připíšeme“ skupině, které má v sobě nejvíc, a ostatní jeho body počítáme jako chybu; body označené za šum se také počítají jako chyba. Míra má smysl, jen když skutečné skupiny známe (tady víme, který bod je z kterého měsíce). Pro k-means vyšlo 221 z 300 bodů, tedy 221/300 ≈ 73,7 %. U dvou stejně velkých skupin dá i zcela náhodné rozdělení kolem 50 % — tohle je výchozí laťka, ne nula.

Prokletí dimenze (Richard Bellman, 1957) není jedna matematická věta, ale souhrnný název pro jevy, kvůli kterým metody fungující v malém počtu rozměrů d selhávají ve velkém. Pro body v jednotkové krychli [0, 1]^d jsou to hlavně tyto tři:

1. objem roste exponenciálně: pokrýt krychli mřížkou s krokem h ∈ (0, 1), 1/h ∈ ℕ, vyžaduje (1/h)^d bodů — pro h = 0,1 a d = 10 deset miliard;
2. data jsou řídká: krychle o hraně r ∈ (0, 1) zabírá podíl r^d objemu, takže okolí, ve kterém leží 10 % rovnoměrně rozložených bodů, musí mít hranu r s r^d = 0,1 — pro d = 10 je to r ≈ 0,79, „okolí“ tedy zabírá skoro celý rozsah každé souřadnice;
3. vzdálenosti se koncentrují: pro dva nezávislé rovnoměrně rozložené body x, y ∈ [0, 1]^d klesá podíl směrodatné odchylky a střední hodnoty ‖x − y‖² k nule jako 1/√d (výpočet níže), takže všechny vzdálenosti jsou si relativně čím dál podobnější.

Rozdíl t = x_r − y_r dvou nezávislých rovnoměrných čísel z [0, 1] má trojúhelníkovou hustotu 1 − |t| na [−1, 1]. Proto E[t²] = 2∫₀¹ t²(1 − t) dt = 2(1/3 − 1/4) = 1/6 a E[t⁴] = 2(1/5 − 1/6) = 1/15. Rozptyl jednoho sčítance t² je E[t⁴] − E[t²]² = 1/15 − 1/36 = 7/180. Střední hodnoty i rozptyly nezávislých sčítanců se sčítají, takže pro d souřadnic dostaneme střední hodnotu d/6 a rozptyl 7d/180, tedy směrodatnou odchylku √7d/180. Poměr (barvy podle vzorce výše):

√7d/180 · 6/d = √36 · 7/180 d = √1,4/d.

Tohle je náš vlastní výpočet pro rovnoměrná data; stejný jev (koncentraci vzdáleností) pro mnohem obecnější rozdělení dokázali Beyer a kol. (1999) v práci „When is ‚nearest neighbor‘ meaningful?“.

Nechť n ≥ 2 a D ≥ 1 jsou přirozená čísla, x₁, …, x_n ∈ ℝ^D body s průměrem x̄ = 1/n ∑_i=1ⁿ x_i a Σ ∈ ℝ^D×D jejich kovarianční matice (definice níže). Matice Σ je symetrická a pozitivně semidefinitní, má proto vlastní čísla λ₁ ≥ λ₂ ≥ … ≥ λ_D ≥ 0 a k nim vlastní vektory w₁, …, w_D ∈ ℝ^D, které jsou ortonormální: Σw_r = λ_rw_r, ‖w_r‖ = 1 a w_r · w_s = 0 pro r ≠ s. Pak

1. r-tá hlavní komponenta (r = 1, …, D) je směr w_r (při navzájem různých vlastních číslech je určen jednoznačně až na znaménko);
2. snížení dimenze na q (q ∈ {1, …, D}) je zobrazení x ↦ (w₁ · (x − x̄), …, w_q · (x − x̄)) ∈ ℝ^q, tj. promítnutí vycentrovaného bodu na prvních q komponent;
3. podíl vysvětleného rozptylu prvními q komponentami (pro Σ ≠ 0) je λ₁ + … + λ_q/λ₁ + … + λ_D.

Vlastní číslo λ_r je přitom rozptyl dat ve směru w_r (věta níže), takže první komponenta je směr, ve kterém se data liší nejvíc.

Nechť n ≥ 2, D ≥ 1 a x₁, …, x_n ∈ ℝ^D jsou body (sloupcové vektory, x_ia je a-tá souřadnice bodu x_i) s průměrem x̄ = 1/n ∑_i=1ⁿ x_i. Jejich (výběrová) kovarianční matice je matice Σ ∈ ℝ^D×D,

Σ = 1/n − 1 ∑_i=1ⁿ (x_i − x̄)(x_i − x̄)^T,   tj.   Σ_ab = 1/n − 1 ∑_i=1ⁿ (x_ia − x̄_a)(x_ib − x̄_b)   pro a, b = 1, …, D.

Prvek Σ_ab je tedy (skoro) průměrný součin odchylek a-té a b-té souřadnice od jejich průměrů. Na diagonále jsou rozptyly jednotlivých souřadnic, mimo ni kovariance: kladná, když obě souřadnice bývají současně nad průměrem (nebo současně pod ním). Dělí se n − 1, ne n, protože průměr byl odhadnut ze stejných dat (Besselova korekce).

Nechť D ≥ 1 a Σ ∈ ℝ^D×D je symetrická matice s největším vlastním číslem λ₁. Pak funkce w ↦ w^TΣw na množině {w ∈ ℝ^D : ‖w‖ = 1} nabývá maxima, toto maximum je rovno λ₁ a nabývá se ve vlastním vektoru Σ k vlastnímu číslu λ₁.

Důkaz. Maximum existuje, protože spojitá funkce na jednotkové sféře (uzavřené a omezené množině) maxima nabývá. Hledáme ho Lagrangeovou funkcí L = w^TΣw − λ(w^Tw − 1). Gradient podle w je 2Σw − 2λw, a ten je nulový, právě když Σw = λw — w je vlastní vektor. Dosazením: w^TΣw = λw^Tw = λ, takže rozptyl ve směru vlastního vektoru je jeho vlastní číslo a maximum dává vlastní vektor s největším λ. ∎ Protože je Σ symetrická, jsou její vlastní vektory navzájem kolmé a součet všech vlastních čísel je stopa Σ, tedy součet rozptylů všech souřadnic — odtud jmenovatel ve vzorci pro vysvětlený rozptyl.

Průměr 150 lidí je 171,9 cm a 71,5 kg; kovarianční matice (v cm², cm · kg a kg²) vyšla

Σ = ( 63,84   54,95 ;   54,95   92,31 ).

Vlastní čísla matice 2 × 2 jsou kořeny λ² − (stopa) λ + det = 0. Stopa = 63,84 + 92,31 = 156,14, det = 63,84 · 92,31 − 54,95² ≈ 2 873,3, a tak λ = 78,07 ± √78,07² − 2 873,3 = 78,07 ± 56,76, tedy λ₁ = 134,83 a λ₂ = 21,31.

Vlastní vektor k λ₁: z prvního řádku (63,84 − λ₁)w₁ + 54,95 w₂ = 0 vyhovuje (54,95; 70,99); vydělením délkou 89,77 dostaneme w₁ = (0,61; 0,79) a w₁^TΣw₁ = λ₁ = 134,83 (barvy podle vzorce výše). Vysvětlený rozptyl 134,83/156,14 = 86,4 %. Směrodatná odchylka podél první komponenty je √134,83 ≈ 11,6, podél druhé √21,31 ≈ 4,6 (smíšené jednotky, protože osy míchají centimetry s kilogramy — v praxi se proto souřadnice s různými jednotkami před PCA často standardizují).

Nechť n ≥ 2, D ≥ 1, x₁, …, x_n ∈ ℝ^D jsou body a u > 0 zvolená perplexita. t-SNE (t-distributed stochastic neighbor embedding; van der Maaten a Hinton, 2008) přiřadí bodům x_i body mapy y₁, …, y_n ∈ ℝ² takto:

1. z bodů x_i a perplexity u spočte podobnosti p_ij ≥ 0 pro i ≠ j se součtem ∑_i≠j p_ij = 1 (další definice) — p_ij vyjadřuje, nakolik jsou i a j sousedé v původním prostoru;
2. pro body mapy definuje podobnosti se Studentovým jádrem (Studentovo rozdělení s jedním stupněm volnosti) q_ij = (1 + ‖y_i − y_j‖²)⁻¹/∑_k≠l (1 + ‖y_k − y_l‖²)⁻¹   pro i ≠ j;
3. body y₁, …, y_n volí tak, aby byla co nejmenší Kullbackova–Leiblerova divergence C(y₁, …, y_n) = KL(P ‖ Q) = ∑_i≠j p_ij ln p_ij/q_ij (sčítají se jen dvojice s p_ij > 0); minimum hledá gradientním sestupem z náhodného startu, takže najde jen lokální minimum.

t-SNE tedy rozmístí body v rovině tak, aby si sousedé z původního prostoru zůstali sousedy i na mapě.

Nechť n ≥ 2, x₁, …, x_n ∈ ℝ^D a pro každé i = 1, …, n nechť σ_i > 0 (šířka jádra bodu i). Pro i ≠ j položme

p_j|i = exp(−‖x_i − x_j‖²/2σ_i²)/∑_k≠i exp(−‖x_i − x_k‖²/2σ_i²)

— pravděpodobnost, že si bod i vybere za souseda bod j; pro pevné i tvoří čísla P_i = (p_j|i)_j≠i rozdělení pravděpodobnosti. Jeho entropie a perplexita jsou

H(P_i) = −∑_j≠i p_j|i log₂ p_j|i,   Perp(P_i) = 2^H(P_i).

Šířka σ_i se pro každý bod zvlášť najde půlením intervalu tak, aby Perp(P_i) = u (zadaná perplexita, u nás 30). Nakonec se podobnosti symetrizují:

p_ij = p_j|i + p_i|j/2n   pro i ≠ j,

takže p_ij = p_ji a ∑_i≠j p_ij = 1 (každé P_i má součet 1, dohromady 2n). Kdyby bod měl přesně 30 stejně pravděpodobných sousedů a ostatní nulové, vyšlo by H = log₂ 30 a perplexita 30 — perplexita je tedy „efektivní počet sousedů“. V hustém místě tak vyjde σ_i malé, v řídkém velké.

Nechť n ≥ 6, y₁, …, y_n ∈ ℝ^q jsou body zobrazení (mapy nebo původního prostoru) a t₁, …, t_n jejich skutečné třídy (u nás číslice 0–9). Pro každé i nechť N₅(i) ⊆ {1, …, n} ∖ {i}, |N₅(i)| = 5, je množina indexů pěti bodů nejbližších k y_i (bez něj samotného; shody vzdáleností se rozhodnou libovolně). Čistota sousedů je číslo

čistota = 1/n ∑_i=1ⁿ |{j ∈ N₅(i) : t_j = t_i}|/5 ∈ [0, 1],

tedy průměrný podíl sousedů se stejnou třídou. Bod, jehož pět sousedů jsou tři stejné číslice a dvě jiné, přispěje 3/5. Míra vůbec nepoužívá vzdálenosti mezi ostrůvky — měří jen to, na čem t-SNE záleží.