# 並列化と GPU：Amdahl・Gustafson の法則、仕事量と依存鎖、roofline モデル

> 並列化で得られる速度向上の上限を Amdahl・Gustafson の法則と Brent の評価で、メモリ帯域による上限を roofline モデルで証明し、部分和の分割による丸めの差の上界と、並列の乱数ストリームの分け方を実行例つきで扱う。
> https://rikai.mugen-giken.com/computer-science/research-computing/parallel-and-gpu-computing

## 0. この記事の要点

計算を速くしたいとき、「コアの数を倍にすれば、時間は半分になるのか」というのは素朴で大事な疑問です。実際には、順番にしかできない仕事が残ることや、データの出し入れで待たされることのせいで、倍にしても倍速くはなりません。この記事は、そうなる理由と限界を、実験の前に紙の上で見積もる道具（Amdahl の法則など）と、並列にしたときに和の末尾の桁や乱数がどうなるかを扱います。読み終えると、並列化や GPU に手を出す前に、どれくらい速くなりそうかを見積もれるようになります。

以下の要点は、記事全体の結果を先にまとめたものです。分からない言葉があっても、各節で定義から順に説明するので、そのまま読み進めて構いません。

- 逐次にしか実行できない部分の割合が $f$ の計算は、何台使っても $1/f$ 倍より速くなりません（**Amdahl の法則**）。問題の規模を台数に比例させて大きくするなら、速度向上は $g + (1-g)p$ と台数にほぼ比例します（**Gustafson の法則**）。2 つは同じ量を別の基準で測ったもので、変換式 $f = g/(g + (1-g)p)$ で移り合います。
- 測定した速度向上から逆算した逐次割合 $e(p)$（**Karp–Flatt 指標**）が台数とともに増えるなら、並列化のオーバーヘッドが増えています。
- 仕事量 $W$、最長の依存鎖 $L$ の計算は、$p$ 台で $\max(W/p, L)$ ステップ以上かかり、貪欲なスケジュール（実行できる作業があれば、空いている台数の分だけ必ず実行する方針）で $W/p + L$ ステップ以下で終わります（**Brent の評価**）。
- 演算強度 $I$（1 バイトあたりの浮動小数演算数）の計算の性能は $\min(\pi, \beta I)$ を超えません（**roofline モデル**）。GPU などの加速器に仕事を移して得をするのは、節約できる計算時間が転送の時間を上回るときだけです。
- 配列を $p$ 個に分けて足した結果は、逐次和と $2\gamma_{n-1}\sum|x_i|$ 以内でしか一致しません。分割をスレッド数によらず固定すれば、結果はビット単位で再現します。
- 並列のワーカーに同じ乱数種を渡すと標本が重複し、平均の分散は $1/p$ 倍になりません。`SeedSequence.spawn`（1 つの種から、重ならない子の種を作る NumPy の関数）で子を作って渡します。

<Aside type="tip" title="この記事で使う事前知識">

