Word Mover's Distance(2)fastWMD のビルド(前編)

Word Mover's Distance (2) Building fastWMD (Part 1)


2021/05/03
藤田昭人


前回、 概要を紹介したので、 続けて fastWMD のビルドについてサクッと手短に説明したいと思います。 実はいろいろ試行錯誤があったのですが、 そのあたりを端折って手順をつらつらと書き出してみたところ、 思いのほか長くなったので前後編に分けました。

実際、fastWMD の実装は論文の評価パートのためのコードのようで、 データ採りに必要最低限の内容しか含まれていない印象です。 README等で示されている手順に補足説明や追加作業を加えておきます。


まずはソースコードのダウンロード

fastWMD のソースコードGithub からダウンロードできます。

github.com

ちなみに Python のラッパーも収録されたバージョンは こちら をご覧ください。

インストール環境について

次に、インストール環境について。fastWMD の README.md によれば、次の環境で動作を確認したとのこと。

  • Ubuntu 16.04 and 18.04
  • Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz,with 8 GB of RAM

ちなみに以降の作業は次の環境で実行しました。

しかし Windows 等、他のプラットフォームでも 概ね同じ手順でインストールできると思います。 たぶん(笑)

インストールの手順を確認すると install_dependencies.shbuild-project.sh を順次実行するだけなのですが、 試してみたところ install_dependencies.sh でズッコケました。 本当に Ubuntu だけ インストールできる仕様になっているようです。


install_dependencies.sh の手順を手作業で実行する

確認したところ install_dependencies.sh では、fastWMD のビルドに必要な Google OR-ToolsEigen3 をインストールしているだけでした。 どうやらディレクトリ構成やセキュリティ・チェックなど 環境の違いがズッコケる原因のようなので、 パッケージを手作業でインストールします。

Google OR-tools

Google OR-tools はその名のとおり オペレーションズ・リサーチ関連、 以前紹介した lpSolve と同じカテゴリのパッケージのようです。 WMDに必要な最適輸送問題を解くアルゴリズムを 呼び出しているのでしょう。

このパッケージ、C++だけでなく、PythonJavaC# もサポートしているようですが、 必要なのは C++ です。実は当初はバイナリ配布のインストールを試みたのですが macOS Big Sur のセキュリティ・チェックで弾かれました。 (Google が配布しているので、今はこの問題を回避できるのかもしれませんが…) やむなくソースコードからインストールをしました。手順は以下のページです。

developers.google.com

この手順説明は macOS での C++ プログラムのビルド環境の 一般的な設定についても丁寧に説明しているので、 もし C++ プログラムを動かしたことがなければ ソースコードからのインストールの実行をお勧めします。 (どのみち fastWMD のビルドの際に設定をしなければなりません)

Eigen3

Eigenは、線形代数、行列およびベクトル演算、幾何学的変換、 数値ソルバーとそれに関連するアルゴリズムのための テンプレート・ヘッダーによる高水準 C++ ライブラリです。 おそらく fastWMD ではベクトル表現や基本演算のために使っているのでしょう。 ソースコードは次のページからダウンロードできます。

eigen.tuxfamily.org

で、ソースツリーを眺めていたら CMake 関連のファイルを見つけたので 「テンプレート・ライブラリなのに?」 と思っていたら Eigen の INSTALL メモには次のような記述がありました。

Eigen consists only of header files, hence there is nothing to compile before you can use it. Moreover, these header files do not depend on your platform, they are the same for everybody.

Eigen はヘッダファイルのみで構成されていますので、使用前にコンパイルする必要はありません。さらに、これらのヘッダファイルはプラットフォームに依存しません。プラットフォームに依存せず、誰でも同じものが使えます。

…ということなので、配布ファイルを展開して しかるべきヘッダーファイルのパスに入れることにしました。 実はちょっと試行錯誤があったのですが、 次の手順でOKなようです。

$ tar xzf eigen-3.3.9.tar.gz
$ sudo mv eigen-3.3.9 /usr/local/include/eigen3
$ 


まとめ

以上、fastWMD のビルドにまでは至ってませんが、 記事的な分量のキリが良いので、本稿はここまで。 間をあけずに後編をアップしますので、 よろしくお願いします。

以上



f:id:Akito_Fujita:20210503123933p:plain

Word Mover's Distance(1)fastWMD ー 高速化の試み

Word Mover's Distance(1) fastWMD


2021/04/29
藤田昭人


本稿では、論文 "Speeding up Word Mover’s Distance and its variants via properties of distances between embeddings" を眺めながら Word Mover's Distance (WMD) の高速化の試みをザックリと紹介したいと思います。


Word2Vec に思うこと

まずは Word2Vec について… 2013年に論文発表されたWord2Vec。 正直に告白すると当時は 「これ、いったい何に使えるんや?」 と思ったものです*1。 冒頭の論文では、Word2Vec の価値(そして意義)について次のように語られています。

数年前までは、意味的な関係(semantic relations)を取得する適切な方法がなかったので、 意味的な関係はほとんど利用されていませんでした。 それ故、研究者たちはこの問題を軽減する方法として、オントロジーを使用する決心をしましたが、 これはアプリケーションが外部の知識ベースに依存することを意味していました。 このシナリオが変わったのは,Word2Vec とその亜種が登場してからです。

この Word2Vec で分析できるとされる 「意味的な関係」が、 人工知能研究が始まって以来の命題である 「知能」に対する理解を 大きく変えつつあると僕は実感しています。

文章の類似性を測る WMD にとって その文章を構成する単語の類似性を測る Word2Vec は 必須のアルゴリズムですが、 同時に WMD は Word2Vec の有用性を 具体的に提示する重要なアプリケーション (のひとつ)でもあるように思います。


WMDの高速化の試み

実際「よく似た内容の文章を見つける」という 類似文書検索は日常的に役立つ アプリケーションでしょう。 その機能的な核となる文章間の類似度を 高い精度で計測できるWMDアルゴリズムは 需要が高い。それ故に登場した当初から 指摘されていた「分析に時間がかかる」という 課題には多くの研究者がチャレンジして来ました。 昨年(2020年)発表された冒頭の論文は、 WMD の Earth Mover's Distance (EMD) の部分、 すなわち以前紹介した最適輸送問題の解法に着目した高速化手法を扱っています。 参考文献には次の4本の先行研究も掲載されていて、 高速化はWMDが登場した2015年以来 毎年のように成果が発表される、 今もホットな研究領域のように見えます。

年号 論文
2015 From word embeddings to document distances
2017 Linear-complexity relaxed word mover’s distance with gpu acceleration
2018 Word mover’s embedding: From word2vec to document embedding
2019 Linear-complexity data-parallel earth movers distance approximations
2020 'Speeding up Word Mover’s Distance and its variants via properties of distances between embeddings'

そもそも高速化問題とは、 最適輸送問題の解法 が輸送元と輸送先の組み合わせを総あたりで計算することに起因します。 WMDに適用した場合には、その計算量は扱う語数の3乗に比例するそうで、 対象とする文書が大きい過ぎると「いつまで経っても終わらない」計算になります。

この問題は2015年のオリジナルのWMD論文でも指摘されていて、 WMDと同じ精度でもっと高速に計算する手法が提案されました。 冒頭の論文は2015年以来の先行研究を網羅した上で、 独自の手法を提案しています。


fastWMD ー C++ による WMD 実装

僕が、冒頭の論文に注目するもうひとつの理由は、 提案の評価に使われた(先行研究も含む) C++ による WMD 実装が Github で公開されていることです。

github.com

もちろん「実装コードがある」というのはありがたいことですが、 特に fastWMD では先行研究の成果であるアルゴリズムも実装されているので コードレベルで直接比較できることは大変ありがたいです。 ザッと調べてみたところ fastWMD に実装されているアルゴリズムは オリジナルの Word Mover ’s Distance (WMD) に加えて次の3つです。

  • Relaxed Word Mover ’s Distance (RWMD)
  • Linear-Complexity Relaxed Word Mover’s Distance (LC-RWMD)
  • Related Relaxed Word Mover ’s Distance (Rel-RWMD)

RWMD はオリジナルの論文で言及されていた高速化手法で、 最適輸送問題の制約条件をカットすることにより計算量を削減するアイデアです。 LC-RWMD は前処理を追加して計算量が語数に比例するように組み替えた高速化手法で、 並列化が容易になるよう配慮されている(?)反面、計算精度が低下するらしい。 Rel-RWMD は与えられた単語に対して 関連性がある単語グループと関連性がない単語グループに分ける 閾値を決めて計算の枝刈りをする高速化手法で、 Word2Vec では対象となるデータセット内での単語埋め込み間の意味的距離は 次のようなグラフを描く特性を活用しているそうです。

f:id:Akito_Fujita:20210429220216p:plain
Amazonのデータセットに含まれる全ての単語埋め込みから、単語 “cat” までの距離

