Pythonの高速化

Question

前提・実現したいこと

Window10のノートPCでAnacondaをインストールして、jupyter notebook上でPython3のコードを書いています。numpyを用いています。
0と1を要素にもつm×n行列Aの列ベクトルをv1,v2,...,vnとします。viに対して、k連続する1をkに置き換えたベクトルをuiとします。u1,u2,...unを列ベクトルとしてもつ行列をBとします。
具体例を示します。
v1=[0,0,1,1,1,0,0,0,1,0,0,1,1,1,1,1,1,1,1,0,1,0]
のとき、
u1=[0,0,3,3,3,0,0,0,1,0,0,8,8,8,8,8,8,8,8,0,1,0]
です。
Aが与えられたとき、Bを得られるようなコードを書きました。しかし例えばm=n=20000だと、計算に40秒ほどかかってしまいます。この計算を数千回繰り返したいと考えておりますので、1回あたりの計算時間をより短くしたいと考えております。自分でもいくつかの方法を試しましたが、計算時間を短縮する有効な方法が分かりません。計算時間を短縮する方法について何かアドバイスを頂けたら幸いです。よろしくお願いいたします。

なお、k連続する1をkに置き換える方法については、
https://maspypy.com/numpy連続同一値の数え上げ-atcoder-abc-129-d
を参考にさせていただきました。

該当のソースコード

コード1

import numpy as np
from tqdm import tqdm

Matrix=np.zeros((2*10**4, 2*10**4), dtype=np.int32)

def f(x):
    temp=np.arange(len(x))
    temp[x>0]=0
    np.maximum.accumulate(temp, out = temp)
    left=temp
    
    x_reversed=x[::-1]
    temp=np.arange(len(x))
    temp[x_reversed>0]=0
    np.maximum.accumulate(temp, out = temp)
    right=len(x)-1-temp[::-1]
    
    y=right
    y -= left + 1
    y[x==0]=0
    return y

Hight=Matrix.shape[0]
Width=Matrix.shape[1]

for i in tqdm(range(Width)):
    Matrix[0,i]=0
    Matrix[-1,i]=0
    Matrix[:,i]=f(Matrix[:,i])

試したこと

①Matrixのデータ型の変更
np.int16にすると計算時間は30秒ほど
np.int8にすると計算時間は22秒ほど
になりましたがデータ型はnp.int32を用いたいと考えています。

②行列に対してまとめて計算
コード1では元の行列の一つ一つの列ベクトルに対して、k連続する1をkに置き換える、という処理をしていますが、行列Aに対して、縦にk連続する1をkに置き換える、という処理をしてみました。

コード2

import numpy as np

Matrix=np.zeros((2*10**4, 2*10**4), dtype=np.int32)

def f(x):
    x[:,0]=np.zeros(len(x), dtype=np.int32)
    x[0]=np.zeros(len(x), dtype=np.int32)
    x[-1]=np.zeros(len(x), dtype=np.int32)
    
    temp=np.arange(len(x)).repeat(len(x)).reshape(len(x),len(x))
    temp[x>0]=0
    np.maximum.accumulate(temp, out = temp)
    left=temp
    
    x_reversed=x[::-1]
    temp=np.arange(len(x)).repeat(len(x)).reshape(len(x),len(x))
    temp[x_reversed>0]=0
    np.maximum.accumulate(temp, out = temp)
    right=len(x)-1-temp[::-1]
    
    y=right
    y -= left + 1
    y[x==0]=0
    return y

%time f(Matrix)

計算時間は45秒ほどとなり、効果はないことが分かりました。

③numbaを用いる
from numba import jit
を加え、
def f(x)
の1行上に
@jit
を付け加えたのですが、エラーメッセージが出たうえで、計算時間はほぼ変わりませんでした。
@jitの代わりに
@jit(nopython=True)
を付け加えたみたところ、エラーが表示され、計算が実行されませんでした。

コード1に対してnumbaを用いる方法をご存じでしたら教えていただけないでしょうか。(→hayataka2049さんの回答により解決しました)

④Cythonを用いる
Cythonというものを用いればPythonの高速化ができる、ということを知りました。そこで次のページを参考にしました。
https://qiita.com/kenmatsu4/items/7c08a85e41741e95b9ba

まず

%load_ext Cython

を実行し、次に

%%cython -n test_cython_code
def fib(int n):
    cdef int i
    cdef double a=0.0, b=1.0

    for i in range(n):
        a, b = a+b, a
    return a

を実行すると、

DistutilsPlatformError: Unable to find vcvarsall.bat

というエラーが出てしまいます。
このエラーを解消する方法をご存じでしたら教えていただけないでしょうか。

