TensorFlow Probability : Tutorials : 確率的 PCA (翻訳/解説)
翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 11/19/2017
* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
TensorFlow Probability : 確率的 PCA
確率的主成分分析 (PCA) は低次元潜在空間を通してデータを解析する次元削減テクニックです (Tipping and Bishop 1999)。それはデータに欠落値があるときやマルチ次元スケーリングのためにしばしば使用されます。
インポート
!pip install -q tensorflow-probability
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import tensorflow as tf
import tensorflow_probability as tfp
from tensorflow_probability import edward2 as ed
import warnings
plt.style.use("ggplot")
warnings.filterwarnings('ignore')
モデル
$N$ データポイントのデータ集合 $\mathbf{X} = \{\mathbf{x}_n\}$ を考えます、そこでは各データポイントは $D$-次元で、$\mathbf{x}_n \in \mathbb{R}^D$ です。各 $\mathbf{x}_n$ をより低次元 $K < D$ である潜在変数 $\mathbf{z}_n \in \mathbb{R}^K$ で表わすことを狙います。主軸の集合 $\mathbf{W}$ は潜在変数をデータに関係づけます。
特に、各潜在変数は正規分布に従うものと仮定します :
$$
\begin{equation*}
\mathbf{z}_n \sim N(\mathbf{0}, \mathbf{I}).
\end{equation*}
$$
対応するデータポイントは射影を通して生成されます :
$$
\begin{equation*}
\mathbf{x}_n \mid \mathbf{z}_n
\sim N(\mathbf{W}\mathbf{z}_n, \sigma^2\mathbf{I}),
\end{equation*}
$$
ここで行列 $\mathbf{W}\in\mathbb{R}^{D\times K}$ は主軸として知られています。確率的 PCA では、典型的には主軸 $\mathbf{W}$ とノイズ項 $\sigma^2$ を推定することに興味があります。
確率的 PCA は古典的 PCA を一般化します。潜在変数を周辺化すると、各データポイントの分布は :
$$
\begin{equation*}
\mathbf{x}_n \sim N(\mathbf{0}, \mathbf{W}\mathbf{W}^\top + \sigma^2\mathbf{I}).
\end{equation*}
$$
ノイズの共分散が無限に小さくなるとき、$\sigma^2 \to 0$、古典的 PCA は確率的 PCA の特定のケースです。
下でモデルをセットアップします。私達の分析では、$\sigma$ は既知と仮定し、モデルパラメータとして $\mathbf{W}$ を点推定する代わりに、主軸に渡る分布を推論するためにそれに渡る事前分布を置きます。
def probabilistic_pca(data_dim, latent_dim, num_datapoints, stddv_datapoints): # (unmodeled) data
w = ed.Normal(loc=tf.zeros([data_dim, latent_dim]),
scale=2.0 * tf.ones([data_dim, latent_dim]),
name="w") # parameter
z = ed.Normal(loc=tf.zeros([latent_dim, num_datapoints]),
scale=tf.ones([latent_dim, num_datapoints]),
name="z") # parameter
x = ed.Normal(loc=tf.matmul(w, z),
scale=stddv_datapoints * tf.ones([data_dim, num_datapoints]),
name="x") # (modeled) data
return x, (w, z)
log_joint = ed.make_log_joint_fn(probabilistic_pca)
データ
Edward2 モデルをデータを生成するために使用できます。
num_datapoints = 5000
data_dim = 2
latent_dim = 1
stddv_datapoints = 0.5
model = probabilistic_pca(data_dim=data_dim,
latent_dim=latent_dim,
num_datapoints=num_datapoints,
stddv_datapoints=stddv_datapoints)
with tf.Session() as sess:
x_train, (actual_w, actual_z) = sess.run(model)
print("Principal axes:")
print(actual_w)
Principal axes: [[-0.10648464] [ 1.3358579 ]]
データセットを可視化します。
plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1)
plt.axis([-20, 20, -20, 20])
plt.title("Data set")
plt.show()
最大事後確率推論
最初に事後確率密度を最大化する潜在変数の点推定を探索します。これは最大事後確率 (MAP) 推論として知られ、事後密度 $p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X}) \propto p(\mathbf{W}, \mathbf{Z}, \mathbf{X})$ を最大化する$\mathbf{W}$ と $\mathbf{Z}$ の値を計算することで成されます。
tf.reset_default_graph()
w = tf.Variable(np.ones([data_dim, latent_dim]), dtype=tf.float32)
z = tf.Variable(np.ones([latent_dim, num_datapoints]), dtype=tf.float32)
def target(w, z):
"""Unnormalized target density as a function of the parameters."""
return log_joint(data_dim=data_dim,
latent_dim=latent_dim,
num_datapoints=num_datapoints,
stddv_datapoints=stddv_datapoints,
w=w, z=z, x=x_train)
energy = -target(w, z)
optimizer = tf.train.AdamOptimizer(learning_rate=0.05)
train = optimizer.minimize(energy)
init = tf.global_variables_initializer()
t = []
num_epochs = 200
with tf.Session() as sess:
sess.run(init)
for i in range(num_epochs):
sess.run(train)
if i % 5 == 0:
cE, cw, cz = sess.run([energy, w, z])
t.append(cE)
w_inferred_map = sess.run(w)
z_inferred_map = sess.run(z)
x = range(1, num_epochs, 5) plt.plot(x, t)
$\mathbf{W}$ と $\mathbf{Z}$ のために推論された値のためのデータをサンプリングするために Edward2 モデルを使用し、そして (その上で条件付けした) 実際のデータセットと比較することができます。
print("MAP-estimated axes:")
print(w_inferred_map)
def replace_latents(w=actual_w, z=actual_z):
def interceptor(rv_constructor, *rv_args, **rv_kwargs):
"""Replaces the priors with actual values to generate samples from."""
name = rv_kwargs.pop("name")
if name == "w":
rv_kwargs["value"] = w
elif name == "z":
rv_kwargs["value"] = z
return rv_constructor(*rv_args, **rv_kwargs)
return interceptor
with ed.interception(replace_latents(w_inferred_map, z_inferred_map)):
generate = probabilistic_pca(
data_dim=data_dim, latent_dim=latent_dim,
num_datapoints=num_datapoints, stddv_datapoints=stddv_datapoints)
with tf.Session() as sess:
x_generated, _ = sess.run(generate)
plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1, label='Actual data')
plt.scatter(x_generated[0, :], x_generated[1, :], color='red', alpha=0.1, label='Simulated data (MAP)')
plt.legend()
plt.axis([-20, 20, -20, 20])
plt.show()
MAP-estimated axes: [[-0.20085734] [ 2.7221382 ]]
変分推論
MAP は事後分布のモード (or モードの一つ) を見つけるために使用できますが、それについて他の洞察は一切提供しません。次に、変分推論を使用します、そこでは事後分布 $p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X})$ は $\boldsymbol{\lambda}$ によりパラメータ化された変分分布 $q(\mathbf{W}, \mathbf{Z})$ を使用して近似されます。目的は、q と事後分布 $\mathrm{KL}(q(\mathbf{W}, \mathbf{Z}) \mid\mid p(\mathbf{W}, \mathbf{Z} \mid \mathbf{X}))$ の間の KL ダイバージェンスを最小化する、あるいは同値にエビデンス下界 $\mathbb{E}_{q(\mathbf{W},\mathbf{Z};\boldsymbol{\lambda})}\left[ \log p(\mathbf{W},\mathbf{Z},\mathbf{X}) – \log q(\mathbf{W},\mathbf{Z}; \boldsymbol{\lambda}) \right]$ を最大化する、変分パラメータ $\boldsymbol{\lambda}$ を見つけることです。
tf.reset_default_graph()
def variational_model(qw_mean, qw_stddv, qz_mean, qz_stddv):
qw = ed.Normal(loc=qw_mean, scale=qw_stddv, name="qw")
qz = ed.Normal(loc=qz_mean, scale=qz_stddv, name="qz")
return qw, qz
log_q = ed.make_log_joint_fn(variational_model)
def target_q(qw, qz):
return log_q(qw_mean=qw_mean, qw_stddv=qw_stddv,
qz_mean=qz_mean, qz_stddv=qz_stddv,
qw=qw, qz=qz)
qw_mean = tf.Variable(np.ones([data_dim, latent_dim]), dtype=tf.float32)
qz_mean = tf.Variable(np.ones([latent_dim, num_datapoints]), dtype=tf.float32)
qw_stddv = tf.nn.softplus(tf.Variable(-4 * np.ones([data_dim, latent_dim]), dtype=tf.float32))
qz_stddv = tf.nn.softplus(tf.Variable(-4 * np.ones([latent_dim, num_datapoints]), dtype=tf.float32))
qw, qz = variational_model(qw_mean=qw_mean, qw_stddv=qw_stddv,
qz_mean=qz_mean, qz_stddv=qz_stddv)
energy = target(qw, qz)
entropy = -target_q(qw, qz)
elbo = energy + entropy
optimizer = tf.train.AdamOptimizer(learning_rate = 0.05)
train = optimizer.minimize(-elbo)
init = tf.global_variables_initializer()
t = []
num_epochs = 100
with tf.Session() as sess:
sess.run(init)
for i in range(num_epochs):
sess.run(train)
if i % 5 == 0:
t.append(sess.run([elbo]))
w_mean_inferred = sess.run(qw_mean)
w_stddv_inferred = sess.run(qw_stddv)
z_mean_inferred = sess.run(qz_mean)
z_stddv_inferred = sess.run(qz_stddv)
print("Inferred axes:")
print(w_mean_inferred)
print("Standard Deviation:")
print(w_stddv_inferred)
plt.plot(range(1, num_epochs, 5), t)
plt.show()
with ed.interception(replace_latents(w_mean_inferred, z_mean_inferred)):
generate = probabilistic_pca(
data_dim=data_dim, latent_dim=latent_dim,
num_datapoints=num_datapoints, stddv_datapoints=stddv_datapoints)
with tf.Session() as sess:
x_generated, _ = sess.run(generate)
plt.scatter(x_train[0, :], x_train[1, :], color='blue', alpha=0.1, label='Actual data')
plt.scatter(x_generated[0, :], x_generated[1, :], color='red', alpha=0.1, label='Simulated data (VI)')
plt.legend()
plt.axis([-20, 20, -20, 20])
plt.show()
References
- [1]: Michael E. Tipping and Christopher M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 61(3): 611-622, 1999.
以上