Robust PCA (Principlal Component Pursuit) の実装

はじめに

多変量データは実世界の様々なところで現れる(e.g., 画像、音声、動画).これらの多変量データの多くは、データ自体がもつ次元 (e.g., ピクセル数)よりも小さい次元(自然画像の多様体的なやつ)で表現することができる. データから低次元の構造(低ランク行列)を取り出す代表的な方法として、PCAがある. しかし、実世界で観測されるデータには外れ値が含まれていることが多く、PCAはこの外れ値の影響を強く受けて低ランク行列を求めてしまうという性質がある. 外れ値にロバストに低ランク行列を推定する方法 (a.k.a., robust PCA)として、Principal Component Pursuitが提案されている. ここではPCPの簡単な説明と実装を行い、PCAと比較する.

低ランク行列

なんらかの低ランクな構造を持っている(日本語あやしい)多変量データ X\in \mathbb{R}^{d \times N}は低ランクな行列 L\in \mathbb{R}^{d \times N}とノイズ E\in \mathbb{R}^{d \times N}の和としてかける.ここで、 Nはデータ数、 dは各データの次元を表す.

 X = L + E

ここでは例として、MNISTのデータを Xとして、低ランクな行列 Lを求める問題を考える.

データのダウンロードなどは以下. 最近scikit-learnのMNIST downloadはうまくいかないらしいので,こちらを参考にしてダウンロードした.(いつのまにかdeprecatedになっていた)

以降では Xの平均は0であるとする. X \boldsymbol{1} = \boldsymbol{0}

import numpy as np
import pylab as plt
from sklearn.decomposition import PCA
from sklearn import datasets
%matplotlib inline
import urllib.request
import os.path

url_base = 'http://yann.lecun.com/exdb/mnist/'
key_file = {
    'train_img':'train-images-idx3-ubyte.gz',
    'train_label':'train-labels-idx1-ubyte.gz',
    'test_img':'t10k-images-idx3-ubyte.gz',
    'test_label':'t10k-labels-idx1-ubyte.gz'
}
dataset_dir = "."

def _download(filename):
    file_path = dataset_dir + '/' + filename
    if os.path.exists(file_path):
        return print('already exist')
    print('Downloading ' + filename + ' ...')
    urllib.request.urlretrieve(url_base + filename, file_path)
    print('Done')

def download_mnist():
    for v in key_file.values():
       _download(v)

download_mnist()
def load_mnist(filename):
    img_size = 28**2
    file_path = dataset_dir + '/' + filename
    with gzip.open(file_path, 'rb') as f:
        data = np.frombuffer(f.read(), np.uint8, offset=16)
    return data.reshape(-1, img_size)
import gzip

nb_samples = 3000

data = load_mnist(key_file['train_img'])[:nb_samples]
data = data/255
img1 = data[0].reshape(28, 28)
import pylab as plt
plt.imshow(img1)

f:id:ksknw:20181117140525p:plain

with gzip.open(key_file["train_label"], 'rb') as f:
    labels = np.frombuffer(f.read(), np.uint8, offset=8)[:nb_samples]
labels.shape
(3000,)
X = data - data.mean(axis=0)
plt.pcolormesh(X)
plt.show()

f:id:ksknw:20181117140750p:plain

PCA

PCAではノイズ Eとして、ガウス分布から発生するノイズを考えて、以下のような最適化問題によって低ランク行列 Lを求める.(実際にPCAすると"分散が最大になる方向"に軸を置き直すが、その部分は今は考えずに、低ランクに近似、つまり次元削減の部分にのみ着目する.)

 \mathrm{argmin}_ L \|X - L \|^2 _F \ \  \mathrm{s.t.} \  \mathrm{rank}(L) \lt p

ここで \| \cdot \|^2 _Fはフロベニウスノルム. つまり、PCAはl2の意味でデータ Xに最も"近い"低ランク行列 Lを求めているといえる. この問題は X特異値分解することで、以下のように Lを求める方法が知られている

 L = U_p \Sigma _p V_p ^\top ここで X = U \Sigma V^\topであり、 \Sigma_pは特異値 \sigma_1, \dots \sigma_dのうち、大きい方から p個を対角要素に並べた行列.  U_p, V_pはそれに対応した行,列を抜き出したもの.

from numpy.linalg import svd
u, Σ , vh = svd(X, full_matrices=False)
L = (u[:,:2] @ np.diag(Σ[:2]) @ vh[:2])
plt.pcolormesh(L)
plt.show()

f:id:ksknw:20181117140806p:plain

ジャギジャギしていたデータがなんか綺麗になっていて,なんとなく二次元ぐらいで表現されているような気がする (適当). ここままだと見にくいので,PCAと同じようにデータを特異ベクトル(特異値が大きいほうから2次元)の方向へ射影して散布図として表現すると以下のようになる.

