【Python3】 numpy.ndarrayの疎視化について速度比較など

要旨

numpy.ndarrayを疎視化する方法を何種類か実験してみました。
ここでいう疎視化とはopencvのresizeのような処理を指します。
ただし以下に記すような疎視化について扱うため、配列の拡大や疎視化した一辺が元の配列の一辺を割り切れない場合などは取り扱っていません。

前書き

格子ガスオートマトンを実装していたところ、結果を描画する段において疎視化をする必要が出ました。
具体的には1000010000の配列を1010ごとに平均して1000*1000の配列を作る処理をしたいなと。
opencvにresizeがあるし、使いたい人も多いはずだから（主観）無いことは無いだろうとググるとなんか違うresizeが出てきました。
じゃあopencvのresizeを使えば良い話ではあるのですが、微妙にやっかいなことに符号なし1byte整数にしか対応しておらずそれ以外の型の配列を入れるとエラーを吐きます。結果の描画だけなのでこれを使ってもまぁそこまで悪いわけではないのですが折角だから精確な値を用いたいので実装することにしました。

実行環境

• Windows 7 Home Premium 64bit
• Intel Core i5-2430M 2.4GHz
• Python 3.6.1 (Anaconda)

試した方法

以下の四種類を試しました。

• 疎視化後の配列を用意して、元の配列をスライスで取り出して平均した値をfor文を用いて代入していく
• 上の手順をリスト内包表記で記述
• 元の配列をスライスして足し合わせる操作を行と列で行う
• 上の手順をnp.rollを用いて記述

