スライスを使わないコードにすると、計算時間がとてもかかります

Question

実現したいこと

「該当のソースコード」は、こちらの途中まで切り出したものです。このコードは正しく動きます。

「該当のソースコード」の31行目以降に度々出てくる[1:-1,1:-1,1:-1] のような書き方にまだ慣れなくて、読みづらいと思い、29～43行目を以下のように直してみました。両者で出力されるファイルの中身（p[HLN]-1.csvとp[HLN]-2.csv）は同じなのですが、p[HLN]-1.csvが出力されるのにはものすごく時間がかかりました。

スライスを使った書き方をした方が圧倒的に計算時間が短くなるので、今後はこのような書き方に慣れた方が良いということでしょうか。もしくは、スライスを使わない書き方で、私のコード（p[HLN]-1.csvを出力するコード）よりも計算時間が早くなるような書き方がありましたら教えてください。スライスに慣れていないので、安全のためスライスを使わない高速コードがあれば嬉しいです。

omega = 1
for I in range(500):
    pd = p.copy()#p[HLN]-2.csvと同じ値になるかを確認するために、あえてコピーする
    
    for i in range(0,nx,1):
        for j in range(0,ny,1):
            for k in range(0,nz,1):
                if k ==0 or k ==nz-1 or j==0 or j==ny-1 or i ==0 or i ==nx-1:
                    p[i][j][k] = 0
                else:
                    p[i][j][k] = (1-omega)*pd[i][j][k]+omega *(pd[i+1][j][k]+pd[i-1][j][k]+pd[i][j+1][k]+pd[i][j-1][k]+pd[i][j][k+1]+pd[i][j][k-1]+delta_L**2*charge[k][j][i]/eps0)/6

filenamecsv = "p[HLN]-1.csv" 
np.savetxt(filenamecsv, p[HLN], delimiter=",")

該当のソースコード

import numpy as np

delta_L=1  
LN=50
HLN=int(LN/2)
nx = LN+1
ny = LN+1
nz = LN+1
xmin = 0
xmax = LN*delta_L
ymin = 0
ymax = LN*delta_L
zmin = 0
zmax = LN*delta_L

p  = np.zeros((nz, ny, nx))
pd = np.zeros((nz, ny, nx))
charge = np.zeros((nz, ny, nx))
x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(ymin, ymax, ny)
z  = np.linspace(zmin, zmax, nz)
print(x)
eps0=1  #とりあえず
Q=1
q=Q/delta_L**3	# 電荷密度に変換

charge[HLN][HLN][HLN]=q

for I in range(500): 
    pd = p.copy()
    p[1:-1,1:-1,1:-1] = (pd[2:,1:-1,1:-1] + pd[:-2,1:-1,1:-1] +
                         pd[1:-1,2:,1:-1] + pd[1:-1,:-2,1:-1] + 
                         pd[1:-1,1:-1,2:] + pd[1:-1,1:-1,:-2] + (delta_L**2)*charge[1:-1, 1:-1, 1:-1]/eps0) /6 

    p[0][:][:]=0
    p[nx-1][:][:]=0
    p[:][0][:]=0
    p[:][ny-1][:]=0
    p[:][:][0]=0
    p[:][:][nz-1]=0

filenamecsv = "p[HLN]-2.csv" 
np.savetxt(filenamecsv, p[HLN], delimiter=",")

試したこと

Googleで「スライスで書くと計算時間が短くなる」で検索すると、以下のように表示されます。

AIによる概要
スライス表記でコードを書くと、計算時間が短縮される場合があります。特に、リストや配列の特定の範囲を扱う際に、スライスを使用することで、元のリストや配列をコピーせずに済むため、処理が高速化されることがあります。ただし、常にスライスが最速とは限りません。データ構造や操作の種類によっては、他の方法がより効率的な場合もあります。

Comments

[melian]

本題とは関係ありませんが、p = np.zeros((nz, ny, nx)) で初期化されていますので以下の処理は不要です。

p[0][:][:]=0
p[nx-1][:][:]=0
p[:][0][:]=0
p[:][ny-1][:]=0
p[:][:][0]=0
p[:][:][nz-1]=0

[meg_]

多重forループの処理は特に遅くなります。単純に処理回数が増えます。
速度が気になるなら配列の使い方を覚える必要があるかと思います。

【追記】もしかして負のインデックスの使い方が分かりにくいという話でしょうか？

[yyicp]

>melian様
コメントありがとうございます。
ご指摘の件、確認しました。

>meg_様
コメントありがとうございます。
ご指摘の通り、負のインデックスが分かりにくく感じるということです。
単純に慣れの問題かもしれません。

[melian]

> 負のインデックスが分かりにくく感じるということです。

負のインデックス -n (n は正の自然数)は arr.size - n と等価なので、以下の様に表記することもできます。