L_d =  X @ np.transpose(vh)
plt.scatter(L_d[:,0], L_d[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140832p:plain

これは符号を除いて,PCAによって求まるものと全く同じになる. 実際にscikit-learnのPCAはSVDで解かれている

pca = PCA()
pcaed = pca.fit_transform(data)
plt.scatter(pcaed[:,0], pcaed[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140850p:plain

外れ値が含まれるデータに対するPCA

ここでデータに外れ値が含まれている場合を考える. 例として、MNISTに含まれるデータのうち、1ピクセルが他のピクセルよりも大きめのランダムな値であるとする.

# data[0] = np.random.randn(784) * 4
data[0, 0] = 100
X = data - data.mean(axis=0)
u, s, vh = svd(X, full_matrices=False)
L_d_with_noise =  X @ np.transpose(vh)
plt.scatter(L_d_with_noise[:,0], L_d_with_noise[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140912p:plain

このようにPCAでは少ない数の外れ値に大きく影響された射影が得られてしまうことがある.

Principal Component Pursuit

外れ値は現実のデータにおいてよく現れるが、データに含まれる構造を見たいという目的で低ランク近似をしていた場合、少ない数の外れ値に大きく影響された結果は望ましくない.

そこで、外れ値の存在を考慮して低ランク行列 Lを求めるアルゴリズム(Principal Component Pursuit)が提案されている.

ここで外れ値を S\in \mathbb{R}^{d\times N}と書く.外れ値は全体の一部であるとして、 Sはスパースな行列であると仮定する. Principal Component Pursuitでは、行列 Xを以下のように表すことができると考える.

 X = L + S

スパース制約をナイーブに考えると Sのl0ノルムを小さくしつつ、 Lを求める必要があるが、この問題はNP-hardなので難しい. そこで、PCPでは以下の緩和した問題を考える.

 \min |L|_* + \lambda |S|_1, \ \ s.t.\ L+S=X

ここで | \cdot |_*はトレースノルム(特異値の和)であり, | \cdot |_1はl1ノルム. この制約付き最適化問題拡張ラグランジュ法を用いて解くことができる. PCPの論文ではこの緩和した問題を解くことで得られる L, Sがいくつかの仮定のもとで、はじめの最適化問題の解と一致することを示している (証明ちゃんと読んでないのでたぶん). 仮定の1つとして(??)、パラメータ \lambdaが以下の値であるとしている.

 \lambda = \frac{1}{\sqrt{n_{(1)}}}、ここで、  n_{(1)} = \max (d, N)である

以上から、PCPを実装すると以下のようになる.

X = data
X.shape
(3000, 784)
S = 0
Y = 0
μ = 2
λ = 1/(max(X.shape))**0.5
print(λ)
nb_epochs = 100
0.018257418583505537
def shrink(x, τ):
    shape = x.shape
    x = x.reshape(-1)
    temp = np.apply_along_axis(lambda a: np.max(np.c_[np.abs(a)-τ, np.zeros(len(a))], axis=1), 0, x)
    return np.sign(x.reshape(shape)) * temp.reshape(shape)
from tqdm import tqdm

hist_normY = []
for i in tqdm(range(nb_epochs)):
    u, Σ, vh = np.linalg.svd(X - S + Y/μ, full_matrices=False)
    L = np.matmul(u * shrink(Σ, 1/μ), vh)
    S = shrink(X - L + Y/μ, λ/μ)
    Y = Y + μ*(X - L - S)
100%|██████████| 100/100 [00:57<00:00,  1.73it/s]
plt.scatter(S[:,0], S[:,1], marker=".", c=labels)
plt.show()

f:id:ksknw:20181117141112p:plain

L_d_PCP =  (X - S - Y)@ np.transpose(vh)

結果

  • PCPで求めた低ランク行列
plt.scatter(L_d_PCP[:,0], L_d_PCP[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117140940p:plain

  • ノイズがない場合のPCAの結果
plt.scatter(L_d[:,0], L_d[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117141007p:plain

  • ノイズがある場合のPCAの結果
plt.scatter(L_d_with_noise[:,0], L_d_with_noise[:,1], marker=".", c=labels, cmap="jet")
plt.show()

f:id:ksknw:20181117141019p:plain

おわりに

外れ値にロバストなPCAとして、Principal Component Pursuitの実装を行った. PCAを行列の低ランク近似だと思うと、スパースな外れ値がのったモデルを考えるのは妥当であるように思える. 実装した結果を見ても、きちんと外れ値を推定することができていた. 今回はMNISTに適当なノイズをのせたデータを用いたが、論文には実際の動画データに適用した例もある(長時間映っている背景部分が低ランク行列で表現され、移動している人が外れ値になる).

この論文はむちゃくちゃ引用されており、Robust Principal Component Pursuitなど(まだロバストにしたいのか)色々発展もありそうなので、チェックしてみたい.

参考

1次元ガウス分布の測地線と双対測地線のプロット

はじめに

情報幾何学の新展開」という本を読んでいる. まだ序盤しか読めてない上に,あまり理解できていないが,自分の理解のために,例として1次元ガウス分布を対象として,以下の導出とプロットをやる.

  • 指数型分布族の標準形および双対座標系
  • ポテンシャル関数,双対ポテンシャル関数
  • 測地線,双対測地線

かなり天下り的にやっている部分が多いので,主にプロットしただけという感もある.

import numpy as np
import pylab as plt
from numpy import pi as π
from numpy import log
from numpy import exp
from mpl_toolkits.mplot3d import Axes3D
e = exp(1)
%matplotlib inline

座標系のプロット

指数型分布の標準形は以下のようなもの.  p({\bf x, \theta}) = \exp\left(\sum_i \theta^{(i)} x_i - \psi({\bf \theta})\right) 1次元ガウス分布 p(x, \mu , \sigma) = \frac{1}{\sqrt{2 \pi \sigma^2} } \exp \left( \frac{(x-\mu)^2}{2 \sigma^2} \right) であるので,これをごちゃごちゃいじって,

 \theta^{(1)} = \frac{\mu}{\sigma^2}, \theta^{(2)} = \frac{1}{2\sigma^2}, x_{(1)}=x, x_{(2)}=-x^2とすると,ごちゃごちゃいじって,

 p(x, \theta) = \exp(\sum_i \theta^{(i)} x_i - \psi({\bf \theta})) - \left( \frac{(\theta^{(1)})^2}{4\theta^{(2)}} + \frac{1}{2}\log \pi -\frac{1}{2} \log \theta^{(2)}\right)

よって,  \psi(\theta) = \frac{(\theta^{(1)})^2}{4\theta^{(2)}} + \frac{1}{2}\log \pi -\frac{1}{2} \log \theta^{(2)} となる.

ここで, \psi(\theta)をプロットすると,以下のように,これが凸関数っぽくなっていることがわかる.(実際にも凸関数である)

def p2θ(μ, σ2):
    θ1 = μ/σ2
    θ2 = 1/(2*σ2)
    return θ1, θ2

def θ2p(θ1, θ2):
    μ = θ1/(2*θ2)
    σ2 = 1/(2*θ2)
    return μ, σ2

def ψ(θ1, θ2):
    return (θ1)**2/ (4*θ2) + 1/2 * log(π) + 1/2*log(θ2) 
def plot_lattice(t1, t2, c, alpha=1):
    for t1_i, t2_i in zip(t1, t2):
        plt.plot(t1_i, t2_i, c=c, alpha=alpha)
    for t1_i, t2_i in zip(t1.transpose(), t2.transpose()):
        plt.plot(t1_i, t2_i, c=c, alpha=alpha)
temp_θ1 = np.linspace(-1,1, 30)
temp_θ2 = np.linspace(1, 3, 30)

θ1, θ2 = np.meshgrid(temp_θ1, temp_θ2)
ψ_θ = ψ(θ1, θ2)
μ, σ2 = θ2p(θ1,θ2)
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_wireframe(θ1, θ2, ψ_θ)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")
ax.set_zlabel("$\psi(\\theta)$")
plt.show()

f:id:ksknw:20180814095351p:plain

%matplotlib inline
plt.figure(figsize=(8,4))

plt.subplot(121)
plot_lattice(μ, σ2, "C0", alpha=0.3)
cont = plt.contour(μ, σ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(122)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
cont = plt.contour(θ1, θ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095418p:plain

ここで \psi(\theta)微分 {\bf \eta} = \nabla \psi (\theta)を考えると,  \eta^{(1)} = \frac{\theta^{(1)}}{2\theta^{(2)}}, \eta^{(2)} = -\frac{(\theta^{(1)})^2}{(\theta^{(2)})^2}-\frac{1}{2\theta^{(2)}}

凸関数の微分は元の座標系に対して1つ求まり,また,同じものはない. このため, \etaを座標系として使うこともできる.これを双対座標という. これをプロットすると以下.

def θ2η(θ1, θ2):
    η1 = θ1/(2*θ2)
    η2 = -θ1**2/(2*θ2)**2 - 1/(2*θ2)
    return η1, η2

def η2θ(η1, η2 ):
    θ1 = -η1/(η2+η1**2)
    θ2 = -1/2 * 1/(η2+η1**2)
    return θ1, θ2
μ, σ2 = θ2p(θ1,θ2)
η1,η2 = θ2η(θ1, θ2)
plt.figure(figsize=(12,4))
plt.subplot(131)
plot_lattice(μ, σ2, "C0")
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")
plt.subplot(132)
plot_lattice(θ1, θ2, "C1")
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")
plt.subplot(133)
plot_lattice(η1, η2, "C2")
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095434p:plain

この座標系でも \psi(\theta)をプロットすると以下.

μ, σ2 = θ2p(θ1,θ2)
η1,η2 = θ2η(θ1, θ2)

plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
cont = plt.contour(μ, σ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
cont = plt.contour(θ1, θ2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
cont = plt.contour(η1, η2, ψ_θ)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095446p:plain

ここで,双対ポテンシャル関数は  \phi(\eta) = \theta \cdot \eta - \psi \left( \theta(\eta)\right )で与えられるので(?),  \phi(\eta) = -\frac{1}{2} \log(2\pi e) -\frac{1}{2} \log{ -(\eta^{(1)})^2 + \eta^{(2)}} これをプロットして以下.

これも双対ポテンシャル関数も凸関数であるので,その微分を座標系にすることもできる. これを求めると,元の \theta座標系になる.

def ϕ(η1, η2):
    return -1/2*log(2*π*e) - 1/2 * log(-(η1**2 + η2))
temp_η1 = np.linspace(-0.4, 0.4, 30)
temp_η2 = np.linspace(-1, -0.2, 30)

η1, η2 = np.meshgrid(temp_η1, temp_η2)


θ1,θ2 = η2θ(η1, η2)
μ, σ2 = θ2p(θ1,θ2)

ϕ_η = ϕ(η1, η2)
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_wireframe(η1, η2, ϕ_η)
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
ax.set_zlabel("$\phi(\eta)$")
plt.show()

f:id:ksknw:20180814095501p:plain

%matplotlib inline

plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
cont = plt.contour(μ, σ2, ϕ_η)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
cont = plt.contour(θ1, θ2, ϕ_η)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
cont = plt.contour(η1, η2, ϕ_η)
cont.clabel(fmt='%1.1f', fontsize=14)
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095517p:plain

測地線と双対測地線

測地線と双対測地線は,それぞれの座標系で直線として与えられる. つまり測地線は  \theta(t) = t \theta_1 + (1-t)\theta_2, 双対測地線は  \eta(t) = t \eta_1 + (1-t)\eta_2, である.

対応する点  \theta \etaについて,それぞれの測地線をプロットすると以下.

θ_1 = np.array([-1, 1])
θ_2 = np.array([1, 3])

list_t = np.arange(0, 1, 0.01)
line_θ = np.transpose([θ_1*t + θ_2*(1-t) for t in list_t])
line_p1 = θ2p(*line_θ)
line_η = θ2η(*line_θ)

temp_θ1 = np.linspace(-1,1, 30)
temp_θ2 = np.linspace(1, 3, 30)

θ1, θ2 = np.meshgrid(temp_θ1, temp_θ2)
μ, σ2 = θ2p(θ1, θ2)
η1,η2 = θ2η(θ1,θ2)


plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
plt.plot(line_p1[0], line_p1[1], "C4")
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
plt.plot(line_θ[0], line_θ[1], "C4")
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
plt.plot(line_η[0], line_η[1], "C4")
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095531p:plain

η_1 = np.array([-0.5, -3/4])
η_2 = np.array([ 1/6, -7/36])

list_t = np.arange(0, 1, 0.01)
line_η= np.transpose([η_1*t + η_2*(1-t) for t in list_t])
line_θ= η2θ(*line_η)
line_p2 = θ2p(*line_θ)
    
temp_η1 = np.linspace(-0.5, 0.5, 30)
temp_η2 = np.linspace(-1, -0.2, 30)

η1, η2 = np.meshgrid(temp_η1, temp_η2)
θ1,θ2 = η2θ(η1,η2)
μ, σ2 = θ2p(θ1, θ2)

    
plt.figure(figsize=(12,4))

plt.subplot(131)
plot_lattice(μ, σ2, "C0", alpha=0.3)
plt.plot(line_p2[0], line_p2[1], "C3")
plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")

plt.subplot(132)
plot_lattice(θ1, θ2, "C1", alpha=0.3)
plt.plot(line_θ[0], line_θ[1], "C3")
plt.xlabel("$\\theta^{(1)}$")
plt.ylabel("$\\theta^{(2)}$")

plt.subplot(133)
plot_lattice(η1, η2, "C2", alpha=0.3)
plt.plot(line_η[0], line_η[1], "C3")
plt.xlabel("$\eta^{(1)}$")
plt.ylabel("$\eta^{(2)}$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095542p:plain

これらを元のパラメータ空間でプロットすると以下のようになる.

temp_μ = np.linspace(-0.5,0.5, 30)
temp_σ2 = np.linspace(0.01, 1.0, 30)

μ, σ2 = np.meshgrid(temp_μ, temp_σ2)

plot_lattice(μ, σ2, "C0", alpha=0.3)
plt.plot(line_p1[0], line_p1[1], "C4")
plt.plot(line_p2[0], line_p2[1], "C3")

plt.text(-0.2, 0.5,"Dual geodesic line", size=14, color="C3")
plt.text(-0.4, 0.2,"Geodesic line", size=14, color="C4")

plt.xlabel("$\mu$")
plt.ylabel("$\sigma^2$")
plt.tight_layout()
plt.show()

f:id:ksknw:20180814095554p:plain

おわりに

とりあえず双対測地線のプロットまでをやった. 理解できてない部分もあるので,何か間違っているかもしれない.

この本は個人的には読みやすく感じるが,数学に詳しい人に聞くと, 先に接続とか平坦とかをやってから,測地線の方程式を求めたりしたほうがいいらしい (微分幾何を専攻したい人生だった). 実際天下り的に感じる部分もあるので,もうちょっと読めたら,ちゃんとそっちの順番で理解していきたいと思う.

参考

情報幾何学の新展開

機械学習とかの実験のパラメータをいい感じにやるコード

機械学習とかやっていると複数のパラメータを,コマンドラインから設定したいことがある. 割といい感じにできるようになった気がするので,備忘録として書いておく.

コマンドラインパーサとしてclickを使う. clickについては,こちらが詳しくていい感じだった.

この記事では以下ができるようにする.

  1. コマンドラインからいくつかのパラメータを設定する.
  2. パラメータを増やしたくなったときに,さくっとできるようにする.
  3. パラメータを名前に含んだディレクトリを作成して,その中に色々入れる.
import click
import datetime
import os
import pickle


def train(batch_size, learning_rate, nb_epochs, **kwargs):
    return 0, 0


def test(test_mode, **kwargs):
    print(test_mode)
    return 0


@click.command()
@click.option("-b", "--batch_size", type=int, default=12, help="mini batch size")
@click.option("-l", "--learning_rate", type=float, default=0.01)
@click.option("-n", "--nb_epochs", type=int, default=10, help="number of epochs")
@click.option("-m", "--test_mode", type=click.Choice(["hoge", "hogehoge"]))
@click.option("-o", "output_prefix", type=str, default="../results")
def main(**kwargs):
    net, loss_hist = train(**kwargs)
    test(**kwargs)

    output_prefix = kwargs.pop("output_prefix")
    output_dir = output_prefix + "/hoge_"
    now = datetime.datetime.now().strftime("%y%m%d%H%M%S")
    output_dir += now
    kwargs["output_dir"] = output_dir

    os.mkdir(output_dir)
    with open(output_dir + "/log.txt", "w") as f:
        for key, param in kwargs.items():
            print("{}, {}".format(key, param), file=f)
    with open(output_dir + "/params.pkl", "wb") as f:
        pickle.dump(kwargs, f)
    with open(output_dir + "/loss_hist.pkl", "wb") as f:
        pickle.dump(loss_hist, f)

    return


if __name__ == '__main__':
    main()

kwargsには辞書型が入ってくるので,あとはお好みで色々できる. 引数の最後に kwargsをdef train(batch_size, learning_rate, nb_epochs, **kwargs):な感じで入れておけば,必要な変数だけとってきて,残りはおいておくことができるので,いい感じに書けた.

独立成分分析の実装

はじめに

統計的因果探索を読んでいる. この中で,独立成分分析が出てくるが,そういえばあんまり理解できていなかったので,実装して理解する. 日本語の本も見たけど,結局 Independent Component Analysis: Algorithms and Applicationsがわかりやすかった.

独立成分分析

確率変数  \textbf{s} = {s_1, s_2, ..., s_N} が線形に"mix"された観測  \bf{x} が得られたとき,元の観測をどうにかして復元したい. 独立成分分析では,元の確率変数  \textbf{s}の要素が独立であるという仮定から,復元を試みる.

"mix"行列を Aとすると,

 \textbf{x} = A\textbf{s}

import pylab as plt
%matplotlib inlinee

import numpy as np
from scipy import signal

from sklearn.decomposition import PCA

import torch
from torch.autograd import Variable

np.random.seed(0)

n_samples = 5000
time = np.linspace(0, 8, n_samples)

s1 = np.sin(2 * time)  # Signal 1 : sinusoidal signal
s2 = np.sign(np.sin(3 * time))  # Signal 2 : square signal

S = np.c_[s1,  s2]

S /= S.std(axis=0)  # Standardize data
A = np.array([[1, 1], 
              [0.5, 2],
             ])  # Mixing matrix

X = np.dot(S, A.T)  # Generate observations

元の信号

plt.plot(S)
plt.title("S")
plt.show()

f:id:ksknw:20180519231649p:plain

観測値

plt.plot(X)
plt.title("X")
plt.show()

f:id:ksknw:20180519231705p:plain

独立成分分析では, \textbf{s}の各要素が独立であるという仮定から,復元行列  Wをいい感じに求めて,復元値 \textbf{y} = W\textbf{x}を決める.

中心極限定理と独立

ここで  \textbf{y} = [y_1, y_2]とすると,

 y_1 = \textbf{w}^\top A \textbf{s} ここで, \textbf{w} Wの行ベクトル

ここで, \textbf{z}=A\textbf{w} とすると  y_1 = \textbf{z}^\top \textbf{s}

ここから復元値  y_1, y_2 は独立な確率変数 \textbf{s}を重み付きで足したものであるといえる.

中心極限定理から独立な確率変数の和は"usually",ガウス分布に近づく. このため,復元がうまく行っていない場合には, \textbf{y} \textbf{s}に比べてガウス分布に近づく. 独立成分分析では,この性質を利用して \textbf{y}がよりガウス分布から離れるように Wを求める.

ガウス分布度合いを測る指標として尖度(kurtosis)がある.  \textbf{y} の分散共分散行列が単位行列のとき,  \mathrm{kurt} ( \textbf{y} ) = \mathbb{ E } [ y^4 ]-3

尖度はガウス分布に近いときに0となる. よって,尖度の絶対値を最大化するように Wを求めることで信号を復元する. 勾配法によって最適化する.ただし, Wの大きさを大きくすると,どこまでの誤差関数を小さくできてしまうので, Wのノルムを1に制約する. 今回は適当に毎ステップでノルムが1になるように割る

pca = PCA(n_components = X.shape[1], whiten=True)
X = pca.fit_transform(X)
def kurt(Y):
    return torch.mean(Y**4, dim=0) - 3
def kurt_loss(Y):
    # Y の分散共分散が単位行列であることが前提
    return -torch.mean(torch.abs(kurt(Y)))
def ICA(X, loss_func = kurt_loss):
    X = Variable(torch.from_numpy(X), requires_grad =False).float()
    W = Variable((torch.randn((2, 2))), requires_grad =True)
    W.data = W.data / torch.norm(W.data)
    hist_loss = []
    
    for i in range(1000):
        Y = torch.mm(X, W)
        loss = loss_func(Y)
        hist_loss.append(loss.data.numpy())
        loss.backward()

        W.data = W.data - W.grad.data * 0.1
        W.grad.data.zero_()
        W.data = W.data / torch.norm(W.data)
    return Y, W, hist_loss
Y, W, hist_loss = ICA(X, loss_func=kurt_loss)
plt.plot(Y.data.numpy())
plt.show()
plt.plot(hist_loss)

f:id:ksknw:20180519231733p:plain

f:id:ksknw:20180519231743p:plain

いい感じに復元することができた.

Negentropyを用いた誤差関数

今回はそれなりにうまく復元できたけど,適当に最適化していることもあってうまく収束しないことも多い. また,尖度は外れ値に対して敏感であることが知られている. そこで、ガウス分布度合いを測る別の指標、Negentropyを用いる.

 J(\textbf{y}) = H(\textbf{y}_\mathrm{gauss}) - H(\textbf{y})

ここで Hエントロピーガウス分布エントロピーが最も高い確率分布であるため,その差を使ってガウス分布度合いを評価する. しかし,結局,この計算は確率密度関数が必要なので難しい. そこで,以下のように近似的にNegentropyを求める.

 J(\textbf{y}) \approx \sum _{i=1}^p k_i { \mathbb{E} [G_i(y) ] - \mathbb{E} [G_i(\nu) ] }^2

ただし  \nu \sim \mathcal{N}(0,\mathbb{I}), k_i>0 . また, Gには色々な関数が設定できるが,例えば, G_1(u) = 1/a_1 \log \mathrm{cosh} a_1 uがある.

def negentropy_loss_logcosh(Y, a=1):
    return -torch.mean(torch.log(torch.cosh(a * Y))/a )**2 
Y, W, hist_loss = ICA(X, loss_func=negentropy_loss_logcosh)
plt.plot(Y.data.numpy())
plt.show()
plt.plot(hist_loss)

f:id:ksknw:20180519231759p:plain

f:id:ksknw:20180519231815p:plain

おわりに

独立成分分析の理解のために,自分で実装してみた. 今回は,天下り的に独立とガウス分布らしさの関係について考えたが,論文では同じ最適化問題相互情報量を最小化するという考えからも導出されている. また,今回は雑に最適化問題を解いたが,FastICAでは不動点反復法(?)で最適化するらしい. 制約付きの最適化問題なので,リーマン多様体上の最適化としても解ける気がする.

参考

リーマン多様体の最適化 (シュティーフェル多様体上の最適化まで)

概要

  • リーマン多様体上での最適化のお勉強をしている.
  • 多様体上での最適化では接空間上での勾配を求めたり,多様体上へレトラクションしたりする.
  • 例として球面上で,自明な最適化とレイリー商(対称行列の固有ベクトルを求める)をやる.
  • レイリー商から複数の固有値を同時に求めようとするとシュティーフェル多様体上での最適化になる.
  • シュティーフェル多様体はわかった気がするけど,グラスマン多様体がわからんのでできてない.
  • 実装は色んなライブラリがあるけど,pytorchを使う.autogradがあればなんでもいい.

はじめに

年末からお正月で多様体上の最適化について勉強している. 今回は自分の理解のために勉強したことをまとめるが,多様体がそもそもよくわかってないレベルなので色々許してください.

やりたいのはユークリッド空間上での制約つき最適化. 制約付きの最適化はラグランジュの未定乗数法とかで解くイメージだが,ニューラルネットとかがくっついているとどうやって解くのかわからない(あるのかもしれないが). 簡単な方法としてペナルティ項として制約条件を表し,制約なし最適化として解く方法が考えられるが,制約条件を満たす保証はなかったり,パラメータ調整が難しかったりする.

ユークリッド空間上での制約付きの最適化は,制約条件を満たす全ての点があるリーマン多様体上にあるとき,リーマン多様体上での制約なし問題として解くことができる. 例えば,以下は後で使うレイリー商(を変形した式)である.この制約付きの最適化は半径1の球面上の制約なし最適化問題として解釈できる.

 \mathrm{min}_{\bf x} \ {\bf x}^\top A {\bf x}\ \ s.t. \ \ ||x||_2=1

ただし{\bf x} \in \mathbb{R}^n ,\ \ A \in \mathbb{R}^{n\times n}

普通の勾配法に加わる操作としては以下の2つである.

  1. ユークリッド空間で求めた勾配 \nabla _x\  f(x)多様体の接空間に射影する( \mathrm{grad}\ f(x)).
  2. 更新後の点(x_t + \xi多様体からはみ出す)を多様体上に戻す(レトラクション  R_x (x_t+\xi)).

図にすると以下のような感じ.

f:id:ksknw:20180103100907p:plain

まだあんまり理解できてないので数学的な議論を避けて,具体例をやっていく. 最適化にもニュートン法とか共役勾配法とか色々あるけど,今回は単純な最急降下法だけを考える.(学習率も固定) pytorchを使ってfloatで計算しているので,数値誤差がそれなりにある.

こちらなどを参考にして勉強した(している). 数学的な話を含めてわかりやすいように思うのでおすすめ.

以下の内容は主に,これらに書かれていることを自分の理解した範疇で書きなおして,理解のためにpython (pytorch) で実装したというもの.

またこれが有名らしい(?)けど,あまり読めてない (英語苦手...)

3次元空間中の球面の例

まずは簡単な例として球面上の最適化を考える. 最適化する関数は以下のようなもの

\mathrm{argmax}_{\bf x \in \mathbb{R}^3} f({\bf x})={\bf x}[2] \ \ s.t.\ || {\bf x} ||_2=1

解は自明で{\bf x}=[0,0,1]^\topだが,これを多様体上の最急降下法で解くことにする.

ペナルティ項としての実装

リーマン多様体上での最適化を行う前に比較手法として,制約をペナルティ項として表すことを考える. つまり以下の関数を制約なし最適化問題として解く.

\mathrm{argmin}_{\bf x \in \mathbb{R}^3} - x[2] + \gamma |(1-||x||_2)|

コードは以下.

求まった解は

([  5.31967416e-05,   3.33191572e-05,   1.02169939e+00])

で0.02ぐらいずれていて,制約を満たせていない.

収束途中の遷移は以下. 青い球が制約条件で,オレンジの点が上に向かって進んでいる.

多様体からはみだしたり戻ったりを繰り返しながらガタガタしつつ進んでいく. 最急降下法で最適解に落とすには学習率を減衰させる必要があるが,今回はやってないので,こんな感じになる.

f:id:ksknw:20180102211837p:plain

これでもこの問題に対してはある程度の解は得られるが,今回はリーマン多様体上によって制約条件を常に(数値的には)厳密に満たしつつ最急降下法を行う.

多様体上の最適化

以下の最適化問題をリーマン多様体上での最適化として解く. つまり球面をS^{n-1}とするとき,目的関数を以下のように変形する.

\mathrm{argmax}_{\bf x \in S^{n-1}} f({\bf x})={\bf x}[2]

そもそもリーマン多様体が何かあまり理解できてないが, 多様体上の任意の点の接空間でリーマン計量\lt \cdot, \cdot> (内積的なもの)が定義できていればいい(たぶん). 内積が定義できれば接ベクトルの長さが定義できるので,測地線の長さがわかるということ(なのかな?)

ということで,まずは球面S^{n-1}とするときの内積を以下のように定義する.

\lt \xi,\eta> =(\xi^\top \eta)\ \ \xi,\eta\in T_xS^{n-1}\ \ x\in S^{n-1}

(これはいわゆる普通の内積だが,後で行列多様体を考える練習として順にやる.)

接空間上での勾配を考える. ユークリッド空間上での勾配は[0,0,1]^\topであるが,これは多様体の接空間上にはない.

リーマン多様体上での点pにおける勾配(\mathrm{grad}\ f)_pは以下を満たすものとして定義される.  \xiは任意の接ベクトル.

\lt (\mathrm{grad}\ f)_p, \xi>_p = \xi^\top \nabla f

なので,多様体上での勾配は以下のようになる.

\mathrm{grad}f\left({\bf x}\right)= \left(\mathbb{I}-{\bf xx}^\top  \right) \nabla f

\xi^\top x=0なのでこの勾配は上の定義を満たす.

次にレトラクションを考える.今回は球面上に点を戻せばいいので,レトラクションは単純にノルムでわるだけにする.

 R_x(\eta)=\frac{x+\eta}{||x+\eta||_2}

以上から,勾配法で解くプログラムは以下のようになる.

得られた解は以下. float32なのでこれぐらいの誤差は許してほしいところ.

array([  1.29403008e-04,  -1.01325495e-04,   1.00000000e+00], dtype=float32)

制約条件は途中の点でも全て同様に満たされている.

f:id:ksknw:20180102223542p:plain

レイリー商の例

レイリー商(を変形したもの)である以下の最適化問題を考える.

 \mathrm{min}_{\bf x} \ {\bf x}^\top A {\bf x}\ \ s.t. \ \ ||x||_2=1

これも上と同様にリーマン多様体(球面)上の最適化問題として解ける. ちなみにこの問題の解は Aの最小固有値に対応する固有ベクトルになる.

このとき勾配,レトラクションは上と同様.

\mathrm{grad}f\left({\bf x}\right)= \left(\mathbb{I}-{\bf xx}^\top  \right) \nabla f

 R_x(\eta)=\frac{x+\eta}{||x+\eta||_2}

プログラムは以下.

図中青がペナルティ法での点の軌跡.オレンジが多様体上で最適化した時の点の軌跡. 紫は解(固有ベクトル)を表している.

オレンジは多様体上を移動しているのに対して,青は制約条件を満たさない点を経由して最適解へと近づいていく. また,(学習率などの問題もあるが)今回の問題では多様体上での最適化のほうが圧倒的に収束が速かった.

f:id:ksknw:20180102230356p:plain

レイリー商(シュティーフェル多様体上での最急降下法)

ここまではベクトル{\bf x}の最適化について考えてきた. 次に行列 X\in \mathbb{R}^{n\times m}の最適化について考える.

解く問題は先ほどと同様に対称行列A固有ベクトルを求めるもので以下の最適化問題である.

 \mathrm{min}_{X\in \mathbb{R}^{n\times m}} \ \mathrm{trace} \ X^\top A X\ \ s.t. \ \ X^\top X=\mathbb{I}_n

ここでA固有値n個あるがここではそのうちm個を求めるときを考えている。つまりXは正方行列に限らない。

これを解きたいがトレースは行列を並び替えても値が変わらないのでこのままでは解けない。 そこで以下のBrockett cost functionを導入する.

 \mathrm{min}_{X\in \mathbb{R}^{n\times m}} \ \mathrm{trace} \ X^\top A X N \ \ s.t. \ \ X^\top X=\mathbb{I}_n

N=\mathrm{diag}(\mu_1,...,\mu_n),\ \ 0\lt\mu_1\lt...\lt \mu_n

Nによって重み付けされるので,固有値の大きさによって解が順序付けられた状態のみが最適解になる.

制約条件 X^\top X=\mathbb{I}_nを満たす行列からなる多様体をシュティーフェル多様体\mathrm{St(m,n)}という。 まずは 接空間上のリーマン計量を以下の標準内積で定義する。 \lt \xi,\eta>_X = \mathrm{trace} (\xi^\top \eta) \ \ \xi,\eta \in T_X\mathrm{St}(m,n)

(あんまりわかってないので中略) このとき,多様体上の勾配は以下(たぶん). \mathrm{grad}\ f(X) = \nabla f(X) - X \mathrm{symm}(X^\top  \nabla f(X) ), \ \ \mathrm{symm}(D) = 0.5 (D+D^\top)

次にレトラクションを考える. シュティーフェル多様体のレトラクションとして

  • 極分解
  • qr分解を使う方法

などがあるが今回はqr分解を使う方法を用いる.

以上よりプログラムは以下.

これにより2つの固有値を同時に求めることができた.

f:id:ksknw:20180102231107p:plain

おわりに

リーマン多様体上での最急降下をやった. ペナルティで実装するよりは良さそう. 接空間を求めてからどうやって接空間上の勾配を求めれば良いのかがまだわかっていない. あとはQR分解とかの計算コストがどれぐらいなのかを考えたほうが良いかもしれない.

今回はpytorchで実装したが,pymanopt (python manifold optimization) というプロジェクトがあり,このライブラリを使うとシュティーフェル多様体上での最適化などが簡単にできるっぽいので良さそう.(Tensorflowなどとも連携できる(?))

そのうち中略をちゃんと埋める. また,Brockett cost functionではなく,グラスマン多様体として最適化する方法もやっていきたい.

参考

JuliaとPythonで擬似コードや数式を動かそうとしたときの比較

この記事はJulia Advent Calendar 2017の17日目の記事です.

普段はpythonばかり書いていて,juliaは最近文法覚えてきたかなレベルなので色々許してください.

コードの全体はここにあります.

github.com

概要

  • この記事では擬似コードや数式を可能な限りそのままプログラムすることを目的とします.
  • unicode文字を使いまくって以下の画像のようなプログラムを作成します.
  • juliaとpythonで実装して,書きやすさと実行速度を比較します.
  • 書きやすさが悪化するので型指定はしません.
  • 結論は以下です.
    • juliaのほうが色んなunicode文字が使えるから,書きやすく可読性が高い.
    • インデックスが1から始まるのがいい.
    • juliaのほうが倍程度速くなることもあるけど,思ったより速くならない (型指定してないから)
    • juliaのeinsumを何も考えずに使うと激遅になる.(僕がやらかしてるかも)

unicode文字は以下の埋め込みではおそらく微妙な感じに表示されますが,ちゃんと設定されたエディタやjupyter notebookでは以下のように表示されます. f:id:ksknw:20171205223939p:plain

Table of Contents

はじめに (式or擬似コードに可能な限り近いプログラム)

数式を見て実装するときや論文に書いてある擬似コードを実装するとき,可能な限りそれらに近いプログラムを書くようにしている. 数式/擬似コードとプログラムを近づけるために,{ \displaystyle \alpha}などunicode文字を多用する. このようなプログラムを書くことで,

という利点がある.

例えばNo-U-Turn Sampler Hoffman+, 2011擬似コード(の一部,論文より引用)は

f:id:ksknw:20171130224045p:plain

これに対して書いたpythonコードは以下.

def BuildTree(θ, r, u, v, j, ε):
    if j == 0:
        θd, rd = Leapfrog(x, θ, r, v * ε)
        if np.log(u) <= (L(*θd) - 0.5 * np.dot(rd, rd)):
            Cd_ = [[θd, rd]]
        else:
            Cd_ = []
        sd = int(np.log(u) < (Δ_max + L(*θd) - 0.5 * np.dot(rd, rd)))
        return θd, rd, θd, rd, Cd_, sd
    else:
        θ_minus, r_minus, θ_plus, r_plus, Cd_, sd = BuildTree(θ, r, u, v, j - 1, ε)
        if v == -1:
            θ_minus, r_minus, _, _, Cdd_, sdd = BuildTree(θ_minus, r_minus, u, v, j - 1, ε)
        else:
            _, _, θ_plus, r_plus, Cdd_, sdd = BuildTree(θ_plus, r_plus, u, v, j - 1, ε)
        sd = sdd * sd * int((np.dot(θ_plus - θ_minus, r_minus) >= 0) and (np.dot(θ_plus - θ_minus, r_plus) >= 0))
        Cd_.extend(Cdd_)

        return θ_minus, r_minus, θ_plus, r_plus, Cd_, sd

pythonではギリシャ文字などの一部のunicode文字を変数に使うことができる. しかし,一部の例えば { \displaystyle
\theta ^+
}{ \displaystyle
\nabla
} などの記号は使うことができないため,微妙に見難い感じになってしまっていた. 探してみるとjuliaで同じことをしている人がいて,こちらのほうがだいぶ良さそうだった.

bicycle1885.hatenablog.com

juliaでは↑であげたような文字に加えて { \displaystyle
\hat a
} みたいな修飾文字や不等号の { \displaystyle
\le
} とかを使える. juliaで使えるunicode文字一覧はこちら

juliaすごい.ということで,以下ではこれまでpythonで実装したコードをjuliaで書きなおし,書きやすさと実行速度を比較する.

NUTSは既にやられているようなので

  • Dynamic Time Warping
  • Stochastic Block Model

について実装する.

juliaは型をちゃんと指定するともっと速くなるが,「そのまま実装する」という目的に反するのでやらない.

ちなみにNUTSのpython実装の話はこっち.

ksknw.hatenablog.com

Dynamic Time Warping

Dynamic Time Warping (DTW) は,2つの連続値系列データのいい感じの距離を求めるアルゴリズム動的計画法で解く. pythonで実装したときの記事はこっち ksknw.hatenablog.com

A global averaging method for dynamic time warping, with applications to clustering (DTW自体の論文ではない) の擬似コードを参考にしてプログラムを書く. 擬似コードは以下.

f:id:ksknw:20171130230716p:plain

f:id:ksknw:20171130230720p:plain

(ただしこの擬似コードは多分間違っているので,m[i,j]の遷移前のインデックスをsecondに入れるように変えた.)

python

δ = lambda a,b: (a - b)**2
first = lambda x: x[0]
second = lambda x: x[1]

def minVal(v1, v2, v3):
    if first(v1) <= min(first(v2), first(v3)):
        return v1, 0
    elif first(v2) <= first(v3):
        return v2, 1
    else:
        return v3, 2 

def dtw(A, B):
    S = len(A)
    T = len(B)

    m = [[0 for j in range(T)] for i in range(S)]
    m[0][0] = (δ(A[0],B[0]), (-1,-1))
    for i in range(1,S):
        m[i][0] = (m[i-1][0][0] + δ(A[i], B[0]), (i-1,0))
    for j in range(1,T):
        m[0][j] = (m[0][j-1][0] + δ(A[0], B[j]), (0,j-1))

    for i in range(1,S):
        for j in range(1,T):
            minimum, index = minVal(m[i-1][j], m[i][j-1], m[i-1][j-1])
            indexes = [(i-1,j), (i,j-1), (i-1,j-1)]
            m[i][j] = (first(minimum)+δ(A[i], B[j]), indexes[index])
    return m

擬似コードや数式ではインデックスを1から始めることが多いが,pythonは0からなので,頭の中でインデックスをずらしながらプログラムを書く必要がある. 他にはv1が添字っぽくない程度で,大体そのまま書けた.

julia

δ(a,b) = (a - b)^2
# first(x) = x[1] firstは元からあるのでいらない
second(x) = x[2]

function minVal(v₁, v₂, v₃)
    if first(v₁) ≤ minimum([first(v₂), first(v₃)])
        return v₁, 1
    elseif first(v₂) ≤ first(v₃)
        return v₂, 2
    else
        return v₃, 3
    end
end

function DTW(A, B)
    S = length(A)
    T = length(B)
    m = Matrix(S, T)
    m[1, 1] = [δ(A[1], B[1]), (0,0)]
    for i in 2:S
        m[i,1] = [m[i-1, 1][1] + δ(A[i], B[1]), [i-1, 1]]
    end
    for j in 2:T
        m[1,j] = [m[1, j-1][1] + δ(A[1], B[j]), [1, j-1]]
    end
    for i in 2:S
        for j in 2:T
            min, index = minVal(m[i-1,j], m[i,j-1], m[i-1,j-1])
            indexes = [[i-1, j], [i, j-1], [i-1, j-1]]
            m[i,j] = first(min) + δ(A[i],B[j]), indexes[index]
        end
    end
    return m
end

endがある分pythonより長い.一方でpythonでは使えない { \displaystyle
v_1}とか { \displaystyle
≤} が使えるので,より忠実な感じになっている.

実際に書いてみるとわかるけど,インデックスが擬似コードと同じなのは結構大事で,全然頭を使わずにそのまま書けた.

実行速度比較

juliaはpythonに比べて速いらしいので実行速度を比較した. シェルのtime を使って実行速度を比較した. コードの全体はここにある.

julia dtw_julia.jl  2.62s user 0.30s system 110% cpu 2.641 total
python dtw_python.py  2.76s user 0.11s system 99% cpu 2.873 total

期待していたよりも全然速くならなかった. 実行時間が短くてコンパイルのオーバヘッドが大きいのかなと思ったから,forで10回実行するようにした結果が以下.

julia dtw_julia.jl  21.97s user 0.66s system 101% cpu 22.355 total
python dtw_python.py  25.81s user 0.78s system 99% cpu 26.591 total

多少速い気がするけど,期待としては数十倍とかだったので,いまいち. 型指定していないとこれぐらいなのかな(?)

Stochastic Block Model

Stochastic Block Model (SBM)は非対称関係データのクラスタリング手法. 隠れ変数の事後分布をサンプリングして近似する(周辺化ギブスサンプラー). 今回の例では特にクラスタ割り当て { \displaystyle
z_1
} の事後分布を求めて,サンプリングする部分をやる.

pythonでの実装したときの記事はこっち. ksknw.hatenablog.com

関係データ学習という本に載っているクラスタzの事後確率に関する数式は以下. { \displaystyle
z_{1,i}}, { \displaystyle
z_{2,j}}をサンプリングする。

他の変数がgivenだとした時の、{ \displaystyle
z_{1,i}}の事後確率は、

ここで、

はガンマ関数で、 はパラメータ

この事後分布を求めてサンプリングする部分をpythonとjuliaで比較する. 全部書くと見づらいので一部だけ,コードの全体は同じくここにある.

python

n_pos = np.einsum("ikjl, ij", np.tensordot(z1, z2, axes=0), X)  # n_pos_kl = n_pos[k][l]
n_neg = np.einsum("ikjl, ij", np.tensordot(z1, z2, axes=0), 1 - X)

m1_hat = lambda i: m1 - z1[i]  # m1_hat_k = m1_hat[k]

n_pos_hat = lambda i: n_pos - np.einsum("kjl, j", np.tensordot(z1, z2, axes=0)[i], X[i])
n_neg_hat = lambda i: n_neg - np.einsum("kjl, j", np.tensordot(z1, z2, axes=0)[i], 1 - X[i])

α_1_hat = lambda i: α1 + m1_hat(i)
a_hat = lambda i: a0 + n_pos_hat(i)
b_hat = lambda i: b0 + n_neg_hat(i)

aᵢhat = a_hat(i)
bᵢhat = b_hat(i)

p_z1ᵢ_left = logΓ(aᵢhat + bᵢhat) - logΓ(aᵢhat) - logΓ(bᵢhat)
p_z1ᵢ_right_upper = logΓ(aᵢhat + np.dot(X[i], z2)) + logΓ(bᵢhat + np.dot((1 - X[i]), z2))
p_z1ᵢ_right_lower = logΓ(aᵢhat + bᵢhat + m2)
p_z1ᵢ = (α_1_hat(i) * exp(p_z1ᵢ_left + p_z1ᵢ_right_upper - p_z1ᵢ_right_lower)).prod(axis=1)
p_z1ᵢ = p_z1ᵢ.real
p_z1ᵢ = p_z1ᵢ / p_z1ᵢ.sum()
new_z1.append(onehot(choice(range(nb_k), p=p_z1ᵢ), nb_k))

数式には{ \displaystyle
\hat{a}
}{ \displaystyle
n^+
} などが頻出するが,pythonではこれらの文字を使うことができない. このため,m1_hatやn_posなどの変数名になってしまっている.

julia

@einsum n⁺[k,l] := X[i,j] * 𝕀z₁[i,k] * 𝕀z₂[j,l]
@einsum n⁻[k,l] := (ones(X)[i,j] - X[i,j]) * 𝕀z₁[i,k] * 𝕀z₂[j,l]

m̂₁ = m(𝕀z₁) - transpose(𝕀z₁[i,:])
@einsum Σ⁺x𝕀z₂[i,l] := X[i,j] * 𝕀z₂[j,l]
@einsum Σ⁻x𝕀z₂[i,l] := (ones(X)[i,j] - X[i,j]) * 𝕀z₂[j,l]
@einsum n̂⁺[k,l] := n⁺[k,l] - 𝕀z₁[i,k] * Σ⁺x𝕀z₂[i,l]
@einsum n̂⁻[k,l] := n⁻[k,l] - 𝕀z₁[i,k] * Σ⁻x𝕀z₂[i,l]

α̂₁ = α₁ + m̂₁
â = a₀ + n̂⁺
b̂ = b₀ + n̂⁻

temp⁺ = zeros(â)
temp⁻ = zeros(â)
temp = zeros(â)
for j in 1:size(temp⁺)[1]
    temp⁺[j,:] = Σ⁺x𝕀z₂[i,:]
    temp⁻[j,:] = Σ⁻x𝕀z₂[i,:]
    temp[j,:] = sum(𝕀z₂,1)
end

@einsum p_z₁[k,l] := exp(logΓ(â + b̂)-logΓ(â)-logΓ(b̂)
    + logΓ(â+temp⁺)+logΓ(b̂+temp⁻)-logΓ(â+b̂+temp))[k,l]
p_z₁ = α̂₁ .* transpose(prod(p_z₁, 2))
p_z₁ /= sum(p_z₁)

𝕀z₁[i,:] = onehot(sample(1:K, Weights(p_z₁[:])), K)

pythonに対してjuliaではα̂₁やn̂⁺などをそのまま変数名に使えて便利.

実行速度比較

python sbm_python.py  8.74s user 0.09s system 354% cpu 2.492 total
julia sbm_julia.jl  たぶん60000sぐらい(終わらない...)

(なにか間違っているかもしれないが,) なんとjuliaのほうが圧倒的に遅くなってしまった. というか同じ量の計算をさせると終わらない... 調べてみるとeinsumがだめっぽいので,以下ではforで書き直す. (とはいえさすがに遅すぎるので何かやらかしている気がする.)

einsumをfor文で書き直す

juliaのマクロ(@のやつ)はコンパイル時に別の式に展開されるらしい. einsumのマクロもそれぞれの@einsumに対して多重のfor文を展開しているらしく,しかも単独でもnumpyのeinsumより遅いこともあるらしい

ということで,普通にfor文による実装に変更して,実行速度を比較し直す. コードは以下.

m̂₁ = m(𝕀z₁) - transpose(𝕀z₁[i,:])
n⁺ = zeros(K, K)
n⁻ = zeros(K, K)

Σ⁺x𝕀z₂ = zeros(K)
Σ⁻x𝕀z₂ = zeros(K)
for l in 1:K
    for j in 1:N₂
        Σ⁺x𝕀z₂[l] += X[i,j] * 𝕀z₂[j,l]
        Σ⁻x𝕀z₂[l] += (1 - X[i,j]) * 𝕀z₂[j,l]
        for k in 1:K
            n⁺[k,l] += 𝕀z₁[i,k] * X[i,j] *  𝕀z₂[j,l]
            n⁻[k,l] += (1 - X[i,j]) * 𝕀z₁[i,k] * 𝕀z₂[j,l]
        end
    end
end

n̂⁺ = zeros(K,K)
n̂⁻ = zeros(K,K)
for k in 1:K
    for l in 1:K
        n̂⁺[k,l] += n⁺[k,l] - 𝕀z₁[i,k] * Σ⁺x𝕀z₂[l]
        n̂⁻[k,l] += n⁻[k,l] - 𝕀z₁[i,k] * Σ⁻x𝕀z₂[l]
    end
end
α̂₁ = α₁ + m̂₁
â = a₀ + n̂⁺
b̂ = b₀ + n̂⁻

p_z₁ = α̂₁
for k in 1:K
    for l in 1:K
        p_z₁[k] *= exp(logΓ(â[k,l] + b̂[k,l])-logΓ(â[k,l])-logΓ(b̂[k,l]) \
            + logΓ(â[k,l]+Σ⁺x𝕀z₂[l])+logΓ(b̂[k,l]+Σ⁻x𝕀z₂[l])-logΓ(â[k,l]+b̂[k,l]+sum(𝕀z₂,1)[l]))
    end
end
p_z₁ /= sum(p_z₁)

𝕀z₁[i,:] = onehot(sample(1:K, Weights(p_z₁[:])), K)

for...endのせいで見難く感じるが,これは普段pythonばかり書いていてfor文に拒絶反応があるからかもしれない. 書いている時の感想としては,行列の向きとか気にしなくていいので,einsumと同じ程度には書きやすい. 実際にfor文を全部とっぱらえば数式と近くなるはず.

実行速度比較

python sbm_python.py  8.74s user 0.09s system 354% cpu 2.492 total
julia sbm_julia_2.jl  3.62s user 0.46s system 114% cpu 3.559 total

無事にpythonの半分以下の時間で処理が終わるようになった. 普段python使っていると無意識的にfor文を避けていたのが, juliaなら何も気にせずに普通にfor文書きまくれば良さそう.

おわりに

DTWとSBMという2つのアルゴリズムについて,pythonとjuliaでの実装を比較した.

書きやすさの改善 実行速度の改善
DTW インデックスが1から.数字の添字 25.81s → 21.97s
SBM (einsum) α̂₁やn̂⁺をそのまま書ける 激遅 (やらかしているかも)
SBM (for) 8.74s → 3.62s

pythonをjuliaで書き直すことで

  • 数式の添字をそのまま書くことができ,書きやすさと可読性が向上した.
  • 素直にfor文を使えば実行速度は改善した.
  • juliaのeinsumは激遅だった.(僕がやらかしてるかも)
  • 型の指定をしていないからか,実行速度の改善は限定的だった.

また,やっていて気づいたこととして

  • jupyter notebook上で実行しつつ書くというやり方だと,コンパイル時間が以外と鬱陶しい
  • a\hat はOKだけど,\hat aはだめとか.\leはOKだけど,\leqqはだめとかunicode関連で微妙につまることがあった.

この記事ではやらなかったこととして以下がある.気が向いたらやろうと思っている.

  • juliaとpythonで型指定してそれぞれどれぐらい速くなるか.
  • pythonのopt-einsumとの比較

はてなだと表示が微妙だけど,jupyter とかでは綺麗に表示されて見やすいよ!(大事なこと)

追記 (2017/12/17)

Qiitaにて,juliaのコードについてのコメントをいただきました.ありがとうございます.

  • DTWについて,こちらに書かれているように(Matrixの要素の型推論とベクトルの代わりにタプルを使用)コードを変更すると,数十倍程度高速化します.
  • SBMについて,グローバル変数になっているパラメータにconstを指定するだけで,実行時間が90%程度になります.

qiita.com

もう少し自分でも高速化のために,あれこれやってみたいと思います.

参考

pytorchで Canonical Correlation Analysis (正準相関分析)の実装

はじめに

  • pytorchの練習も兼ねて,Canonical Correlation Analysis (正準相関分析)をpytorchを使って実装する.
  • 本当は分散共分散行列からなる行列の一般化固有値問題を解くが,今回は勾配法で解を求める.
  • pytorchのプログラムが間違っていないことを確認するためにscikit-learnでもやる.

Canonical Correlation Analysis (正準相関分析)

こちらの資料を参考にプログラムを書く.

主成分分析が多次元の値に対して,分散が大きい方向に射影するアルゴリズムなのに対して,正準相関分析では2つの多次元変数を射影先で相関が大きくなるように射影するアルゴリズムである.

多次元のデータ x\in\mathbb{R}^{T\times K_x}y\in\mathbb{R}^{T\times K_y} 間の正準相関を考える. ここでTはデータ数(系列データのときはデータ長),K_x, K_yx_i,y_iの次元を表す.

x,yの射影ベクトルをそれぞれ a\in\mathbb{R}^{K_x} , b\in\mathbb{R}^{K_y} とする. x,yの平均がそれぞれ0のとき,射影後の分散 S_{xx},S_{yy} と共分散 S_{xy} は, S_{xx} = (a^\top x) (a^\top x)^\top=a^\top V_{xx}a , S_{yy} = (b^\top y) (b^\top y)^\top=b^\top V_{yy}bS_{xy} = (a^\top x) (b^\top y)^\top=a^\top V_{xy}b

よって,射影後の相関\rho\rho = \frac{a^\top V_{xy}b}{\sqrt{a^\top V_{xx}a}\sqrt{b^\top V_{yy}b}}になる. ここで,aとかbを単純にk倍しても相関は変化しない. 最適化のときに解が複数あるのは面倒なので,制約条件a^\top V_{xx}a = b^\top V_{yy}b = 1を加える.

以上から以下を解くことで正準相関が求まる. \max_{a,b} a^\top V_{xy}b \ \mathrm{s.t.}\ a^\top V_{xx}a = b^\top V_{yy}b = 1

これはラグランジュの未定乗数法で解くことができて,一般化固有値問題になるが,今回はこの式のまま勾配法で解く.

生成データでのテスト

以下のように適当に生成したデータを使って,正準相関分析のプログラムをテストする.

以下では x,y,\mathrm{dummy}\in\mathbb{R}^{200 \times 2}の3つの信号を使う. x,yはもともと同じ信号に大きめのノイズが乗ったもの.dummyは適当なノイズだけのもの. x,y間では正準相関の値が大きくなり,x,dummy間では正準相関の値が小さくなってほしい.

import pylab as plt
import numpy as np
import torch 
from torch.autograd import Variable
import logging
import coloredlogs
from tqdm import tqdm

logger = logging.getLogger(__name__)
coloredlogs.install(
    #fmt="%(asctime)s %(levelname)s %(message)s"
    fmt="%(levelname)s %(message)s"
)

%matplotlib inline

logger.setLevel("DEBUG")
if logger.level >= logging.WARNING:
    tqdm = lambda x:x
u1 = np.r_[np.zeros(40), np.ones(30), np.ones(50)*2, np.ones(30), np.zeros(20), np.ones(30) ]
u2 = np.r_[np.ones(25), np.zeros(25), np.ones(25)*2, np.ones(55), np.zeros(35), np.ones(35) ]
u = np.c_[u1, u2]
plt.plot(u)
plt.show()

f:id:ksknw:20171213230955p:plain

sx = np.random.random(len(u)) * 10
sy = np.random.random(len(u)) * 10
sd = np.random.random(len(u)) * 10
x = np.c_[
    u1 + sx + np.random.random(len(u)),
    u2 - sx + np.random.random(len(u)),   
]
plt.plot(x)
plt.show()

f:id:ksknw:20171213231026p:plain

y = np.c_[
    u1 + sy + np.random.random(len(u)),
    u2 - sy + np.random.random(len(u)),   
]
plt.plot(y)
plt.show()

f:id:ksknw:20171213231039p:plain

dummy = np.c_[
    sd + np.random.random(len(u)),
    -sd + np.random.random(len(u)),   
]
plt.plot(dummy)
plt.show()

f:id:ksknw:20171213231050p:plain

以上のように,見た目では全部ぐちゃぐちゃなノイズにしか見えない3つの信号を生成した. 見た目では全部バラバラだが,x,yに関しては"いい感じの方向"に写像できればノイズを除去できる. x,yそれぞれだけでは"いい感じの方向"がわからないのでノイズと信号を区別できないが,両方を使って相関が大きくなる方向をいい感じの方向と定義することで, 相関する信号を抽出できる.

Scikit-Learnによる実装

sciikit-learnには既に実装されているので,以下のように簡単に正準相関と射影後の変数を求めることができる.

from sklearn.cross_decomposition import CCA
cca = CCA(n_components=1)

x,y間の正準相関

t_x, t_y = cca.fit_transform(x,y)
plt.scatter(t_x, t_y)
print(np.corrcoef(t_x.T, t_y.T )[0,1])
plt.show()

正準相関の値

0.906462806701

射影した後のデータの分布

f:id:ksknw:20171213231110p:plain

x,\mathrm{dummy}間の正準相関

t_x, t_d = cca.fit_transform(x,dummy)
plt.scatter(t_x, t_d)
print(np.corrcoef(t_x.T, t_y.T )[0,1])
plt.show()

正準相関の値

0.141422242479

射影した後のデータの分布

f:id:ksknw:20171213231121p:plain

これと同じ結果を得ることを目標にpytorchで実装する.

PyTorchによる実装

pytorchにはtorch.eigという固有値を求める関数が用意されているので,これを頑張って使えば解けそうではある. が,諸々の事情からペナルティ関数法で制約条件を加えた誤差関数を設定し,autogradで微分を求めて勾配法で解く.

元の最適化問題はこれだった. \max_{a,b} a^\top V_{xy}b \ \mathrm{s.t.}\ a^\top V_{xx}a = b^\top V_{yy}b = 1

制約条件をペナルティ関数にして,以下のような誤差関数を最適化する. \mathrm{argmin}_{a,b} \mathcal{L} = -a^\top V_{xy}b + \gamma \left| 1 - a^\top V_{xx}a \right| + \gamma \left| 1 - b^\top V_{yy}b \right|

\gammaは適当に0.5にした.小さい値だと制約があまり満たされず,逆に大きい値とかにするとうまく収束しないことが多い.この辺がよくわからないので,「機械学習のための連続最適化」という本を買ったけどまだ読めていない.

optimizerは何でもいいけど,単純にSGDにした.pytorch等で実装してると,なんとなくでoptimizerを選べるのでいいなと思う.

x -= np.mean(x, axis=0)
y -= np.mean(y, axis=0)
dummy -= np.mean(dummy, axis=0)
N = len(x)

def CCA(x,y,
        nb_steps = 100000,
        lr = 0.001,
        γ = 0.5):
    vxx = np.dot(x.transpose(), x)/N
    vyy = np.dot(y.transpose(), y)/N
    vxy = np.dot(x.transpose(), y)/N
    
    X = Variable(torch.from_numpy(x)).float()
    Y = Variable(torch.from_numpy(y)).float()
    
    a = Variable(torch.randn(2,1), requires_grad=True)
    b = Variable(torch.randn(2,1), requires_grad=True)
    
    Vxx = Variable(torch.from_numpy(vxx)).float()
    Vyy = Variable(torch.from_numpy(vyy)).float()
    Vxy = Variable(torch.from_numpy(vxy)).float()
    
    hist = []
    from torch import optim
    optimizer = optim.SGD([a, b], lr = 0.0001)

    for i in tqdm(range(nb_steps)):
        Loss = - torch.mm(torch.mm(torch.t(a), Vxy), b)  \
             + γ * torch.abs(1 - torch.mm(torch.mm(torch.t(a),Vxx),a))  \
             + γ * torch.abs(1 - torch.mm(torch.mm(torch.t(b),Vyy),b)) 
        hist.append(Loss.data.numpy()[0,0])
        optimizer.zero_grad()
        Loss.backward()
        optimizer.step()
        
    x_t = np.dot(a.data.numpy().transpose(), x.transpose())
    y_t = np.dot(b.data.numpy().transpose(), y.transpose())
    cca_score = torch.mm(torch.mm(torch.t(a), Vxy), b)
    
    logger.info("a Vxx a: %s"%torch.mm(torch.mm(torch.t(b),Vyy),b).data.numpy())
    logger.info("b Vyy b: %s"%torch.mm(torch.mm(torch.t(a),Vxx),a).data.numpy())
    return hist, x_t, y_t, cca_score

x,y間の正準相関

hist, x_t, y_t, cca_score = CCA(x,y)
plt.plot(hist)
plt.show()
100%|██████████| 100000/100000 [00:19<00:00, 5232.01it/s]
INFO a Vxx a: [[ 1.00017715]]
INFO b Vyy b: [[ 1.00021267]]

lossの履歴

f:id:ksknw:20171213231145p:plain

plt.scatter(x_t, y_t)
plt.show()
print(cca_score)

正準相関の値

Variable containing:
 0.9066
[torch.FloatTensor of size 1x1]

射影した後のデータの分布

f:id:ksknw:20171213231156p:plain

hist, x_t, y_t, cca_score = CCA(x,dummy)
plt.plot(hist)
plt.show()
100%|██████████| 100000/100000 [00:19<00:00, 5199.69it/s]
INFO a Vxx a: [[ 1.00273335]]
INFO b Vyy b: [[ 0.99813366]]

lossの履歴

f:id:ksknw:20171213231211p:plain

x,\mathrm{dummy}間の正準相関

plt.scatter(x_t, y_t)
plt.show()
cca_score

正準相関の値

Variable containing:
1.00000e-02 *
  9.6734
[torch.FloatTensor of size 1x1]

射影した後のデータの分布

f:id:ksknw:20171213231220p:plain

無事にpytorchによる実装で,ある程度正しそうな結果を得ることができた.

おわりに

pytorchで勾配法を使ったCCAの実装をやった. 制約付き最適化を勾配法で解く良いやり方がいまいちわかってないので,勉強したい.

参考