Word2Vecを説明した記事 では distance コマンドを紹介しましたが、 その表示でグラフを描くとこのようなカーブを描くようですね。 Rel-RWMD の弱点は閾値は単語ごとに異なる可能性があることでしょうか?


まとめ

本稿では冒頭の論文を軸に fastWMD がサポートする WMDの高速化手法についてザックリと紹介しました。 この他にもWMDの高速化手法はたくさん提案されているようですが、 やはり実装がある fastWMD から手を付けるべきでしょう。

僕は対話システムへの応用を考えているので、 必ずしも大規模な文書が想定される訳ではないように思います。 いずれの手法が有効なのか?は対話システムの設計次第でしょう。 個々の高速化手法の詳細が気になるところですが、 以後、まずは fastWMD のビルド方法から説明することにします。

以上



*1:その僕がまさか 記事 を書くとはね。

Word Mover's Distance(WMD)

Word Mover's Distance


2021/04/27
藤田昭人


僕は一度気が削がれると なかなか元の状態に戻れないヤツなのですが…
とある キッカケ を掴んでブログ執筆を再開することができました。

ここからの数回は Word Mover's Distance (WMD) について取り上げます。

BookBot の紹介をして以来、 Word2VecEarth Mover's Distance と対話機能の説明からドンドン離れて行くので 「迷走している」 と感じている方もいらっしゃるかもしれませんが、 この WMD の説明で「なるほど」と納得してもらえる…つもりです(笑)


WMD って何?

WMDは、類似文書検索(あるいは 概念検索 ともいう)、つまり 「文書検索で任意の文を与えると、 その文と意味的に近い文を検索する」 ためのアルゴリズムなんだそうです。

言い換えれば 「2つの文の意味的な距離を測る」 アルゴリズムなんですが…

  1. 文を単語単位に分解し、
  2. Word2Vec を使って、各々の単語毎に「単語の意味的な距離」を測り、
  3. EMD を使って、単語間の意味的距離を合算する

…のがWMD、つまり直感的には WMD = Word2Vec + EMD と説明すると分かり易いでしょう。

実は、コーパスに書籍を使うBookBotでは、応答する場合に 「質問文と意味的に近い文を書籍の中から探す」 必要があるため、このアルゴリズムに着目しました。


WMDに関する「論文」的な話

WMD は2015年に発表された論文 "From Word Embeddings To Document Distances" が初出なようで、ごく最近登場したアルゴリズムのようです。

この論文を読むには類似文書検索などの 統計的自然言語処理の基礎知識」 が必要で、特に「テキストを行列に変換する」 ベクトル空間モデル単語埋め込み、 単語分散表現など (これ全部同じことを意味する用語のようなんですが) が基礎になっています。

WMDをザッと理解したい場合には 前述の論文の日本語による概要説明が役に立ちます。

yubessy.hatenablog.com

この概要説明でも次の定番のオバマ元大統領に関する例文が登場します。

Sym Sentance
D0 'The President greets the press in Chicago.'
D1 'Obama speaks to the media in Illinois.'
D2 'The band gave a concert in Japan.'


この3つの例文の間の意味的距離の計算方法を図示したのが次の図です。

f:id:Akito_Fujita:20210423195546p:plain

この図を見ると Word2Vecで求めた単語どうしの意味的距離輸送問題の解法で集約して 文の間の意味的距離を計算していることが直感的に理解できます。


WMDをちょっと試してみる

WMD を手軽に試してみるには Python を使うのが一番楽です。
次のブログ記事では Python を対話的に使って WMD を計算する手順を説明しています。

hironsan.hatenablog.com

Python3 を使ってこの手順どおりにWMDを計算してみました。

