はじめに
多次元空間のデータを二次元平面で可視化した結果を評価する指標である、k3n error (k-nearest neighbor normalized error for visualization and reconstruction) というものを見つけたので試してみた。
k3n errorとは
こちらによると、以下の2つを同時に評価する指標である。
① 多次元空間において近いサンプル同士は、2次元平面においても近いか？
② 2次元平面において近いサンプル同士は、多次元空間においても近いか？
プログラムは、こちらで公開されている。
k3n errorの使い方
以下のようにインポートして実行する。
from k3nerror import k3nerror
X = .. # 可視化前の多次元空間のデータ
Z = .. # 可視化後の2次元平面のデータ
k3n_Z_error = k3nerror(X, Z, k_in_k3nerror)
k3n_X_error = k3nerror(Z, X, k_in_k3nerror)
ソース内の解説を読んだところ、k3nerror(X, Z, k_in_k3nerror)で①、knerror(Z, X, k_in_k3nerror)で②を評価するための計算が行われ、最終的に①＋②の値が、k3nerrorとなるようである。
やってみよう
今回以下のデータを、PCA, MDS, tSNE, UMAP, GTMで可視化した際の可視化指標を求めてみる。
- ケモインフォマティクスのデータということで、いつものごとくRDKitに付属のこちらのデータを利用
- RDkitを利用し、化学構造から説明変数を生成(記述子計算)
- 記述子計算して得られた123個の説明変数から、計算エラーが含まれる説明変数、分散が0になる説明変数、相関係数が0.9以上の説明変数(どちらかを削除)を削除した81の説明変数を利用
ソース
可視化および指標の計算を行うプログラムは以下の通りである。実行には、各可視化ライブラリに加え、GPy, GPyOptが必要となる。GTMはこちらのプログラムを利用し、ハイパーパラメータチューニングは、 こちらを流用させてもらった。GTM以外はパラメータチューニングはしなかった。
deduct_dmension.py
import numpy as np
import pandas as pd
import argparse
import csv
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE, MDS
import umap
from k3nerror import k3nerror
from gtm import GTM
import GPy, GPyOpt
def output_csv(path, ids, xs, ys):
with open(path, 'w') as f:
writer = csv.writer(f, delimiter=",", lineterminator="\n")
for id, x, y in zip(ids, xs, ys):
id = id.replace('\uffb4', '')
writer.writerow([id, x, y])
f.close()
def main():
parser = argparse.ArgumentParser()
parser.add_argument("-input", type=str, required=True, help="input csv file.")
parser.add_argument("-method", type=str, required=True, choices=["PCA", "tSNE", "MDS", "UMAP", "GTM"])
parser.add_argument("-output", type=str, required=True, help="output csv file.")
args = parser.parse_args()
# データの読み込み
df = pd.read_csv(args.input, index_col=0)
# オートスケーリング
ss = StandardScaler()
X = ss.fit_transform(df.values)
k_in_k3nerror = 10
Z = None
if "PCA" in args.method:
pca = PCA(n_components=2)
Z = pca.fit_transform(X)
if "MDS" in args.method:
mds = MDS(n_jobs=4)
Z = mds.fit_transform(X)
if "tSNE" in args.method:
tsne_model = TSNE(n_components=2)
Z = tsne_model.fit_transform(X)
if "UMAP" in args.method:
Z = umap.UMAP().fit_transform(X)
if "GTM" in args.method:
# 変数探索範囲
bounds = [
#{'name': 'shape_of_map_grid', 'type': 'discrete', 'domain': np.arange(30, 31, dtype=int)},
{'name': 'shape_of_map_grid', 'type': 'discrete', 'domain': np.arange(50, 51, dtype=int)},
{'name': 'shape_of_rbf_centers_grid', 'type': 'discrete', 'domain': np.arange(2, 22, 2, dtype=int)},
{'name': 'variance_of_rbfs_grid', 'type': 'discrete', 'domain': np.arange(-5, 4, 2, dtype=float)},
{'name': 'lambda_in_em_algorithm_grid', 'type': 'discrete', 'domain': np.arange(-4, 0, dtype=float)},
{'name': 'number_of_iterations', 'type': 'discrete', 'domain': np.arange(300, 301, dtype=float)},
]
k_in_k3nerror = 10
def gtmf(x):
shape_of_map_grid = int(x[:, 0][0])
shape_of_rbf_centers_grid = int(x[:, 1][0])
variance_of_rbfs_grid = int(x[:, 2][0])
lambda_in_em_algorithm_grid = int(x[:, 3][0])
number_of_iterations = int(x[:, 4][0])
display_flag = 0
model = GTM([shape_of_map_grid, shape_of_map_grid],
[shape_of_rbf_centers_grid, shape_of_rbf_centers_grid],
2 ** variance_of_rbfs_grid, 2 ** lambda_in_em_algorithm_grid, number_of_iterations,
display_flag)
model.fit(X)
if model.success_flag:
# calculate of responsibilities
responsibilities = model.responsibility(X)
# calculate the mean of responsibilities
means = responsibilities.dot(model.map_grids)
# calculate k3n-error
k3nerror_of_gtm = k3nerror(X, means, k_in_k3nerror)
# 評価指標を出力
k3n_Z_error = k3nerror(X, means, k_in_k3nerror)
k3n_X_error = k3nerror(means, X, k_in_k3nerror)
print("k3n-Z-error={0}".format(k3n_Z_error))
print("k3n-X-error={0}".format(k3n_X_error))
print("k3n-error ={0}".format((k3n_Z_error + k3n_X_error)))
else:
k3nerror_of_gtm = 10 ** 100
return k3nerror_of_gtm
myBopt = GPyOpt.methods.BayesianOptimization(f=gtmf, domain=bounds)
#myBopt.run_optimization(max_iter=300)
myBopt.run_optimization(max_iter=30)
print(myBopt.x_opt)
print(-myBopt.fx_opt)
shape_of_map = [int(myBopt.x_opt[0]), int(myBopt.x_opt[0])]
shape_of_rbf_centers = [int(myBopt.x_opt[1]), int(myBopt.x_opt[1])]
variance_of_rbfs = 2 ** myBopt.x_opt[2]
lambda_in_em_algorithm = 2 ** myBopt.x_opt[3]
number_of_iterations = int(myBopt.x_opt[4])
display_flag = 0
# construct GTM model
model = GTM(shape_of_map, shape_of_rbf_centers, variance_of_rbfs, lambda_in_em_algorithm, number_of_iterations, display_flag)
model.fit(X)
# calculate of responsibilities
responsibilities = model.responsibility(X)
Z = responsibilities.dot(model.map_grids)
output_csv(args.output, df.index.values, Z[:, 0], Z[:, 1])
k3n_Z_error = k3nerror(X, Z, k_in_k3nerror)
k3n_X_error = k3nerror(Z, X, k_in_k3nerror)
print("k3n-Z-error={0}".format(k3n_Z_error))
print("k3n-X-error={0}".format(k3n_X_error))
print("k3n-error ={0}".format((k3n_Z_error + k3n_X_error)))
if __name__ == "__main__":
main()
結果
結果を以下に示す。k3n-Z-errorが①、k3n-X-errorは②の値である。数値が小さい程よい値である。
計算時間は掲載していないが、GTMはパラメータチューニングで30回実行しており、また1回あたりの計算時間も他手法より時間がかかった。
PCA
k3n-Z-error=2.6223228035570054
k3n-X-error=0.578405109666366
k3n-error =3.2007279132233712
MDS
k3n-Z-error=1.7886900303710693
k3n-X-error=0.44521214922051466
k3n-error =2.233902179591584
tSNE
k3n-Z-error=1.033589128650684
k3n-X-error=0.12006930439514711
k3n-error =1.153658433045831
UMAP
k3n-Z-error=1.970052284296805
k3n-X-error=0.17142884248639179
k3n-error =2.1414811267831966
GTM
k3n-Z-error=2.2293570682113013
k3n-X-error=0.40085975080186265
k3n-error =2.630216819013164
考察
- tSNEが①、②、総合でダントツで良い結果になった。
- GTMはパラメータチューニングの余地があるのか、思わしくない結果に。
ちなみにGTMの可視化結果はこんな図になりました。
今後
同じ画面内で複数の手法による可視化結果を比較できるようなインターフェースを作ってみたい。
参考
- k3nerror
- k3nerror k3nerrorのソース
- GTM 今回使ったGTM
- GTM (Generative Topographic Mapping) のハイパーパラメータチューニングでベイズ最適化を使った GTMパラメータチューニングで参考にした