p[1:-1,1:-1,1:-1] => p[1:nz-1,1:ny-1,1:nx-1]
pd[:-2,1:-1,1:-1] => pd[:nz-2,1:ny-1,1:nx-1]

[yyicp]

>melian様
ありがとうございます。

pd[:-2,1:-1,1:-1] の最初の0も省略されているようなので、なかなか慣れませんでした。

[yyicp]

>負のインデックス -n (n は正の自然数)は arr.size - n と等価
負のインデックス -n (n は正の自然数)は len(arr) - n と等価、と理解すれば良いですよね？

[melian]

はい、その通りです。(numpy の ndarray には size という属性(attribute)があって、今回の場合は nx == x.size, ny == y.size, nz == z.size になります)

[yyicp]

どうもありがとうございました。

Answer 1

pythonではfor文を使用すると遅い
は、現時点ではかなり怪しい言説だと思います。
（逆に言うと、10年くらい前はかなり正しくて、内包表記などを駆使してfor文を避けるのが速く動かすテクニックでした。現在は「愚直なforループ」を対象にした最適化が重ねられていて、昔ほどの違いはない気がします。読みやすさを犠牲にしてまで「愚直なforループ」を避ける理由はないと思います）

Pythonにおいては、「Pythonのコードを動かすこと」全般が（他のネイティブ実行コードを動かす言語に比べると）遅いのです。

スライスを使った書き方をしたら速い
もそんなに正しくはないです。
いまスライスを使った操作をしている対象はnumpyのndarrayになっていると思います。
ndarrayへのスライス操作は Python のコード実行ではなくて、numpy が提供するC言語のコードが実行されます。だから速いのです。

---

✗ Python のforループが遅い
◯ Python のコード実行が遅い

✗ スライス操作で書くと速い
◯ numpy（やCPython）のC言語側のコードが実行されるように書くと速い かつ numpy.ndarray ではスライス操作で書くとPuthonのコード実行ではなくてC言語ライブラリの実行になる から速い

が精確かなと思います。

---

質問者さんに必要なことは
自分がPythonの標準ではないサードパーティライブラリであるnumpyを使っているという自覚
numpy では速いコード実行を実現するためにどういう操作を提供しているのかという視点
かと思います。

Comments

[yyicp]

ご回答ありがとうございます。

>ndarrayへのスライス操作は Python のコード実行ではなくて、
>numpy が提供するC言語のコードが実行されます。
が勉強になりました。

Answer 2

ループの処理の速度を上げるにはなるべくループの数を減らすことが重要だと思います。
なので

    for i in range(0,nx,1):
        for j in range(0,ny,1):
            for k in range(0,nz,1):
                if k ==0 or k ==nz-1 or j==0 or j==ny-1 or i ==0 or i ==nx-1:
                    p[i][j][k] = 0
                else:
                    p[i][j][k] = (1-omega)*pd[i][j][k]+omega *(pd[i+1][j][k]+pd[i-1][j][k]+pd[i][j+1][k]+pd[i][j-1][k]+pd[i][j][k+1]+pd[i][j][k-1]+delta_L**2*charge[k][j][i]/eps0)/6

このループの処理を確認します。
例えばfor k in range(0,nz,1):このループですが、k＝0の時とk=nz-1の時は０をセットする形になっています。
ループ前に初期化コードがあるので０をセットする必要はなくループを1～nz-2までと回数を減らすことが可能になります。
その上位のjのループやiのループも同様ですので

    for i in range(1,nx-1,1):
        for j in range(1,ny-1,1):
            for k in range(1,nz-1,1):
                p[i][j][k] = (1-omega)*pd[i][j][k]+omega *(pd[i+1][j][k]+pd[i-1][j][k]+pd[i][j+1][k]+pd[i][j-1][k]+pd[i][j][k+1]+pd[i][j][k-1]+delta_L**2*charge[k][j][i]/eps0)/6

と書き換えることが可能になるかと思います。
ループ回数はnx*ny*nzから(nx-2)*(ny-2)*(nz-2)に減ることとなるのでその分高速化されることとなります。

しかしpythonではfor文を使用すると遅いです。これを速くするとなるとNumba等を利用するなどコンパイルを行って実行すると高速化できます。
ただ全ての機能がコンパイルして実行できるわけではないので注意が必要です。

私の実行環境では下記のようになりました

	実行時間
提示されたコード	123.54秒
私が提示したループ回数を減らしたコード	116.07秒
提示されているスライス表記のコード	0.74秒
提示されたコード（Numba使用）	0.73秒

Comments

[yyicp]

どうもありがとうございました。

時間を計測すると、定量的に分かって、スライスに慣れるべきだと思いました。
「pythonではfor文を使用すると遅いです」が勉強になりました。