$ python3
Python 3.9.1 (v3.9.1:1e5d33e9b9, Dec  7 2020, 12:10:52) 
[Clang 6.0 (clang-600.0.57)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from gensim.models.keyedvectors import KeyedVectors
from gensim.models.keyedvectors import KeyedVectors
>>> model = KeyedVectors.load_word2vec_format('./GoogleNews-vectors-negative300.bin', binary=True)
model = KeyedVectors.load_word2vec_format('./GoogleNews-vectors-negative300.bin', binary=True)
>>> model.most_similar(positive=['japanese'])
model.most_similar(positive=['japanese'])
[('japan', 0.6607722043991089), ('chinese', 0.6502295732498169), ('Japanese', 0.6149078607559204), ('korean', 0.6051568984985352), ('german', 0.5999272465705872), ('american', 0.5906798243522644), ('asian', 0.5839767456054688), ('san', 0.5834755897521973), ('jap', 0.5764404535293579), ('swedish', 0.5720360279083252)]
>>> sent1 = 'But other sources close to the sale said Vivendi was keeping the door open to further bids and hoped to see bidders interested in individual assets team up.'.split()
sent1 = 'But other sources close to the sale said Vivendi was keeping the door open to further bids and hoped to see bidders interested in individual assets team up.'.split()
>>> sent2 = 'But other sources close to the sale said Vivendi was keeping the door open for further bids in the next day or two.'.split()
sent2 = 'But other sources close to the sale said Vivendi was keeping the door open for further bids in the next day or two.'.split()
>>> distance = model.wmdistance(sent1, sent2)
distance = model.wmdistance(sent1, sent2)
>>> print(distance)
print(distance)
0.8738126733213613
>>> d0 = 'The President greets the press in Chicago.'.split()
d0 = 'The President greets the press in Chicago.'.split()
>>> d1 = 'Obama speaks to the media in Illinois.'.split()
d1 = 'Obama speaks to the media in Illinois.'.split()
>>> d2 = 'The band gave a concert in Japan.'.split()
d2 = 'The band gave a concert in Japan.'.split()
>>> dist01 = model.wmdistance(d0, d1)
dist01 = model.wmdistance(d0, d1)
>>> dist02 = model.wmdistance(d0, d2)
dist02 = model.wmdistance(d0, d2)
>>> print(dist01)
print(dist01)
1.968269276774589
>>> print(dist02)
print(dist02)
2.2929445610887638
>>> quit()
quit()
$ 

ここではブログ記事と同じ手順で、 前述の3つの例文の意味的距離も計算しています。 D0-D1 間の距離が 1.968269276774589 に対し、 D0-D2 間の距離は 2.2929445610887638 となりました。 Word2Vecの学習済みデータが異なるので、 前述の図とは値が異なりますが、 D0-D1 間の距離の方が D0-D2 間の距離よりも近いことがわかります。


まとめ

本稿ではWMDについて大変ザックリと説明しました*1

やはり Python は Gensim などパッケージが充実しているので、 今どきの統計的自然言語処理のちょっとした実験には大変便利な言語です。 日本語による解説記事もたくさんありますしね。

でも、さまざまなパッケージが高度にインテグレートされているので、 アルゴリズムの中身を調べようとすると骨が折れると言うのが僕の印象です。 とにかく覚えなければならないことが多いので、正直ウザい(笑) そもそも "Simple is Beautiful" の世代の僕には、 C言語のような「必要最小限」アプローチがフィットしているようです。

以降の解説では敢えて Python を使わないアプローチで WMD に迫ってみたいと考えてます。

以上

*1:前回 語ったように「週3回更新」が可能か?試してみる気になりました。
10分〜15分で読み切れるように 1回の記事の分量を調整したつもりですが、 いかがだったでしょうか?

10万というマジックナンバー

100,000 as the magic number.


2021/04/20
藤田昭人


約1ヶ月ぶりのブログです。
実は締め切りの3日前にどうにか確定申告しまして、 その後はボーッと過ごしてました。
ブログを再開するにあたり、 マジックナンバーの話を書きたいと思います*1


顔芸ベーシスト、ぴんはげ氏の話

このところの僕のマイブームはベーシスト YouTuber の ぴんはげ 氏でして、ほぼ毎晩、寝る前に彼が制作した映像を見てます。 この方、ベースの腕前はもちろんのこと、 超絶技巧をふんだんに盛り込んだ課題曲を連発する作曲能力、 ルーパーやDTMツール(?)を駆使したカバー曲の演奏で見せるアレンジ能力、 それにハイトーン・ボーカル*2… 音楽的才能の塊のような人なのですが、 音楽的高いスキルと独特の顔芸とのギャップ感が人気の秘密なんでしょうねぇ。 やはり一芸で YouTube でブレークするには プロを凌駕する圧倒的なスキルが必要なんでしょうねぇ。

…なもんで、連日ぴんはげ氏の映像アーカイブを漁っているのですが、
いつもとはテイストの違う次の映像を見つけました。

www.youtube.com

この映像の 4:30 あたりからを観て欲しいのですが、
「本アカウントの登録者数が10万を超えたあたりから状況が変わり始めた」
とぴんはげ氏は仰っています*3


10万というマジックナンバー

情報屋にとってマジックナンバーとは 「理由はよくわかんないけど、 意味があるだろうとなんとなく確信してしまう 魔法の定数」なんですけども、 10万という数字には 幾つか身に覚えがあります。

僕の一番古い記憶では、 かつてパソコン通信が華やかだった頃の Nifty Serve で 「登録者数が10万を超えたあたりから状況が劇的に変化した」 との運営者のコメントを聞いたことがあります。 もっともその時は 「10万とは興行や広告・宣伝の業界人が媒体として注目する人為的な数値」 と理解したのですがね。レコードが10万枚売れるとか、 コンサートツアーで10万人動員とか… 90年代は1日で10万人を集めるコンサートが よく話題になりましたよね。


ソーシャルデータの分析での「10万」

ところが…

その後、この10万という数字には 別の意味がありそうなことに気づく機会がありました。

古くから僕を知っている方は 僕がIIJ時代に「Wikipedia ランキング」なるサービスを 運営していたことを覚えてらっしゃる方もいるかもしれません。 Wikipedia は各ページ毎の1時間単位のページビュー数をデータとして公開しているのですが、 毎時、このデータをソートしてランキングで表示するサービスを運営していたのでした。

その延長でWikipedia日本語版の民放のドラマページのページビュー数と視聴率との相関性を 調べ始めました。運の良いことに今世紀最高視聴率のドラマ『半沢直樹』が放送された直後だったので、 そのページビュー数と視聴率を調べて得た公式を使って、 その他の民放ドラマに適用していったのです*4。 確か、ドラマページの1週間の累積ページビュー数が10万を超えると 視聴率とのシンクロ率が急激に上がる知見を得ました。 が、全てのドラマページの累積ページビュー数が10万を超える訳ではなく… その当時の社長に「で、結局ドラマページのページビュー数から視聴率は予測できるのか?」 と聞かれた時に「うまくいく場合とそうでない場合があります」と答えたら、 「バカヤロー!!」と怒られましたが(笑)

ともあれ…

この経験から 「10万」という数字は 興行や広告・宣伝の業界人の直感だけではない なんらかの意味が内包されているマジックナンバーだ と僕は考えるようになりました。


対話システム開発における「10万」

その後、対話システムでも このマジックナンバーに遭遇しました。

LINEに在籍している頃は 用途が限定されたスマートスピーカーのスキル (スマートスピーカー用アプリ) の開発を手がけていました。 この用途に有効な対話コーパスは存在しなかったので、 昔ながらの ELIZA スタイルの応答生成エンジンを使っていたのですが、 スキルが自然な会話ができるよう 機械学習を使って一般的な対話コーパスを取り込む方法を 模索したことがあります。

対話コーパスとは質問文とそれに対する回答文のペアが たくさん記録された会話録のようなデータなのですが、 同様の研究をしている文献をしらみつぶしに探して観たところ、 どの論文をみても10万組の質問文と回答文のペアを学習すると 対話システムがもっともらしい応答をすることがわかりました。

もっとも、これには反例があります。 以前書いたブログ記事 で紹介しましたが、 チューリング・テストのコンテストであるロブナー賞の2000年の覇者でもある Artificial Linguistic Internet Computer Entity (通称: AliceBot)の開発者である Richard Wallace は AliceBot の開発において、オフィスの同僚に対話をしてもらい、 AliceBot が答えられなかった質問について AliceBot になったつもりで考えた応答文を随時ルールに追加していったそうです。 その知見によれば、概ね4万組の質問文・回答文ペアがルールに登録できれば、 ボットは人間のどんな質問にも答えられるようになると見積もっていたそうです。

この時のざっくりした調査では、 無作為な文例(例えばSNSなどのログデータから質問&応答ペアをピックアップする) では10万件以上、 作為的な文例(質問文に対する回答を人間が考える) では4万件以上の文例が必要と結論づけたのですが、 その時点で手元にある文例は数百のオーダーだったので 機械学習の導入は「もう少し文例が溜まってから…」と棚上げにしました。


ふたたび、ぴんはげ氏の話

このように僕は「10万」にまつわる様々な経験をしてきたことから、 ぴんはげ氏の口から「10万」という数字を聞いた時、 僕は「やっぱりねぇ…」と思ったのでした。

もちろん、この「10万」というマジックナンバージップの法則 などと同様のソーシャルな経験則に基づくルールなので 論理的な根拠を示すのが難しい代物です。 でも、いわゆるビッグデータを扱うときの目安値としては重宝します。

さてさて…冒頭のぴんはげ氏の話に戻ります。

冒頭の引っ越しビデオをよく聞いてると「10万を超えるための戦略」は 「去年(2019年)の12月から『週3本投稿します』と宣言して…」 ってことでした。ぴんはげ氏のどちらかというと陽気なビデオの内容とは 裏腹に過酷な生活を送ってるんだろうなぁ…と想像しているところ。

僕もいわゆる「ブログを書いている人」なので、 「アクセス数が10万を突破して状況が劇的に変化する」 には大変に関心があるのですが、 今の内容でブログ記事をアップするには2週間に1本が精一杯。 週3でアップするとなると月刊連載時に培った執筆スタイルを捨て、 ブログ向けシフトをしなくちゃなぁ…と思い悩んでいるところです。

以上

*1:って言うか、このネタを思いついた時 「これは書いておかなければと…」 と思ったのでした。これがなければ ブログを再開するキッカケを掴めなかったかも。

*2:僕が彼にハマったのは次の映像を見たからです。

www.youtube.com

これ、例の劇場版『鬼滅の刃』の主題曲なんだそうで…
作詞・作曲は Kalafina プロデューサーの梶浦由記で、 元々彼女の楽曲は僕のツボに刺さることが多いのですが、 この映像でのぴんはげアレンジは物凄くて 「ベース1本でこんなことができるのかぁ…」 と衝撃を受けた次第。

実はこの曲をフルコーラスで聞いたのはこの映像が初めてで 「オリジナルよりこっちの方が絶対良い」 と僕は信じて疑いません(笑)

*3:しかし、 楽器メーカーからシグネチャモデルの提供を受けるって、 プロでもなかなかないことなんじゃないかと…

*4:この試みはその時の僕のデータサイエンスの師匠との連名で特許も取得したんですがね。

Earth Mover's Distance(2)lpSolveを使ったC実装

Earth Mover's Distance (2)
EMD with lpSolve


2021/03/11
藤田昭人


だいぶん間が空いてしまいましたが…

実は Earth Mover's Distance の実装を巡って悪戦苦闘してました。 その顛末をダラダラを書き連ねた記事が予想外に長くなってしまったので、 要点のみをピックアップして再構成したものが本稿です。

今日 Earth Mover's Distance ついてはさまざまな実装が出回っていますが、 調べてみたところ Earth Mover's Distance を最初に提案した論文 "A Metric for Distributions with Applications to Image Databases" の著者 Yossi Rubner が作成した(と思われる) Code for the Earth Movers Distance (EMD) に掲載されているC実装がオリジナルのようです。 ですが、何分にも20年以上前の1998年に作成されたコードなので、 現在のCコンパイル環境では正しく動いてくれません。

やむなく、代わりとして使える リファレンス実装のコードを探すところから 本稿の作業を始めました。


Rによる EMD の実装

さまざまな条件でググってみたところ…

前述のオリジナルのC実装にも言及している 次のブログ記事を見つけました。

aidiary.hatenablog.com

この記事では EMD が (前回紹介した) 線形計画法の輸送問題を応用した指標であることを説明したのちに、 前述のオリジナル実装に付属するデモコード example1.c を例題として Rと Python の実装例を紹介しています。

記事によれば EMD は次の公式で表されますが…


{\rm EMD} (P, Q) = \frac{\displaystyle \sum_{i=1}^m \displaystyle \sum_{j=1}^n d_{ij} f_{ij}^*}{\displaystyle \sum_{i=1}^m \displaystyle \sum_{j=1}^n f_{ij}^*}

この式の分母部分の  \displaystyle \sum_{i=1}^m \displaystyle \sum_{j=1}^n f_{ij} が輸送される荷物の重さの総和、 分子部分の  \displaystyle \sum_{i=1}^m \displaystyle \sum_{j=1}^n d_{ij} f_{ij} が輸送問題の解なんだそうです。 つまり、この式は輸送問題を解いた後、 移動させる荷物の総重量で割って 単位重量あたりの作業量を計算していることになります。

f:id:Akito_Fujita:20210309195233p:plain

上記は記事から拝借してきた オリジナル実装の example1 の条件を示した図です。 この条件で輸送問題を解くと、 荷物は次のように移動するとコストは最小となります。

移動元 移動先 荷物量
 P_{1}  Q_{1} 0.4
 P_{2}  Q_{1} 0.1
 P_{2}  Q_{2} 0.2
 P_{3}  Q_{3} 0.2
 P_{4}  Q_{2} 0.1

この場合、移動する荷物量の総和は 1.0 となるので、 輸送問題の解である 160.542763 が EMD の解になります。


Rパッケージの lpSolve

Rでの EMD 実装に使われている輸送問題のソルバー lp.transport は lpSolve パッケージに収録されています。詳細は lpSolve.pdf を見ていただくとして、インターフェースのみを次に抜粋します。

名前 説明
lp lp_solve線形・整数計画法システムへのインターフェース
lp.assign 割当問題を解くために特化した lp_solve へのインタフェース
lp.transport 輸送問題を解くために特化した lp_solve へのインターフェース
make.q8 8クイーン問題のための疎な制約行列の生成
lp.object lpオブジェクトの構造
print.lp lpオブジェクトをプリントするメソッド

線形計画法のソルバーとしては 汎用の lp割当問題 に特化した lp.assign輸送問題 に特化した lp.transport の 3つのインターフェースが定義されています。 いずれのインターフェースも lp.object で パラメータの受け渡しをしていますが、 その内訳を次に示します。

メンバー 説明
direction 入力された最適化の方向
x.count 目的関数の変数数
objective 目的関数の係数のベクトル
const.count 入力された制約の数
constraints 制約|入力されたとおりの制約行列(lp.assignやlp.transportでは返されません)
int.count 整数変数の数
int.vec 整数変数のインデックスのベクトル
objval 最適時の目的関数の値
solution 最適な係数のベクトル
num.bin.solns 返された解の数を数値で表示します
status 数字のインジケータです。0 = 成功, 2 = 実現可能な解決策がない

このパッケージ、 実はオープンソース線形計画法ソルバーCライブラリである lpSolve そのものです。 ですが、このライブラリのAPIはエントリーが多くて少々煩雑なので、 複数のエントリーを集約した3つのRコードを パッケージのインターフェースとして定義されているようです。

ちなみにライブラリAPIをそのままインタフェース化した lpSolveAPI も存在します。詳細は lpSolveAPI.pdf を確認してください。両方のパッケージともCコードは同じです。


線形計画法ソルバーライブラリ lpSolve

Rのパッケージが元はCライブラリだとわかって 「それじゃ直接コールすればいいじゃん」 と考えるのは今どきのプログラマの流儀ではないのかもしれませんが(笑)

ソルバーライブラリを探してみたところ、 なんと sourceforge(懐かしい)で、 今もメンテナンスが続いているようです。 ドキュメントを次に示します。

lpsolve.sourceforge.net

で、ライブリラリのビルドの手順を 書こうかと考えたのですが、 素敵なことに下記の Qiita の記事に 手順が書いてあったので、 こちらを参考にしてください。

qiita.com

この記事は lpSolve を Python に 取り込む方法を紹介していますが、 本家の lpSolve API バインディングを 使用しているようです。 元のブログ記事では Rパッケージを Python から使う事例が 紹介されてたのは 輸送問題のAPIである lp.transport が 使いたかったからですかねぇ。 「それ、むっちゃ効率が悪くないか?」 と思うのは僕が老人だからでしょうかねぇ?


Cライブラリの lpSolve を使って輸送問題を解く

ライブラリもビルドできたところで、 輸送問題を解くCプログラムを物色し始めました。 「わりと一般的だろうなぁ…」と想像して lpSolve を呼び出すCコードのサンプルを ググりまくったのですが、 日本語はもちろんのこと、英語でも見つからない。 「今どきのプログラマはもうCは書かないのかな?」 などと思いつつ、 輸送問題の解説 と先程の図をつきわせて 「lpSolve を使って輸送問題を解く」 勉強をしました。 (大学の情報演習みたいだ)

そもそも lpSolve で輸送問題 (というか最適化問題)を解くためには 目的関数と制約条件を設定する必要があります。 輸送問題の目的関数と制約条件を簡単に説明します。

f:id:Akito_Fujita:20210309195233p:plain

前述の図を再掲します。この図を見ると どうしても分布Pと分布Qの 都合7つのマルに目がいってしまうのですが、 輸送問題で注目するのはマルを繋ぐ都合12本の矢印だったりします。

つまり、輸送前は分布Pには図に示すとおり荷物が積まれており、分布Qには荷物は全くない状態だと考えます。 そして、輸送後は分布Pには荷物は全くなく、分布Qには図に示すとおりに荷物が積まれている状態になります。 では、輸送中は?というと都合12本の矢印に荷物が乗っている、つまり荷物が輸送されている状態です。 この時にコストを最小にするためにはどの矢印に荷物をどれくらい載せるべきか?を考えるのが 前回 紹介した「ヒッチコックの輸送問題」なんだそうです。

ここでコストと呼んでるのは分布内の各所の間の距離が違うからです。 例えばP1ーQ1間の距離とP1ーQ2間の距離は異なりますが、 この距離というのはどんな場合にも一定なので定数として扱えます。 各矢印毎のコストは距離と運ぶ荷物量を掛けた値となり、 輸送全体のコストは全ての矢印毎のコストを足しあわせたものと考えます。 前述の EMD の公式の分子部分はこの計算していますが、 これが輸送問題の「目的関数」でもあり、 輸送全体のコストが最小になるように 各矢印の間での運ぶ荷物量の調整を行います。

次に輸送問題の「制約条件」をですが、 これは分布Pと分布Qの各所が扱える荷物量に着目します。 例えば、P1はQ1、Q2、Q3に荷物を発送できますが、 その3箇所に発送する荷物量の合計は 0.4 以下に制約されます。 同様にP2、P3、P4も保有する荷物量で制約されます。 一方、Q1に着目するとP1、P2、P3、P4から荷物を受け取れますが、 その合計が 0.5 以上でなければなりません。 同様にQ2、Q3も受け取る荷物量の下限が制約されます。 このように輸送問題では分布の構成要素の総数分、 この例題では7つの「制約条件」が設定されます。

以上のように「目的関数」と「制約条件」が決められれば、 あとは lpSolve がこの問題を解いてくれます。 例題にそって条件を設定して lpSolve を呼び出す サンプルプログラム を末尾に掲載します*1。 先程ビルドした lpSolve からライブラリとヘッダーをコピーすれば、 次の手順をコンパイルできます。

$ ls inc
lp_Hash.h   lp_lib.h    lp_mipbb.h  lp_utils.h
lp_SOS.h    lp_matrix.h lp_types.h
$ ls lib
liblpsolve55.a
$ c++ -I./inc -g -Wno-macro-redefined -Wno-format -Wno-c++11-compat-deprecated-writable-strings -c example.cpp
$ c++ -L./lib -llpsolve55 -o example example.o
$ ./example
obj_val: 160.542763
variables:
0.400000
0.000000
0.000000
0.100000
0.200000
0.000000
0.000000
0.000000
0.200000
0.000000
0.100000
0.000000
$ 

ようやく冒頭のブログ記事と同じ EMD の計算値 160.542763 が確認できました。 続く variables として表示されてる12の数値は各矢印毎の輸送する荷物量、 つまり分布Pから分布Qへ荷物を輸送する差配を示しています。

しかし、この解を見てると、なんだか人間の代わりに コストが最小になるように lpSolve が考えてくれたように僕には見えてしまいます。 ある意味では人工知能の正体を覗いている感覚になります。 やはり「人工知能」という言葉やそのイメージは、 こういった数学的な解法を巧みにカモフラージュしてしまう 側面が否めないなぁと感じてしまいます。


ここまでのまとめ

lpSolve は1990年代から存在する歴史あるCライブラリで、 さまざまなプログラミング言語へのバインディングを始め、 機能が充実しています。例えば example.cpp の85行目をコメントアウトしてもらうと、 次のように目的関数と制約条件を書いた lp フォーマットのファイルが生成されます。

/* Objective function */
min: +109.927248669 C1 +97.2830920561 C2 +352.90083593 C3 +211.955183942 C4 +195.971936766 C5 +348.09481467 C6
 +244.180261283 C7 +115.429632244 C8 +254.909787964 C9 +141.435497666 C10 +52 C11 +334.752147118 C12;

/* Constraints */
+C1 +C2 +C3 <= 0.4;
+C4 +C5 +C6 <= 0.3;
+C7 +C8 +C9 <= 0.2;
+C10 +C11 +C12 <= 0.1;
+C1 +C4 +C7 +C10 >= 0.5;
+C2 +C5 +C8 +C11 >= 0.3;
+C3 +C6 +C9 +C12 >= 0.2;

このファイルはライブラリに付属する lp_solve コマンドの入力として使えます。 なので JavaScript と接続するには、 この lp フォーマットのファイルを生成して、 子プロセスで lp_solve コマンドを呼び出す方法が 最もお手軽かな?などと考えています。

ともあれ…

Rによる EMD 実装をなぞることにより Earth Mover's Distance は「輸送問題の解を正規化した指標」である ことがわかったことは収穫でした。 もっとも、今どき「輸送問題」をCで解こうとするとある穴全部にハマる… まぁ僕個人の良い勉強になったということにしておきましょう。

しかし、今どきのプログラミング言語のトレンドである 過去に実装されたコードをブラックボックスのまま取り込む風潮には、 正直「これで良いのか?」などと考えていたら、 次の書籍があと数日で出版されることを知りました。

www.hanmoto.com

僕の心配は杞憂で、この3週間は「くたびれ儲け」だった というオチが付いたところで本稿を一旦締めます。

以上

付録: example.cpp

#include <stdio.h>
#include <stdlib.h>
#include "lp_lib.h"

double f1[4][3] = { { 100,  40,  22},
            { 211,  20,   2},
            {  32, 190, 150},
            {   2, 100, 100} };
double f2[3][3] = { {   0,   0,   0},
            {  50, 100,  80},
            { 255, 255, 255} };
double w1[4] = { 0.4, 0.3, 0.2, 0.1 };
double w2[3] = { 0.5, 0.3, 0.2 };

#include <math.h>

double dist(double *F1, double *F2)
{
  double dx = F1[0] - F2[0];
  double dy = F1[1] - F2[1];
  double dz = F1[2] - F2[2];
  return(sqrt(dx*dx + dy*dy + dz*dz));
}

double a[13];

double b[13] = { 0,  //  <= 0.4
         1,  1,  1,
         0,  0,  0,
         0,  0,  0,
         0,  0,  0 };
double c[13] = { 0,  //  <= 0.3
         0,  0,  0,
         1,  1,  1,
         0,  0,  0,
         0,  0,  0 };
double d[13] = { 0,  //  <= 0.2
         0,  0,  0,
         0,  0,  0,
         1,  1,  1,
         0,  0,  0 };
double e[13] = { 0,  //  <= 0.1
         0,  0,  0,
         0,  0,  0,
         0,  0,  0,
         1,  1,  1 };


double f[13] = { 0,  //  <= 0.5
         1,  0,  0,
         1,  0,  0,
         1,  0,  0,
         1,  0,  0 };
double g[13] = { 0,  //  <= 0.3
         0,  1,  0,
         0,  1,  0,
         0,  1,  0,
         0,  1,  0 };
double h[13] = { 0,  //  <= 0.2
         0,  0,  1,
         0,  0,  1,
         0,  0,  1,
         0,  0,  1 };

void
init()
{
  a[0] = 0;
  for (int i = 0; i < 4; i++) {
    for (int j = 0; j < 3; j++) {
      int n = i*3 + j + 1;
      a[n] = dist(f1[i], f2[j]);
    }
  }
}

int
main()
{
  init();

  lprec *lp;
  int ret;

  lp = make_lp(0, 12);
  if (lp == NULL) exit(1);

  ret = set_obj_fn(lp, a);
  if (ret == 0) exit(1);

  add_constraint(lp, b, LE, 0.4);
  add_constraint(lp, c, LE, 0.3);
  add_constraint(lp, d, LE, 0.2);
  add_constraint(lp, e, LE, 0.1);
  add_constraint(lp, f, GE, 0.5);
  add_constraint(lp, g, GE, 0.3);
  add_constraint(lp, h, GE, 0.2);

  //write_lp(lp, "example.lp");

  set_verbose(lp, 1); // CRITICAL; NORMAL = 4; FULL = 6;
  ret = solve(lp);
  if (ret != 0) {
    printf("status: %d\n", ret);
    print_lp(lp);
    exit(1);
  }

  double obj_val = get_objective(lp);
  printf("obj_val: %f\n", obj_val);

  double var[100];
  get_variables(lp, var);
  printf("variables:\n");
  for (int i = 0; i < 12; i++) {
    printf("%f\n", var[i]);
  }

  delete_lp(lp);
  exit(0);
}

*1:.cpp となってますが、 これはコードの途中で変数宣言をしたかったからだけで、 頭から最後までCで記述しています。

Earth Mover's Distance(1)古くて新しいアルゴリズム

Earth Mover's Distance (1)
Old but New Algorithms


2021/02/08
藤田昭人


本稿では Earth Mover's Distance というアルゴリズムを紹介します。 このアルゴリズムWikipedia 日本語版にはページがないので英語版から引用しますと…

In statistics, the earth mover's distance (EMD) is a measure of the distance between two probability distributions over a region D. In mathematics, this is known as the Wasserstein metric. Informally, if the distributions are interpreted as two different ways of piling up a certain amount of dirt over the region D, the EMD is the minimum cost of turning one pile into the other; where the cost is assumed to be amount of dirt moved times the distance by which it is moved.

統計学では、Earth mover's distance (EMD)は、領域D上の2つの確率分布間の距離の尺度である。 数学では、これはワッサーシュタイン計量法として知られている。 直感的に説明すると、分布が領域D上にある量の土を積み上げる2つの異なる方法であると解釈される場合、 EMDは、一方の土をもう一方の土に変えるための最小コストであり、 コストは土の移動量に移動距離をかけたものと仮定される。

確かに "Informally, " とは書いてありますが、数学がダメダメな僕の直感にはちっとも響かない。 こういう時は "History" の項目に目を移すのが僕の常套手段です。


「モンジュの問題」

元々は「1781年に ガスパール・モンジュGaspard Monge) が(数学の)輸送理論の文脈で初めて導入した」とのこと。 そこで「ガスパール・モンジュ」「 輸送理論」でググってみると 「モンジュの最適輸送問題をめぐる話題について」 というスライドがひっかかりました。そこから「モンジュの最適輸送問題」で調べまくったところ次の記事に行き当たりました。

