Multidimensional Scaling (多次元尺度構成法, MDS) の計算方法と実装

Multidimensional Scaling (MDS) について少し勉強したのでメモしておきます。

Multidimensional Scaling とは

はじめに

日本語では「多次元尺度構成法」と呼ばれる統計テクニックの一つです。英語版ウィキペディアの記事が詳細です。

多次元尺度構成法 - Wikipedia
Multidimensional scaling - Wikipedia

MDS には色々な派生がありますが、ここでは Classical MDS (古典的多次元尺度構成法) と呼ばれる最も基礎的で古典的なものについてのみ記述します。

どういうときに使えるか

いくつか要素があり、要素間の距離がなんらかの形で定義できる場合に、その要素を多次元空間に配置する場合に使えます。

具体例

例えば、いま目の前に日本酒が10銘柄(鳳凰美田、花陽浴、...)あるとし、これらを一つずつ味見し「どの銘柄とどの銘柄がどのぐらい似ているか(似ていないか)」という指標によって、以下のように任意の2銘柄間の距離が定義できるとします。


要素間の距離の定義のイメージ*1d(\cdot, \cdot)は要素間の距離を表すとする。
d(農口, 村祐) = 2 d(農口, 鳳凰美田) = 4 d(農口, 花陽浴) = 7 ...
d(十四代, 村祐) = 2 d(十四代, 呼友) = 2 d(十四代, 花陽浴) = 5 ...
... ... ... ...


MDS を用いると、これらの距離の情報から、例えば2次元の座標系に銘柄の「マップ」を可視化することができたりします。ここで、各点の座標は MDS によって計算され、その際に各点間の距離は先ほど定義した距離が反映されます。



MDS の使用イメージ*2。MDS によって各点の座標が計算される。

入力と出力

MDS の入力は要素間の距離と、対象とする空間の次元です。一般的には、距離は距離行列 (distance matrix) として表現しておきます。要素数n であった場合、D_{i,j}i 番目と j 番目の要素の距離の二乗*3であるような正方行列 \mathbf{D} \in \mathbb{R}^{n \times n} が入力です。このとき、距離の対称性と非退化性*4から、 \mathbf{D} は対角成分が 0 であるような対称行列となります。

MDS の出力は、各要素の座標 (\mathbf{x}_1, \ldots, \mathbf{x}_n) です。ここで、指定された空間の次元が m だとすると、\mathbf{x}_i \in \mathbb{R}^m です。これらは i 番目の行を \mathbf{x}_i とする行列 \mathbf{X} \in \mathbb{R}^{m \times n} として表現したりします*5

次元圧縮 (Low-Dimensional Embedding)

高次元の点群データから距離行列 \mathbf{D} を計算し、元の次元よりも小さな次元 m を指定し MDS を実行することで、次元圧縮 (Low-Dimensional Embedding) を実現することができたりもします*6

非線形な次元圧縮手法としてよく使われる ISOMAP も MDS を拡張したものです。

MDS の計算方法

行列 \mathbf{X} を求めるためには、まず

\mathbf{K} = \mathbf{X}^T\mathbf{X} \in \mathbb{R}^{n \times n}

を満たす行列 \mathbf{K}*7*8を、入力として与えられた距離行列 \mathbf{D} から計算できるとします。すると、この行列 \mathbf{K}固有値分解 (eigenvalue decomposition) によって
\mathbf{K} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^T

と分解することができます*9*10。ここで \mathbf{V}固有ベクトルを並べた行列、\mathbf{\Lambda}固有値を対角成分に持つ行列です。これらを用いて、
\mathbf{X} = \mathrm{sqrt}(\mathbf{\Lambda})\mathbf{V}^T

とすることで、所望の \mathbf{X} を計算することができます*11。ここで \mathrm{sqrt}(\cdot) は要素毎に累乗根をとることを意味しています*12。また、このままでは空間の次元は n ですが、これをより小さい次元 m に削減するには、固有値の上位 m 個と、それに対応する固有ベクトルのみを使って前式を計算することで実現できます*13

与えられた距離行列 \mathbf{D} からカーネル行列 \mathbf{K} を計算する方法ですが、これは

