2. 線形変換，白色化，同時対角化【Math2Code】

研究している中で，パターン認識の数式を自分でアルゴリズムに直して記述する機会があったので，メモ．


# 線形変換
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/e6241ffa-4d01-a273-d61e-45026a7beca8.png)

```
import numpy as np
import matplotlib.pyplot as plt

cov=np.matrix([[1, 0], [0, 3]])
icov=np.linalg.inv(cov)
r=np.matrix([[np.cos(np.pi/3), -np.sin(np.pi/3)], [np.sin(np.pi/3), np.cos(np.pi/3)]])
xx,yy = np.meshgrid(np.linspace(-5,5),np.linspace(-5,5))
plt.figure()
plt.axis('equal')
p=1./(2.*np.pi*np.sqrt(np.linalg.det(cov))) * np.exp(-1./2.*(icov[0,0]*xx*xx+(icov[0,1]+icov[1,0])*xx*yy+icov[1,1]*yy*yy))
plt.contour(xx,yy,p)
#plt.savefig('linearnorm-1.eps')
plt.figure()
plt.axis('equal')
cov2=r.T.dot(cov).dot(r)
icov2=np.linalg.inv(cov2)
pp=1./(2.*np.pi*np.sqrt(np.linalg.det(cov2))) * np.exp(-1./2.*(icov2[0,0]*xx*xx+(icov2[0,1]+icov2[1,0])*xx*yy+icov2[1,1]*yy*yy))
plt.contour(xx,yy,pp)
#plt.savefig('linearnorm-2.eps')
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/ef2a7ad2-57aa-67b6-2260-c2a9cdb7392f.png)


# 白色化
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/a7a5b4c0-402a-e3db-c180-1e12201dd5fc.png)


```
import numpy as np
import matplotlib.pyplot as plt

d=2; n=1000
x=np.matrix(np.random.randn(n,d))
s=np.matrix([[.7, 0.],[0., 2.]])
r=np.matrix([[np.cos(np.pi/3), -np.sin(np.pi/3)],[np.sin(np.pi/3), np.cos(np.pi/3)]])
t=np.matrix([0.5, -1.])
x=x.dot(s).dot(r)+np.ones([n,1]).dot(t)
m=np.mean(x,axis=0);
cov=(x-np.ones([n,1]).dot(m)).T.dot(x-np.ones([n,1]).dot(m))/n
icov=np.linalg.inv(cov)
xx,yy=np.meshgrid(np.linspace(-5,5),np.linspace(-5,5))
xt=xx-m[0,0]
yt=yy-m[0,1]
p=1./(2.*np.pi*np.sqrt(np.linalg.det(cov))) * np.exp(-1./2.*(icov[0,0]*xt*xt+(icov[0,1]+icov[1,0])*xt*yt+icov[1,1]*yt*yt))
plt.figure()
plt.axis('equal')
plt.scatter([x[:,0]],[x[:,1]])
plt.contour(xx,yy,p,cmap='hsv')
#plt.savefig('whitening-before.eps')

ed,ev=np.linalg.eig(cov)
tr=np.diag(ed**(-1./2.)).dot(ev.T)
x1=x.dot(tr.T)
m1=np.mean(x1,axis=0)
xt1=xx-m1[0,0]
yt1=yy-m1[0,1]
a=tr.dot(cov).dot(tr.T)
ia=np.linalg.inv(a)
p1=1./(2.*np.pi*np.sqrt(np.linalg.det(a))) * np.exp(-1./2.*(ia[0,0]*xt1*xt1+(ia[0,1]+ia[1,0])*xt1*yt1+ia[1,1]*yt1*yt1))
plt.figure()
plt.axis('equal')
plt.scatter([x1[:,0]],[x1[:,1]])
plt.contour(xx,yy,p1,cmap='hsv')
#plt.savefig('whitening-after.eps')
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/3394a867-c0ea-549e-78b4-01f35ee32699.png)


###  マハラノビス距離
ちなみに，上記の等高線はマハラノビス距離で描いている．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/d70f307e-1ac6-5e04-73ec-1c6451427e5e.png)


# 同時対角化
2つの正規分布を同時に対角化する．1つ目を白色化して，2つ目を対角化するだけ．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/32a7b498-7a24-2cfa-fe37-3d939c52a7b1.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/d1743e1a-dd43-214b-52d2-82b9f20b233c.png)




```
from scipy.io import loadmat
data = loadmat('data.mat')
x = data['x']
cov1 = data['cov1']
cov2 = data['cov2']
m1 = data['m1']
m2 = data['m2']

import numpy as np
import matplotlib.pyplot as plt

plt.figure()
plt.scatter([x[:, 0]], [x[:, 1]], s=10)

xx1,yy1=np.meshgrid(np.linspace(-10,10),np.linspace(-10,10))
xt1=xx1-m1[0,0]
yt1=yy1-m1[0,1]
icov1 = np.linalg.inv(cov1)
p1=1./(2.*np.pi*np.sqrt(np.linalg.det(cov1))) * np.exp(-1./2.*(icov1[0,0]*xt1*xt1+(icov1[0,1]+icov1[1,0])*xt1*yt1+icov1[1,1]*yt1*yt1))
plt.contour(xx1,yy1,p1,cmap='hsv')

xx2,yy2=np.meshgrid(np.linspace(-10,10),np.linspace(-10,10))
xt2=xx2-m2[0,0]
yt2=yy2-m2[0,1]
icov2 = np.linalg.inv(cov2)
p2=1./(2.*np.pi*np.sqrt(np.linalg.det(cov2))) * np.exp(-1./2.*(icov2[0,0]*xt2*xt2+(icov2[0,1]+icov2[1,0])*xt2*yt2+icov2[1,1]*yt2*yt2))
plt.contour(xx2,yy2,p2,cmap='hsv')

ed,ev=np.linalg.eig(cov1)
tr=np.diag(ed**(-1./2.)).dot(ev.T)
y = x.dot(tr.T)
y = x.dot(ev).dot(np.diag(ed**(-1./2.)))
plt.figure()
plt.axis('equal')
plt.scatter([x1[:,0]],[x1[:,1]], s=1)
plt.show()

i = tr.dot(cov1).dot(tr.T)
k = tr.dot(cov2).dot(tr.T)
print(f'I: {i}')
print(f'K: {k}')
ed_k, ev_k = np.linalg.eig(k)
z = y.dot(ev_k.T)

print(ev_k.T.dot(k).dot(ev_k))

plt.figure()
plt.scatter([z[:,0]],[z[:,1]], s=10)

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/3bcb4fd0-bf92-c1f1-f360-040598549df9.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/711751/70a02386-2e73-8106-60e5-768d1c7f02be.png)