fukuoka-u.repo.nii.ac.jp

この記事は福岡大学の桑江先生がお書きになった「研究雑談」なのですが、 何が素敵かというと「数学者が書いた文章なのに数式が全く出てこない」ってところです。

本項では、記事の中ほどにある「モンジュの問題」の説明を杖に以降の話を進めていきます。


元々は「砂山を移動させる」問題だった

桑江先生の説明では 「モンジュはラプラスフーリエと並んで ナポレオンに仕えた数学者の一人であり、 真面目で正義感の強い人物でしたが、 それが災いしてか、 ナポレオンに最後まで信奉したために 最後には悲惨な末路を迎えた逸話が残っています」 とのことなんですが、Wikipedia でのモンジュの経歴をみると数学に秀でた軍人だったように見えます。

ja.wikipedia.org

そのように理解すると、 当時ヨーロッパで最強の陸軍だった ナポレオンの軍団の上級将校だったとの想像から、 問題1「ある砂山をそれと同じ体積の穴に移したい。 砂粒の移動には移動距離に依存したコストがかかる時、最適な移動のさせ方は何か?」 という問いも非常にイメージしやすくなります。 陣地を構築するために大量の砂を 輸送する必要があったとか、 騎馬を進めるために行く手にある 大量の穴ぼこを全部埋める必要があるとか… と状況は腑に落ちます。

