numpyで対称行列初期化速度を比較

はじめに

数値計算では対称行列を扱うことが多くあります。
今回はnumpyで対称行列の初期化方法の速度比較を行いました。

検証

比較内容

以下の2方法を比較しました。

方法1
配列の全要素にアクセスし値を格納する。

方法2
計算したい対称行列の上三角部分について要素ごとにアクセスして値を格納し、上三角行列Aに対して次式の計算で対称行列を計算する。

A + A^T - diag(A)

比較時は行列サイズを振って2方法を比較しました。

実行環境

• CPU: 2.5 GHz Intel Core i7
• python 3.6
• numpy 1.17.0
• matplotlib 3.1.1

コード

import time
import numpy as np
import matplotlib.pyplot as plt

def calc1(n_dim: int) -> tuple:
    """ 全配列要素を計算して対称行列を作成する。

    Args:
        n_dim (int): 対称行列のサイズ

    Returns:
        tuple : 計算時間[sec], 対称行列
    """

    t_start = time.time()

    mat = np.zeros((n_dim, n_dim))
    for i_row in range(n_dim):
        for i_col in range(n_dim):
            mat[i_row, i_col] = i_row + i_col

    return time.time() - t_start, mat

def calc2(n_dim: int) -> tuple:
    """ 上三角部分だけ計算して転置行列との和から対角成分を引いて対称行列を作成する。

    Args:
        n_dim (int): 対称行列のサイズ

    Returns:
        tuple : 計算時間[sec], 対称行列
    """

    t_start = time.time()

    mat = np.zeros((n_dim, n_dim))
    for i_row in range(n_dim):
        for i_col in range(i_row, n_dim):
            mat[i_row, i_col] = i_row + i_col

    return time.time() - t_start, mat + mat.T - np.diag(mat.diagonal())


if __name__ == "__main__":

    # 対角行列計算速度比較
    n_try = 10
    n_dim_list = [5, 20, 100, 500, 2000]
    t_calc1_list = []
    t_calc2_list = []

    # 計算比較
    for n_dim in n_dim_list:

        print("n_dim: {}".format(n_dim))

        # 方法1
        t_calc = 0
        for i_try in range(n_try):

            t_calc_, mat1 = calc1(n_dim)
            t_calc += t_calc_

        t_calc1_list.append(t_calc / n_try)
        print("elapsed_time:{} [sec]".format(t_calc / n_try))

        # 方法2
        t_calc = 0
        for i_try in range(n_try):

            t_calc_, mat2 = calc2(n_dim)
            t_calc += t_calc_

        t_calc2_list.append(t_calc / n_try)
        print("elapsed_time:{} [sec]".format(t_calc / n_try))
        print()

    # グラフの作図
    x = np.array(range(len(n_dim_list)))
    y = np.array(t_calc2_list) / np.array(t_calc1_list)
    label = [str(n_dim) for n_dim in n_dim_list]
    plt.bar(x, y, color='b', tick_label=label, label='t_calc2 / t_calc1', align='center')
    plt.show()

結果

計算時間の結果は下のようになりました。

n_dim: 5
elapsed_time:1.0991096496582032e-05 [sec]
elapsed_time:7.677078247070313e-06 [sec]

n_dim: 20
elapsed_time:6.659030914306641e-05 [sec]
elapsed_time:5.0449371337890624e-05 [sec]

n_dim: 100
elapsed_time:0.0015201091766357422 [sec]
elapsed_time:0.001001739501953125 [sec]

n_dim: 500
elapsed_time:0.03668978214263916 [sec]
elapsed_time:0.018210577964782714 [sec]

n_dim: 2000
elapsed_time:0.5522727251052857 [sec]
elapsed_time:0.28664069175720214 [sec]

また、方法2の計算時間を方法1の計算時間で割ったものをグラフ化したのが下図です。
横軸は行列のサイズです。
この結果から方法2のほうが1の約2倍の速度で計算できていることがわかります。
行列サイズが小さい場合(~100)は試行によって、(方法2計算時間/方法１計算時間)は0.6~1の間で変動しますが、行列サイズが大きい場合(500~)は値が安定していました。

結論

numpyで対称行列の作成を行う場合は行列の上三角部分を計算し、次式

A + A^T - diag(A)

を用いて計算するほうが全成分の計算を行うよりも速い。