NumPyレシピ | 複数のベクトル同士の内積は成分計算するのが速い

Ciao！みなさんこんにちは！このブログでは主に
(1)pythonデータ解析,
(2)DTM音楽作成,
(3)お料理,
(4)博士転職
の4つのトピックについて発信しています。

今回はPythonデータ解析、複数のベクトル同士の内積を高速に計算する方法です。一組のベクトルの内積はNumPyのnp.dot()で計算できます。ところが、np.dot()は複数のベクトルの組み合わせについての内積を一度に計算することができないため、計算したいベクトルの数だけループを回す必要があります。このような計算を高速化するには成分計算を行う必要があります。

この記事を読めば、複数のベクトルの内積を高速計算する方法がわかります。

Kaiko

この記事はこんな人におすすめ

Pythonでデータ解析をしている
Numpyでベクトルの内積を高速に計算したい

Abstract | 成分計算とループ方向の工夫で高速化！
Background | やりたいこと
Issue | ループとnp.dot()を使う方法では時間がかかる
1. np.dot()はひとつの内積しか計算できない
2. ループとnp.dot()を使って1000個の内積を計算
Data | np.tensordot()を使うとループは排除できるが遅い
Method | 成分計算で実装するのが速い
Result | 3つの方法の比較
1. 複数のベクトルの内積を計算するなら成分計算で実装するのが吉
2. 今回作成したPythonコード
Conclusion | まとめ

Abstract | 成分計算とループ方向の工夫で高速化！

複数のベクトル同士の内積を計算するときには、成分計算で実装することでループを避けて高速化することができます。NumPyにはnp.dot()といったベクトル計算用の便利な関数が用意されています。ところが、np.dot()は一つのベクトル同士の内積の計算しか行なえません。複数のベクトル同士の内積を計算するには計算したいベクトルの数だけループを回す必要があり、速度が遅くなります。また、np.tensordot()といった一度に複数のベクトルの内積を計算できる関数もありますが、これも速度が遅くメモリ負荷も高いため高速化にはあまり役に立ちません。

そこで、成分計算で実装することで、ループを排除し複数のベクトル同士の内積の計算を高速化できます。予めベクトルの次元が定まっている場合にはループを完全に排除することができます。ベクトルの次元が定まっていない場合でも、ループをベクトルの数の方向ではなく成分方向に設定した上で各成分の成分計算を行うことでループ回数を最小限にとどめて高速化することができます。

Pythonのベクトルの内積計算の高速化でお悩みの方はぜひご参考に最後までご覧になってください！

Background | やりたいこと

ここでは、$n$個のベクトル$\vec{u}_i$と$\vec{v}_i$の内積を求めるということを行います。例えば、以下のように1000個の2次元ベクトル$\vec{u}$と$\vec{v}$の内積を求めることを考えます。すなわち1000個のベクトルペアについて内積を計算したいという状況です。

import numpy as np

# fix random seed
np.random.seed(0)

# get 5 vectors
u = np.random.normal(size=[1000,2]) # 1000 2-dim vectors
v = np.random.normal(size=[1000,2]) # 1000 2-dim vectors

# print first element of array
print('u:', u[0]) # 2-dim vector
print('v:', v[0]) # 2-dim vector

Issue | ループとnp.dot()を使う方法では時間がかかる

まずはループとnp.dot()を組み合わせる方法でやってみます。

np.dot()はひとつの内積しか計算できない

np.dotはNumPyのリファレンスにあるように、2つの引数にはそれぞれベクトルまたは行列を取り、一つの内積を計算します。以下のように実行することができます。

print(np.dot(u[0], v[0])) # inner prod. of two vectors

一方で、1000個の内積を一度に計算することはできません。以下のようにするとエラーになります。

print(np.dot(u, v)) #　error

ループとnp.dot()を使って1000個の内積を計算

np.dot()を使うなら、ループを使って1000回の計算を行う必要があります。以下のように行います。

import time

# time measurement befort calculation
t1 = time.process_time()

# prepare output array for inner product
inner_prod = np.zeros([len(u)]) # 1000 scalar array

# loop
for i in range(len(u)):
    inner_prod[i] = np.dot(u[i], v[i])
#