桑江先生によれば、 この問題「密に詰まった無限の各砂山から 密に詰まった無限の各穴に総コストが最小となるような 1対1上への写像を決める問題」と理解できるそうで 「たいへんな難問」なんだそうです。 この説明、僕にはピンと来ないのですが、 それでも、山と積み上っている砂粒を 一粒ずつ行き先を決める関数を書いてコンピュータで 実行するとしても、いつ計算が終わるかわからない 無謀なプログラムが出来上がることぐらいは 想像が付きます。

ともあれ…

この問題が今から200年以上前の ナポレオンの時代に提示されていたことに、 まずはビックリっと言ったところでしょう。


問題を「荷物の移動」に変更すると…

そこで、写像の問題(らしい)「モンジュの問題」を 解くために、問題の方を置き換えるそうです。 問題2「n 個の工場と n 個の店舗があり、 各工場から各店舗にそれそれ一個の製品のみを移動させ距離に応じてコストがかかるとき、 総コストを最小にする移動のさせ方はなにか?」

これで「工場の散らばり(確率分布)と 店舗の散らばり(確率分布)が与えられているときに 移動の総コストが最小となるような 工場と店舗の配置の結合分布を求めよ」と 置き換えられたそうで、 確かに工場や店舗に収まっている荷物は 砂粒よりは(物理的に)ずっと大きいので 現実的なプログラムが書けそうです。

本稿の冒頭にある「領域D上の2つの確率分布」とは こういうことかぁ…と思った次第。 つまり Earth mover's distance (EMD) とは 「モンジュ・カントロヴィッチ問題(MK問題)」 を解くアルゴリズムという訳ですね。 ちなみにカントロヴィッチはこの方です。