なお以下のソースでは元の配列をarray（サイズ10000*10000）、疎視化する一辺の長さをmean_size(以下の説明では10）として説明していきます。

疎視化後の配列を用意して、元の配列をスライスで取り出して平均した値をfor文を用いて代入していく

以下がソース

mean_loop

def mean_loop(array, mean_size):
    meaned_size = np.shape(array)[0] // mean_size
    return_array = np.empty((meaned_size, meaned_size))

    for y in range(meaned_size):
        for x in range(meaned_size):
            return_array[y, x] = array[mean_size*y:mean_size*(y+1), mean_size*x:mean_size*(x+1)].mean()

    return return_array

やりたいことをそのまま素直に実装したものになります。
Cで実装するならこうするのではないでしょうか。
ただPythonで二重ループという時点で嫌な雰囲気がします。
しかも今回は1000010000の配列を10001000にするので実行時間はひどいことになるでしょう。

上の手順をリスト内包表記で記述

mean_list_comprehension

def mean_list_comprehenion(array, mean_size):
    meaned_size = np.shape(array)[0] // mean_size
    return np.array([[array[mean_size*y:mean_size*(y+1), mean_size*x:mean_size*(x+1)].mean() for x in range(meaned_size)] for y in range(meaned_size)])

Pythonと言えばリスト内包表記、リスト内包表記といえばPython。
速くなるとも言われている（主観）ので期待が持てます。
ただしこちらもループは1000*1000です。
しかし書いといてなんですが二重のリスト内包表記はクソ長いうえに読み辛いですね。

元の配列をスライスして足し合わせる操作を行と列で行う

mean_slice

def mean_slice(array, mean_size):
    sum_row_array = np.array(array[::mean_size])

    for row in range(1, mean_size):
        sum_row_array += array[row::mean_size]

    sum_column_array = sum_row_array[:,::mean_size]

    for column in range(1, mean_size):
        sum_column_array += sum_row_array[:,column::mean_size]

    return sum_column_array / mean_size**2

少々泥臭めの方法。
10行足し合わせて10列足し合わせるため何とループはたったの20回。
どう考えても上の二つよりは速いはずです。

上の手順をnp.rollを用いて記述

mean_np_roll

def mean_np_roll(array, mean_size):
    sum_row_array = np.array(array)
    for row in range(1, mean_size):
        sum_row_array += np.roll(array, -row, axis=0)

    sum_column_array = np.array(sum_row_array)
    for column in range(1, mean_size):
        sum_column_array += np.roll(sum_row_array, -column, axis=1)

    return sum_column_array[::mean_size, ::mean_size] / mean_size ** 2

やっていることが若干異なりますが、思想としては同じもの。
ndarrayは要素のアクセスが遅いと聞いたことがあるのでもしかしたらsliceは遅いのではと思ったので作ってみました。
とはいえsliceの方法では足し合わせる領域が小さくなるのに対し、こちらでは元の配列と同サイズのままなのでこの点で遅くなりそうではあります。

結果

いろんな実験設定で見てみたかったので、元の配列の一辺と平均するサイズを(3000, 10), (6000, 10), (9000, 10), (3000, 100), (6000, 100), (3000, 1000)に変化させて計算してみました。
以下が実行結果です。

sliceが上の二つより速いのは想像通りでしたが、np.rollを用いた手法が予想以上に遅かったです。
特に平均する領域が大きいときの遅さがなかなか凄いことになっています。
それとリスト内包表記が別にそんなに速いわけでもない結果になりました。
こうした結果を考えると可読性の点からも無理して内包表記にするのはマズイのではと思います。
とりあえずsliceを用いた上記の関数がよさそうなのでこれを用いるのが良いのではないかと考えています。

（おまけ？）Cとの比較

一番上の手法をCで実装してみました。拡張子はcppですが内容は普通のCです。
なんだかんだで速いと聞くので生ポを用いています。
一次元配列を二次に見立てている部分など細部が異なりますが、Cでの実装を行うとなるとこんなコードになるはずですのでそこまでの公平性は考えませんでした。
Visual Studio 2017 Communityにてx64 Releaseビルドを行いました。

mean_loop.cpp

#include "stdafx.h"
#include <stdio.h>
#include <time.h>
#include <Windows.h>

double* mean_loop(int mean_size, int array_size, double* array) {
	double *return_array;
	int meaned_size = array_size / mean_size;

	return_array = new double[meaned_size * meaned_size];


	for (size_t i = 0; i < meaned_size; i++)
	{
		for (size_t j = 0; j < meaned_size; j++)
		{
			double temp = 0;
			for (size_t k = 0; k < mean_size; k++)
			{
				for (size_t l = 0; l < mean_size; l++)
				{
					temp += array[(i + k) * meaned_size * mean_size + j*mean_size + l];
				}
			}
			return_array[i * meaned_size+ j] = temp / (mean_size * mean_size);
		}
	}

	return return_array;
};

int main()
{
	double *array, *return_array;
	int mean_size = 10, array_size = 3000;
	clock_t start, end, span;
	int array_sizes[] = {3000, 6000, 9000};
	int mean_sizes[] = { 10, 100, 1000 };

	for (size_t as = 0; as < 3; as++)
	{
		for (size_t ms = 0; ms < 3; ms++)
		{
			array_size = array_sizes[as];
			mean_size = mean_sizes[ms];

			array = new double[array_size * array_size];

			for (size_t i = 0; i < array_size * array_size; i++)
			{
				array[i] = i;
			}

			start = clock();

			return_array = mean_loop(mean_size, array_size, array);

			end = clock();

			span = end - start;
			printf("%d %4d took: %lf\n", array_size, mean_size, (double)span / CLOCKS_PER_SEC);

			delete return_array, array;

		}
		
	}
	
	system("pause");

	return 0;
}

この計算の結果です。

やっぱり速いですね！
生ポなのでメモリリークが怖いですがそこは努力ですし（極論）、数値計算では大規模な計算を行うことが多いのでやはりC/C++を用いるべきではないでしょうか。
とはいえそうした恐怖と引き換えにこの速度であるのと、Pythonの上記の比較的安全でまぁまぁ速いコードを比較しますと適材適所という気もします。

まとめ

数値計算にスクリプト言語を使うのはよしましょう。