⑤CuPyを用いる
CuPyというものを使って、GPUを搭載したコンピューターで計算を行うと高速化できる可能性がある、ということを知りました。(自分のPCにはGPUは搭載されていません)
Google Colaboratoryを使うと、CuPyの性能を確認することができるそうです。以下のようなページを読んで、CuPyについて調べているところです。
https://qiita.com/samacoba/items/d18e6cf09f544477aff4

⑥並列処理を行う
複数のCPUがあるコンピュータでは、複数のCPUを用いることによって計算時間は短くなると思います。そこで並列処理について調べてみて、まずは以下のページを参考にしました。

http://iatlex.com/python/parallel_first

このページ内の以下のコードを実行してみました。

######## 並列計算を使えるように #########
from multiprocessing import Pool
 
##### 並列計算させる関数(処理):引数1つ ###
##### この場合は，引数の二乗を返す関数 ###
def nijou(x):
    print( x*x )
 
###### 並列計算させてみる #########
if __name__ == "__main__":
    p = Pool(4)
    p.map( nijou, range(10) )#nijouに0,1,..のそれぞれを与えて並列演算

しかしjupyter notebookのセルの数字が*になったままになり、結果が出力されません。

コード1に対して複数のCPUを用いた並列処理をする方法をご存じでしたら教えていただけないでしょうか。

補足情報

その他に高速化する方法をご存じでしたら教えていただきたいです。

Comments

[meg_]

・リンクは「リンクの挿入」で記入してください。

[s_miyazaki]

修正しました

[hayataka2049]

解決済みの質問ってそんなに見られないので、トラブルであれば別途立てたほうが良いかもしれません。cythonのvcvarsall.batのエラーは検索すれば色々出てくるのと、環境が絡む質問は回答する側としては再現するのも回答するのも難しいので、できるだけ自力で解決していただいた方が良いとは思いますが……。

[s_miyazaki]

修正しました。
そうなのですね。色々とありがとうございます。

Answer 1

numbaで書いてみました。手元でいくつかのテストケースでは確認しましたが、絶対に正しいとは言い切れないので、ちゃんと動くかはご自身で確認してください。

方針としては、とにかくpython側で処理してしまうと遅いので、配列ごとnumbaに投げます。

numba側のプログラムの書き方は見ての如くです。CとかFORTRANのノリで書いてください。numpyの関数を呼び出すよりそっちの方が速いのです（新たなnumpy配列を返す関数はメモリ上に新たな配列を作るのですべて本質的には遅いのです）。

import numpy as np
from numba import jit

@jit("i4[:,:](i4[:,:])", nopython=True)
def f(A):
    B = A.copy()
    for i in range(B.shape[1]): # 列のループ
        before = 0
        start_pos = 0
        cnt = 0
        for j in range(B.shape[0]): # 行のループ
            if B[j, i] == 1:  # 1のとき
                if before == 0:  # 直前の状態が0なら1にして数え始める
                    start_pos = j
                    before = 1
                cnt += 1  # 数える
            else:  # 0のとき
                if before != 1:  # 直前の状態が0なら無視して続ける
                    continue
                else:  # 直前の状態が1のとき
                    for k in range(start_pos, j):  # 1が出た範囲をcntで埋める
                        B[k, i] = cnt
                    # 状態をリセットする
                    before = 0
                    cnt = 0
        # 行が終わって状態が1のとき
        if before == 1:
            for k in range(start_pos, B.shape[0]):
                B[k, i] = cnt
            
    return B

2万の正方行列で10秒くらいでした。

---

列方向に見ていくのはキャッシュ効率の観点からするとあまり好ましくありません。ということで、同じロジックで列ベクトルではなく行ベクトルを扱うバージョンの関数も作ってみました。

入力を転置して与えてください。結果も転置されたものが出てきます。
（転置がビューかコピーかで変わるかな？　とも思ったのですが、これで速くなったのでたぶんいいのでしょう）

@jit("i4[:,:](i4[:,:])", nopython=True)
def f_trans(A):
    B = A.copy()
    for i in range(B.shape[0]): # 行のループ
        before = 0
        start_pos = 0
        cnt = 0
        for j in range(B.shape[1]): # 列のループ
            if B[i, j] == 1:  # 1のとき
                if before == 0:  # 直前の状態が0なら1にして数え始める
                    start_pos = j
                    before = 1
                cnt += 1  # 数える
            else:  # 0のとき
                if before != 1:  # 直前の状態が0なら無視して続ける
                    continue
                else:  # 直前の状態が1のとき
                    for k in range(start_pos, j):  # 1が出た範囲をcntで埋める
                        B[i, k] = cnt
                    # 状態をリセットする
                    before = 0
                    cnt = 0
        # 行が終わって状態が1のとき
        if before == 1:
            for k in range(start_pos, B.shape[1]):
                B[i, k] = cnt
            
    return B