ja.wikipedia.org

でもね、この問題2は「ヒッチコック型輸送問題」って言いませんでしたっけ? このヒッチコックって誰?カントロヴィッチとはどういう関係?


ヒッチコック型輸送問題」とは?

普通、ヒッチコックといえば映画監督の アルフレッド・ヒッチコック を思い浮かべる方が多いと思うのですが、 ここで登場するヒッチコックはこの方です。

en.wikipedia.org

どうやらMITの先生でいらっしゃったようですね。 彼自身がモンジュ由来の輸送問題に言及したのは、1941年の論文 "The Distribution of a Product from Several Sources to Numerous Localities(複数の供給元から複数の地域への製品の流通)" です。この論文が問題2が紹介された初めてのケースだったようです。 また「ヒッチコック型輸送問題」は "HITCHCOCK TRANPORTATION PROBLEM" という報告書でも紹介されています。

The transportation problem was first formulated by F. L. Hitchcock in 1941; he also gave a computational procedure, much akin to general simplex method, for solving the problem. Independently, during World War II, T. C. Koopmans arrived at the same problem in connection with his work as a member of the Joint Shipping Board. The Problem is thus frequently referred to as the Hitchcock-Koopmans problem.

輸送問題は1941年に F. L. Hitchcock によって最初に定式化され、問題を解くための一般的なシンプレックス法に近い計算手順を与えた。 それとは独立して、第二次世界大戦中、T. C. Koopmans は英米共同海運委員会のメンバーとしての仕事に関連して、同じ問題に到達した。 このような経緯から、この問題はヒッチコック・クープマンズ問題と呼ばれることが多い。

新たに登場したクープマンズとはこちら方です。

ja.wikipedia.org

ザッと経歴をみてみると… 1940年にオランダからアメリカへ移民した人、 名前から察するにユダヤ系のようなので ナチス・ドイツの迫害を逃れてきたようですね。 オランダでは理論物理学者として 研究をしていたようですが、 移民を契機に統計学に転向し、 その後、経済学者として活躍した人のようです。

しかし、次から次へと新たな人物が登場するのは、 数学に関わるトピックだからでしょうか? ちょっとキリがない気がしてきました。 ここで探索のアプローチを変えて、 最後の報告書に着目します。 というのも、この文書が保管されていたサイトが手がかりになりそうだから。

このサイトはアメリカの国防総省が管轄した過去の委託研究の報告書のアーカイバーです。主にARPA (国防高等研究計画局) の委託研究、例えばインターネットの技術基盤となった ARPANETの開発時の研究報告などが収蔵されているわけなのですが、 ARPAが設立されたのが1958年ですので、 この報告書はそれ以前のものです。

つまり「ヒッチコック型輸送問題」とは 国防総省が関心を持つような研究だったことになります。


3人の接点

そこでカントロヴィッチ、ヒッチコック、 クープマンズの3人に何かしらの接点がないかと探したところ、Wikipedia線形計画法 の英語版 Linear programming のページで次の記述を見つけました。

In 1939 a linear programming formulation of a problem that is equivalent to the general linear programming problem was given by the Soviet mathematician and economist Leonid Kantorovich, who also proposed a method for solving it. It is a way he developed, during World War II, to plan expenditures and returns in order to reduce costs of the army and to increase losses imposed on the enemy. Kantorovich's work was initially neglected in the USSR. About the same time as Kantorovich, the Dutch-American economist T. C. Koopmans formulated classical economic problems as linear programs. Kantorovich and Koopmans later shared the 1975 Nobel prize in economics. In 1941, Frank Lauren Hitchcock also formulated transportation problems as linear programs and gave a solution very similar to the later simplex method. Hitchcock had died in 1957 and the Nobel prize is not awarded posthumously.

1939年、ソ連の数学者で経済学者のレオニード・カントロヴィッチが 一般的な線形計画問題に相当する問題を線形計画法で定式化し、 彼もその解法を提案した。 それは、第二次世界大戦中に、軍隊のコストを削減し、敵に課せられる損失を増やすために、 支出とリターンを計画するために彼が開発した方法である。 カントロヴィッチの研究は、当初ソ連では軽視されていた。 カントロヴィッチとほぼ同時期に、 オランダ系アメリカ人の経済学者T.C.クープマンスが、古典的な経済問題を線形プログラムとして定式化した。 カントロヴィッチとクープマンスは、後に1975年のノーベル経済学賞を共同受賞した。 1941年には、フランク・ローレン・ヒッチコックもまた、交通問題を線形プログラムとして定式化し、 後のシンプレックス法と非常によく似た解を与えました。 ヒッチコックは1957年に死去しており、ノーベル賞は死後に授与されていない。

どうやらカントロヴィッチ、ヒッチコック、 クープマンズの3人は互いに独立して、 後に線形計画法と呼ばれる研究に 先鞭をつけていたようですね。

線形計画法といえば ジョージ・ダンツィーグGeorge Dantzig) の シンプレックス法 を思い出します。 ダンツィーグがこの技法を発表したのは1947年のことですが、 カントロヴィッチ、ヒッチコック、クープマンズの3人の取り組みは、 概ねその10年近く前の第2次世界大戦が勃発した前後のことだったようです。

つまり、今日オペレーションズ・リサーチ(OR)と呼ばれる研究分野を 切り開いたのはこの3人ということになるわけですね。 なかでも、1939年に「線形計画法」を体系的に語った ”The Mathematical Method of Production Planning and Organization(生産計画と組織の数理的方法)" を発表したカントロヴィチは、 今では「線形計画法」の創始者として 広く認知されてようです。

いや、しかし、 ORがロシア起源だったとは知らなかったなぁ…


オペレーションズ・リサーチの歴史

オペレーションズ・リサーチOperations Research) は「数学的・統計的モデル、 アルゴリズムの利用などによって、 さまざまな計画に際して 最も効率的になるよう決定する科学的技法」 と定義されています。 僕は第2次世界大戦後に アメリカで生まれた学問だと 思い込んでいたのですが、 その歴史 を辿ってみると第2次世界大戦中の 「軍事作戦のための兵站計画の立案」 を起源としていたようです。 そもそも「モンジュの問題」自体も この兵站計画の範疇に収まりますからねぇ。

第2次世界大戦後、 この技法は「線形計画法」という名前を冠し、 オペレーションズ・リサーチの技法のひとつとして 非軍事的な用途へも急速に拡大していった訳です。 でも、そこは何分にも冷戦真っ只中の時代のことですから、 技法の由来や起源については 各者が都合よく解釈した…それが、 数学、統計学、経済学の間で史観が微妙に異なり、 同じ問題に(関わったとされる)さまざまな人物の名前が 付けられることになった理由なんじゃないかと 僕は想像しています。

桑江先生も紹介しておられるように、 カントロヴィチは1975年にノーベル賞を受賞しています。 が、それはクープマンズと共同でノーベル経済学賞の受賞でした。

www.nobelprize.org

その受賞理由は次のとおり…

The Sveriges Riksbank Prize in Economic Sciences in Memory of Alfred Nobel 1975 was awarded jointly to Leonid Vitaliyevich Kantorovich and Tjalling C. Koopmans "for their contributions to the theory of optimum allocation of resources."

1975年のアルフレッド・ノーベルを記念したスヴェリヒス・リクスバンク経済学賞は、 「資源の最適配分理論への貢献に対して」 レオニード・ヴィタリェヴィチ・カントロヴィッチとティヤリング・クープマンズに共同で授与されました。

一般にノーベル賞は「表彰されるまで待たされる」賞であることで有名ですが、 1975年は冷戦にくたびれていた米ソ両陣営が デタント に向かっていた時期でもあり、 両陣営を代表する学者が共同受賞することで米ソ融和をアピールする狙いがあったのかも? それが彼らが35年間も待たされた 理由かもしれません。


「Earth Mover's Distance」という新しい器

ここまで長々と説明して来た 「モンジュ・カントロヴィッチ問題」
ある意味では前世紀に決着済みのはずのこの問題には、 実はその後があります。 そもそも「砂山を運ぶためのコスト」を計算するために 考え出されたはずのこのアルゴリズムですが、 全く別の思いもよらない使い方があったのでした。

その新しい使い方を紹介しているのが論文 "The Earth Mover’s Distance as a Metric for Image Retrieval " (画像検索のための指標としてのEarth Mover’s Distance ) です。2000年に発表されたこの論文では、 なんと類似画像検索での類似度判定にEMDを使うこと提案してます。

そもそも画像検索には「2つの画像がどれくらい似ているかを判定する」アルゴリズムが必要になりますが、 この「画像が似ている」を判定するのは難しいですよね? 画像の形に着目するのか?色に着目するのか? はたまた「画像の一部分だけがそっくり」といった例もあるでしょう。 人間が「似ている」と感じることを判定するという意味ではAI的ですし、 近年では機械学習が大きな成果を上げている研究分野でもあります。

