地点相互の距離

Numpyではループを使うと遅くなるので，様々な技法が使われます。ベクトルや行列形式で書けば良いものは簡単なのですが，ブロードキャストを用いるものがあります。その例として，距離の行列の生成があります。複数の場所の位置がxに格納されているとして，相互の距離を計算するにはどのようにしたら良いでしょうか。

とりあえずループ

筆者はFortranに長年親しんでいるので，ループを使って書きたくなります。

do j = 1, n
  do i = 1, n
    r(i, j) = sqrt((x(i) - x(j))**2)
  end do
end do

Pythonに直訳すると

for i in range(len(x)):
    for j in range(len(x)):
        r[i, j] = math.sqrt((x[i] - x[j])**2)

となります。ndarrayにしても長さ10程度では差はなく，むしろリストのままの方が速くなりました。

ブロードキャスト

Numpyの（積も含む）配列演算は，基本的に要素毎です。4x3の2次元配列など同じ形状の配列の演算は，それぞれの要素に対して行われます。4x3の2次元配列と1x3の1次元配列のようなサイズが違うものに対しては，サイズが同一になるように1次元配列の行が繰り返されているとして計算が行われます。この変換を「ブロードキャスト」と呼びます。

位置を格納したxについては，その転置との演算を行えば配列が返ってきそうですが，そうはなりません。同一形状とみなされ0の羅列になってしまいます。

次元を増やす

ブロードキャストを使いループなしに求める距離行列を生成するには，片方の次元を増やします。次元を増やすには，np.newaxisを使います。Noneでも同様です。np.newaxisまたはNoneを指定した次元の長さは1になります。なお，筆者は使ったことがないのですが，Fortran95にはspreadという組込関数があります。

np.sqrt((x - x[:, np.newaxis])**2)

nx1と1xnの演算にしても良いでしょう。

np.sqrt((x[np.newaxis, :] - x[:, np.newaxis])**2)

筆者の環境では10地点に対して%timeitで計測したところ，ndarrayのループが69.1 μs，リストのループが54.1 μs，ループを回避して行列の次元を増やしたものは3.04 μsでした。1000地点だとndarrayのループが660 ms，行列の次元を増やしたものは2.47 msとなりました。

n次元化

...を使うと次元の数が配列の次元になるように:を並べたのと同じになります。n次元空間の点のユークリッド距離は次のように計算できます。

np.sqrt(((x[..., np.newaxis, :] - x[..., :, np.newaxis])**2).sum(axis=0))

参考

Scipy lecture notes 1.3.2.3 Broadcasting 和訳
scipy.interpolate.Rbfのソース