Gauss hatalı gözlem modeli

$A$ bir matris, $y$ ve $x$ de uygun boyutlarda birer vektör olsun. $y$'yi ekteki gibi saklı $x$ bileşenlerinin doğrusal dönüşümüne Gauss hata eklenmiş şekilde gözlemlediğimizi varsayalım:
\begin{align}
y = Ax + \eta
\end{align}

$\eta$ için Gauss dağılsın demiştik. Hatanın neden genellikle Gauss dağıldığını varsaydığımızı merak ediyorsanız merkezi limit teoremi yazısının sondan üçüncü paragrafını okuyabilirsiniz. $I$ birim matrisi göstermek üzere $\eta \sim \mathcal{N}(0,I)$ olsun. Yani hata ortalaması $0$, kovaryansı birim matris olan bir Gauss dağılımından geliyor. $A$ ve $x$'i bildiğimiz taktirde $y$'nin olabilirliği nedir?

\begin{align*}
p(y|x) = \mathcal{N}(y;Ax,I)
\end{align*}

Fakat biz $x$'i bilmiyoruz, bulmak istediğimiz zaten o. O halde $y$'nin, yani gözlemlenen verinin olabilirliği en yüksek yapan $x$'i arıyoruz. Aradığımız $x$ yıldızlı $x$ olsun, $x^*$ ile gösterelim:
\begin{align*}
x^* = \argmax_x p(y|x)
\end{align*}

Bir değeri en büyüklerken kendisi yerine logunu kullanabiliriz:
\begin{align}
x^* &= \argmax_x p(y|x) \nonumber \\
&= \argmax_x \log p(y|x) \nonumber \\
&= \argmax_x \log \mathcal{N}(y;Ax,I) \label{loglik}
\end{align}

Biliyoruz ki,
\begin{align*}
\mathcal{N}(y;Ax,I) = \frac{1}{\sqrt{(2\pi)^d}} \exp\left(-\frac{1}{2} (y-Ax)^T (y-Ax) \right)
\end{align*}
$\log \mathcal{N}(y;Ax,I)$ değerini en büyükleyen $x$'i aradığımızı hatırlayalım. Denklem \eqref{loglik}'yi yeniden yazalım:
\begin{align*}
x^* &= \argmin_x \frac{1}{2} \|y - Ax\|_2^2
\end{align*}

Bu sonuca nasıl vardık? Tanım gereği, $\|y - Ax\|_2 = \sqrt{(y-Ax)^T (y-Ax)}$. Öyleyse, $(y-Ax)^T (y-Ax)$ ile $\|y - Ax\|_2^2$ aynı. Denklemdeki diğer terimler $x$'den bağımsız. Bu sebeple optimizasyon sonucunu etkilemiyorlar.

En küçük karelerden Moore-Penrose tersine

Hatayı Gauss kabul edip olabilirliği en büyüklemek ile en küçük kareler yaklaşımının aynı şey olduğunu gördük. Burada bırakmayalım, çözüme devam edelim:

\begin{align*}
x^* &= \argmin_x \frac{1}{2} \|y - Ax\|_2^2 \\
&= \argmin_x (y-Ax)^T (y-Ax) \\
&= \argmin_x y^T y - y^T A x - x^T A^T y + x^T A^T A x \\
&= \argmin_x - y^T A x - x^T A^T y + x^T A^T A x
\end{align*}

Görülüyor ki aslında en küçük değerini bulmaya çalıştığımız basit bir masraf fonksiyonu var:

\begin{align*}
J(x) = - y^T A x - x^T A^T y + x^T A^T A x
\end{align*}

Tam bu noktada matematiğin sihirli değneğini kullanalım. $J(x)$'daki tüm terimler skaler. Dolayısıyla bu terimlerin izleri (trace) de skaler olacak. Unuttuysak hatırlayalım; $\operatorname{Tr}(B)$ işlevi herhangi kare bir $B$ matrisinin köşegenlerinin toplamına eşit. İz işlevinin sihirli değneğinin dokunduğu masraf fonksiyonuna dönelim:

\begin{align*}
J(x) = \operatorname{Tr}(x^T A^T A x) - \operatorname{Tr}(y^T A x) - \operatorname{Tr}(x^T A^T y)
\end{align*}