論文の要約は次のとおりです。

We investigate the properties of a metric between two distributions, the Earth Mover’s Distance (EMD), for content-based image retrieval. The EMD is based on the minimal cost that must be paid to transform one distribution into the other, in a precise sense, and was first proposed for certain vision problems by Peleg, Werman, and Rom. For image retrieval, we combine this idea with a representation scheme for distributions that is based on vector quantization. This combination leads to an image comparison framework that often accounts for perceptual similarity better than other previously proposed methods. The EMD is based on a solution to the transportation problem from linear optimization, for which efficient algorithms are available, and also allows naturally for partial matching. It is more robust than histogram matching techniques, in that it can operate on variable-length representations of the distributions that avoid quantization and other binning problems typical of histograms. When used to compare distributions with the same overall mass, the EMD is a true metric. In this paper we focus on applications to color and texture, and we compare the retrieval performance of the EMD with that of other distances.

本研究では、コンテンツベースの画像検索における2つの分布間のメトリックであるEarth Mover's Distance (EMD)の特性を調べる。 EMDは、正確な意味で、一方の分布を他方の分布に変換するために支払わなければならない最小コストに基づいており、 Peleg, Werman, Romによって、ある種の視覚問題に対して最初に提案された。 画像検索のために、我々はこの考えを、ベクトル量子化に基づいた分布の表現スキームと組み合わせる。 この組み合わせにより、以前に提案された他の手法よりも知覚的類似度をよく考慮した画像比較のフレームワークが得られる。 EMDは、効率的なアルゴリズムが利用可能な線形最適化からの輸送問題の解に基づいており、 部分的なマッチングも自然に行うことができる。 EMDはヒストグラムマッチング手法よりもロバストであり, ヒストグラムに典型的な量子化やその他のビニングの問題を回避した分布の可変長表現で動作することができる。 同じ全体の質量を持つ分布を比較するために使用される場合、EMDは真のメトリックである。 この論文では、色とテクスチャへの応用に焦点を当て、EMDの検索性能を他の距離と比較する。

どうやら、この提案は比較する2つの画像の各々から特徴を抽出した後、 画像がどれくらい似ているかを数値として表現するためにEMDを使うようです。 さらに詳しい説明は論文を読んでもらうとして…

ちなみに、この論文では Earth Mover’s Distance という名前の由来も紹介されています。

We give the name of Earth Mover’s Distance (EMD), suggested by Stolfi (1994), to this metric in this new context. The transportation problem is to find the minimal cost that must be paid to transform one distribution into the other. The EMD is based on a solution to the transportation problem for which efficient algorithms are available, and it has many desirable properties for image retrieval, as we will see.

Stolfi (1994)によって提案された Earth Mover's Distance (EMD)という名前を、この新しいコンテキストにおけるこのメトリックに与える。 輸送問題は、ある分布を他の分布に変換するために支払わなければならない最小コストを見つけることである。 EMDは、効率的なアルゴリズムが利用可能な輸送問題の解決策に基づいており、画像検索に多くの望ましい特性を持っている。

「Stolfi (1994)によって提案された」とありますが、参考文献には ”Stolfi*1 , J. 1994. Personal communication. ” とあるので、 コンピュータ・ビジョンの専門家の間での 仲間うちでのトークで登場した呼び名なのでしょうか?

しかし、モンジュの問題を踏まえた洒落た名前ですよねぇ。


とりあえず中締めを…

本稿では、ガスパール・モンジュが提起した「モンジュの問題」から "Earth Mover’s Distance" までを (本職の数学者やコンピュータ・ビジョンの専門家に突っ込まれないよう気をつけて) 紹介して来ました。高校数学も怪しい僕にはこのあたりが限界なのですが、 同じテーマで数学者が書いた書籍も出ているようです。

www.springer.com

それから EMD といえばコンピュータ・ビジョンの専門家にはよく知られているアルゴリズムだそうですが、 画像認識に関わる基礎技術でもあるので機械学習でも頻繁に登場するようです。

とはいえ、これ以上数式は勘弁して欲しい僕としては、 以下のC言語のコードでアプローチしたい思ってます。

users.cs.duke.edu

以上

*1:Stofi とはこの人です。

en.wikipedia.org

アフィン演算 の考案者として有名なんだとか…

Word2Vec(2)distance.js

JavaScript implementation of Word2Vec


2021/01/27
藤田昭人


本稿は 前回 の続編です。

Googleのオリジナル実装を使うと比較的お手軽に Word2Vec が使えることがわかりました。が、 BookBotJavaScript 専用 PaaS である Glitch で動いているので(コーパスの学習は word2vec コマンドを使うとしても)word2vec の ご利益を得るには、 学習済みデータへアクセスする JavaScript のコードが必要になります。

で、最近では JavaScript にも慣れてきたので 「サクッと作れるだろう…」と 当初は甘くみていたのですが、さにあらず。 UNIXプログラミングでは 使い慣れてたはずのストリームの扱いに ハマるハマる。 本稿ではそのドタバタの結果を 紹介したいと思います。


学習済みデータが取り込めない

まずは nodejs の常套句のような次の2行で、前回紹介した 東北大学の日本語 word2vec 学習済みデータ を取り込もうとしたのですが…

const fs = require('fs');
var buf = fs.readFileSync('entity_vector/entity_vector.model.txt', 'utf8');

これを nodejs で実行すると…

$ node a.js
buffer.js:608
    slice: (buf, start, end) => buf.utf8Slice(start, end),
                                    ^

Error: Cannot create a string longer than 0x1fffffe8 characters
    at Object.slice (buffer.js:608:37)
    at Buffer.toString (buffer.js:805:14)
    at Object.readFileSync (fs.js:421:41)
    at Object.<anonymous> (/Users/fujita/xtr/BookBot/WikiEntVec/a.js:2:14)
    at Module._compile (internal/modules/cjs/loader.js:1063:30)
    at Object.Module._extensions..js (internal/modules/cjs/loader.js:1092:10)
    at Module.load (internal/modules/cjs/loader.js:928:32)
    at Function.Module._load (internal/modules/cjs/loader.js:769:14)
    at Function.executeUserEntryPoint [as runMain] (internal/modules/run_main.js:72:12)
    at internal/main/run_main_module.js:17:47 {
  code: 'ERR_STRING_TOO_LONG'
}
$ 

…ズッコケます。

使ったことのある方はご存知でしょうが、 nodejs では fs.readFileSync の第2引数に文字コードを指定しておくと、 ファイルから取り込んだデータを指定した文字コードの文字列に変換してくれるのです。 その変換途中でどうやらズッコケているようです。

そこで読み込んでいるテキストファイルのサイズを確認すると…

$ ls -l entity_vector/
total 5475416
-rw-r--r--@ 1 fujita  staff   830354941  2 17  2017 entity_vector.model.bin
-rw-r--r--@ 1 fujita  staff  1951880197  2 17  2017 entity_vector.model.txt

…2GB弱の巨大なテキストファイルでした。 (Wikipedia日本語版の全文だがら当然か)

ちなみに nodejs の場合 readFileSync など インメモリに取り込むファイルのサイズの上限は2GBです。 このファイルの場合は2GBを超えてはないので、 文字コードを指定しなければ(readFileSyncを第1引数だけで呼ぶ)エラーは出ませんが、 その後文字コードの変換ができないので「万事休す」ことに変わりはありません。


nodejs でのストリーム操作

大容量ファイルへの対処方法は「ストリームによる処理」が定石と言われてます。 もちろん javascript もストリーム処理をサポートしており「nodejs ストリーム」で ググるとそのための情報をたくさん集められます。 例えば、次のページは Stream API を網羅的に紹介してくれているので便利…

qiita.com

…なんですが、javascript の Stream API は非同期で実行されるハンドラを 引き渡す仕様になっていて、そのハンドラの書き方の解説を見つけるの ひと苦労します。

そこでハンドラを書いてみました。 次のソースは2GB以上のファイルを1行ずつリードしてコンソールに表示する、 UNIXの cat コマンドみたいな javascript コードです。

const fs = require('fs');
const readline = require('readline');

var file = 'entity_vector/entity_vector.model.txt';

const reader = async function(rl) {
  for await (const chunk of rl) {
    console.log(chunk);
  }
}

async function main() {
  const rs = fs.createReadStream(file, 'utf8');
  const rl = readline.createInterface(rs, {});
  await reader(rl);
}

main();

このプログラムでは関数 reader が非同期実行できるハンドラなんですが、 関数 main では reader が全てのストリームを受け取って処理し終えるまで 待ってくれます。javascript の非同期実行といえば 昔からある promise/then を使ったコードをよく見かけますが、 このプログラムでは async/await を使ってます。

