💨
結局 Python は遅いのか
結局 Python は遅いのか
はじめに
「Python は遅い」，「Python で for 使うな！」 と言われることは多いですが，実際にどれくらい遅いのか，正直わかっていませんでした．
個人的に，最近は処理速度が必要な場面も増え，「C++ で書くべきか…」と悩むことが多くなってきました．
そこで今回は，Python で素朴な for ループを書いた場合や，他の高速化手法を使った場合，そして C++ で書いた場合とで，どれくらい速度や精度に違いが出るのかを実験し，実際にどれほど差があるのかを確かめてみます.
問題設定
今回は，計算量がそれなりにあって実装が簡単な モンテカルロ積分 を問題として選択しました．
目的
- Python の書き方による速度の差を確認する
- C++ で同様の処理を書いた時の差を確認する
問題
モンテカルロ法を用いて，次の定積分を近似する：
真値（比較用）
アルゴリズム定義
モンテカルロ積分の基本形は次の通りです：
where:
-
x_i \sim U(0, 1) [0,1] -
N
条件
- 乱数は事前生成
- 単一スレッド
-
N 1 \times 10^5 1 \times 10^6 1 \times 10^7 - 各
N - 中央値を代表値とする
比較指標
- 実行時間（ms）
- 相対速度（C++ = 1）
- 推定値の誤差（一回のみ計測）
比較対象
| Lv | 実装 | 説明 |
|---|---|---|
| 1 | Python 素朴for | Python で単純に実装 |
| 2 | Python list 内包表記 | 高速なイメージのある list 内包表記 |
| 3 | Python generator | list を作らない generator |
| 4 | NumPy ベクトル化 | 内部の多くの処理は C++ |
| 5 | Numba JIT（fastmath=False） | Just-In-Time コンパイルで機械語変換 |
| 6 | C++ (-O3) | C++ で単純に実装 |
まず，Python で 1. シンプルな for-loop，2. list comprehension，3. generator，4. Numpy，5. Numba（Just-In-Time コンパイル）を使って比較します．
最後に，C++ でシンプルな for-loop の速度を算出します．
この時使用したコードは Appendix と github に記載します．
実験結果
絶対誤差
まず絶対誤差の結果です．以下に，真値と近似値の差を示します．
| Exp | 実装 | |||
|---|---|---|---|---|
| 1 | Python 素朴for | 3.82e-04 | 2.95e-04 | 5.03e-06 |
| 2 | Python 内包表記 | 3.82e-04 | 2.95e-04 | 5.03e-06 |
| 3 | Python generator | 3.82e-04 | 2.95e-04 | 5.03e-06 |
| 4 | NumPy ベクトル化 | 3.82e-04 | 2.95e-04 | 5.03e-06 |
| 5 | Numba JIT | 3.82e-04 | 2.95e-04 | 5.03e-06 |
| 6 | C++ (-O3) | 8.48e-04 | 5.57e-05 | 7.41e-05 |
当たり前ですが，
Python の各実装では，同一の乱数列を使用しているため，数値誤差を除き同一の結果となっています．
一方，C++ の場合は Python よりも誤差が大きいケースも見られました．これは，使用している乱数生成器が C++ と Python で異なることに加え，浮動小数点演算の順序や sin 関数のライブラリ実装差によるものであり，アルゴリズム上の問題ではありません．（Python はどの実装も seed 固定で同じ）
ここで，モンテカルロ法による推定値の標準誤差（
where:
上記より，
-
N=1 \times 10^5 SE \approx \sqrt{(0.0614/1\times 10^5)} \approx 7.8 \times 10^{-4} -
N=1 \times 10^6 SE \approx \sqrt{(0.0614/1\times 10^6)} \approx 2.5 \times 10^{-4} -
N=1 \times 10^7 SE \approx \sqrt{(0.0614/1\times 10^7)} \approx 7.8 \times 10^{-5}
C++ の誤差（8.48e-04, 5.57e-05, 7.41e-05）はいずれもこの標準誤差のオーダーに収まっており，モンテカルロ法の統計的な誤差と考えられます．
計算速度
以下に，200 回実行した時の処理時間の中央値（単位：ms）を示します．
| Exp | 実装 | |||
|---|---|---|---|---|
| 1 | Python 素朴for | 14.425 | 144.264 | 1432.188 |
| 2 | Python 内包表記 | 14.895 | 148.467 | 1640.013 |
| 3 | Python generator | 15.411 | 155.766 | 1548.547 |
| 4 | NumPy ベクトル化 | 0.553 | 6.894 | 69.083 |
| 5 | Numba JIT | 0.604 | 6.280 | 65.199 |
| 6 | C++ (-O3) | 0.589295 | 5.89161 | 59.7527 |
同じ処理をしていても実装方法によって大きな差が見られます．
特に，Exp 1，Exp 2，Exp 3 は見てられません...．
今回のタスクでは圧倒的に，Numpy・Numba が高速です．
まず
また，Exp 4，Exp 5はわずかに増加しています．
Exp 2 で遅くなる理由（と思われること）
list 内包表記では，以下の処理が発生します：
-
sin(x)の計算 - Python float オブジェクトの生成
- list への append
-
sum()による再走査
つまり，
- 2
N -
N
が発生します．
Python の list はポインタ配列なので，各要素は
[ptr][ptr][ptr]...
↓
Python float object
となっています．
そのため，
- キャッシュミスの増加
- メモリ書き込みの増加
が発生し，CPU 計算よりもメモリアクセスがボトルネックになります．
for-loop （Exp 1）の場合
- list 生成しない
-
N
などから，追加のオブジェクト生成やメモリ書き込みが発生しないため，余計なオブジェクト生成がなく，最も効率的に実行されます．
NumPy・Numba（Exp 4, 5）がわずかに線形からズレる理由
Exp 4，Exp 5 が若干線形からズレるのは，関数呼び出しやメモリ帯域，ライブラリコストの影響だと考えられます．
相対速度（c++=1.00）
以下に，C++ での速度を 1.00 としたときの Python での速度を示します．
| Exp | 実装 | |||
|---|---|---|---|---|
| 1 | Python 素朴for | 24.478402 | 24.486346 | 23.968591 |
| 2 | Python 内包表記 | 25.275965 | 25.199733 | 27.446676 |
| 3 | Python generator | 26.151588 | 26.438614 | 25.915934 |
| 4 | NumPy ベクトル化 | 0.938409 | 1.170139 | 1.156149 |
| 5 | Numba JIT | 1.024954 | 1.065923 | 1.091147 |
| 6 | C++ (-O3) | 1.00 | 1.00 | 1.00 |
C++ の処理速度と比較すると，Python for-loop・list 内包表記・generator 表記は 23 倍〜27 倍ほどの差が出ました．
Numpy・Numba はほぼ C++ と同等で肉薄しています．
※NOTE: なお，NumPy や Numba が高速に見えるのは，Python のループが高速になっているわけではなく，内部でネイティブコードが実行されているためです．
NumPy のベクトル演算は C で実装された連続メモリアクセスのループとして実行され，Numba は Python コードを LLVM を用いて機械語に JIT コンパイルします．
そのため，実行時には Python のインタプリタを介さず，ネイティブコードとして動作しており，純粋な Python 実行とは本質的に異なります．
結果
- Python の素朴な for-loop やリスト内包表記は，C++ と比べて大幅に遅い（約25倍前後）
- NumPy や Numba を活用すれば，C++ とほぼ同じレベルの高速化が可能
- Python で速度が重要な場合は，NumPy や Numba の利用が非常に有効
今回の結果を見ると，（想定通り）Python をそのまま書くとかなり遅いですが，実装次第では C++ と比較してそこまで遅くならないという結果でした．
さいごに
今回のタスクでは
Python で実装する場合，NumPy や Numba を用いると、ネイティブコードとして実行されるため，単純な C++ 実装と同等レベルの処理速度が得られました．
今回はシンプルなタスクだったので NumPy で実装ができましたが，タスクによっては NumPy が使えないこともあると思うので，Numba の方が良い選択肢に見えます．
一方，Numba は非常に強力ですが、
- インストールがやや面倒
- Python の最新バージョンに追従が遅い
- サポートされない構文がある
といった制約もあります。
なので，結局 Python・C++・他の言語...うまく使い分けるのが良いですね...．
というわけで，今回は Python での実装による処理速度差と C++ での処理速度差を算出してみました．
Python はそのまま書くと遅くなりがちですが，NumPy や Numba を適切に活用することで，ネイティブ言語に近い性能を引き出せることが分かりました．
Python は遅いのではなく，「Python のループ」が遅い
Appendix
config.py
import time
import numpy as np
Ns = [100_000, 1_000_000, 10_000_000]
trials = 200
REAL = 0.45969769413186023
def measure(xs, func):
start = time.perf_counter()
result = func(xs)
end = time.perf_counter()
return (end - start) * 1000 # ms
def bootstrap_ci(data, n_boot=1000, ci=95):
medians = []
n = len(data)
for _ in range(n_boot):
sample = np.random.choice(data, n, replace=True)
medians.append(np.median(sample))
low = np.percentile(medians, (100 - ci) / 2)
high = np.percentile(medians, 100 - (100 - ci) / 2)
return low, high
01_for_loop.py
import numpy as np
import math
from config import Ns, trials, REAL, measure, bootstrap_ci
def compute(xs):
s = 0.0
for x in xs:
s += math.sin(x)
return s / len(xs)
for N in Ns:
print(f"N = {N}")
np.random.seed(42)
xs = np.random.rand(N)
# ---- Warmup ----
compute(xs)
warm = compute(xs)
print(f"Warmup Result: {warm}, Abs Error: {abs(warm - REAL)}")
times = []
for _ in range(trials):
times.append(measure(xs, compute))
times = np.array(times)
median = np.median(times)
min_v = np.min(times)
p90 = np.percentile(times, 90)
p95 = np.percentile(times, 95)
low, high = bootstrap_ci(times)
print(f"Median (ms): {median:.3f}")
print(f"Min (ms): {min_v:.3f}")
print(f"P90 (ms): {p90:.3f}")
print(f"P95 (ms): {p95:.3f}\n")
print(f"Bootstrap CI (95%) (ms): [{low:.3f}, {high:.3f}]\n")
02_list_comprehension.py
import numpy as np
import math
from config import Ns, trials, REAL, measure, bootstrap_ci
def compute_pythonic(xs):
return sum([math.sin(x) for x in xs]) / len(xs)
for N in Ns:
print(f"N = {N}")
np.random.seed(42)
xs = np.random.rand(N)
compute_pythonic(xs)
warm = compute_pythonic(xs)
print(f"Warmup Result: {warm}, Abs Error: {abs(warm - REAL)}")
times = []
for _ in range(trials):
times.append(measure(xs, compute_pythonic))
times = np.array(times)
median = np.median(times)
min_v = np.min(times)
p90 = np.percentile(times, 90)
p95 = np.percentile(times, 95)
low, high = bootstrap_ci(times)
print(f"Median (ms): {median:.3f}")
print(f"Min (ms): {min_v:.3f}")
print(f"P90 (ms): {p90:.3f}")
print(f"P95 (ms): {p95:.3f}\n")
print(f"Bootstrap CI (95%) (ms): [{low:.3f}, {high:.3f}]\n")
03_generator.py
import numpy as np
import math
from config import Ns, trials, REAL, measure, bootstrap_ci
def compute_pythonic(xs):
return sum(math.sin(x) for x in xs) / len(xs)
for N in Ns:
print(f"N = {N}")
np.random.seed(42)
xs = np.random.rand(N)
compute_pythonic(xs)
warm = compute_pythonic(xs)
print(f"Warmup Result: {warm}, Abs Error: {abs(warm - REAL)}")
times = []
for _ in range(trials):
times.append(measure(xs, compute_pythonic))
times = np.array(times)
median = np.median(times)
min_v = np.min(times)
p90 = np.percentile(times, 90)
p95 = np.percentile(times, 95)
low, high = bootstrap_ci(times)
print(f"Median (ms): {median:.3f}")
print(f"Min (ms): {min_v:.3f}")
print(f"P90 (ms): {p90:.3f}")
print(f"P95 (ms): {p95:.3f}\n")
print(f"Bootstrap CI (95%) (ms): [{low:.3f}, {high:.3f}]\n")
04_numpy.py
import numpy as np
from config import Ns, trials, REAL, measure, bootstrap_ci
def compute_numpy(xs):
return np.sin(xs).mean()
for N in Ns:
print(f"N = {N}")
np.random.seed(42)
xs = np.random.rand(N)
compute_numpy(xs)
warm = compute_numpy(xs)
print(f"Warmup Result: {warm}, Abs Error: {abs(warm - REAL)}")
times = []
for _ in range(trials):
times.append(measure(xs, compute_numpy))
times = np.array(times)
median = np.median(times)
min_v = np.min(times)
p90 = np.percentile(times, 90)
p95 = np.percentile(times, 95)
low, high = bootstrap_ci(times)
print(f"Median (ms): {median:.3f}")
print(f"Min (ms): {min_v:.3f}")
print(f"P90 (ms): {p90:.3f}")
print(f"P95 (ms): {p95:.3f}\n")
print(f"Bootstrap CI (95%) (ms): [{low:.3f}, {high:.3f}]\n")
05_numba.py
import numpy as np
import math
from numba import njit
from config import Ns, trials, REAL, measure, bootstrap_ci
@njit
def compute_numba(xs):
s = 0.0
for x in xs:
s += math.sin(x)
return s / len(xs)
for N in Ns:
print(f"N = {N}")
np.random.seed(42)
xs = np.random.rand(N)
# ---- JIT Warmup (IMPORTANT) ----
compute_numba(xs) # Compile and warmup
warm = compute_numba(xs)
print(f"Warmup Result: {warm}, Abs Error: {abs(warm - REAL)}")
times = []
for _ in range(trials):
times.append(measure(xs, compute_numba))
times = np.array(times)
median = np.median(times)
min_v = np.min(times)
p90 = np.percentile(times, 90)
p95 = np.percentile(times, 95)
low, high = bootstrap_ci(times)
print(f"Median (ms): {median:.3f}")
print(f"Min (ms): {min_v:.3f}")
print(f"P90 (ms): {p90:.3f}")
print(f"P95 (ms): {p95:.3f}\n")
print(f"Bootstrap CI (95%) (ms): [{low:.3f}, {high:.3f}]\n")
calc_mc.cpp
#include <vector>
#include <cmath>
#include <chrono>
#include <iostream>
#include <numeric>
#include <random>
#include <algorithm> // For std::sort
const double real = 0.45969769413186023;
double compute(const std::vector<double> &xs)
{
double s = 0.0;
for (size_t i = 0; i < xs.size(); ++i)
{
s += std::sin(xs[i]);
}
return s / xs.size();
}
double measure(const std::vector<double> &xs)
{
auto start = std::chrono::high_resolution_clock::now();
volatile double result = compute(xs); // Avoid optimizing away the compute call
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> diff = end - start;
return diff.count(); // ms
}
std::pair<double, double> bootstrap_ci(
const std::vector<double> ×,
int n_boot = 1000)
{
std::mt19937 rng(123);
std::uniform_int_distribution<int> dist(0, times.size() - 1);
std::vector<double> medians;
medians.reserve(n_boot);
std::vector<double> sample(times.size());
for (int b = 0; b < n_boot; ++b)
{
// Resampling with replacement
for (size_t i = 0; i < times.size(); ++i)
sample[i] = times[dist(rng)];
// Compute median of the sample
std::sort(sample.begin(), sample.end());
medians.push_back(sample[sample.size() / 2]);
}
std::sort(medians.begin(), medians.end());
double lo = medians[n_boot * 0.025];
double hi = medians[n_boot * 0.975];
return {lo, hi};
}
int main()
{
std::vector<size_t> Ns = {100000, 1000000, 10000000};
const int trials = 200; // Same number of trials with Python for consistency
for (size_t N : Ns)
{
std::mt19937 gen(42);
std::uniform_real_distribution<double> dist(0.0, 1.0);
std::cout << "N = " << N << "\n";
// ---- Generate random data ----
std::vector<double> xs(N);
for (size_t i = 0; i < N; ++i)
xs[i] = dist(gen);
// ---- Warmup ----
compute(xs);
double warm = compute(xs);
std::cout << "Warmup Result: " << warm
<< ", Abs Error: " << std::abs(warm - real) << "\n";
// ---- Measurement ----
std::vector<double> times;
times.reserve(trials);
for (int t = 0; t < trials; ++t)
{
times.push_back(measure(xs));
}
// ---- Statistics (median-based) ----
std::sort(times.begin(), times.end());
double median;
if (trials % 2 == 0)
median = (times[trials / 2 - 1] + times[trials / 2]) / 2.0;
else
median = times[trials / 2];
double min_v = times.front();
// Percentiles (simple implementation)
auto idx90 = static_cast<size_t>(trials * 0.90);
auto idx95 = static_cast<size_t>(trials * 0.95);
double p90 = times[idx90];
double p95 = times[idx95];
auto [ci_lo, ci_hi] = bootstrap_ci(times);
std::cout << "Median (ms): " << median << "\n";
std::cout << "Min (ms): " << min_v << "\n";
std::cout << "P90 (ms): " << p90 << "\n";
std::cout << "P95 (ms): " << p95 << "\n";
std::cout << "95% CI Median (ms): ["
<< ci_lo << ", " << ci_hi << "]\n\n";
}
return 0;
}
Discussion