İz işlevinin şahane özellikleri var, birkaçını kullanacağız
\begin{align*}
\operatorname{Tr}(A) &= \operatorname{Tr}(A^T) \\
\operatorname{Tr}(AB) &= \operatorname{Tr}(BA)
\end{align*}

Son özelliğin tüm dairesel permütasyonlar için genellenebilir olduğunu görün. Örneğin $\operatorname{Tr}(ABC) = \operatorname{Tr}(BCA) = \operatorname{Tr}(CAB)$.

$J(x)$'in en küçük değerini bulmak için türevini almalıyız. İz işlevinin türevine ihtiyaç duyacağız. Doğrusal cebirden hatırlayalım:

\begin{align*}
\frac{\partial \operatorname{Tr}(Ax)}{\partial x} &= A^T \\
\frac{\partial \operatorname{Tr}(Ax^T)}{\partial x} &= A \\
\frac{\partial \operatorname{Tr}(A x x^T)}{\partial x} &= Ax + A^T x
\end{align*}

Bu kuralları uygulayarak $x$'e göre masraf fonksiyonunun kısmi türevini alalım
\begin{align*}
\frac{\partial J(x)}{\partial x} &= \frac{\partial (-\operatorname{Tr}(y^T A x) - \operatorname{Tr}(x^T A^T y) + \operatorname{Tr}(x^T A^T A x))}{\partial x} \\
&= 2A^T A x - 2 A^T y
\end{align*}
ve türevi sıfıra eşitleyelim,
\begin{align*}
2A^T A x = 2 A^T y
\end{align*}

Dolayısıyla,
\begin{align*}
x^* = (A^T A)^{-1} A^T y
\end{align*}

Bulduğumuz $(A^T A)^{-1} A^T$ terimine $A$'nın Moore-Penrose tersi (pseudoinverse) diyoruz.

Ve yazıyı sonlandırırken Deniz'e bu paylaşım için teşekkür ediyoruz.

Soru 1: PCA'de kullandığınız matlab kodlarına, princomp komutu alternatif olabilir mi? Farkı nedir?

Cevap: Sonuç olarak aralarında pozitif/negatif işareti dışında bir fark yok. princomp fonksiyonu daha okunabilir, anlaşılabilir bir fonksiyon. Eğer MATLAB'ınızda yoksa svd ile aynısını kolaylıkla yapabilirsiniz. Örnek üstünden bakalım:

load hald;
X = ingredients; % Veri X olsun. Her satırda bir örnek var.
X

X =
     7    26     6    60
     1    29    15    52
    11    56     8    20
    11    31     8    47
     7    52     6    33
    11    55     9    22
     3    71    17     6
     1    31    22    44
     2    54    18    22
    21    47     4    26
     1    40    23    34
    11    66     9    12
    10    68     8    12

Öncelikle halihazırdaki fonksiyonu kullanalım:

[pc1,score1] = princomp(X)

pc1 =
   -0.0678   -0.6460    0.5673    0.5062
   -0.6785   -0.0200   -0.5440    0.4933
    0.0290    0.7553    0.4036    0.5156
    0.7309   -0.1085   -0.4684    0.4844
score1 =
   36.8218   -6.8709   -4.5909    0.3967
   29.6073    4.6109   -2.2476   -0.3958
  -12.9818   -4.2049    0.9022   -1.1261
   23.7147   -6.6341    1.8547   -0.3786
   -0.5532   -4.4617   -6.0874    0.1424
  -10.8125   -3.6466    0.9130   -0.1350
  -32.5882    8.9798   -1.6063    0.0818
   22.6064   10.7259    3.2365    0.3243
   -9.2626    8.9854   -0.0169   -0.5437
   -3.2840  -14.1573    7.0465    0.3405
    9.2200   12.3861    3.4283    0.4352
  -25.5849   -2.7817   -0.3867    0.4468
  -26.9032   -2.9310   -2.4455    0.4116

Şimdi kendimiz yazalım. Önce her satırdan ortalamayı çıkaralım. Ardından svd yapalım. Dikkat, burada pc2 değişkeni aslında V'ye denk gelen yerde, çünkü örnekler satırlarda. Eğer örnekler sütunlarda olsaydı, [pc2, ~, ~] biçiminde olacaktı.