こちらは6秒でした。

どちらにしても、処理時間のかなりの割合は1.6GBもある配列の生成とコピーで食っています。ということで、これで実用上問題にはならないでしょう。Cythonで同じロジックをやるともう少し速い可能性はあるのですが、numbaのjitコンパイルだって優秀です。

（ベンチマーク取ってるページがありましたが、これを見るとnumbaでいいじゃんとなる・・・
Python を高速化する Numba, Cython 等を使って Julia Micro-Benchmarks してみた - Qiita）

int16やint8にすれば多少速くなりますが、8bitだと128個1が続いただけでオーバーフローなので都合が悪いでしょう。

---

行方向走査で、配列のコピーをやめてインプレース処理にすると関数自体の速度は1秒を切ります（0.6くらい）。使いたいかどうかはわかりませんが、一応載せておきます。

@jit("i4[:,:](i4[:,:])", nopython=True)
def f2_i(B):
    for i in range(B.shape[0]): # 行のループ
        before = 0
        start_pos = 0
        cnt = 0
        for j in range(B.shape[1]): # 列のループ
            if B[i, j] == 1:  # 1のとき
                if before == 0:  # 直前の状態が0なら1にして数え始める
                    start_pos = j
                    before = 1
                cnt += 1  # 数える
            else:  # 0のとき
                if before != 1:  # 直前の状態が0なら無視して続ける
                    continue
                else:  # 直前の状態が1のとき
                    # 1が出た範囲をcntで埋める
                    B[i, start_pos:j] = cnt

                    # 状態をリセットする
                    before = 0
                    cnt = 0
        # 行が終わって状態が1のとき
        if before == 1:
            B[i, start_pos:B.shape[1]] = cnt            
    return B

逆に言うとコピーに時間がかかると思うので、新しい配列を返すつもりであれば本体の処理の高速化で受けられる恩恵はそんなにないのかもしれません。

Comments

[s_miyazaki]

回答ありがとうございます！
説明もコードも大変分かりやすく、助かりました！
入力する行列は、関数を通した後は使わないのでin place処理で大丈夫です。
自分のノートPCでin place処理のコードを実行してみたところ、0.7sほどとなっていて、大幅な高速化ができていました。
これから細部の確認をしていきたいと思います。

[hayataka2049]

速いのは転置して行ベクトルにした場合のものです。列ベクトル版も同様に作れますが、そこそこ遅くなるかと思います。その点はご留意ください。

Answer 2

for文を使わないようにしたらちょっとはやくなりました

Matrix = np.random.randint(0, 2, (2*10**4, 2*10**4), dtype=int)
Matrix[[0, -1], :] = 0

# %%
Matrix1d = Matrix.ravel('F')
repeat = np.ediff1d(np.concatenate(
    ([True], Matrix1d[1:] != Matrix1d[:-1], [True])).nonzero()[0])
base_array = repeat.copy()
base_array[::2] = 0
result_Matrix = base_array.repeat(repeat).reshape(Matrix.shape).T

Comments

[s_miyazaki]

回答ありがとうございます！
①元の行列を1つのベクトルVにする
②Vと、Vの要素を1つずらしたベクトルV1を利用してVにおける0と1の切り替わる場所(I)を検出する
③Iを利用して、0または1が連続する回数を求める(repeat)
④Vの1つ目の要素は0と決まっているのでrepeatを使えば、Vのk連続の1をkに置換したベクトル(U)が得られる
ということですね！
実行時間は40秒ほどとなりましたが、k連続の1をkに置換する別の考え方を知ることができました。ありがとうございました。

Answer 3

• Cythonについては解説記事がたくさんありますので、それらを読めばCの知識がなくても大丈夫かと思います。（複雑な処理をする場合は必要になるかもしれませんが）

• 並列処理については「引数の二乗を返す関数」とありますが、値を返す処理が書かれていません。セルの処理が終わらない原因だと思います。（※windowsのjupyter notebookにおいて）

---

【追記】
windowsにおいてjupyter notebook では並列処理中のprint()が表示されませんね。処理自体はmap()の戻り値を受け取れば可能です。print()したい場合はターミナルで実行するなど他の方法でスクリプトを実行してみてください。

Comments

[s_miyazaki]

回答ありがとうございます！
Cythonの取り入れ方について、調べてみようと思います。