それからストリームから1行ずつデータを取り出す readline も使ってます。 関数 reader の for await of ループは UNIX の gets を使ったループと等価で console.log(chunk); で1行ごとコンソールに出力してます。 この部分を任意に書き換えて UNIX の フィルター系コマンドと 等価なコードを書くことができます。


コサイン類似度

さて、オリジナル実装に含まれている distance コマンド。 このコマンドは任意の単語を指定すると、 その単語と意味的に距離が近い単語を距離の近い順にソートして、上位40個表示してくれます。 この「意味的な距離が近い」計算とはコサイン類似度の計算で、 実は「2つの単語ベクトル間の角度」を求めています。 計算結果は1〜−1の範囲の値をとり、 1なら角度ゼロ(完全一致)、 0なら角度は90°(意味的には無関係) −1なら角度は180°(意味的には正反対?) を意味します。

Wikipedia日本語版をチェックしてみたのですが、 「コサイン類似度」というページはない( 「ベクトルのなす角」 というページ転送されます)ので、 代わりに次のページを参考にしてください。

mathtrain.jp

ちなみに、Wikipedia英語版では "Cosine similarity” というページで 「情報検索やテキストマイニングでは2つの文書がその主題に関してどれだけ類似しているかを示す有用な尺度となる」 と紹介されています。

さて、肝心のコサイン類似度の計算方法。 当初は「オリジナル実装のコードを見れば」 と考えていたのですが……… とても読解できない「古き悪しき実装」でした。 やむなく上記のページの公式を見ながら 次のコードを書きました。

function Absolute(m) {
  var ret = 0;
  for (var i = 0; i < m.length; i++) {
    ret += m[i] * m[i];
  }
  return(Math.sqrt(ret));
}

function DotProduct(m1, m2) {
  var ret = 0;
  for (var i = 0; i < m1.length; i++) {
    ret += m1[i] * m2[i];
  }
  return(ret);
}

function CosineSimilarity(dot, nrm1, nrm2) {
  return(dot/(nrm1*nrm2));
}

コサイン類似度は2つのベクトルの 内積ドット積) を各ベクトルの 絶対値 で割ることで求めます。この絶対値は各々のベクトルの長さを表し ノルム と呼ばれることもあります。 ベクトルの絶対値、内積、コサイン類似度に対応する AbsoluteDotProductCosineSimilarity を定義しました。


JavaScript 版 distance

以上のコードを組み合わせてJavaScript版 distance.js を作成しました。
ソースコードは末尾に添付しておきます。

オリジナル実装との違いは 次のようにターゲットとなる単語を引数で指定するところです。

$ ./distance.js entity_vector/entity_vector.model.txt '[アラン・チューリング]'

Word: [アラン・チューリング]
1015474 times
                                                Word       Cosine distance
--------------------------------------------------------------------------
[ジョン・フォン・ノイマン]                                      0.8238200737420581
[クロード・シャノン]                                          0.7706864375586782
[ジョン・マッカーシー]                                        0.7634118658850984
チューリング                                               0.7550014937957042
[クルト・ゲーデル]                                          0.7505430963870695
[ハーバート・サイモン]                                       0.7503034618503429
[ダフィット・ヒルベルト]                                       0.7381630672368353
[アルベルト・アインシュタイン]                                  0.7370111227611059
[スタニスワフ・ウラム]                                        0.7355724275046345
[ヴォルフガング・パウリ]                                       0.7320032112916035
[エンリコ・フェルミ]                                           0.729639429509676
[マックス・ボルン]                                          0.7238478243360493
[アロンゾ・チャーチ]                                        0.7203194651116861
[ロバート・オッペンハイマー]                                   0.7152767971441714
[ポール・ディラック]                                         0.7129907170231481
[エミー・ネーター]                                          0.7093252080731395
[ヘルマン・ワイル]                                          0.708357522309429
[エルヴィン・シュレーディンガー]                                 0.707707742898606
[ノーバート・ウィーナー]                                       0.7072144624267143
[バートランド・ラッセル]                                       0.7030623602860854
[ロジャー・ペンローズ]                                        0.7027514788704126
[ゴッドフレイ・ハロルド・ハーディ]                                0.7015682037971648
[レフ・ランダウ]                                             0.6995820039570446
アインシュタイン                                            0.6988088557431229
[アンドレ・ヴェイユ]                                         0.6978330851832386
[ヴェルナー・ハイゼンベルク]                                   0.6968193353952211
[ライナス・ポーリング]                                        0.696506621225639
数学者                                                   0.6874893585915669
[マービン・ミンスキー]                                        0.6872785430574049
[フリーマン・ダイソン]                                        0.6863152690167025
[アンリ・ポアンカレ]                                           0.6851716370245848
[ユージン・ウィグナー]                                        0.6812997448505773
[アレン・ニューウェル]                                        0.6774566693465804
[ハンス・ベーテ]                                             0.6767724362750837
[アルバート・アインシュタイン]                                    0.6750565951905274
[フリッツ・ロンドン]                                           0.6739032952575766
[リチャード・P・ファインマン]                                      0.6737536291870052
[ニールス・ボーア]                                          0.6737340516366506
[朝永振一郎]                                               0.6729520531327631
[ゴットロープ・フレーゲ]                                         0.6717680358048904
$ 

コンソール表示を見る限り オリジナル実装と遜色ない結果がでるのですが… とにかく遅い。これは教育済みデータを (ターゲットとなる 単語ベクトルを見つけるために1回、 ターゲットの単語とその他の単語との コサイン類似度を計算するためにもう1回の) 都合2回走査しているためです。


まとめ

JavaScriptスクリプト言語のなかでも高速な部類なんですが、 そのパワーを持ってしてもビッグデータを扱うことが難しいことを実感しました。 これが JavaScript機械学習系のコード実装が少ない理由なのかもしれません。

東北大学の日本語 word2vec 学習済みデータは十分にビッグデータで、 何か工夫をして扱うデータのサイズを削減しないと、 JavaScript本来のスピードで処理できないと言わざる得ません。

新しい技術的課題を見つけちゃったなぁ…

以上

PS 前回と今回のソースコードGithub にアップしました。

github.com



付録1.distance.js

#!/usr/bin/env node

var file = process.argv[2];
var keyword = process.argv[3];

const fs = require('fs');
const readline = require('readline');

function Absolute(m) {
  var ret = 0;
  for (var i = 0; i < m.length; i++) {
    ret += m[i] * m[i];
  }
  return(Math.sqrt(ret));
}

function DotProduct(m1, m2) {
  var ret = 0;
  for (var i = 0; i < m1.length; i++) {
    ret += m1[i] * m2[i];
  }
  return(ret);
}

function CosineSimilarity(dot, nrm1, nrm2) {
  return(dot/(nrm1*nrm2));
}

var n1 = 0;
var a = {};

const reader1 = async function(rs1) {
  for await (const chunk of rs1) {
    var elm = chunk.split(' ');
    if (elm[0] == keyword) {
      a = {};
      a.key = elm[0];
      a.mtx = [];
      for (var i = 1; i < elm.length; i++) {
        a.mtx.push(parseFloat(elm[i]));
      }
      a.nrm = Absolute(a.mtx);
      break;
    }
    n1++;
  }
};

var n2 = 0;
var distance = [];
const reader2 = async function(rs2) {
  for await (const chunk of rs2) {
    var elm = chunk.split(' ');
    if (elm.length > 2 && elm[0] != keyword) {
      var b = {};
      b.key = elm[0];
      b.mtx = [];
      for (var i = 1; i < elm.length; i++) {
        b.mtx.push(parseFloat(elm[i]));
      }
      b.nrm = Absolute(b.mtx);
      dot = DotProduct(a.mtx, b.mtx);
      sim = CosineSimilarity(dot, a.nrm, b.nrm);
      distance.push({ 'key': b.key, 'sim': sim });
    }
    process.stderr.write(n2+" times\r");
    n2++;
  }
};

const N = 40; // number of closest words that will be shown

function Compare(a, b) {
  return(b.sim - a.sim);
}

async function main() {
  const rs1 = fs.createReadStream(file, { encoding: 'utf8' });
  const rl1 = readline.createInterface(rs1, {});
  await reader1(rl1);
  if (!a.key) {
    console.log("Out of dictionary word!");
    return;
  }
  console.log("\nWord: %s", a.key);
  const rs2 = fs.createReadStream(file, { encoding: 'utf8' });
  const rl2 = readline.createInterface(rs2, {});
  await reader2(rl2);
  distance.sort(Compare);
  console.log("\n                                                Word       Cosine distance\n--------------------------------------------------------------------------");
  for (var i = 0; i < N; i++) {
    var pad = '                                                  ';
    var str = (distance[i].key+pad).slice(0, 50);
    console.log("%s\t%f", str, distance[i].sim);
  }
}

main();