\mathbf{K} = - \frac{1}{2} \mathbf{H} \mathbf{D} \mathbf{H}

によって算出するようです。ただし
\mathbf{H} = \mathbf{I} - \frac{1}{n}\mathbf{1}\mathbf{1}^T

は中心化行列*14と呼ばれ、得られる結果の原点を重心に移す効果があるそうです*15*16

MDS の実装

C++ によって MDS を実装した例を以下に示します。行列演算ライブラリとして Eigen を用いています。

#include <Eigen/Core>
#include <vector>
#include <utility>
using namespace std;
using namespace Eigen;

// This function computes low-dimensional embedding by using multi-dimensional scaling (MDS)
// - Input:  A distance (dissimilarity) matrix and a target dimension for embedding
// - Output: A coordinate matrix whose i-th column corresponds to the embedded coordinates of the i-th entry
MatrixXd computeMDS(const MatrixXd &D, unsigned dim)
{
    assert(D.rows() == D.cols());
    assert(D.rows() >= dim);
    const unsigned n = D.rows();
    const MatrixXd H = MatrixXd::Identity(n, n) - (1.0 / static_cast<double>(n)) * VectorXd::Ones(n) * VectorXd::Ones(n).transpose();
    const MatrixXd K = - 0.5 * H * D * H;
    const EigenSolver<MatrixXd> solver(K);
    VectorXd S = solver.eigenvalues().real();
    MatrixXd V = solver.eigenvectors().real();
    extractNLargestEigen(dim, S, V);
    const MatrixXd X = DiagonalMatrix<double, Dynamic>(sqrt(S)) * V.transpose();
    return X;
}

以下にテストコード等を含めた git repository 公開しました:
github.com

*1:これらの銘柄は実在するものですが、その距離として与えられている値には全く根拠はなく、極めて乱数に近いものであると考えてください。

*2:この図はあくまでイメージを掴むためのもので、実際のデータから計算されたものではありません。

*3:この部分について明示していない文献も多く私は途中まで勘違いしていたのですが、距離行列は距離の二乗を要素として持つようです。

*4:距離空間 - Wikipedia

*5:文献によっては i 番目の\mathbf{x}_{i}^T とする行列 \mathbf{X} \in \mathbb{R}^{n \times m} として表現したりすることもあります。

*6:この操作は Principal Component Analysis (PCA) と同値?

*7:グラム行列PRML 6.2章)

*8:文献によって \mathbf{B} と表記されることもあります。

*9:\mathbf{K} は対称行列であるため、特異値分解 (SVD) を用いても同様の計算ができます。

*10:\mathbf{D}ユークリッド距離に基づいているとき、まとそのようなときに限り \mathbf{K} は半正定値行列となる?

*11:\mathbf{X}^T\mathbf{X} = (\mathrm{sqrt}(\mathbf{\Lambda})\mathbf{V}^T)^T (\mathrm{sqrt}(\mathbf{\Lambda})\mathbf{V}^T) = \mathbf{V} \mathrm{sqrt}(\mathbf{\Lambda}) \mathrm{sqrt}(\mathbf{\Lambda}) \mathbf{V}^T = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^T

*12:ただし、\mathbf{K} が半正定値行列でない場合には固有値に負の値が現れます。その場合はその固有値を 0 だとみなして以降の計算を行います。これはある種の近似を行っていると考えられるそうです

*13:そのようにして良い根拠はこの資料で説明されています。

*14:文献によって \mathbf{G} と表記されることもあります。

*15:このあたりはよく理解できていませんが、距離を一定に保った点群には剛体移動(平行移動と回転移動)の自由度があるので、そのうち平行移動については重心が原点にくるように制約をかけるということだと思います。

*16:この変換は「ダブル・センタリング変換」と呼ばれる?

日本人研究者によるSIGGRAPH論文

概要

日本人研究者によるSIGGRAPH論文 (2014--2015年分) の情報を以下のページにまとめました。せっかくなので公開します:

SIGGRAPH Papers by Japanese Researchers

YouTube動画の埋め込みや著者の個人ページへのリンクも含んでおり、日本人研究者の活躍ぶりをまとめて把握する上で役に立つかもしれません。

