静かなる名辞

pythonとプログラミングのこと


【python】scipyのpdistとsquareformの使い方と仕組み

python scipy Tips

はじめに

　scipyで距離行列を扱うときはscipy.spatial.distanceのpdist, squareformなどを主に使いますが、長年よくわからないままに使っていたので、整理してまとめておきます。

　なお、以下のドキュメントを参考にします。
scipy.spatial.distance.pdist — SciPy v1.2.1 Reference Guide
scipy.spatial.distance.squareform — SciPy v1.2.1 Reference Guide

pdist

　pdistを使うと様々な距離の距離行列を計算できます。デフォルトはユークリッド距離ですが、他にもいろいろ指定できる距離の種類があります。詳細はドキュメントを見てください。

　そしてpdistは「condensed distance matrix」という表現方式の距離行列を返します。この方式にはあまり馴染みがありませんが、このように使えます。

>>> import numpy as np
>>> a = np.random.random(size=(5,3))
>>> from scipy.spatial.distance import pdist
>>> pdist(a)
array([0.58791701, 0.71055164, 0.46750397, 0.3389748 , 0.39126501,
       0.28675164, 0.79462635, 0.39543688, 0.93094202, 0.59666923])

　結果は常に1次元配列で帰ります。5つのベクトルの距離を計算してみた結果、要素数は10になっています。5*4/2=10なので、距離行列の上下どちらか半分の三角行列がflattenされて返っているのだとわかります。上下の三角行列は対称ですし、対角要素（自分同士の距離）は常に0なので、最低限これだけあれば表現できる訳です。

　結果として省メモリな表現になりますが、人間にとってはわかりづらくなります。また、関数やライブラリによってはこの形式を受け付けないこともままあります。そこで、次のsquareformが登場します。

squareform

　squareformはpdistで生成された配列と、馴染みのある距離行列（「square-form distance matrix」と称しています）の相互変換を担います。

　上の実行例から続けてこう打つと、距離行列が得られます。

>>> from scipy.spatial.distance import squareform
>>> squareform(pdist(a))
array([[0.        , 0.58791701, 0.71055164, 0.46750397, 0.3389748 ],
       [0.58791701, 0.        , 0.39126501, 0.28675164, 0.79462635],
       [0.71055164, 0.39126501, 0.        , 0.39543688, 0.93094202],
       [0.46750397, 0.28675164, 0.39543688, 0.        , 0.59666923],
       [0.3389748 , 0.79462635, 0.93094202, 0.59666923, 0.        ]])

　なお、squareformは1次元配列か2次元配列を引数として受け取ります。1次元配列を受け取ると「condensed distance matrix」→「square-form distance matrix」の変換、2次元配列を受け取ると「square-form distance matrix」→「condensed distance matrix」の変換として動くという仕組みです。要するに何が言いたいかというと、二回通すと元に戻ります。

>>> squareform(squareform(pdist(a)))
array([0.58791701, 0.71055164, 0.46750397, 0.3389748 , 0.39126501,
       0.28675164, 0.79462635, 0.39543688, 0.93094202, 0.59666923])

　便利かどうかはちょっと即答しかねますが、けっこう面白い気はします。

インデックスの変換

　さて、ここまで読んできた方は「squareformのインデックスとpdistのインデックスはどう対応付けられるの？」という興味を持ったかと思います。

　pdistのドキュメントには意味深な一文があります。

Returns a condensed distance matrix Y. For each and (where ),where m is the number of original observations. The metric dist(u=X[i], v=X[j]) is computed and stored in entry ij.

　結論から先にいうと、この文に深い意味はないので、無視して構いません。真に受けるとかえって混乱します。

　日本人で英語が苦手な私には上の一文の意図するところはよくわかりませんが、英語圏の方もわからないようで、stackoverflowでこんなスレッドを見つけました。

python - How does condensed distance matrix work? (pdist) - Stack Overflow

　基本的にはここの回答に書いてあることが参考になります。ただし、python2時代のスレッドなので、整数の割り算が/で記載されていることには注意する必要があります（python3では/を使うとfloatの結果になり、整数の結果を得たければ//を使う必要があります）。

　とはいえ、stackoverflow見てねというのも不親切なので、この記事でも説明しておきます。

　とりあえず4つの要素の距離の「square-form distance matrix」の側の三角行列のインデックスについて考えることにすると、以下のようになります。

上三角行列のインデックス
上三角行列のインデックス

　pdistで作る「condensed distance matrix」の場合は以下のようなものです。

pdistで作った配列のインデックス
pdistで作った配列のインデックス

　対応をまとめると、

0,1→0
0,2→1
0,3→2
1,2→3
1,3→4
2,3→5

　となります。

　以下は上記stackoverflowの回答に基づいた記述ですが、pythonの数式で変換するとすると、

n*j - j*(j+1)//2 + i - 1 - j

　こんなですが、覚える必要はないです。

　また、itertools.combinationsの結果と同じ順で返ったりします。

>>> from itertools import combinations
>>> list(combinations(range(4), 2))  # 4つから2つ取り出す
[(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]

　概念的にはこちらの方がわかりやすいですね。

まとめ

　scipyで距離行列の絡む処理をしようとすると避けて通れないpdistとsquareformですが、実はやっていることは単純です。慣れればそれほど困惑することもないので、気楽に使えるようになると良いと思います。