#はじめに
多次元空間のデータを二次元平面で可視化した結果を評価する指標である、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パラメータチューニングで参考にした