# time measurement after calculation
t2 = time.process_time()

# elapsed time
delta_t = t2 - t1 # sec
print('elapsed [ms]: ', delta_t * 10**3)

# print first 5 elements
print('inner prod:', inner_prod[:5])

後で比較するために、timeのprocess_time()というメソッドを使って時間を測定しておきます。私の環境では5から10ミリ秒程度となりました。

Data | np.tensordot()を使うとループは排除できるが遅い

1000個の内積を計算したい場合、np.tensordot()を使うとループを排除できます。NumPyのリファレンスにあるように、2つの多次元配列を引数にし、キーワード引数のaxesで和を取る方向を指定することで内積を作ることができます。今回の場合は1次元目の方向に和を取れば内積を作ることができます。

np.tensordot()をaxes=(1,1)で実行したときの振る舞い

np.tensordot()をaxes=(1,1)で実行すると、全部の組み合わせについての内積が計算されます。すなわち、

np.tensordot(u, v, axes = (1,1))

では以下が出力されます。
$$
\left[
\begin{array}[cccc]\\
\vec{u}[0] \cdot \vec{v}[0] & \vec{u}[0] \cdot \vec{v}[1] & \dots & \vec{u}[0] \cdot \vec{v}[n-1] \\
\vec{u}[1] \cdot \vec{v}[0] & \vec{u}[1] \cdot \vec{v}[1] & \dots & \vec{u}[1] \cdot \vec{v}[n-1] \\
\dots & \dots & \dots & \dots \\
\vec{u}[n-1] \cdot \vec{v}[0] & \vec{u}[n-1] \cdot \vec{v}[1] & \dots & \vec{u}[n-1] \cdot \vec{v}[n-1]
\end{array}
\right]
$$
したがって、np.tensordot(u, v, axes = (1,1))の対角成分を取れば、1000個のベクトル同士の内積を取り出すことができます。対角成分の取り出しにはnp.diag()を使うことができます。

np.tensordot()を使って複数のベクトル同士の内積を計算する方法

np.tensordot()とnp.diag()を使って複数のベクトル同士の内積を計算してみます。今回の1000個のベクトルの例では以下のようにします。

# time measurement befort calculation
t1 = time.process_time()

# get inner product by tensordot and diag
inner_prod1 = np.diag(np.tensordot(u, v, axes = (1, 1)))

# time measurement after calculation
t2 = time.process_time()

# elapsed time
delta_t = t2 - t1 # sec
print('elapsed [ms]: ', delta_t * 10**3)

# print first 5 elements
print('inner prod:', inner_prod1[:5])

私の環境では100,000回の計算に10から20ミリ秒程度かかりました。
念の為、ループとnp.dot()で計算したものと答えが一致していることを確かめておきます。以下のようにします。

print(np.max(np.fabs(inner_prod - inner_prod1)))

結果は

8.881784197001252e-16

となり、数値誤差の範囲で一致していることがわかります。

np.tensordot()の弱点 | メモリを食う&遅い

np.tensordot()を使うデメリットはメモリを余計に食うこととループを書かない割には遅いことです。np.tensordot()ではすべての組み合わせの内積（$\vec{u}[0] \cdot \vec{v}[1]$など）を計算してしまうので、メモリを余計に食う上に動作が遅いです。

今回は1000個のベクトルの例なので実行できましたが、例えば、最初の「n = 1000」のところを「n = 100000」（10万）にすると、私の環境ではメモリが確保できずにクラッシュします。興味があれば試しにやってみてください。

ループとnp.dot()を使う方法では、変数inner_prodの分、今回の例ではn = 1000要素分のメモリを確保すれば済みましたが、np.tensordot()を使う場合は1000$\times$1000で100万要素分のメモリを確保する必要があります。仮に、n = 100000（10万）ならば、np.tensordot()で確保されるメモリ領域は10万$\times$10万で100億要素分となります。

Method | 成分計算で実装するのが速い

複数のベクトルの内積は成分計算で実装するのが速いです。ただし、ベクトルの次元数が予めわかっていることが条件です。今回は2次元とわかっているので、比較的簡単に実装できます。

成分計算で内積を求める関数の準備

まずは成分計算で内積を求める関数を用意しておきます。