SIGGRAPH論文

当該ページの冒頭にも書きましたが、SIGGRAPHはコンピュータグラフィクスの分野のトップカンファレンスとして知られ、またSIGGRAPHに論文が採択されることは非常に難しいことで有名です。SIGGRAPHに論文が採択されると、分野の研究者だけでなくメディアなどを通じて様々な人の目に触れることになります。

ちなみに、SIGGRAPH AsiaACM Transactions on Graphicsも"uniform high quality"であると明示的に宣言されており、まとめて「SIGGRAPH論文」と表現することも多いです。SIGGRAPHSIGGRAPH Asiaで発表された論文は自動的にACM Transactions on Graphicsに掲載され、またACM Transactions on Graphicsに掲載された論文はSIGGRAPHSIGGRAPH Asiaで発表する権利が与えられます。

また、ここでの論文とは正式にはTechnical Papersと表記されるものであり、Emerging TechnologiesやPostersとは別のものです。

なお、Eigenfactor scoreによると、ACM Transactions on Graphicsはコンピュータグラフィクス分野だけでなくComputer Science, Software Engineering分野においてもトップジャーナルなのだそうです。

日本人研究者

例年、SIGGRAPHに採択される日本人研究者による論文は片手で数えられる程度、0件であることも珍しいことではないようです。例えば、直近だとSIGGRAPH Asia 2014では日本人研究者による論文は採択されなかったようです。ただし今年のSIGGRAPHには非常に多くの論文が日本から採択されたようです。

なお、日本人がSIGGRAPHに論文が採択されると、コンピュータグラフィクスの国内シンポジウムの一つVisual Computing/グラフィクスとCAD合同シンポジウム (VC/GCAD) にて招待講演を依頼されるという制度があります。しかしながら、ACM Transactions on GraphicsやSIGGRAPH Asiaの論文は取りこぼしがあるようで、完全ではありません。

そこで、日本人研究者(海外に在住する人や、外国籍でも日本に深い縁のある人を含む)によるSIGGRAPH論文を、YouTube動画や著者の個人ページへのリンクも含めて集めてみました。収集手法は目grep法を用いたため、精度に問題があるかもしれません。情報に欠けを見つけた場合はご連絡ください。当該ページは、今後順次更新していくかもしれないし、現状のまま放置するかもしれません。

論文読み: Solid Simulation with Oriented Particles (SIGGRAPH 2011)

今回は古い論文の紹介です。

2年前にこの手法を実装したことがあるのですが、そのときにキャプチャした動画もいくつか紹介します。

Matthias Müller and Nuttapong Chentanez
Solid Simulation with Oriented Particles
SIGGRAPH 2011

研究概要

Position-Based Dynamics の上で動く、ロバストで高速な変形可能物体の物理シミュレーション手法「Oriented Particles 法」を提案しています。

これは Shape Matching 法と呼ばれる手法の拡張になっていて、従来は三次元的な構造(ソリッド)しか扱うことができなかったのに対して、布などの二次元的な構造(シェル)やひもなどの一次元的な構造(ロッド)など、様々な構造を統一的に扱うことができるようになった点が特徴的です。

f:id:yuki-koyama:20150130131454p:plain

Position-Based Dynamics については以下を参照してください。


実装してみた

2年ほど前、ちょうど私が情報処理推進機構の未踏IT人材発掘・育成事業でのプロジェクトで物理エンジンを開発していたころ、実際にこの Oriented Particles 法を実装したことがあります。

以下がそのときにキャプチャした動画です。いずれも当時の普通の MacBook Air でリアルタイムに計算しています。C++OpenGL を用いて実装しています。

三次元的な構造 (ソリッド)

Shape Matching 法と同様、ソリッドな構造をシミュレーションできます。パーティクルのサイズを小さく設定した場合には Shape Matching 法とほぼ同等の挙動となります。

二次元的な構造 (シェル)

布のシミュレーションをしています。非伸縮制約 (inextensibility) をいれていないため、引っ張るとよく伸びます。非伸縮制約を入れた場合は、これが伸びの挙動を主に決定して、Oriented Particles 法による制約が曲げの挙動を主に決定することになります。

一次元的な構造 (ロッド)

