確率的PCA

TensorFlow.orgで表示 GoogleColabで実行 GitHubでソースを表示ノートブックをダウンロードする

確率的主成分分析(PCA)は、低次元の潜在空間を介してデータを分析する次元削減手法です( Tipping and Bishop 1999 )。これは、データに欠落値がある場合、または多次元尺度構成法でよく使用されます。

輸入

import functools
import warnings

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

import tensorflow.compat.v2 as tf
import tensorflow_probability as tfp

from tensorflow_probability import bijectors as tfb
from tensorflow_probability import distributions as tfd

tf.enable_v2_behavior()

plt.style.use("ggplot")
warnings.filterwarnings('ignore')

モデル

$ N $データポイントのデータセット$ \ mathbf {X} = {\ mathbf {x} _n} $について考えます。ここで、各データポイントは$ D $次元、$ \ mathbf {x} _n \ in \ mathbb { R} ^ D $。各$ \ mathbf {x} _n $を、潜在変数$ \ mathbf {z} _n \ in \ mathbb {R} ^ K $の下で、より低い次元、$ K <D $で表すことを目指しています。主軸のセット$ \ 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*} $$

古典的PCAは、ノイズの共分散が非常に小さくなり、$ \ sigma ^ 2 \ to 0 $になる確率的PCAの特定のケースです。

以下にモデルを設定します。私たちの分析では、$ \ sigma $が既知であると仮定し、モデルパラメーターとして$ \ mathbf {W} $を点推定する代わりに、主軸上の分布を推測するために、その上に事前分布を配置します。モデルをTFPJointDistributionとして表現します。具体的には、 JointDistributionCoroutineAutoBatchedを使用します

def probabilistic_pca(data_dim, latent_dim, num_datapoints, stddv_datapoints):
  w = yield tfd.Normal(loc=tf.zeros([data_dim, latent_dim]),
                 scale=2.0 * tf.ones([data_dim, latent_dim]),
                 name="w")
  z = yield tfd.Normal(loc=tf.zeros([latent_dim, num_datapoints]),
                 scale=tf.ones([latent_dim, num_datapoints]),
                 name="z")
  x = yield tfd.Normal(loc=tf.matmul(w, z),
                       scale=stddv_datapoints,
                       name="x")
num_datapoints = 5000
data_dim = 2
latent_dim = 1
stddv_datapoints = 0.5

concrete_ppca_model = functools.partial(probabilistic_pca,
    data_dim=data_dim,
    latent_dim=latent_dim,
    num_datapoints=num_datapoints,
    stddv_datapoints=stddv_datapoints)

model = tfd.JointDistributionCoroutineAutoBatched(concrete_ppca_model)

データ

モデルを使用して、共同事前分布からサンプリングすることでデータを生成できます。

actual_w, actual_z, x_train = model.sample()

print("Principal axes:")
print(actual_w)
Principal axes:
tf.Tensor(
[[ 2.2801023]
 [-1.1619819]], shape=(2, 1), dtype=float32)

データセットを視覚化します。

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()

png

最大事後推論

最初に、事後確率密度を最大化する潜在変数の点推定を検索します。これは最大事後(MAP)推論として知られており、事後密度$ p(\ mathbf {W}、\を最大化する$ \ mathbf {W} $と$ \ mathbf {Z} $の値を計算することによって行われます。 mathbf {Z} \ mid \ mathbf {X})\ property p(\ mathbf {W}、\ mathbf {Z}、\ mathbf {X})$。

w = tf.Variable(tf.random.normal([data_dim, latent_dim]))
z = tf.Variable(tf.random.normal([latent_dim, num_datapoints]))

target_log_prob_fn = lambda w, z: model.log_prob((w, z, x_train))
losses = tfp.math.minimize(
    lambda: -target_log_prob_fn(w, z),
    optimizer=tf.optimizers.Adam(learning_rate=0.05),
    num_steps=200)
plt.plot(losses)
[<matplotlib.lines.Line2D at 0x7f19897a42e8>]

png

モデルを使用して、$ \ mathbf {W} $と$ \ mathbf {Z} $の推定値のデータをサンプリングし、条件付けした実際のデータセットと比較できます。

print("MAP-estimated axes:")
print(w)

_, _, x_generated = model.sample(value=(w, z, None))

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:
<tf.Variable 'Variable:0' shape=(2, 1) dtype=float32, numpy=
array([[ 2.9135954],
       [-1.4826864]], dtype=float32)>

png

変分推論

MAPは、事後分布のモード(またはモードの1つ)を見つけるために使用できますが、それに関する他の洞察は提供しません。次に、変分推論を使用します。ここで、事後分布$ p(\ mathbf {W}、\ mathbf {Z} \ mid \ mathbf {X})$は、変分分布$ q(\ mathbf {W}、\ mathbf {Z})$は$ \ boldsymbol {\ lambda} $によってパラメータ化されます。目的は、qと後部の間のKL発散を最小化する変分パラメーター$ \ boldsymbol {\ lambda} $を見つけることです。 \ミッドP(\ mathbf {W}、\ mathbf {Z} \中間\ mathbf {X})を最大証拠が下限こと)$、または等価的に、$ \ 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] $。

qw_mean = tf.Variable(tf.random.normal([data_dim, latent_dim]))
qz_mean = tf.Variable(tf.random.normal([latent_dim, num_datapoints]))
qw_stddv = tfp.util.TransformedVariable(1e-4 * tf.ones([data_dim, latent_dim]),
                                        bijector=tfb.Softplus())
qz_stddv = tfp.util.TransformedVariable(
    1e-4 * tf.ones([latent_dim, num_datapoints]),
    bijector=tfb.Softplus())
def factored_normal_variational_model():
  qw = yield tfd.Normal(loc=qw_mean, scale=qw_stddv, name="qw")
  qz = yield tfd.Normal(loc=qz_mean, scale=qz_stddv, name="qz")

surrogate_posterior = tfd.JointDistributionCoroutineAutoBatched(
    factored_normal_variational_model)

losses = tfp.vi.fit_surrogate_posterior(
    target_log_prob_fn,
    surrogate_posterior=surrogate_posterior,
    optimizer=tf.optimizers.Adam(learning_rate=0.05),
    num_steps=200)
print("Inferred axes:")
print(qw_mean)
print("Standard Deviation:")
print(qw_stddv)

plt.plot(losses)
plt.show()
Inferred axes:
<tf.Variable 'Variable:0' shape=(2, 1) dtype=float32, numpy=
array([[ 2.4168603],
       [-1.2236133]], dtype=float32)>
Standard Deviation:
<TransformedVariable: dtype=float32, shape=[2, 1], fn="softplus", numpy=
array([[0.0042499 ],
       [0.00598824]], dtype=float32)>

png

posterior_samples = surrogate_posterior.sample(50)
_, _, x_generated = model.sample(value=(posterior_samples))

# It's a pain to plot all 5000 points for each of our 50 posterior samples, so
# let's subsample to get the gist of the distribution.
x_generated = tf.reshape(tf.transpose(x_generated, [1, 0, 2]), (2, -1))[:, ::47]

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()

png

謝辞

このチュートリアルは、もともとEdward 1.0(ソース)で書かれました。そのバージョンの作成と改訂に貢献してくれたすべての人に感謝します。

参考文献

[1]:Michael E.TippingとChristopherM.Bishop。確率的主成分分析。 Journal of the Royal Statistics Society:Series B(Statistical Methodology) 、61(3):611-622、1999。