- **単位丸め**: 浮動小数の 1 回の演算で、結果が真の値から相対的にずれる幅の上限のこと（記号は $u$）。$\gamma_{n-1}$ は $n-1$ 回の演算が重なったときの上限を表す量。→ [浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)
- **和の木**: 数の和を、葉が数で内部の頂点が 1 回の加算である二分木として表したもの。括弧の付け方がそのまま木の形になる。→ [和の木](/computer-science/research-computing/reproducible-experiments#def-summation-tree)
- **擬似乱数生成器**: 種と呼ぶ初期値から、決まった乱数の列を作る仕組み。同じ種なら同じ列が出る。→ [擬似乱数生成器](/computer-science/research-computing/reproducible-experiments#def-prng)
- **有向グラフと有向道**: 辺に向きのある図と、辺の向きに沿ってたどった頂点の並び。閉路がないとは、たどって元の頂点に戻れないこと。
- **キャッシュと主記憶**: プロセッサの近くにある小さくて速い記憶と、大きいが遅い本体のメモリ。
- **分散と共分散**: 分散は平均からのばらつきの大きさ、共分散は 2 つの量が一緒に動く度合い。無相関とは共分散が $0$ のこと。

</Aside>

## 1. 動機

数値実験が 1 日で終わらなくなると、多くの人はまずコア数の多い計算機や GPU を探します。ところが、16 コアの計算機に移しても 3 倍しか速くならなかった、GPU に移したらかえって遅くなった、スレッド数を変えたら結果の末尾の桁が変わった、という経験もまたよくあります。これらは実装の巧拙以前に、計算の構造から決まる上限と、浮動小数の加法の性質から説明できる現象です。

この章の目的は、並列化の**前に**紙の上で見積もりをするための道具をそろえることです。逐次部分の割合が分かれば Amdahl の法則が速度向上の上限を与え、依存関係の深さが分かれば Brent の評価が必要なステップ数を上下から挟み、演算数とメモリの読み書きの量が分かれば roofline モデルが「コアを増やしても速くならない」計算を見分けます。どれも証明は短く、仮定がはっきりしています。仮定がはっきりしていれば、測定がモデルから外れたときに、どの仮定が破れたのかを問うことができます。

並列化はまた、[再現できる数値実験](/computer-science/research-computing/reproducible-experiments)の章で扱った決定性を壊しやすい操作でもあります。和の順序が変わることと、乱数の消費の順序が変わることの 2 点について、何が保証され何が保証されないかを命題として述べます。

この章では時間の計測値を載せません。計測値は計算機ごとに違い、読者の手元で再現しないからです。速度向上はモデルの式と演算数で論じ、コードの出力は、命題が予言する挙動を見るための例にとどめます。

## 2. 速度向上・効率・スケーリング

<Definition id="def-speedup" title="速度向上と効率">
ある計算を $p$ 台（$p \ge 1$ は整数。コア・プロセス・計算ノードのいずれでもよい）で実行したときの所要時間を $T(p) \gt 0$ とする。
$$
S(p) = \frac{T(1)}{T(p)}, \qquad E(p) = \frac{S(p)}{p}
$$
をそれぞれ**速度向上**（speedup）、**効率**（efficiency）という。$T(1)$ のうち、台数を増やしても短くならない部分の時間 $T_{\mathrm{s}}$ の割合 $f = T_{\mathrm{s}}/T(1)$ を**逐次部分の割合**という。
</Definition>

<Definition id="def-scaling" title="強スケーリングと弱スケーリング">
問題の規模を固定して台数 $p$ を増やし、$T(p)$ がどう減るかを調べることを**強スケーリング**という。台数に比例して問題の規模を大きくし（1 台あたりの仕事を一定に保ち）、所要時間がどれだけ一定に保たれるかを調べることを**弱スケーリング**という。
</Definition>

並列化の形にも名前を付けておきます。同じ操作を配列の多数の要素に同時に適用する形を**データ並列**、互いに異なる仕事を同時に進める形を**タスク並列**といいます。仕事が互いに独立で、結果を最後に集めるだけで済む計算（パラメータの掃引、乱数種を変えた試行、画像ごとの処理）は「恥ずかしいほど並列」（embarrassingly parallel）と呼ばれ、逐次部分は入力の配布と結果の集計だけです。

実行の単位はプロセス（<Ref to="computer-science/cs-basics/operating-systems#def-process" />）とスレッドです。プロセスはアドレス空間を別々に持ち、スレッドは 1 つのプロセスの中でアドレス空間を共有します。標準の CPython には、Python のバイトコードを同時に 1 スレッドしか実行させないロック（GIL: global interpreter lock）があるので、純粋な Python のループをスレッドで分けても速くなりません。NumPy の配列演算の多くは計算の間 GIL を手放し、内部で呼ぶ BLAS（行列やベクトルの基本演算をまとめた、高速な標準ライブラリ）は自前のスレッドを使うので、配列演算はスレッドでも並列に走ります。Python のループそのものを並列にしたいときは、プロセスを分ける（`concurrent.futures.ProcessPoolExecutor`）か、ループを機械語に翻訳して GIL の外で実行する（Numba などの実行時翻訳、JIT）ことになります。

## 3. Amdahl の法則と Gustafson の法則

### 3.1. Amdahl の法則

[AI 時代のエンジニアの生存戦略](/computer-science/ai-era/engineer-survival-strategy)の章の一般化アムダール則（<Ref to="computer-science/ai-era/engineer-survival-strategy#prop-amdahl" />）は、作業の一部分だけが速くなるときの全体の速度向上を与えました。ここでは部分の数を増やした形で述べ直し、並列計算の法則をその系として導きます。

<Theorem id="thm-amdahl" title="Amdahl の法則">
1 台での所要時間 $T(1)$ が $m$ 個の部分の時間の和 $T(1) = T_1 + \cdots + T_m$（$T_i \ge 0$）に分かれ、$p$ 台では第 $i$ の部分が $s_i \gt 0$ 倍速くなって $T(p) = \sum_{i=1}^m T_i/s_i$ になるとする。$f_i = T_i/T(1)$ とおくと
$$
S(p) = \frac{1}{\displaystyle\sum_{i=1}^{m}\frac{f_i}{s_i}}
$$
である（$s_i$ の、重み $f_i$ による調和平均）。とくに、逐次部分（割合 $f$、$s_1 = 1$）と並列化できる部分（割合 $1 - f$、$1 \le s_2 \le p$）の 2 つに分かれるとき
$$
S(p) \le \frac{1}{f + \dfrac{1-f}{p}} \le \frac{1}{f}
$$
が成り立つ（$f = 0$ のとき右端は $+\infty$ と読む）。左の等号は $s_2 = p$ のときに成り立つ。
</Theorem>

<Proof of="thm-amdahl">
$T_i = f_iT(1)$ なので $T(p) = \sum_i T_i/s_i = T(1)\sum_i f_i/s_i$ で、$S(p) = T(1)/T(p) = 1/\sum_i(f_i/s_i)$ です。

2 つの部分の場合は $S(p) = 1/\bigl(f + (1-f)/s_2\bigr)$ で、これは <Ref to="computer-science/ai-era/engineer-survival-strategy#prop-amdahl" /> で高速化される割合を $1 - f$、速度倍率を $s_2 \ge 1$ とした式そのものです。$s_2 \le p$ と $1 - f \ge 0$ から $(1-f)/s_2 \ge (1-f)/p$ なので、分母は $f + (1-f)/p$ 以上で、$S(p) \le 1/\bigl(f + (1-f)/p\bigr)$ です。等号は $s_2 = p$（または $f = 1$）のときです。さらに $(1-f)/p \ge 0$ から $f + (1-f)/p \ge f$ で、$f \gt 0$ なら $1/\bigl(f + (1-f)/p\bigr) \le 1/f$ です。
</Proof>

仮定 $s_2 \le p$ は「$p$ 台に分けた仕事は、1 台のときの $p$ 倍より速くはならない」というモデルの仮定です。実際の計算機では、分割によって 1 台あたりのデータがキャッシュに収まるようになり、$p$ 倍を超える速度向上（超線形の速度向上）が測定されることがあります。これは定理の反例ではなく、仮定が破れた場合です。

定理の読み方は単純です。$f = 0.05$ なら上限は $20$ 倍で、$p = 16$ でも $S \le 1/(0.05 + 0.95/16) \approx 9.1$ にとどまります。効率は $E(16) \le 0.57$ です。並列化できる部分をどれだけ速くしても、残りの $5\%$ が全体を支配するようになります。[クラウドの章](/computer-science/software-engineering/cloud-computing)の拡張則（<Ref to="computer-science/software-engineering/cloud-computing#prop-usl" />）は、これに台数間の調整コストの項を足して、台数を増やすとかえって遅くなる領域があることを示したものでした。命令のパイプライン化で段数倍に達しない（<Ref to="computer-science/cs-basics/computer-architecture#prop-pipeline-speedup" />）のも、段ごとの固定の遅延が「逐次部分」として残るという同じ構造です。

### 3.2. Gustafson の法則

Amdahl の法則は問題の規模を固定しています。しかし台数が増えたら、同じ問題を速く解くのではなく、同じ時間でより大きい問題（細かい格子、多い標本）を解く、という使い方のほうが普通です。このとき基準にするのは 1 台での時間ではなく、並列機での時間です。

<Theorem id="thm-gustafson" title="Gustafson の法則">
$p$ 台の並列機である計算を実行した所要時間が、逐次部分の時間 $t_{\mathrm{s}} \ge 0$ と並列部分の時間 $t_{\mathrm{p}} \ge 0$ の和 $T(p) = t_{\mathrm{s}} + t_{\mathrm{p}} \gt 0$ で、並列部分では $p$ 台がそれぞれ時間 $t_{\mathrm{p}}$ ずつ仕事をしているとする。$g = t_{\mathrm{s}}/T(p)$（並列機での実行時間に占める逐次部分の割合）とおく。同じ計算を 1 台で実行すると、並列部分の仕事を順に行うので $T(1) = t_{\mathrm{s}} + p\,t_{\mathrm{p}}$ かかり、
$$
S(p) = \frac{T(1)}{T(p)} = g + (1 - g)\,p
$$
である（**拡大速度向上**）。この計算の <Ref to="def-speedup" /> の意味の逐次部分の割合 $f = t_{\mathrm{s}}/T(1)$ は
$$
f = \frac{g}{g + (1-g)p}, \qquad g = \frac{fp}{1 - f + fp}
$$
を満たし、この $f$ を <Ref to="thm-amdahl" /> の式に入れた値 $1/\bigl(f + (1-f)/p\bigr)$ は $g + (1-g)p$ に一致する。
</Theorem>

<Proof of="thm-gustafson">
$t_{\mathrm{s}} = gT(p)$、$t_{\mathrm{p}} = (1-g)T(p)$ なので $T(1) = t_{\mathrm{s}} + p\,t_{\mathrm{p}} = \bigl(g + (1-g)p\bigr)T(p)$ で、最初の式を得ます。

$f = t_{\mathrm{s}}/T(1) = gT(p)/\bigl((g + (1-g)p)T(p)\bigr) = g/(g + (1-g)p)$ です。これを $g$ について解きます。分母を払うと $fg + f(1-g)p = g$、すなわち $fp = g(1 - f + fp)$ で、$1 - f + fp \gt 0$（$0 \le f \le 1$、$p \ge 1$ から $1 - f + fp \ge \min(1, p) \gt 0$）なので $g = fp/(1 - f + fp)$ です。

最後に $1 - f = (1-g)p/(g + (1-g)p)$ から
$$
f + \frac{1-f}{p} = \frac{g + (1-g)}{g + (1-g)p} = \frac{1}{g + (1-g)p}
$$
で、逆数をとると $g + (1-g)p$ です。
</Proof>

2 つの法則は矛盾しません。違いは、逐次部分の割合を何に対して測るかだけです。Amdahl の $f$ は「1 台で走らせたときの時間」に対する割合、Gustafson の $g$ は「$p$ 台で走らせたときの時間」に対する割合で、同じ実行でも $p$ 台のほうが全体の時間が短いぶん、逐次部分の割合は大きく見えます。たとえば $p = 16$ 台で時間 $1$ のうち逐次部分が $0.05$ なら $g = 0.05$ で、1 台なら $0.05 + 16 \times 0.95 = 15.25$ かかるので $f = 0.05/15.25 \approx 0.0033$ です。

同じ点を、もう少し詳しく述べます。同じ 1 回の実行について、逐次部分の割合を 1 台での時間に対して測ったのが $f$、$p$ 台での時間に対して測ったのが $g$ です。$g$ を固定して $p$ を増やすと $f = g/(g + (1-g)p)$ は $0$ に近づきます。つまり Gustafson の法則が台数に比例する速度向上を与えるのは、「問題を大きくすると並列部分だけが増え、逐次部分は増えない」という仮定のおかげです。格子を細かくすると逐次の前処理も重くなる計算では、この仮定が成り立ちません。

### 3.3. 測定から逐次割合を逆算する

$f$ は実測しにくい量です。逆に、測定した $S(p)$ から「Amdahl の法則が等号で成り立つとしたら $f$ はいくらか」を逆算することはできます。

<Proposition id="prop-karp-flatt" title="Karp–Flatt 指標">
$p \ge 2$ と $S \gt 0$ に対し
$$
e(p) = \frac{\dfrac{1}{S(p)} - \dfrac{1}{p}}{1 - \dfrac{1}{p}}
$$
とおく。所要時間が、逐次割合 $f \in [0,1]$ とオーバーヘッドの割合 $\omega(p) \ge 0$ を使って
$$
T(p) = T(1)\Bigl(f + \frac{1-f}{p} + \omega(p)\Bigr)
$$
と書けるとき、次が成り立つ。

1. $e(p) = f + \dfrac{p}{p-1}\,\omega(p)$。
2. $e(p) = f$ であることと、$\omega(p) = 0$、すなわち $S(p) = 1/\bigl(f + (1-f)/p\bigr)$ であることは同値。
3. $2 \le p \lt p'$ で $e(p) \lt e(p')$ ならば $\omega(p) \lt \omega(p')$。
</Proposition>

証明は式変形だけです。$\omega(p)$ は、仕事を分けたことで新しく生じる余分な時間（通信、同期の待ち、負荷の偏りなど）の、1 台での時間に対する割合です。$e(p)$ の式から $f$ の部分を取り出すと、残りが $\omega(p)$ の定数倍になることを示します。

<Proof of="prop-karp-flatt">
1. $1/S(p) = T(p)/T(1) = f + (1-f)/p + \omega(p)$ なので
   $$
   \frac{1}{S(p)} - \frac{1}{p} = f - \frac{f}{p} + \omega(p) = f\Bigl(1 - \frac1p\Bigr) + \omega(p)
   $$
   です。$1 - 1/p = (p-1)/p \gt 0$ で割ると $e(p) = f + \frac{p}{p-1}\omega(p)$ です。

2. 1 から $e(p) = f \iff \omega(p) = 0$ で、$\omega(p) = 0$ は $1/S(p) = f + (1-f)/p$ と同じです。

3. 1 から $\frac{p}{p-1}\omega(p) \lt \frac{p'}{p'-1}\omega(p')$ です。$\frac{p}{p-1} = 1 + \frac{1}{p-1}$ は $p$ について減少するので $\frac{p'}{p'-1} \lt \frac{p}{p-1}$ で、$\omega(p) \ge 0$ から $\frac{p'}{p'-1}\omega(p) \le \frac{p}{p-1}\omega(p) \lt \frac{p'}{p'-1}\omega(p')$ です。両辺を $\frac{p'}{p'-1} \gt 0$ で割ると $\omega(p) \lt \omega(p')$ です。
</Proof>

速度向上や効率の表だけを見ていると、効率が落ちる原因が逐次部分（台数によらない）なのか、通信・同期・負荷の偏りなどのオーバーヘッド（台数とともに増える）なのかが区別できません。$e(p)$ を並べると、一定なら前者、増えていれば後者と読めます。対策は違います。前者なら逐次部分そのものを速くするしかなく、後者なら分割の粒度や通信の回数を見直します。

<Example id="ex-amdahl-table" title="Amdahl・Gustafson の数値表と Karp–Flatt 指標">
3 つの式を表にします。Karp–Flatt 指標は、$T(1) = 1$、$f = 0.05$ として、時間のモデル $T(p) = f + (1-f)/p + c\,(p-1)$（台数に比例するオーバーヘッド）から計算した $S(p)$ に適用します。測定値ではなく、モデルの式の値です。

```python
# Python 3.12.14（標準ライブラリのみ）
def amdahl(f, p):            # S(p) = 1/(f + (1-f)/p)
    return 1.0 / (f + (1.0 - f) / p)

def gustafson(g, p):         # 拡大速度向上 g + (1-g)p
    return g + (1.0 - g) * p

def karp_flatt(S, p):        # e(p) = (1/S - 1/p)/(1 - 1/p)
    return (1.0 / S - 1.0 / p) / (1.0 - 1.0 / p)

ps = [1, 2, 4, 8, 16, 64, 256, 1024]
print("Amdahl  S(p) = 1/(f + (1-f)/p)")
print("     f " + "".join(f"{p:>8d}" for p in ps) + "     1/f")
for f in (0.01, 0.05, 0.2):
    print(f"  {f:4.2f} " + "".join(f"{amdahl(f, p):8.2f}" for p in ps) + f"{1/f:8.1f}")
print("Gustafson  g + (1-g)p")
for g in (0.01, 0.05, 0.2):
    print(f"  {g:4.2f} " + "".join(f"{gustafson(g, p):8.2f}" for p in ps))

# Karp–Flatt: 時間のモデル T(p) = f + (1-f)/p + c(p-1)（T(1) = 1）
print("Karp-Flatt  e(p)   （f = 0.05）")
print("     c " + "".join(f"{p:>8d}" for p in ps[1:]))
for c in (0.0, 0.001):
    row = []
    for p in ps[1:]:
        T = 0.05 + 0.95 / p + c * (p - 1)
        row.append(karp_flatt(1.0 / T, p))
    print(f" {c:5.3f} " + "".join(f"{e:8.3f}" for e in row))
```

```text
Amdahl  S(p) = 1/(f + (1-f)/p)
     f        1       2       4       8      16      64     256    1024     1/f
  0.01     1.00    1.98    3.88    7.48   13.91   39.26   72.11   91.18   100.0
  0.05     1.00    1.90    3.48    5.93    9.14   15.42   18.62   19.64    20.0
  0.20     1.00    1.67    2.50    3.33    4.00    4.71    4.92    4.98     5.0
Gustafson  g + (1-g)p
  0.01     1.00    1.99    3.97    7.93   15.85   63.37  253.45 1013.77
  0.05     1.00    1.95    3.85    7.65   15.25   60.85  243.25  972.85
  0.20     1.00    1.80    3.40    6.60   13.00   51.40  205.00  819.40
Karp-Flatt  e(p)   （f = 0.05）
     c        2       4       8      16      64     256    1024
 0.000    0.050   0.050   0.050   0.050   0.050   0.050   0.050
 0.001    0.052   0.054   0.058   0.066   0.114   0.306   1.074
```

**検算.** (a) $f = 0.05$、$p = 16$ は手計算で $1/(0.05 + 0.95/16) = 1/0.109375 \approx 9.14$ で、表と一致します。どの行も $1/f$ を超えていません（<Ref to="thm-amdahl" />）。(b) $g = 0.05$、$p = 16$ は $0.05 + 0.95 \times 16 = 15.25$ です。(c) $c = 0$ の行は $e(p) = f = 0.05$ で一定（<Ref to="prop-karp-flatt" /> の 2）、$c = 0.001$ の行は 1 の式 $0.05 + \frac{p}{p-1}\cdot0.001(p-1) = 0.05 + 0.001p$ に一致し（$p = 64$ で $0.114$）、台数とともに増えています。$p = 1024$ で $e \gt 1$ になるのは、$T(1024) \approx 1.07 \gt T(1)$、つまり 1 台より遅くなっていることを表します。
</Example>

## 4. 仕事量と依存鎖：Brent の評価

Amdahl の法則の $f$ は、計算のどこが「逐次」なのかを問いません。依存関係の構造から速度向上の上限を導くには、計算を有向グラフで表します。

<Definition id="def-work-span" title="計算グラフ・仕事量・依存鎖の長さ">
有限個の単位時間の作業を頂点とし、作業 $u$ の結果を作業 $v$ が使うときに辺 $u \to v$ を引いた、閉路のない有向グラフを**計算グラフ**という。頂点の数 $W$ を**仕事量**（work）、有向道に含まれる頂点の数の最大値 $L$ を**最長の依存鎖の長さ**（span、深さ）という。$p$ 台での**スケジュール**とは、各頂点 $v$ に実行するステップの番号 $t(v) \in \{1, 2, \ldots\}$ を割り当てる写像で、辺 $u \to v$ があれば $t(u) \lt t(v)$、同じ番号の頂点は $p$ 個以下、という条件を満たすものをいう。使うステップの数の最小値を $T_p$ と書く。
</Definition>

<Proposition id="prop-work-span" title="Brent の評価">
仕事量 $W \ge 1$、最長の依存鎖の長さ $L$ の計算グラフと整数 $p \ge 1$ について
$$
\max\Bigl(\frac{W}{p},\ L\Bigr) \le T_p \le \Bigl\lfloor\frac{W}{p}\Bigr\rfloor + L
$$
が成り立つ。右の不等式は、各ステップで「前提の作業がすべて済んでいる未実行の作業」を、$p$ 個以上あれば $p$ 個、なければ全部実行する貪欲なスケジュールで達成される。
</Proposition>

下界は、「1 ステップで $p$ 個までしか進められない」ことと「依存の鎖は順番にしか進められない」ことから直ちに出ます。上界の証明の方針は次のとおりです。貪欲なスケジュールのステップを、$p$ 個すべてを使えたステップ（満杯）と、使い切れなかったステップ（不足）に分けます。満杯のステップは多くても $W/p$ 回です。不足のステップでは、実行できる作業をすべて実行するので、残りの作業の依存鎖が必ず 1 つ短くなり、不足のステップは $L$ 回までしかありません。両者を足したものが上界です。以下の証明に出てくる「源点」は、まだ実行していない作業のうち、前提となる作業がすべて済んでいるもの、つまり今すぐ実行できる作業のことです。

<Proof of="prop-work-span">
**下界.** 1 ステップで実行できる頂点は $p$ 個以下なので、$T_p$ ステップで実行できる頂点は $pT_p$ 個以下で、$W \le pT_p$ です。また、$L$ 個の頂点を含む有向道 $v_1 \to v_2 \to \cdots \to v_L$ をとると、スケジュールの条件から $t(v_1) \lt t(v_2) \lt \cdots \lt t(v_L)$ で、これらは相異なる正の整数なので $t(v_L) \ge L$ です。よって $T_p \ge L$ です。

**上界.** 貪欲なスケジュールを考えます。ステップ $t$ の開始時点で未実行の頂点が誘導する部分グラフを $G_t$、$G_t$ の有向道に含まれる頂点の数の最大値を $\ell_t$ とします（$G_1$ はもとのグラフで $\ell_1 = L$）。頂点 $v$ の前提がすべて済んでいることは、$v$ に入る辺の始点が $G_t$ に無いこと、つまり $v$ が $G_t$ の中で入る辺をもたないこと（$G_t$ の**源点**であること）と同じです。$G_t$ が空でなければ源点があります。頂点から入る辺を逆向きにたどり続けると、閉路がないので同じ頂点に戻らず、頂点は有限個なので、入る辺をもたない頂点で止まるからです。したがって未実行の頂点がある限り、各ステップで少なくとも 1 個が実行され、スケジュールは有限ステップで終わります。辺 $u \to v$ があるとき、$v$ は $u$ が実行された後のステップで初めて源点になるので、$t(u) \lt t(v)$ も満たされます。

ステップを 2 種類に分けます。$p$ 個を実行するステップを**満杯のステップ**、$p$ 個未満（源点の全部）を実行するステップを**不足のステップ**と呼びます。

満杯のステップは 1 回につき $p$ 個の頂点を実行するので、その回数は $\lfloor W/p\rfloor$ 以下です。

不足のステップでは $\ell_{t+1} \le \ell_t - 1$ を示します。$G_{t+1}$ は $G_t$ から源点をすべて除いたグラフです。$G_{t+1}$ の有向道 $w_1 \to \cdots \to w_k$ を任意にとると、$w_1$ は $G_t$ の源点ではないので、$G_t$ の中に辺 $w_0 \to w_1$ をもつ頂点 $w_0$ があります。閉路がないので $w_0$ は $w_1, \ldots, w_k$ のどれとも異なり、$w_0 \to w_1 \to \cdots \to w_k$ は $G_t$ の有向道で $k + 1$ 個の頂点を含みます。よって $k + 1 \le \ell_t$ で、$\ell_{t+1} \le \ell_t - 1$ です。満杯のステップでは $G_{t+1}$ は $G_t$ の部分グラフなので $\ell_{t+1} \le \ell_t$ です。$\ell_1 = L$、$\ell_t \ge 0$ なので、不足のステップは $L$ 回以下です。

合わせて、貪欲なスケジュールのステップ数は $\lfloor W/p\rfloor + L$ 以下で、$T_p$ はそれ以下です。
</Proof>

速度向上に直すと、$T_1 = W$ なので $S(p) = W/T_p \le \min(p, W/L)$ です。$W/L$ は台数をいくら増やしても超えられない上限で、計算の**並列度**と呼ばれます。上界の側は、台数が並列度より十分小さければ（$p \ll W/L$）、$T_p \le W/p + L \approx W/p$ で、ほぼ台数倍の速度向上が達成できることを言っています。

<Example id="ex-sum-span" title="逐次和と対ごとの和の依存鎖">
$n$ 個の数の和を考えます。どの計算方法でも加算は $n - 1$ 回なので $W = n - 1$ です。

逐次和 $((x_1 + x_2) + x_3) + \cdots$ では、各加算が直前の加算の結果を使うので $L = n - 1 = W$ で、並列度は $1$ です。何台あっても速くなりません。

対ごとの和（配列を半分に分けて、それぞれの和を足す）では、$n$ 個の和の依存鎖の長さ $D(n)$ は $D(1) = 0$、$D(n) = 1 + D(\lceil n/2\rceil)$ を満たします（大きいほうの半分が長い）。$2^{k-1} \lt n \le 2^k$ なら $2^{k-2} \lt \lceil n/2\rceil \le 2^{k-1}$ なので、$k$ についての帰納法で $D(n) = k = \lceil\log_2 n\rceil$ です。<Ref to="prop-work-span" /> から、$p$ 台で $\lfloor(n-1)/p\rfloor + \lceil\log_2 n\rceil$ ステップ以下で終わります。$n = 10^6$、$p = 16$ なら $62499 + 20 = 62519$ 以下で、下界 $\lceil(10^6 - 1)/16\rceil = 62500$ とほとんど差がありません。

同じ $n - 1$ 回の加算でも、括弧の付け方で並列度が $1$ から $n/\log_2 n$ の程度まで変わります。その代わり、括弧を付け替えると丸めの入り方が変わり、結果のビット列は一般に変わります（§6）。
</Example>

## 5. メモリ帯域と加速器：roofline モデル

### 5.1. 演算強度と性能の上限

コアを増やしても速くならない計算のもう 1 つの典型は、演算ではなくメモリの読み書きが律速している計算です。

<Definition id="def-arithmetic-intensity" title="演算強度">
ある計算が行う浮動小数演算の回数を $F$、その間に主記憶とプロセッサ（キャッシュ）の間で転送されるデータの量を $Q \gt 0$ バイトとする。$I = F/Q$ を**演算強度**（arithmetic intensity）という。単位は演算数／バイトである。
</Definition>

<Proposition id="prop-roofline" title="roofline モデル">
1 秒あたり高々 $\pi \gt 0$ 回の浮動小数演算を行え（ピーク演算性能）、主記憶との間で 1 秒あたり高々 $\beta \gt 0$ バイトを転送できる（メモリ帯域）計算機で、演算数 $F \gt 0$、転送量 $Q \gt 0$、演算強度 $I = F/Q$ の計算を時間 $T$ で実行したとする。このとき性能 $P = F/T$ は
$$
P \le \min(\pi,\ \beta I)
$$
を満たす。2 つの上限が一致する演算強度 $I^{*} = \pi/\beta$ より $I$ が小さい計算では、上限は $\beta I$ で、$\pi$ によらない。
</Proposition>

<Proof of="prop-roofline">
時間 $T$ の間に行える演算は $\pi T$ 回以下なので $F \le \pi T$、転送できるデータは $\beta T$ バイト以下なので $Q \le \beta T$ です。前者から $P = F/T \le \pi$、後者から $T \ge Q/\beta$ で $P = F/T \le F\beta/Q = \beta I$ です。$I \lt I^{*}$ なら $\beta I \lt \beta I^{*} = \pi$ なので $\min(\pi, \beta I) = \beta I$ です。
</Proof>

<Figure caption="roofline モデル（両対数）。横軸は演算強度 I、縦軸は性能の上限。斜めの部分はメモリ帯域 β、水平の部分はピーク演算性能 π で決まる。数値は仮の計算機（π = 100 GFLOP/s、β = 20 GB/s）のもの。">
<svg viewBox="0 0 640 310" width="100%" role="img" aria-label="roofline モデルの図。演算強度に対する性能の上限">
  <line x1="90" y1="250" x2="610" y2="250" stroke="currentColor" stroke-width="1.5" />
  <line x1="90" y1="250" x2="90" y2="25" stroke="currentColor" stroke-width="1.5" />
  <polyline points="90,240.3 406.1,50.7 590,50.7" fill="none" stroke="var(--sl-color-accent)" stroke-width="2.5" />
  <line x1="406.1" y1="250" x2="406.1" y2="50.7" stroke="currentColor" stroke-width="1" stroke-dasharray="5 4" opacity="0.6" />
  <circle cx="190" cy="180.3" r="5" fill="var(--sl-color-accent)" />
  <circle cx="440" cy="50.7" r="5" fill="var(--sl-color-accent)" />
  <text x="200" y="196" fill="currentColor" font-size="12">行列ベクトル積（I = 1/4）</text>
  <text x="448" y="42" fill="currentColor" font-size="12">ブロック化した行列積（I = 8）</text>
  <text x="180" y="120" fill="currentColor" font-size="12">β I（帯域で決まる）</text>
  <text x="500" y="70" fill="currentColor" font-size="12">π（演算で決まる）</text>
  <text x="406" y="268" fill="currentColor" font-size="12" text-anchor="middle">I* = π/β = 5</text>
  <text x="90" y="268" fill="currentColor" font-size="12" text-anchor="middle">1/16</text>
  <text x="190" y="268" fill="currentColor" font-size="12" text-anchor="middle">1/4</text>
  <text x="290" y="268" fill="currentColor" font-size="12" text-anchor="middle">1</text>
  <text x="490" y="268" fill="currentColor" font-size="12" text-anchor="middle">16</text>
  <text x="590" y="268" fill="currentColor" font-size="12" text-anchor="middle">64</text>
  <text x="82" y="254" fill="currentColor" font-size="12" text-anchor="end">1</text>
  <text x="82" y="154" fill="currentColor" font-size="12" text-anchor="end">10</text>
  <text x="82" y="55" fill="currentColor" font-size="12" text-anchor="end">100</text>
  <text x="350" y="295" fill="currentColor" font-size="12" text-anchor="middle">演算強度 I（演算数／バイト）</text>
  <text x="20" y="20" fill="currentColor" font-size="12">性能の上限（GFLOP/s）</text>
</svg>
</Figure>

図や例で使う GFLOP/s は、1 秒あたり $10^9$ 回の浮動小数演算を表す単位です。

命題は自明に近いのですが、使い道ははっきりしています。$I \lt I^{*}$ の計算では、コアを増やして $\pi$ を上げても上限 $\beta I$ は動きません。メモリ帯域は同じソケットのコアが共有するからです。速くするには $I$ を上げる、つまり同じ演算数をより少ないメモリ転送で行うように計算の順序を組み替えるしかありません。

<Example id="ex-intensity" title="行列ベクトル積と行列積の演算強度">
倍精度（1 要素 8 バイト）の $n \times n$ 行列について、演算数と転送量を数えます。転送量は「キャッシュに収まるデータは 1 回読めば使い回せ、収まらないデータは使うたびに主記憶から読む」という単純なモデルで数えます。$n$ は行列全体（$8n^2$ バイト）がキャッシュに収まらないほど大きく、長さ $n$ のベクトル数本（$8n$ バイトずつ）は収まるとします。

**行列ベクトル積** $\boldsymbol{y} = A\boldsymbol{x}$。乗算 $n^2$ 回と加算 $n(n-1)$ 回で $F = 2n^2 - n$ です。$A$ の各要素を 1 回読み（$8n^2$）、$\boldsymbol{x}$ を読み $\boldsymbol{y}$ を読み書きする（$24n$）ので $Q = 8n^2 + 24n$、$I = (2n^2 - n)/(8n^2 + 24n) \to 1/4$ です。$A$ の要素は 1 回ずつしか使わないので、順序をどう工夫しても $I$ は $1/4$ を超えません。

**行列積 $C = AB$、ブロック化なし。** $c_{ij} = \sum_k a_{ik}b_{kj}$ を $i, j, k$ の三重ループで計算します。$F = n^2(2n - 1)$ です。$A$ の第 $i$ 行は $j$ を動かす間キャッシュに残りますが、$B$ の第 $j$ 列は $j$ ごとに別の $8n$ バイトで、$B$ 全体は収まらないので、$(i, j)$ の組ごとに読み直します。$Q \ge 8n^3$ で、$I \le (2n^3 - n^2)/(8n^3) \lt 1/4$ です。

**行列積、ブロック化あり。** $n$ を $b$ の倍数とし、行列を $b \times b$ のブロックに分けて $C_{IJ} \mathrel{+}= A_{IK}B_{KJ}$ をブロック単位で行います。$b$ は 3 つのブロック（$24b^2$ バイト）がキャッシュに収まるように選びます。$(I, J, K)$ の組は $(n/b)^3$ 個で、1 組につき $A_{IK}$ と $B_{KJ}$ を読み（$16b^2$ バイト）、乗算 $b^3$ 回と加算 $b^3$ 回を行います。$C_{IJ}$ は $K$ を動かす間キャッシュに置き、$(I, J)$ ごとに 1 回読んで 1 回書きます（$16b^2$ バイト）。合計は
$$
F = 2n^3, \qquad Q = \Bigl(\frac nb\Bigr)^3 16b^2 + \Bigl(\frac nb\Bigr)^2 16b^2 = \frac{16n^3}{b} + 16n^2, \qquad I = \frac{2n^3}{16n^3/b + 16n^2} \approx \frac{b}{8}
$$
です。演算強度は $O(1)$ から $O(b)$ に上がります。$b = 64$ なら $I \approx 8$ です。

図の仮の計算機（$\pi = 100$ GFLOP/s、$\beta = 20$ GB/s、$I^{*} = 5$）では、<Ref to="prop-roofline" /> の上限は、行列ベクトル積とブロック化なしの行列積で $20 \times 1/4 = 5$ GFLOP/s（ピークの $5\%$）、ブロック化した行列積で $\min(100, 160) = 100$ GFLOP/s です。最適化された BLAS の行列積が、自分で書いた三重ループより桁違いに速い主な理由がこれです。行列積を自分で書かない、というのが実務上の結論です。
</Example>

### 5.2. 加速器に仕事を移す損得

GPU は多数の演算器と専用のメモリをもつ加速器です。NumPy と同じ書き方で配列を GPU の上に置いて計算するライブラリ（CuPy など）を使うと、配列演算を GPU に移せます。ただし配列は、計算機本体（ホスト）のメモリと GPU（デバイス）のメモリの間を転送しなければなりません。

<Proposition id="prop-offload-breakeven" title="加速器に移す損益分岐">
ホストで時間 $W \gt 0$ かかる計算を、同じ計算を $s \gt 0$ 倍の速さで行う加速器に移すとする。移すにはデータを $D \ge 0$ 単位だけ転送する必要があり、転送には 1 単位あたり時間 $\tau \ge 0$ と、データ量によらない固定の時間 $c_0 \ge 0$（呼び出しと同期の待ち）がかかるとする。このとき、加速器に移したほうが所要時間が短くなることは
$$
W\Bigl(1 - \frac1s\Bigr) \gt D\tau + c_0
$$
と同値である。とくに $s \le 1$ なら決して得にならず、$s$ がどれだけ大きくても $W \le D\tau + c_0$ なら得にならない。
</Proposition>

<Proof of="prop-offload-breakeven">
加速器に移したときの所要時間は $W/s + D\tau + c_0$ です。これが $W$ より小さいことは、$W/s$ を移項して $W - W/s \gt D\tau + c_0$ と同じです。$s \le 1$ なら左辺は $W(1 - 1/s) \le 0 \le D\tau + c_0$ で不等式は成り立ちません。$s \gt 0$ のとき $1 - 1/s \lt 1$ なので左辺は $W$ より小さく、$W \le D\tau + c_0$ なら不等式は成り立ちません。
</Proof>

<Figure caption="加速器に移す損益分岐の模式図（単位は任意）。ホストの時間は仕事量に比例して増える。加速器の時間は傾きが小さいが、固定の時間 c0 から始まる。2 本の線が交わる点より小さい計算では、ホストのほうが速い。">
<svg viewBox="0 0 640 300" width="100%" role="img" aria-label="ホストでの時間と、加速器での時間（転送を含む）を、計算の大きさに対して比べた模式図">
  <line x1="80" y1="250" x2="580" y2="250" stroke="currentColor" stroke-width="1.5" />
  <line x1="80" y1="250" x2="80" y2="30" stroke="currentColor" stroke-width="1.5" />
  <line x1="80" y1="250" x2="560" y2="40" stroke="currentColor" stroke-width="2.5" />
  <line x1="80" y1="208" x2="560" y2="124" stroke="var(--sl-color-accent)" stroke-width="2.5" />
  <line x1="240" y1="250" x2="240" y2="180" stroke="currentColor" stroke-width="1" stroke-dasharray="5 4" opacity="0.6" />
  <circle cx="240" cy="180" r="5" fill="var(--sl-color-accent)" />
  <text x="440" y="86" fill="currentColor" font-size="12">ホストでの時間 W</text>
  <text x="340" y="150" fill="var(--sl-color-accent)" font-size="12">加速器 W/s + Dτ + c0</text>
  <text x="68" y="212" fill="currentColor" font-size="12" text-anchor="end">c0</text>
  <text x="240" y="268" fill="currentColor" font-size="12" text-anchor="middle">損益分岐</text>
  <text x="330" y="290" fill="currentColor" font-size="12" text-anchor="middle">計算の大きさ（データ量 D に比例）</text>
  <text x="20" y="22" fill="currentColor" font-size="12">所要時間</text>
</svg>
</Figure>

小さい配列で GPU が遅い理由は、この式から読み取れます。要素ごとの演算（$W = an$、$D = dn$、$n$ は要素数）では、条件は $n\bigl(a(1 - 1/s) - d\tau\bigr) \gt c_0$ です。括弧の中が正でも、$n$ が $c_0/\bigl(a(1-1/s) - d\tau\bigr)$ より小さければ固定の時間 $c_0$ を取り返せません。括弧の中が負（1 要素の計算より 1 要素の転送のほうが高くつく）なら、$n$ をいくら大きくしても損です。一方、行列積のように $W$ が $n^3$、$D$ が $n^2$ に比例する計算は、$n$ が大きくなるほど有利になります。

実務での対処は、転送の回数を減らすことです。配列を一度デバイスに置いたら、一連の計算をデバイスの上で済ませ、最後の結果（多くはスカラーか小さな配列）だけをホストに戻します。ループの中でホストとデバイスを往復させると、$D\tau + c_0$ を反復のたびに払うことになります。

## 6. 並列の和はどこまで一致するか

配列の和を $p$ 個のスレッドで計算すると、各スレッドが担当分の部分和を作り、最後に $p$ 個の部分和を足します。これは逐次和とは別の括弧づけなので、[再現性の章](/computer-science/research-computing/reproducible-experiments)で示したとおり結果は一般に変わります。変わる幅には上界があります。

次の定理のねらいは 2 つです。1 つめは、スレッド数を変えて結果が変わったとき、それが丸めの差で説明できる範囲かを判定する基準を作ることです。上界の $\gamma_{n-1}$ は §0 の事前知識で触れた、$n-1$ 回の加算で生じる相対誤差の上限を表す量で、再現性の章の結果をそのまま使います。2 つめは、分割の仕方を決めておけば結果がビット単位で同じになることを示すことです。

<Theorem id="thm-parallel-sum-bound" title="部分和に分けた和の差の上界と、分割の固定による再現">
$x_1, \ldots, x_n$（$n \ge 2$）を浮動小数点数、$u$ を単位丸めとし、$(n-1)u \lt 1$、$\gamma_{n-1} = (n-1)u/(1 - (n-1)u)$ とする。添字の集合 $\{1, \ldots, n\}$ を空でない $p$ 個の組に分け、各組の和を浮動小数の加算で（任意の順序・任意の括弧づけで）計算し、得られた $p$ 個の部分和を浮動小数の加算で（任意の順序・任意の括弧づけで）足した結果を $\hat s_{\mathrm{par}}$ とする。逐次和の結果を $\hat s_{\mathrm{seq}}$ とする。どの計算でもオーバーフローとアンダーフローが起きないならば、次が成り立つ。

1. $|\hat s_{\mathrm{par}} - \hat s_{\mathrm{seq}}| \le 2\gamma_{n-1}\displaystyle\sum_{i=1}^n|x_i|$。
2. 組への分け方、各組の中の加算の順序と括弧づけ、部分和を足す順序と括弧づけがすべて同じ 2 回の実行は、加算が同じ丸めの方向で正しく丸められる限り、どの加算をどのスレッドがいつ実行したかによらず、同じ浮動小数点数を返す。
</Theorem>

<Proof of="thm-parallel-sum-bound">
1. 各組の和の計算は、その組の $x_i$ を葉とし、内部の頂点が浮動小数の加算である二分木で表されます。部分和を足す計算は、$p$ 個の部分和を葉とする二分木です。後者の葉に前者の木を差し込むと、$x_1, \ldots, x_n$ を葉とし、$n - 1$ 個の加算を内部の頂点とする 1 つの二分木になります。これは再現性の章の意味での $x_1, \ldots, x_n$ の和の木（<Ref to="computer-science/research-computing/reproducible-experiments#def-summation-tree" />）で、$\hat s_{\mathrm{par}}$ はその計算結果です。$\hat s_{\mathrm{seq}}$ も、添字の順に左から括弧を付けた和の木（高さ $n - 1$ の一本道の木）の計算結果です。同じ $x_1, \ldots, x_n$ の 2 つの和の木の計算結果の差は $2\gamma_{n-1}\sum_i|x_i|$ 以下（<Ref to="computer-science/research-computing/reproducible-experiments#thm-summation-order-bound" />）なので、主張を得ます。

2. 正しく丸められる加算の結果は、厳密な和 $a + b$ を決まった方向に丸めた値で（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#def-correct-rounding" />）、2 つの入力 $a, b$ だけで決まります。1 の二分木の各頂点の値が、2 回の実行で一致することを、頂点の高さ（葉までの最長の道の長さ）についての帰納法で示します。高さ $0$ の頂点は葉で、値は入力 $x_i$ そのものなので一致します。高さ $h \ge 1$ の頂点の値は、2 つの子（高さ $h - 1$ 以下）の値の加算の結果で、子の値は帰納法の仮定から 2 回の実行で一致するので、この頂点の値も一致します。仮定から 2 回の実行の二分木は同じなので、根の値、すなわち結果が一致します。どのスレッドがいつ計算したかは、二分木にも入力にも現れません。
</Proof>

1 は「スレッド数を変えたら結果が変わった」ときの判定規則を与えます。差がこの上界の内側なら丸めの差として説明がつきます。上界を超えていたら、和の順序ではない別の原因（競合状態、初期化されていない変数、担当範囲の重複や抜け）を疑います。上界 $2\gamma_{n-1}\sum|x_i|$ は最悪の場合の評価で、実際の差はずっと小さいのが普通です。対ごとの和どうしなら、<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#thm-pairwise-sum" /> の $\gamma_{\lceil\log_2 n\rceil}$ を使ったより小さい評価ができます。

2 は対処法を与えます。多くの並列ライブラリは、既定ではスレッド数に応じて分割を決めるので、スレッド数を変えると二分木が変わります。分割を「長さ $B$ のブロックごと」のようにスレッド数によらない規則で決め、部分和を番号順に足せば、何スレッドで実行しても同じビット列になります。その分、スレッド数に合わせた最適な分割はできなくなります。

<Example id="ex-parallel-sum" title="分割数を変えた和と、分割を固定した和">
乱数種を固定して作った $10^6$ 個の binary64 の数（大きさが $10^{-8}$ から $10^{8}$ に散らばる）を、$p$ 個の連続した区間に分けて部分和を作り、部分和を番号順に足します。プロセスやスレッドは起こさず、分割だけを再現しています（<Ref to="thm-parallel-sum-bound" /> の 2 から、結果は実際に並列で実行した場合と同じです）。正しく丸めた和は `math.fsum` で求めます。

```python
# Python 3.12.14 / NumPy 2.5.2
import math
import numpy as np

rng = np.random.default_rng(20261008)
n = 10**6
x = rng.standard_normal(n) * 10.0 ** rng.integers(-8, 9, size=n)

def seq_sum(a):                      # 逐次和（左から 1 つずつ）
    s = 0.0
    for v in a.tolist():
        s += v
    return s

exact = math.fsum(x)                 # 正しく丸めた和
sum_abs = math.fsum(np.abs(x))
u = 2.0 ** -53
gamma = (n - 1) * u / (1 - (n - 1) * u)
bound = 2 * gamma * sum_abs
s1 = seq_sum(x)
print(f"n = {n},  Σ|x_i| = {sum_abs:.6e},  2γ_(n-1)Σ|x_i| = {bound:.3e}")
print("  p  p 個の部分和を足した結果      逐次和との差  fsum との差")
for p in (1, 2, 4, 8, 16):
    edges = [n * k // p for k in range(p + 1)]
    parts = [seq_sum(x[edges[k]:edges[k + 1]]) for k in range(p)]
    s = seq_sum(np.array(parts))
    print(f"{p:3d}  {s:.17e}  {abs(s - s1):.3e}     {abs(s - exact):.3e}")

# 分割を固定する: ブロック長 4096（ワーカー数によらない）、部分和は番号順に足す
B = 4096
nblocks = (n + B - 1) // B
def fixed_blocks(workers):
    parts = [None] * nblocks
    for w in range(workers):                     # ワーカー w の担当は j ≡ w (mod workers)
        for j in range(w, nblocks, workers):
            parts[j] = seq_sum(x[j * B:(j + 1) * B])
    return seq_sum(np.array(parts))
vals = [fixed_blocks(w) for w in (1, 2, 4, 8, 16)]
print("固定ブロック:", [v.hex() for v in vals])
print("すべて同じビット列:", len({v.hex() for v in vals}) == 1)
```

```text
n = 1000000,  Σ|x_i| = 5.249892e+12,  2γ_(n-1)Σ|x_i| = 1.166e+03
  p  p 個の部分和を足した結果      逐次和との差  fsum との差
  1  -1.97246295270851440e+10  0.000e+00     5.150e-04
  2  -1.97246295270854721e+10  3.281e-04     1.869e-04
  4  -1.97246295270855408e+10  3.967e-04     1.183e-04
  8  -1.97246295270857315e+10  5.875e-04     7.248e-05
 16  -1.97246295270857353e+10  5.913e-04     7.629e-05
固定ブロック: ['-0x1.25eb7d85c57b7p+34', '-0x1.25eb7d85c57b7p+34', '-0x1.25eb7d85c57b7p+34', '-0x1.25eb7d85c57b7p+34', '-0x1.25eb7d85c57b7p+34']
すべて同じビット列: True
```

**検算.** (a) 分割数 $p$ を変えると結果は有効数字の 15 桁目から変わり、逐次和との差は $3 \times 10^{-4}$ から $6 \times 10^{-4}$ です。上界 $1.2 \times 10^{3}$ の内側にあり、上界より 6 桁以上小さい値です。(b) 上界の値は手計算でも確かめられます。$\gamma_{n-1} \approx 10^6 \times 2^{-53} \approx 1.11 \times 10^{-10}$、これに $2 \times 5.25 \times 10^{12}$ を掛けて $1.17 \times 10^{3}$ です。(c) `fsum` との差（真の丸め誤差）も同じ程度の大きさで、この例では $p$ が大きいほうがおおむね小さくなっています。部分和に分けると、足し込まれる側の値が小さくなり、1 回の丸めの絶対的な大きさが減る、と説明できます。(d) ブロック長を固定した和は、担当の割り振りを 1・2・4・8・16 通りに変えても同じビット列でした。
</Example>

## 7. 並列の乱数：種の渡し方

Monte Carlo 法は「恥ずかしいほど並列」な計算の代表です。$p$ 個のワーカーがそれぞれ $m$ 個の標本を作り、平均を集めます。このとき各ワーカーに渡す乱数種を誤ると、台数を増やした意味がなくなります。

次の命題は、平均の分散（ばらつき）が、ワーカー間の相関によってどう変わるかを計算します。ワーカーどうしが無相関なら分散は台数に反比例して減り、同じ列を使うと台数を増やしても減らない、というのが結論です。証明は、和の分散を各項の分散と共分散の和に分けて数えるだけです。

<Proposition id="prop-stream-independence-usage" title="同じ種を渡すと分散は減らない">
$X_{j,i}$（$j = 1, \ldots, p$、$i = 1, \ldots, m$）を、平均 $\mu$、分散 $\sigma^2 \lt \infty$ の確率変数とし、各 $j$ について $X_{j,1}, \ldots, X_{j,m}$ は互いに無相関とする。$\bar Y_j = \frac1m\sum_{i=1}^m X_{j,i}$（ワーカー $j$ の標本平均）、$\bar Y = \frac1p\sum_{j=1}^p\bar Y_j$（全体の標本平均）とおく。

1. 定数 $\rho$ があって、すべての $j \ne k$ で $\mathrm{Cov}(\bar Y_j, \bar Y_k) = \rho\,\sigma^2/m$ ならば
   $$
   \mathrm{Var}(\bar Y) = \frac{\sigma^2}{pm}\bigl(1 + (p-1)\rho\bigr).
   $$
2. ワーカーどうしが無相関（$\rho = 0$）なら $\mathrm{Var}(\bar Y) = \sigma^2/(pm)$。
3. すべてのワーカーが同じ列を使う（すべての $j, i$ で $X_{j,i} = X_{1,i}$）なら $\mathrm{Var}(\bar Y) = \sigma^2/m$ で、2 の $p$ 倍である。
</Proposition>

<Proof of="prop-stream-independence-usage">
分散の基本性質（<Ref to="mathematics/probability/random-variables#prop-variance-properties" />）の 2 と 4 を使います。まず各 $j$ について、$X_{j,i}$ が互いに無相関なので 4 の共分散の項は消え、
$$
\mathrm{Var}(\bar Y_j) = \frac{1}{m^2}\mathrm{Var}\Bigl(\sum_{i=1}^m X_{j,i}\Bigr) = \frac{1}{m^2}\cdot m\sigma^2 = \frac{\sigma^2}{m}
$$
です。

1. 同じ性質の 4 を $\bar Y_1, \ldots, \bar Y_p$ に使うと、$j \lt k$ の組は $p(p-1)/2$ 個あるので
   $$
   \mathrm{Var}\Bigl(\sum_{j=1}^p\bar Y_j\Bigr) = p\,\frac{\sigma^2}{m} + 2\cdot\frac{p(p-1)}{2}\cdot\frac{\rho\,\sigma^2}{m} = \frac{p\,\sigma^2}{m}\bigl(1 + (p-1)\rho\bigr)
   $$
   で、$\bar Y$ はこれを $p$ で割ったものなので、2 から $p^2$ で割って主張の式を得ます。

2. 1 で $\rho = 0$ とします。

3. 仮定から $\bar Y_j = \bar Y_1$ がすべての $j$ で成り立つので $\bar Y = \bar Y_1$ で、$\mathrm{Var}(\bar Y) = \mathrm{Var}(\bar Y_1) = \sigma^2/m$ です（このとき $\mathrm{Cov}(\bar Y_j, \bar Y_k) = \mathrm{Var}(\bar Y_1) = \sigma^2/m$ なので、1 で $\rho = 1$ とした式とも一致します）。
</Proof>

[擬似乱数生成器](/computer-science/research-computing/reproducible-experiments#def-prng)の出力列は種で一意に決まります（<Ref to="computer-science/research-computing/reproducible-experiments#prop-seed-determinism" />）。したがって、すべてのワーカーに同じ種を渡すと 3 の状況がそのまま起こります。$pm$ 個の標本を作ったつもりで、実質は $m$ 個です。さらに悪いことに、$pm$ 個が独立だとして標準誤差を $s/\sqrt{pm}$ と計算すると、真の標準誤差 $\sigma/\sqrt m$ の $1/\sqrt p$ 倍に過小評価します。結果は「誤差棒の小さい、精度の高い推定」に見えます。

プロセスを fork（実行中のプロセスのメモリをそのまま複製して、子プロセスを作る方式）で複製する方式では、親プロセスで作った生成器の状態がそのまま子に複製されるので、種を明示的に渡していなくても同じことが起こります。種を `seed + j` のように連番で作るのも、生成器によっては列の間に相関が残るので勧められません。NumPy では `SeedSequence.spawn` で、1 つの種から互いに独立とみなせる子の種を作って渡します。子の列がなぜ重ならないとみなせるのか（ストリームの分け方の理論）は、擬似乱数生成器そのものを扱う章に譲り、ここでは使い方だけを示します。

<Example id="ex-spawn" title="SeedSequence.spawn の子を渡す雛形と、同じ種を渡した場合">
単位正方形の一様乱数が 4 分円に入る割合から $\pi$ を推定します。8 個の仕事（それぞれ標本 $10^5$ 個）をプロセスプールに渡します。

```python
# Python 3.12.14 / NumPy 2.5.2
from concurrent.futures import ProcessPoolExecutor
import numpy as np

M = 100_000                                   # ワーカー 1 つあたりの標本数

def worker(seed):
    """単位正方形の一様乱数で π/4 を推定する。戻り値は (的中数, 最初の乱数)"""
    rng = np.random.default_rng(seed)
    xy = rng.random((M, 2))
    hits = int(np.count_nonzero((xy ** 2).sum(axis=1) < 1.0))
    return hits, float(xy[0, 0])

def run(seeds, workers):
    with ProcessPoolExecutor(max_workers=workers) as ex:
        return list(ex.map(worker, seeds))    # map は入力の順に結果を返す

if __name__ == "__main__":
    p = 8
    good = run(np.random.SeedSequence(20261008).spawn(p), workers=4)
    again = run(np.random.SeedSequence(20261008).spawn(p), workers=2)
    bad = run([20261008] * p, workers=4)      # 誤り: 全員に同じ種
    for name, res in (("spawn の子", good), ("同じ種", bad)):
        hits = [r[0] for r in res]
        print(f"{name}: 的中数 = {hits}")
        print(f"   異なる的中数 = {len(set(hits))} 通り,  π の推定 = {4 * sum(hits) / (p * M):.5f}")
    print("ワーカー数 4 と 2 で同じ結果:", good == again)
    print("同じ種のときの最初の乱数:", sorted({r[1] for r in bad}))
```

```text
spawn の子: 的中数 = [78330, 78433, 78746, 78457, 78764, 78605, 78959, 78272]
   異なる的中数 = 8 通り,  π の推定 = 3.14283
同じ種: 的中数 = [78552, 78552, 78552, 78552, 78552, 78552, 78552, 78552]
   異なる的中数 = 1 通り,  π の推定 = 3.14208
ワーカー数 4 と 2 で同じ結果: True
同じ種のときの最初の乱数: [0.12563646004385953]
```

**検算.** (a) 同じ種を渡した 8 個の仕事は、的中数も最初の乱数もすべて同じで、同じ標本を 8 回数えただけです。(b) 1 個の標本（的中なら $4$、外れなら $0$）の標準偏差は $4\sqrt{q(1-q)}$（$q = \pi/4$）$\approx 1.642$ なので、標準誤差は、独立な $8 \times 10^5$ 個なら $1.642/\sqrt{8 \times 10^5} \approx 0.0018$、実質 $10^5$ 個なら $1.642/\sqrt{10^5} \approx 0.0052$ です（<Ref to="prop-stream-independence-usage" /> の 2 と 3。比は $\sqrt8$）。spawn の子の推定の誤差は $|3.14283 - \pi| \approx 0.0012$ で、標準誤差 $0.0018$ と矛盾しません。(c) 同じ種の推定 $3.14208$ はこの 1 回ではたまたま $\pi$ に近い値ですが、誤差の大きさを $0.0018$ と報告すれば、$\sqrt8 \approx 2.8$ 倍の過小評価です。1 回の結果の近さから方法の良し悪しは判断できません。(d) 子の種を仕事の番号に結びつけ、`map` が入力の順に結果を返すので、ワーカー数を 4 から 2 に変えても結果のリストは一致します。乱数の消費の順序がプロセスの実行順に左右されないからです。
</Example>

弱スケーリングの観点からも、Monte Carlo 法は並列に向いています。標準誤差を半分にするには標本を 4 倍にする必要があり、それには台数を 4 倍にして同じ時間だけ走らせればよく、逐次部分は種の配布と結果の集計だけです。<Ref to="thm-gustafson" /> の $g$ がほぼ $0$ の計算です。

## 8. 実務の作法

**まず測る。** 並列化の前に、プロファイラ（標準ライブラリの `cProfile` など）で時間のかかっている箇所を特定します。<Ref to="thm-amdahl" /> から、全体の $30\%$ を占める部分をいくら並列にしても $1/0.7 \approx 1.4$ 倍が上限です。逐次部分を先に速くする（Python のループを配列演算に直す、同じ計算の繰り返しをやめる）ほうが、たいてい効果が大きく、再現性も損ないません。

**スレッドの掛け算を避ける。** NumPy が呼ぶ BLAS は既定でコア数だけスレッドを使います。その上でプロセスを $p$ 個起こすと、スレッドの総数は $p \times$ コア数になり、コアの取り合いで文脈の切り替え（<Ref to="computer-science/cs-basics/operating-systems#def-context-switch" />）が増えて、1 プロセスのときより遅くなることがあります。<Ref to="prop-karp-flatt" /> の $e(p)$ が増える典型的な原因です。プロセスで並列にするときは、環境変数 `OMP_NUM_THREADS`、`OPENBLAS_NUM_THREADS`、`MKL_NUM_THREADS` を `1` にして BLAS を 1 スレッドにし、実行記録にスレッド数を残します。

**大きな配列をプロセス間で渡さない。** プロセスはアドレス空間が別なので、引数と戻り値は直列化（データをバイト列に変換すること）して複製されます。<Ref to="prop-offload-breakeven" /> と同じ損得勘定（$D\tau + c_0$）がここでも成り立ちます。ワーカーには配列ではなく「種と設定」や「ファイル名と範囲」を渡し、配列はワーカーの中で作るか読みます。戻り値は集計に必要な小さな量にします。どうしても共有するなら共有メモリ（`multiprocessing.shared_memory`）かメモリマップしたファイルを使います。

**JIT が効く形と効かない形。** Numba のような実行時翻訳は、数値の配列に対する単純なループ（添字の計算と四則演算）を機械語にします。効くのは、配列演算に直しにくい逐次的な漸化式や、要素ごとに条件分岐があるループです。Python のオブジェクト（辞書、文字列、任意のクラス）を扱うループには効きません。並列化の指定で和の順序が変わる点は §6 と同じです。

**GPU と精度。** CuPy の配列の既定の型は NumPy と同じ倍精度ですが、GPU 向けに書かれたコード（とくに深層学習の枠組み）は単精度 float32 を既定にしていることが多く、単精度のほうが大幅に速い GPU もあります。単精度の単位丸めは $2^{-24} \approx 6 \times 10^{-8}$ なので、必要な桁数が単精度で足りるかどうかは[精度の選び方](/computer-science/research-computing/precision-and-error-budget)の章の手順で先に見積もります。GPU の縮約演算（和、内積）は並列に行われるので、ホストの逐次和とは <Ref to="thm-parallel-sum-bound" /> の範囲で異なります。

**チャンクの大きさと負荷の偏り。** 仕事を細かく分けすぎると 1 件ごとの受け渡しのコストが効き、粗すぎると、重い仕事を引いたワーカーだけが遅れて他が待つことになります。仕事の所要時間にばらつきがあるときは、細かめのチャンクを動的に配るほうが偏りを吸収できます。ただし動的に配ると、どのワーカーがどの仕事を処理したかが実行ごとに変わります。

**結果を集める順序を固定する。** 完了した順に結果を足すと、和の順序が実行ごとに変わります。結果は仕事の番号をキーにして保存し、番号順に集計します（<Ref to="thm-parallel-sum-bound" /> の 2）。乱数の種も仕事の番号から決めます（<Ref to="ex-spawn" />）。これで「どのワーカーが、いつ処理したか」が結果に入らなくなります。数日かかる計算では、これに中断と再開が加わります。[長時間ジョブの管理](/computer-science/research-computing/long-running-jobs)の章で扱います。

## 9. 演習

<Exercise id="exr-amdahl-nine" difficulty="易">
逐次部分の割合が $f = 0.1$ の計算について、<Ref to="thm-amdahl" /> の等号が成り立つとする。(i) 速度向上 $5$ を得るのに必要な台数を求めてください。(ii) 速度向上 $10$ は有限の台数では得られないことを示してください。
<Solution>
(i) $1/(0.1 + 0.9/p) = 5$ から $0.1 + 0.9/p = 0.2$、$p = 9$ です。効率は $5/9 \approx 0.56$ です。

(ii) 有限の $p$ では $0.9/p \gt 0$ なので $0.1 + 0.9/p \gt 0.1$ で、$S(p) \lt 10$ です。$S(p) = 9.9$ でも $0.9/p = 1/9.9 - 0.1 = 0.1/99$、$p = 891$ 台が必要です。
</Solution>
</Exercise>

<Exercise id="exr-karp-flatt-data" difficulty="標準">
ある計算の速度向上を測ったところ $S(4) = 3.2$、$S(8) = 5.0$、$S(16) = 6.4$ でした。Karp–Flatt 指標 $e(p)$ を計算し、効率の低下の原因について何が言えるかを述べてください。
<Solution>
$e(4) = (1/3.2 - 1/4)/(1 - 1/4) = (0.3125 - 0.25)/0.75 \approx 0.0833$、$e(8) = (0.2 - 0.125)/0.875 \approx 0.0857$、$e(16) = (0.15625 - 0.0625)/0.9375 = 0.1$ です。

$e(p)$ は $p$ とともに増えています。所要時間を <Ref to="prop-karp-flatt" /> の形に書くと、3 からオーバーヘッドの割合は $\omega(4) \lt \omega(8) \lt \omega(16)$ です。効率の低下は逐次部分だけでは説明できず、台数とともに増えるコスト（通信、同期、スレッドの取り合い、負荷の偏り）があります。もし $e$ が $0.083$ で一定だったなら、$S(16)$ は $1/(0.083 + 0.917/16) \approx 7.1$ になっていたはずです。
</Solution>
</Exercise>

<Exercise id="exr-matvec-span" difficulty="標準">
$n \times n$ 行列とベクトルの積を、各行の内積の和を対ごとの和で計算するとします。乗算と加算をそれぞれ単位時間の作業として、仕事量 $W$ と最長の依存鎖の長さ $L$ を求め、$n = 1024$、$p = 64$ のときのステップ数の下界と上界を <Ref to="prop-work-span" /> から求めてください。
<Solution>
乗算は $n^2$ 回、加算は各行 $n - 1$ 回で $n(n-1)$ 回なので $W = 2n^2 - n$ です。乗算どうしは依存せず、各行の加算は深さ $\lceil\log_2 n\rceil$ の木（<Ref to="ex-sum-span" />）で、木の葉が乗算の結果なので、最長の有向道は乗算 1 個と加算 $\lceil\log_2 n\rceil$ 個からなり $L = 1 + \lceil\log_2 n\rceil$ です。異なる行の間に依存はありません。

$n = 1024$ では $W = 2 \times 1048576 - 1024 = 2096128$、$L = 11$ です。$W/p = 32752$（割り切れる）なので、下界は $\max(32752, 11) = 32752$、上界は $32752 + 11 = 32763$ です。速度向上は $2096128/32763 \approx 63.98$ 以上で、ほぼ台数倍です。並列度 $W/L \approx 1.9 \times 10^5$ が $p$ よりずっと大きいからです。ただしこれは演算の依存関係だけのモデルで、<Ref to="ex-intensity" /> のとおり、実際の性能はメモリ帯域で頭打ちになります。
</Solution>
</Exercise>

<Exercise id="exr-offload-matmul" difficulty="標準">
仮の数値で考えます。ホストは $n \times n$ の行列積（演算数 $2n^3$）を毎秒 $5 \times 10^{10}$ 演算の速さで計算し、加速器はその $20$ 倍の速さです。転送するのは 3 つの行列（$24n^2$ バイト）で、転送は 1 バイトあたり $10^{-10}$ 秒、固定の時間は $c_0 = 10^{-3}$ 秒とします。$n = 100$ と $n = 1000$ のそれぞれで、加速器に移すと得かどうかを <Ref to="prop-offload-breakeven" /> で判定してください。
<Solution>
$W = 2n^3/(5 \times 10^{10}) = 4 \times 10^{-11}n^3$ 秒、$s = 20$ なので左辺は $W(1 - 1/20) = 3.8 \times 10^{-11}n^3$ 秒、右辺は $D\tau + c_0 = 2.4 \times 10^{-9}n^2 + 10^{-3}$ 秒です。

$n = 100$: 左辺 $3.8 \times 10^{-5}$、右辺 $2.4 \times 10^{-5} + 10^{-3} = 1.024 \times 10^{-3}$。左辺のほうが小さく、損です。ホストで $4 \times 10^{-5}$ 秒で終わる計算に、固定の時間だけで $10^{-3}$ 秒を払うことになります。

$n = 1000$: 左辺 $3.8 \times 10^{-2}$、右辺 $2.4 \times 10^{-3} + 10^{-3} = 3.4 \times 10^{-3}$。左辺のほうが大きく、得です。所要時間は $4 \times 10^{-2}$ 秒から $2 \times 10^{-3} + 3.4 \times 10^{-3} = 5.4 \times 10^{-3}$ 秒になり、約 $7.4$ 倍です。加速器そのものは $20$ 倍速いのに、全体では $7.4$ 倍にとどまります。転送が <Ref to="thm-amdahl" /> の逐次部分の役をしています。
</Solution>
</Exercise>

## 参考文献

- G. M. Amdahl, "Validity of the single processor approach to achieving large scale computing capabilities", *AFIPS Conference Proceedings* 30 (1967), 483–485.
- J. L. Gustafson, "Reevaluating Amdahl's law", *Communications of the ACM* 31 (1988), 532–533. [doi:10.1145/42411.42415](https://doi.org/10.1145/42411.42415)
- A. H. Karp, H. P. Flatt, "Measuring parallel processor performance", *Communications of the ACM* 33 (1990), 539–543. [doi:10.1145/78607.78614](https://doi.org/10.1145/78607.78614)
- R. P. Brent, "The parallel evaluation of general arithmetic expressions", *Journal of the ACM* 21 (1974), 201–206. [doi:10.1145/321812.321815](https://doi.org/10.1145/321812.321815)
- S. Williams, A. Waterman, D. Patterson, "Roofline: an insightful visual performance model for multicore architectures", *Communications of the ACM* 52 (2009), 65–76. [doi:10.1145/1498765.1498785](https://doi.org/10.1145/1498765.1498785)
- N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — Chapter 4（和の計算の誤差解析）。
- NumPy のドキュメント「Parallel random number generation」 [https://numpy.org/doc/stable/reference/random/parallel.html](https://numpy.org/doc/stable/reference/random/parallel.html)