ひものシミュレーションをしています。伸びや曲げについては良い具合にシミュレートできていますが、ねじったときの複雑な挙動は再現することができません。

Unity に移植してみた

当時急速に人気になりつつあった Unity に移植してみました。外部ライブラリとしてではなく、C# で直接書いています。

Unity 上で物理エンジンを実装しようとして試行錯誤している最中に書いた記事はこちらです。

ちなみにこの Unity 移植版を元に別の手法の研究開発を行い、その後学術論文として発表したりしました。

Yuki Koyama and Takeo Igarashi
View-Dependent Control of Elastic Rod Simulation for 3D Character Animation
SCA '13
http://www-ui.is.s.u-tokyo.ac.jp/~koyama/project/ViewDependentRodSimulation/index.html

またこの成果を「こだわり物理エンジン」という名称でパッケージにまとめたものが以下になります。

Web ページ:こだわり物理エンジン

その他:最近の研究動向と私見

Oriented Particles 法は 2011 年の夏に発表された研究です。それから3年半も経っていますから、当然リアルタイム物理シミュレーションの動向もだいぶ変わってきています。

ソリッドな弾性体をシミュレーションしたい場合

基本的な選択肢としては Shape Matching 法か Oriented Particles 法の2択になりそうです。
Oriented Particles 法の方が

  • データ構造の拡張が必要
  • 計算がほんの少しだけ重い

という小さなデメリットがある反面、

  • 極端な変形に対してもロバストに動く
  • 各パーティクルが回転情報を保持しているため、スキニングの際に利用できる

といったメリットがあります。特にゲーム等での利用を考えると2個目のメリットは大きそうです。

また Position-Based Dynamics の拡張として以下のような研究もありますが、Oriented Particles 法に比べて巨大な行列を保持しておく必要があったり、実行時速度が不十分だったりします。

布をシミュレーションしたい場合

Oriented Particles 法は使うべきではないと思います。

Oriented Particles 法は3行3列の行列の特異値分解 (SVD) または極分解 (polar decomposition) をパーティクルの数だけ毎フレーム行う必要があります。それに対し、特異値分解を必要としないより高速な布の表現手法が多数存在します。

Tae-Yong Kim, Nuttapong Chentanez, Matthias Müller
Long Range Attachments - A Method to Simulate Inextensible Clothing in Computer Games
SCA '12

ひもをシミュレーションしたい場合

これも、Oriented Particles 法は使うべきではないと思います。

再び、Oriented Particles 法は特異値分解を必要とするため、効率的ではありません。それに対し、特異値分解を必要とせず、さらにねじりに対する非線形なひもの振る舞いをモデル化できている手法が最近提案されました。

Nobuyuki Umetani, Ryan Schmidt, Jos Stam
Position-based Elastic Rods
SCA '14

Oriented Particles 法が役に立つシナリオ

複合的な構造

ソリッド・シェル・ロッドが組み合わさっているような物体をシミュレーションするときには Oriented Particles 法は最もシンプルで良い解決法だと思います。例えばこの記事の冒頭で掲載した動画ではタコをシミュレーションしていますが、胴体はソリッド、足の部分はロッドのような構造になっています。Oriented Particles 法ではソリッドもロッドも特に区別せず統一的に扱うことができます。

メッシュを用意するのが難しい場合

Oriented Particles 法ではパーティクル同士が形成しているネットワーク (パーティクルを頂点とするグラフ構造) を対象にシミュレーションを走らせることができる点も利点の一つです。通常のシミュレーションではメッシュ (三角形メッシュや四面体メッシュ) を用意する必要があるのに対して、グラフ構造は生成が簡単です (例えば各パーティクルをある一定の距離以下に存在するパーティクルと接続すると簡単にグラフ構造を作れる)。メッシュを用意するのが大変な状況においても Oriented Particles 法は役に立つといえます。なおこのような考え方を ``mesh-less'' だとか ``mesh-free'' と表現したりするようです。

情報検索のための Learning to Rank の概要

この記事では以下の文献の第一章の最初の方をざっと読んで得た知識を簡単にメモしています。事前に断っておくと、私は情報検索 (Information Retrieval, IR) の専門家でもありませんし、機械学習の専門家でもありません。誤り等ありましたらコメント欄にてご教示ください。