X = X - repmat(mean(X,1), size(X,1), 1);
[~, ~, pc2] = svd(X, 'econ');
score2 = X * pc2

pc2 =
    0.0678    0.6460   -0.5673    0.5062
    0.6785    0.0200    0.5440    0.4933
   -0.0290   -0.7553   -0.4036    0.5156
   -0.7309    0.1085    0.4684    0.4844
score2 =
  -36.8218    6.8709    4.5909    0.3967
  -29.6073   -4.6109    2.2476   -0.3958
   12.9818    4.2049   -0.9022   -1.1261
  -23.7147    6.6341   -1.8547   -0.3786
    0.5532    4.4617    6.0874    0.1424
   10.8125    3.6466   -0.9130   -0.1350
   32.5882   -8.9798    1.6063    0.0818
  -22.6064  -10.7259   -3.2365    0.3243
    9.2626   -8.9854    0.0169   -0.5437
    3.2840   14.1573   -7.0465    0.3405
   -9.2200  -12.3861   -3.4283    0.4352
   25.5849    2.7817    0.3867    0.4468
   26.9032    2.9310    2.4455    0.4116

Dikkat ederseniz işaretler dışında sonuçlar aynı. Bu durum uygulamanız için sorun teşkil etmez.

Soru 2: PCA ile örüntü analizi nasıl yapılıyor? Yani elimizde bir veri seti var diyelim ki. Örüntü analizi için kullanılacak. Verinin bir kısmını PCA, bir kısmını da test için mi kullanmak gerekiyor. Test ve sınıflandırma işlemi nasıl gerçekleştiriliyor?

Cevap: PCA boyut indirgemeye yarar. Bu bağlamda, diğer uygulamalar için genellikle önişlem olarak kullanılıyor. Örneğin bağlanım (regression) yapacaksanız çok fazla boyutla uğraşmak yerine öncelikle PCA uygulayıp verinin boyutunu düşürüyorsunuz. Aynı durum sınıflama için de geçerli. Yukarıdaki örnek üstünden bakacak olursak X'teki her örnek normalde 4 boyutlu. score1'in ilk iki sütununu alır, gerisini atarsak her örneğin yeni uzaydaki boyutu 2'ye düşer. Yarı yarıya kâr! Yani sizin yeni örnekleriniz score1 ile ifade edilen koordinatlara sahip örnekler. Peki yeni boyut sayısı kaç olmalı. Literatürde genel olarak özdeğerlerin birikimsel olarak %95'ini alacak şekilde boyut sayısı seçiliyor:

[pc1,score1, eigvalues] = princomp(X);
cumsum(eigvalues)./sum(eigvalues)

ans =
    0.8660
    0.9789
    0.9996
    1.0000

Görüldüğü üzere, 2 tane koordinatı seçtiğimizde verinin %97.89'unu ifade edebiliyoruz. Bunun yeterli olacağını, gerisinin gürültü olduğunu varsayabiliriz. İllâ ki böyle yapacaksınız diye bir şart yok, isterseniz limiti %98 seçebilirsiniz.

Şimdi eğitim ve test aşamalarında ne yapacağımıza bakalım. PCA'i ifade etmek için temel bileşenlere ve ortalama vektöre ihtiyaç var. Bunları eğitim kümesini kullanarak elde edeceksiniz. Ardından test kümesindeki her örnekten bu ortalama değeri çıkarıp temel bileşenler ile çarpacak ve bu örneklerin indirgenmiş uzaydaki hallerini bulacaksınız. Sonra bağlanım mı yaparsınız, sınıflandırma mı, kullanacağınız hep bu yeni değerler olacak.

Alıştırma:

Burada okura bir sınıflandırma alıştırması vereyim. PCA ile temel bir uygulama alıştırması:

1. MNIST elle yazılmış rakam veritabanını indirin.
2. Normalde eğitim kümesinde 50000, test kümesinde ise 10000 adet örnek var. Her örneği düzleştirin ve yeni örnekler 20x20 boyutlu yerine 400x1 boyutlu olsun.
3. 50000 örnek ek teknik kullanmadan işlem yapabilmek için çok fazla. Her bir rakam için 500 adet rastgele örneği alın ve böylece 5000 elemanlık bir eğitim kümesi oluşturun.
4. PCA uygulayın ve temel bileşenler ile ortalama vektörleri bulun. Sadece %95 birikimsel özdeğeri saklayacak kadar temel bileşeni tutun, gerisini atın.
5. Eğitim kümesindeki her örneğin indirgenmiş boyutlarını hesaplayın ve saklayın.
6. 10000 test örneğinden her biri için
   1. Eğitim kümesindeki örneklere yaptığınız gibi ortalamayı çıkarın ve temel bileşenler ile çarparak boyut indirgeyin.
   2. Eğitim kümesindeki her örneğe olan Öklit uzaklığını hesaplayın. Vektörler 400x1 boyuttan daha az olacağı için işlemler hızlı olacak.
   3. En yakın vektörü tespit edin ve bu test örneğinin etiketini o vektörün etiketine atayın. Yani en yakın vektör 5 rakamına ait  ise sizin gözlemlediğiniz vektör de 5 rakamı olacak.
7. 10000 vektörün kaç tanesinin doğru hesaplandığını bulun.
8. Aynı işlemi PCA uygulamadan deneyin, farkı görün.

$$\log\sum\limits_k e^{x_k}$$

$e^{x_k}$ değerlerinin -1000lerde seyrettiğini düşünün. $e^{1000}$ gibi bir sayıyı açıkça hesaplamaya çalışırsak sıfır elde ederiz.

Probleme tane tane değil, tümden bakalım. Elimizde 3 tane değer olsun: -1000, -1005 ve -1010. Bunları işleyelim. $e^{1005}$,  $e^{1000}$'e göre oldukça küçük.  $e^{1010}$ ise daha da küçük. Burada baskın olan  $e^{1000}$ değeri. Yani, en büyük değer baskın görünüyor. En büyük değere $y$ diyelim ve bunu hesaba katarak sihirli değneğimizi kullanalım:

\begin{align}
\log\sum\limits_k e^{x_k}
& = \log\left[\left(\sum\limits_k e^{x_k}\right)e^{-y} e^y \right] \\
& = \log\left[\left(\sum\limits_k e^{x_k}\right)e^{-y}\right] + y \\
& = \log\left[\sum\limits_k e^{x_k - y}\right] + y
\end{align}

Yani öncelikle tüm değerlerden en büyüğünü çıkartacağız (yeni dizilimin en büyük değeri 0 olacak) ve işlemi yaptıktan sonra bu değeri geri ekleyeceğiz.

Fonksiyonumuz çok basit:

import numpy as np

def logSumExp(x):
    maxX = np.max(x)
    r = maxX + np.log(np.sum(np.exp(x-maxX)))
    return r

Deneyelim:

r = logSumExp(np.array([-1213, -1214]))  # sonuç: -1212.6867383124818
r2 =  np.log(np.sum(np.exp(np.array([-1213, -1214]))))  # sonuç: -Inf

Yöntemimiz sonucu başarıyla hesapladı. Açıkça yaptığımızda ise underflow oldu ve log(0) hesabı yüzünden -sonsuz sonucu döndü.

Not: Kaynak olarak şu videoya bakabilirsiniz.

Önceden Lilypond'u öğrendiğimde düşüncelerimi aktarmıştım. Hatta Oynaya Oynaya Gelin Çocuklar'ı eklemiştim. Musescore, Lilypond'a göre çocuk oyuncağı. Üstelik notaları dinlemek de çok eğlenceli. Bu işi yapan Sibelius, Finale gibi ücretli programlar olduğunu biliyorum. Musescore özgür ve ücretsiz bir yazılım. Üstelik musescore.com adresinde kullanıcıların paylaştığı tonla hazır eser mevcut. Ben de bugün can sıkıntısına oturdum, Haluk Levent'in en sevdiğim şarkılarından biri olan Beni Biraz Anlasana'nın notalarını ve akorlarını ekledim. Alttaki pencerede hem notaları görebilir hem de dinleyebilirsiniz. Eğer düzgün görünmezse şarkı bu adreste. İsterseniz indirebilir, sesinize/zevkinize göre transpose edebilir veya dilediğinizce değiştirebilirsiniz.

Bir de akordeon ile çalıp Youtube'a ekledim. Eğer düzgün görüntülenmezse o da bu adreste.

Not: Akor düzenlemesi, N. Gül'e aittir. Ezgiyi kulaktan yazdım. Hatalar olabilir.