# function to get innder product of 2d vectors
def get_inner_prod_2d(u_s, v_s):
    """
    arguments:
        u_s <2d array [n,2]>: vectors
        v_s <2d array [n,2]>: vectors
    returns:
        inner_prod <1d array [n]>: inner product
    """
    return u_s[:,0] * v_s[:,0] + u_s[:,1] * v_s[:,1]

上記の関数ではn個の2次元ベクトルを2つ引数として与えると、n個の内積が返されます。ちなみに引数を「u_s, v_s」のように「_s」とつけているのは、これらが複数のベクトルであることを示すためです。

複数ベクトルの内積計算の実行

上記の関数を使って内積を計算します。以下のように行います。

# time measurement befort calculation
t1 = time.process_time()

# get inner product by tensordot and diag
inner_prod2 = get_inner_prod_2d(u, v)

# time measurement after calculation
t2 = time.process_time()

# elapsed time
delta_t = t2 - t1 # sec
print('elapsed [ms]: ', delta_t * 10**3)

# print first 5 elements
print('inner prod:', inner_prod2[:5])

私の環境では0.1から0.9ミリ秒程度と他の2つの方法の1/10程度に短縮されました。

念の為、ループとnp.dot()で求めたものと答えが一致していることを確認しておきます。

print(np.max(np.fabs(inner_prod - inner_prod2)))

私の環境では0.0となり、一致していることが確認できました。

ベクトルの次元が決まっていない場合への応用¶

ベクトルの次元が予め定まっていない場合、成分計算で高速化するには工夫が必要です。この場合、完全にループを排除することはできません。しかし、成分計算方向にループを回すことで、ベクトルの数が多い場合には、ベクトルの数方向にループを回してnp.dot()を使うよりも高速に計算できます。

次元未定の複数ベクトルの内積を成分計算する関数

以下のように成分計算方向にループを回しながら内積を計算する関数を作っておきます。

# function to get innder product of m-d vectors
def get_inner_prod_md(u_s, v_s):
    """
    arguments:
        u_s <2d array [n,m]>: m-dim vectors
        v_s <2d array [n,m]>: m-dim vectors
    returns:
        inner_prod <1d array [n]>: inner product
    """
    # get parameters
    n, m = np.shape(u_s) # number and dimension of vector
    
    # prepare output array
    inner_prod = np.zeros([n]) # <1d array [n]>
    
    # loop for dimension
    for j in range(m):
        inner_prod += u_s[:,j] * v_s[:,j]
    
    # end 
    return inner_prod

次元未定の複数ベクトルの内積を計算する方法

上記の関数を使って、次元が予め定まっていない複数ベクトルについての内積を計算することができます。今回はこれまでと同様に1000個の2次元ベクトルを使いますが、3次元以上でも同様です。以下のように実行します。

# time measurement befort calculation
t1 = time.process_time()

# get inner product by tensordot and diag
inner_prod3 = get_inner_prod_md(u, v)

# time measurement after calculation
t2 = time.process_time()

# elapsed time
delta_t = t2 - t1 # sec
print('elapsed [ms]: ', delta_t * 10**3)

# print first 5 elements
print('inner prod:', inner_prod3[:5])

私の環境では0.1から0.9ミリ秒程度でした。一つ前のループを全く使わない手法に遜色ない速度が得られました。

念の為、ループとnp.dot()で求めたものと答えが一致していることを確認しておきます。

print(np.max(np.fabs(inner_prod - inner_prod3)))

私の環境では0.0となり、一致していることが確認できました。

Result | 3つの方法の比較

今回見てきた3つの方法(+1つの応用)を比較しておきます。

#	方法	ループ	速度	メモリ負荷	柔軟性	メリット	デメリット
1	ループ+np.dot	多い (ベクトル数方向)	✗	○	○	単純	遅い
2	np.tensordot	なし	✗	✗	○	コードがシンプル	遅いベクトル数が多いとメモリがクラッシュする
3	成分計算 (次元既知)	なし	◎	○	✗	速い	次元が既知である必要がある
4	成分計算 (次元未定)	少ない (成分方向)	○	○	○	そこそこ速い	次元が高いと遅くなる可能性がある