Tie-Yan Liu
Learning to Rank for Information Retrieval
Foundations and Trends in Information Retrieval: Vol. 3: No. 3, pp 225-331. (2009)
http://www.nowpublishers.com/article/Details/INR-016

導入

情報検索とは情報のデータ群の中からクエリ (query) に合う要素を取り出すことだそうです。
情報検索 - Wikipedia

情報検索には

  • 協調フィルタリング (collaborative filtering)
  • 質問応答システム (question answering)
  • マルチメディア検索 (multimedia retrieval)
  • 文書要約 (text summarization)

など様々な応用があり得ますが、この文献では特にウェブページの検索エンジンを想定して議論を進めています。

この情報検索について、近年急速に流行っているのが Learning to Rank というアプローチだそうです。特に 2005 年以降

  • SIGIR
  • ICML
  • NIPS
  • WWW

などの国際会議において流行っていることが丁寧に説明されています。

ちなみにこの Learning to Rank はどうやら日本語ではランキング学習などと訳されるようです。

従来の情報検索

Learning to Rank の定義について説明する前に、それ以前の情報検索手法について簡単に説明を加えています。

Relevance Ranking Models (関連度によるモデル)

このモデルではクエリとデータの関連度を計算し、関連度の高いものを検索結果として提示するというアプローチをとります。

関連度の計算としては、ベクトル空間モデル (Vector Space Model, VSM) がまず取り上げられています。VSM ではクエリとデータの両方がベクトルとして表されており、その関連度は例えば両者の内積などによって計算するというものです。

VSM の具体的な実装として、ここでは TF-IDF モデルが紹介されていました。
tf-idf - Wikipedia

別のアプローチとして Latent Semantic Indexing (LSI) や Probabilistic Ranking Models についても言及されていましたが、ここでは省略します。

Importance Ranking Models (重要度によるモデル)

このモデルではデータそのものの重要度を計算し、重要度の高いものを優先的に検索結果として採用するというアプローチをとります。

その具体例として、PageRank が紹介されていました。
ページランク - Wikipedia

ちなみに私の理解では以下の論文が PageRank を提案しているということになっているのですが、この文献の中では何も引用されていませんでした。PageRank はもはや引用なしで議論できるほど有名になっているということでしょうか。

Brin, S. and Page, L.
The Anatomy of a Large-Scale Hypertextual Web Search Engine
WWW '98
The Anatomy of a Large-Scale Hypertextual Web Search Engine - Stanford InfoLab Publication Server

Learning to Rank による情報検索

さて、従来の情報検索手法における問題点を挙げるとするなら、それはパラメタのチューニングが必要であるという点でした。例えば PageRank では damping factor と呼ばれるパラメタを適切にチューニングしてやる必要があります。

そこで登場するのが、機械学習を活用した情報検索という流れです。機械学習によって、パラメタの自動チューニングや、複数の手法を適切な重みで組み合わせるといったことが可能となります。

Learning to Rank の定義

一般的には、情報検索の問題を機械学習の手法によって解決することを目指すような一連の手法を Learning to Rank と呼んでしまうことが多いことが指摘されています。例えば機械学習によるパラメタの自動チューニングなども Learning to Rank に含まれることがあります。

しかしながら、この文献ではより狭い定義で Learning to Rank という言葉を使おうとしています。より具体的には、以下の 2 つの特徴を持つ手法を Learning to Rank であると定義しています:

  1. Feature Based: データが全て特徴量を用いて表現されていること。
  2. Discriminative Training: 学習のプロセスが一般的な (教師付き) 機械学習フレームワークによく沿っていること。

この辺りの読み込みが甘く、厳密な理解はできていませんが、簡単に言って、Learning to Rank は「(教師付き) 機械学習によって特徴量の最適な組み合わせ方を発見する一連の手法群」と考えることができそうです。

Learning to Rank のフレームワーク

ここから先はよく読み込めていないので省略しますが、多くの手法は以下の 3 種類に大別できるそうです:

  • The Pointwise Approach
  • The Pairwise Approach
  • The Listwise Approach