# NumPy で書く配列計算：ブロードキャスト・メモリ配置・縮約の順序 > NumPy の配列をデータ領域・shape・strides・offset の組として定式化し、ブロードキャストの結合が順序によらないこと、ビューが番地の式の変換であること、重なりのある更新、行列の鎖の最適な縮約順序、固定幅整数の剰余演算を証明つきで扱う。 > https://rikai.mugen-giken.com/computer-science/research-computing/numpy-array-computing ## 0. この記事の要点 NumPy で配列どうしを足したり、配列の一部を取り出したりするとき、中でどんなことが起きているかを知らなくても、たいていは動きます。けれども「形のちがう配列がなぜ足せるのか」「取り出した配列に書き込むと元が変わることがあるのはなぜか」「整数の計算結果がなぜ負になるのか」といった疑問には、仕組みを知らないと答えられません。この記事は、配列を「メモリ上の並びと、その読み方の規則」として書き下し、これらの疑問に規則から答えられるようにします。 - NumPy の配列は「データ領域・要素の型・shape・strides・offset」の組です。添字 $(i_1,\ldots,i_d)$ の要素の番地は $\text{offset}+\sum_ki_k\,\text{strides}_k$ で、転置・スライス・逆順・軸の追加はデータを複製せず、この式の係数を書き換えるだけです（**ビュー**）。strides は、各軸で添字を 1 つ進めると番地がいくつ増えるかを軸ごとに並べた組です。 - **ブロードキャスト**は、shape を右詰めにそろえ、各軸で「長さが等しいか一方が 1」のときに両立するという規則です。両立する shape の結合は可換・結合的・冪等（同じものどうしを結合しても変わらない）なので、何個の配列をどの順に組み合わせても結果の shape は同じです。 - 出力と入力が同じ番地を共有する更新（`a[1:] += a[:-1]`）は、要素を 1 つずつ書き換える実行と、いったん全部を計算してから書く実行で結果が違います。番地が交わらなければ、どちらでも同じです。 - 行列の鎖 $A_1\cdots A_n$ の積は、括弧づけで乗算回数が大きく変わります。最小の回数は区間についての漸化式で $O(n^3)$ 時間で求まります。一般のテンソルの縮約では、順序で演算量の**指数**が変わります。 - ベクトル化は要素ごとの演算の回数を変えません。変わるのは 1 回あたりの定数と、**一時配列の大きさ**です。全対の距離行列は一時配列が $n^2d$ 個になり、ブロック化で $nbd$ 個に抑えられます。 - 固定幅の整数の加減乗は $2^w$ を法とする演算です。真の値が表現範囲に入っていれば結果は正しく、出ていれば黙って別の数になります。 ## 1. 動機 [数値実験](/computer-science/research-computing/reproducible-experiments#def-experiment)のコードの大半は、配列に対する演算です。NumPy では、要素ごとのループを自分で書く代わりに、配列全体に対する演算を 1 行で書きます。速く、短く、読みやすくなりますが、引き換えに「この 1 行が実際に何をしているか」が見えにくくなります。 たとえば次の疑問に、規則から答えられるでしょうか。 - shape が `(4, 1, 1)`, `(2, 1)`, `(3,)` の 3 つの配列を足すと、結果の shape は何か。足す順序で変わらないか。 - `b = a.T` のあとで `b[0, 1] = 99` とすると、`a` は変わるか。`a.T.reshape(-1)` ではどうか。 - `a[1:] += a[:-1]` は累積和になるか。 - 全対の距離を 1 行で計算したら、メモリが足りなくなった。演算量は同じはずなのに、なぜか。 - `int32` の配列の内積が負になった。 どれも、配列の構造と演算の規則を定義として書き下せば、命題として答えが出ます。この記事ではそれを行います。計算量は一様コスト RAM モデル（）で数え、$O$ 記法は のものです。コード例の約束事（版と乱数種の明記）は[再現できる数値実験](/computer-science/research-computing/reproducible-experiments)に従います。 ## 2. 配列の構造：番地の式 ### 2.1. 定義 計算機の記憶装置は、番地で指定するセルの列です（、）。多次元の配列は、この 1 次元の列の上に、添字から番地への対応を決めることで作ります。 **多次元配列**（$d$ 次元）とは、次の組である。 - **データ領域**: 連続した番地 $0,1,\ldots,B-1$ のバイトの列。 - **dtype**: 1 要素のバイト数 $w$ と、$w$ バイトを値（整数・浮動小数点数など）として読む方法。 - **shape**: $0$ 以上の整数の組 $(n_1,\ldots,n_d)$。$n_k$ を第 $k$ 軸の**長さ**という。 - **strides**: 整数の組 $(s_1,\ldots,s_d)$。 - **offset**: 整数 $o$。 添字 $(i_1,\ldots,i_d)$（$0\le i_k\lt n_k$）の要素は、番地 $$ \mathrm{addr}(i_1,\ldots,i_d)=o+\sum_{k=1}^di_k\,s_k $$ から始まる $w$ バイトである（これがデータ領域に収まるものとする）。 $s_k=w\prod_{l\gt k}n_l$ のとき **C 順**（行優先。最後の添字が最も速く動く）、$s_k=w\prod_{l\lt k}n_l$ のとき **Fortran 順**（列優先）であるといい、どちらかであるとき配列は**連続**であるという（このとき要素は、番地 $o$ から $w$ バイトずつ隙間なく並ぶ）。 配列 $a$ から作った配列 $b$ が、$a$ と同じデータ領域を使い、shape・strides・offset だけが異なるとき、$b$ を $a$ の**ビュー**という。新しいデータ領域を確保して値を書き写したものを**コピー**という。 たとえば 3×4 の `int64`（$w=8$）の C 順では、$s=(8\cdot4,\,8)=(32,8)$ です。行を 1 つ進むと 1 行分の 4 要素（32 バイト）飛び、列を 1 つ進むと 1 要素（8 バイト）進みます。下の図の `a` がこの配列です。 ビューを通じて書き込むと、元の配列の対応する要素も変わります。同じ番地だからです。
0123 4567 891011 番地 032648 a[0,0] → a[0,1]：+8　　a.T[0,0] → a.T[1,0]：+8 a[0,0] → a[1,0]：+32　　a.T[0,0] → a.T[0,1]：+32
### 2.2. ビューは番地の式の変換 配列 $a$ の shape を $(n_1,\ldots,n_d)$、strides を $(s_1,\ldots,s_d)$、offset を $o$ とする。次の各操作の結果 $b$ は、$a$ と同じデータ領域の上で、右に示す shape・strides・offset をもつ配列として表せる。すなわち、$b$ の各要素の番地は、対応する $a$ の要素の番地に等しい。 1. **転置**（軸の置換 $\pi$。$b[j_1,\ldots,j_d]=a[i_1,\ldots,i_d]$, $i_{\pi(k)}=j_k$）: shape $(n_{\pi(1)},\ldots,n_{\pi(d)})$、strides $(s_{\pi(1)},\ldots,s_{\pi(d)})$、offset $o$。 2. **スライス**（第 $k$ 軸で始点 $\alpha$、歩幅 $t\ne0$、個数 $m$。$0\le j\lt m$ で $0\le\alpha+jt\lt n_k$ とし、$b[\ldots,j,\ldots]=a[\ldots,\alpha+jt,\ldots]$）: 第 $k$ 軸の長さを $m$、stride を $ts_k$ に変え、offset を $o+\alpha s_k$ にする。**逆順**は $\alpha=n_k-1$, $t=-1$, $m=n_k$ の場合である。 3. **整数による添字**（第 $k$ 軸を $\alpha$ に固定。$b$ は $d-1$ 次元）: 第 $k$ 軸を shape と strides から取り除き、offset を $o+\alpha s_k$ にする。 4. **軸の追加**（長さ $1$ の軸を挿入。`None` や `np.newaxis`）: その位置に長さ $1$、stride $0$ を挿入する。offset は $o$。 5. **長さ 1 の軸の引き伸ばし**（第 $k$ 軸が $n_k=1$ のとき、長さ $m$ にして $b[\ldots,j,\ldots]=a[\ldots,0,\ldots]$）: 第 $k$ 軸の長さを $m$、stride を $0$ にする。offset は $o$。 どれも、$b$ の番地の式（）に右の値を代入して、対応する $a$ の要素の番地と一致することを確かめます。 1. $b$ の添字 $(j_1,\ldots,j_d)$ の番地は $o+\sum_kj_k\,s_{\pi(k)}$ です。$j_k=i_{\pi(k)}$ を代入すると $o+\sum_ki_{\pi(k)}s_{\pi(k)}$ で、$\pi$ は $\{1,\ldots,d\}$ の置換なので、和の添字を $l=\pi(k)$ に取り替えて $o+\sum_li_ls_l=\mathrm{addr}_a(i_1,\ldots,i_d)$ です。 2. $b$ の番地は $(o+\alpha s_k)+j\cdot ts_k+\sum_{l\ne k}i_ls_l=o+(\alpha+jt)s_k+\sum_{l\ne k}i_ls_l$ で、これは $a$ の第 $k$ 添字を $\alpha+jt$ とした要素の番地です。逆順では $\alpha+jt=n_k-1-j$ で、$j=0,\ldots,n_k-1$ のとき $0\le n_k-1-j\lt n_k$ を満たします。 3. $b$ の番地は $(o+\alpha s_k)+\sum_{l\ne k}i_ls_l$ で、$a$ の第 $k$ 添字を $\alpha$ とした要素の番地です。 4. 挿入した軸の添字は $0$ しかとらないので、その項 $0\cdot0$ は和に寄与せず、番地は $a$ の対応する要素と同じです。 5. $b$ の番地は $o+j\cdot0+\sum_{l\ne k}i_ls_l$ で、$j$ によりません。これは $a$ の第 $k$ 添字を $0$ とした要素の番地 $o+0\cdot s_k+\sum_{l\ne k}i_ls_l$ に等しいです。 5 は次の節のブロードキャストの実装そのものです。長さ 1 の軸を長さ $m$ に「複製」するとき、NumPy はデータを書き写さず、stride を $0$ にして同じ要素を $m$ 回読みます。 3×4 の `int64` 配列（C 順、$w=8$）について、いくつかの操作の結果の shape・strides・offset を表示します。`np.shares_memory` はデータ領域を共有しているかを厳密に調べます。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64） import numpy as np a = np.arange(12, dtype=np.int64).reshape(3, 4) # C 順、1 要素 8 バイト base = a.__array_interface__["data"][0] def info(name, v): view = np.shares_memory(a, v) off = v.__array_interface__["data"][0] - base if view else "-" print(f"{name:<11} shape={v.shape!s:<9} strides={v.strides!s:<10}" f" offset={off!s:>2} {'ビュー' if view else 'コピー'}") info("a", a) info("a.T", a.T) info("a[1:, ::2]", a[1:, ::2]) info("a[::-1]", a[::-1]) info("a[:, None]", a[:, None]) info("a[2]", a[2]) info("a[[0, 2]]", a[[0, 2]]) # 整数配列による添字 info("a[a > 5]", a[a > 5]) # マスク r1 = a.reshape(12) # 連続な配列の reshape r2 = a.T.reshape(12) # 連続でないビューの reshape print(np.shares_memory(a, r1), np.shares_memory(a, r2), r2.tolist()) ``` ```text a shape=(3, 4) strides=(32, 8) offset= 0 ビュー a.T shape=(4, 3) strides=(8, 32) offset= 0 ビュー a[1:, ::2] shape=(2, 2) strides=(32, 16) offset=32 ビュー a[::-1] shape=(3, 4) strides=(-32, 8) offset=64 ビュー a[:, None] shape=(3, 1, 4) strides=(32, 0, 8) offset= 0 ビュー a[2] shape=(4,) strides=(8,) offset=64 ビュー a[[0, 2]] shape=(2, 4) strides=(32, 8) offset= - コピー a[a > 5] shape=(6,) strides=(8,) offset= - コピー True False [0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11] ``` **検算.** (a) 命題との照合: C 順の strides は $(w\cdot4,w)=(32,8)$。転置は strides の入れ替え（ の 1）、`a[1:, ::2]` は第 1 軸が $\alpha=1,t=1$、第 2 軸が $\alpha=0,t=2$ で strides $(32,16)$、offset $1\cdot32=32$（同 2）、逆順は stride $-32$、offset $(3-1)\cdot32=64$、`a[:, None]` は stride $0$ の挿入（同 4）、`a[2]` は offset $2\cdot32=64$（同 3）で、すべて出力と一致します。(b) 別の確かめ方: `np.shares_memory` の結果は、命題が「複製なし」とした 6 つでビュー、整数配列とマスクによる添字でコピーです。後者は取り出す要素の番地が等間隔とは限らないので、1 組の strides では表せません。(c) `reshape`: 連続な `a` の `reshape(12)` はビュー、連続でない `a.T` の `reshape(12)` はコピーです。理由は次の例で示します。 上の例の `a.T.reshape(12)` は、要素を $0,4,8,1,5,9,\ldots$ の順に並べた 1 次元配列です。これらの要素の番地（元のデータ領域での）は $0,32,64,8,40,72,\ldots$ です。1 次元の配列の番地は $o'+js'$（$j=0,1,2,\ldots$）という等差数列でなければなりません（）。ところが $0,32,64,8$ は、はじめの 2 つの差が $32$、3 番目と 4 番目の差が $-56$ で等差数列ではないので、どんな $(o',s')$ を選んでも同じデータ領域の上のビューとしては表せません。NumPy は新しいデータ領域に書き写します。 実務上の帰結は、`reshape` の結果に書き込んだとき、元の配列が変わるかどうかが**元の配列の strides による**ことです。書き込みを元に反映させたいのか、させたくないのかが大事な場面では、`np.shares_memory` で確かめるか、`.copy()` で明示的にコピーを作ります。 ## 3. ブロードキャスト ### 3.1. 規則 $L=\{0,1,2,\ldots\}\cup\{\top\}$ とする（$\top$ は「両立しない」を表す記号）。$a,b\in L$ の**結合** $a\vee b\in L$ を $$ a\vee b=\begin{cases}a&(a=b\ \text{または}\ b=1)\\ b&(a=1)\\ \top&(\text{それ以外})\end{cases} $$ と定める。$0$ 以上の整数 $a,b$ が $a\vee b\ne\top$ を満たす、すなわち $a=b$ または $a=1$ または $b=1$ のとき、長さ $a$ と $b$ は**両立する**という。 shape $n=(n_1,\ldots,n_d)$ と整数 $D\ge d$ に対し、先頭に $1$ を $D-d$ 個補った $(1,\ldots,1,n_1,\ldots,n_d)$ を $\mathrm{pad}_D(n)$ と書く（右詰め）。2 つの shape $n,n'$ の結合を、$D$ を両者の次元の大きい方として、$\mathrm{pad}_D(n)$ と $\mathrm{pad}_D(n')$ の成分ごとの結合と定める。結合のどの成分も $\top$ でないとき $n$ と $n'$ は**両立する**といい、結合を**ブロードキャストした shape** という。 両立する shape の配列どうしの要素ごとの演算では、それぞれの配列を、軸の追加と長さ 1 の軸の引き伸ばし（ の 4 と 5）によってブロードキャストした shape の配列とみなしてから、同じ添字の要素どうしを演算する。 長さが正のとき、両立する $a,b$ の結合は $\max(a,b)$ です。長さ $0$ の軸も許すので（空の配列）、「大きい方」ではなく「$1$ でない方」と定義しています。$0\vee1=0$ で、NumPy も同じ結果を返します。 ### 3.2. 結合は順序によらない 1. の $\vee$ は $L$ の上で可換（$a\vee b=b\vee a$）、結合的（$(a\vee b)\vee c=a\vee(b\vee c)$）、冪等（$a\vee a=a$）で、$1$ が単位元（$a\vee1=a$）である。 2. shape $n^{(1)},\ldots,n^{(m)}$ を、どんな順序・どんな括弧づけで 2 つずつ結合していっても、結果（$\top$ を含む場合も含めて）は同じである。したがって、複数の配列のブロードキャストが成功するかどうかも、成功したときの shape も、演算を適用する順序によらない。 3. 両立するという関係は推移的でない。shape $(2,)$ と $(1,)$、$(1,)$ と $(3,)$ はそれぞれ両立するが、$(2,)$ と $(3,)$ は両立しない。 次の証明のねらいは、「結合が順序によらない」ことを、2 つずつの結合の性質（可換・結合的）に帰着させることです。そのために、まず $L$ に大小関係 $\preceq$ を入れて、$a\vee b$ が $a,b$ の上限（どちらよりも小さくない元のうち最小のもの）であることを確かめます。上限は $\{a,b\}$ という集合だけで決まるので、可換性と結合性はそこから出ます。反対称律・推移律の場合分けは $\preceq$ が大小関係として振る舞うことの確認なので、細部を飛ばして「上限であること」以降から読んでも、2 の議論は追えます。 **1.** $L$ の上の関係 $\preceq$ を「$a\preceq b$ とは、$a=b$ または $a=1$ または $b=\top$」と定めます。これが半順序で、$a\vee b$ が $\{a,b\}$ の上限（最小上界）であることを示します。 *反射律*は定義から成り立ちます。*反対称律*: $a\preceq b$, $b\preceq a$, $a\ne b$ とすると、「$a=1$ または $b=\top$」かつ「$b=1$ または $a=\top$」です。4 つの組合せ $a=1=b$、$a=1,a=\top$、$b=\top,b=1$、$b=\top=a$ は、どれも $a\ne b$ か $1\ne\top$ に反します。よって $a=b$ です。*推移律*: $a\preceq b\preceq c$ とします。$a=1$ か $c=\top$ なら $a\preceq c$ です。$a=b$ なら $b\preceq c$ から、$b=c$ なら $a\preceq b$ から $a\preceq c$ です。残るのは $a\ne1$, $a\ne b$ の場合で、$a\preceq b$ から $b=\top$ です。すると $b\preceq c$ は「$\top=c$ または $\top=1$ または $c=\top$」なので $c=\top$ となり、$a\preceq c$ です。 *上限であること*: $a=b$ なら $a\vee b=a$ で、これは $\{a\}$ の上限です。$a=1$ なら $a\vee b=b$ で、$1\preceq b$, $b\preceq b$ から $b$ は上界、任意の上界 $c$ は $b\preceq c$ を満たすので $b$ が最小です。$b=1$ のときも同様に $a$ が上限です。それ以外（$a\ne b$, $a\ne1$, $b\ne1$）では、上界 $c$ は $a\preceq c$ から「$a=c$ または $c=\top$」、$b\preceq c$ から「$b=c$ または $c=\top$」を満たします。$c\ne\top$ なら $a=c=b$ となり $a\ne b$ に反するので、上界は $\top$ だけで、上限は $\top=a\vee b$ です。 上限は $\{a,b\}$ という集合で決まるので、可換性と冪等性が出ます。$a\vee1=a$ は定義です。結合性: $x=(a\vee b)\vee c$ とおきます。$x\succeq a\vee b\succeq a,b$ かつ $x\succeq c$ なので、$x$ は $\{a,b,c\}$ の上界です。$y$ を $\{a,b,c\}$ の任意の上界とすると、$y$ は $\{a,b\}$ の上界なので $y\succeq a\vee b$、また $y\succeq c$ なので、$y$ は $\{a\vee b,c\}$ の上界で $y\succeq x$ です。よって $x$ は $\{a,b,c\}$ の上限です。同じ議論で $a\vee(b\vee c)$ も $\{a,b,c\}$ の上限で、上限は反対称律からただ 1 つなので、両者は等しいです。 **2.** $D$ を $n^{(1)},\ldots,n^{(m)}$ の次元の最大値とします。2 つの shape の結合は、両者の次元の大きい方 $D'\le D$ にそろえて成分ごとにとりますが、これをさらに $\mathrm{pad}_D$ で補ったものは、はじめから $\mathrm{pad}_D$ でそろえて成分ごとに結合したものに等しいです（補った先頭の成分は $1\vee1=1$ だからです）。したがって、すべての shape を $\mathrm{pad}_D$ でそろえたうえでの成分ごとの演算として考えてよく、各成分で 1 の可換性と結合性が成り立つので、結果は順序と括弧づけによりません。 成功・失敗について。最終的な結合 $N$ がどの成分にも $\top$ を含まないとします。途中に現れる結合 $N'$ は、一部の shape の成分ごとの上限なので、成分ごとに $N'_k\preceq N_k$ です。$N_k\ne\top$ かつ $N'_k\preceq N_k$ なら、$\preceq$ の定義（$N'_k=N_k$ または $N'_k=1$）から $N'_k\ne\top$ です。よって途中のどの段階でも失敗しません。逆に $N$ が $\top$ を含むなら、どの順序で進めても最後の結果は同じ $N$ なので、どこかの段階で $\top$ が現れて失敗します。 **3.** $2\vee1=2$, $1\vee3=3$ はどちらも $\top$ でなく、$2\vee3=\top$ です。 3 は、「$a$ と $b$ を足せる、$b$ と $c$ を足せる」から「$a$ と $c$ を足せる」は出ないことを言っています。3 つを結合すると $2\vee1\vee3=\top$ で、2 のとおりどの順に足しても失敗します。 をそのまま書いた関数を、NumPy の `np.broadcast_shapes` と照合します。長さ $0$ を含む 11 個の shape から、重複を許して 2 個または 3 個を選ぶすべての組（$11^2+11^3=1452$ 通り）で比べます。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64） import itertools import numpy as np def broadcast_shape(*shapes): d = max(len(s) for s in shapes) out = [1] * d for s in shapes: s = (1,) * (d - len(s)) + tuple(s) # 右詰めにそろえる for k in range(d): if s[k] != 1: if out[k] not in (1, s[k]): raise ValueError(f"軸 {k}: {out[k]} と {s[k]} は両立しない") out[k] = s[k] return tuple(out) pool = [(), (1,), (3,), (4,), (1, 3), (2, 1), (2, 3), (4, 1, 1), (4, 2, 3), (0,), (1, 0)] ok = bad = 0 for k in (2, 3): for shapes in itertools.product(pool, repeat=k): try: mine = broadcast_shape(*shapes) except ValueError: mine = None try: ref = np.broadcast_shapes(*shapes) except ValueError: ref = None assert mine == ref, shapes ok += mine is not None bad += mine is None print("照合した組", ok + bad, "両立", ok, "両立しない", bad) print(broadcast_shape((4, 1, 1), (2, 1), (3,))) print(broadcast_shape((0,), (1,)), np.broadcast_shapes((0,), (1,))) for s in [((2,), (1,)), ((1,), (3,)), ((2,), (3,))]: try: print(s, broadcast_shape(*s)) except ValueError as e: print(s, e) ``` ```text 照合した組 1452 両立 818 両立しない 634 (4, 2, 3) (0,) (0,) ((2,), (1,)) (2,) ((1,), (3,)) (3,) ((2,), (3,)) 軸 0: 2 と 3 は両立しない ``` **検算.** (a) 別実装: 1452 通りすべてで、成功・失敗と結果の shape が `np.broadcast_shapes` と一致しました。(b) 命題との対応: 3 個の組は順序を入れ替えたものもすべて含むので、順序によらないこと（ の 2）の確認にもなっています。動機で挙げた $(4,1,1)$, $(2,1)$, $(3,)$ の結合は $(4,2,3)$ です。(c) 境界: $0\vee1=0$ で、長さ $0$ の軸は $1$ と両立して $0$ になります。最後の 3 行は命題の 3 の例です。 ## 4. 要素ごとの演算と、重なりのある更新 スカラーの関数 $f$ を、ブロードキャストした配列の同じ添字の要素どうしに適用して、同じ shape の配列を返す関数を **ufunc**（universal function）という。`np.add`, `np.multiply`, `np.exp` などで、`a + b` は `np.add(a, b)` である。 2 変数の ufunc $f$ による、第 $k$ 軸に沿った**縮約**（reduce）とは、shape $(n_1,\ldots,n_d)$ の配列 $a$ から、第 $k$ 軸を取り除いた shape の配列 $b$ を $$ b[i_1,\ldots,\widehat{i_k},\ldots,i_d]=f(\cdots f(f(a[\ldots,0,\ldots],a[\ldots,1,\ldots]),a[\ldots,2,\ldots])\cdots) $$ で作る操作である（$\widehat{\ }$ は取り除くことを表す）。`np.add` の縮約が `sum`、`np.multiply` の縮約が `prod`、`np.maximum` の縮約が `max` である。 浮動小数の和の縮約では、NumPy は上の式の括弧づけどおりではなく、対ごとの和に近い順序で足すことがあります。順序を変えたときの結果の差は、[再現できる数値実験](/computer-science/research-computing/reproducible-experiments)の で抑えられます。 ufunc の出力先を、入力と同じデータ領域のビューにすることができます（`a += b` や `out=` の指定）。このとき、出力と入力の番地が重なると、書き込みの順序が結果に影響しえます。 メモリを、番地から値への写像 $\mu$ とみる（1 要素を 1 つの番地で代表させる）。有限集合 $J$ の各 $j$ に、出力の番地 $\omega(j)$、入力の番地 $\rho_1(j),\ldots,\rho_m(j)$ が定まっていて、$\omega$ は単射とする。関数 $f$ による更新の 2 通りの実行を考える。 - **一時配列を使う実行**: すべての $j$ について $\tau_j=f(\mu(\rho_1(j)),\ldots,\mu(\rho_m(j)))$ を先に計算し、そのあとで $\mu(\omega(j))$ を $\tau_j$ に書き換える。 - **逐次実行**: $J$ の元をある順序 $j_1,j_2,\ldots$ に並べ、順に「その時点のメモリから $\rho_1(j_t),\ldots,\rho_m(j_t)$ を読み、$f$ の値を $\omega(j_t)$ に書く」を行う。 このとき次が成り立つ。 1. $j\ne j'$ ならばすべての $l$ で $\omega(j)\ne\rho_l(j')$ である（ある要素の出力先を、別の要素の計算が読まない）とき、逐次実行の結果は順序によらず、一時配列を使う実行の結果と一致する。とくに、出力の番地の集合と入力の番地の集合が交わらないとき、および出力と入力がまったく同じ配置（$\omega(j)=\rho_1(j)$ で、ほかに重なりがない）のときにそうである。 2. 長さ $n\ge3$ の配列 $a$ に対する `a[1:] += a[:-1]`（$J=\{1,\ldots,n-1\}$, $\omega(j)=\rho_1(j)=j$, $\rho_2(j)=j-1$, $f(x,y)=x+y$）は 1 の条件を満たさない。一時配列を使う実行の結果は $a'_j=a_j+a_{j-1}$、添字の小さい順の逐次実行の結果は $a'_j=a_0+a_1+\cdots+a_j$（累積和）で、$a_0\ne0$ ならば $j=2$ で異なる。 直感的には、一時配列を使う実行は「全部読んでから書く」、逐次実行は「1 つ書くたびに、次の読みがその影響を受ける」実行です。`a = [1, 1, 1, 1, 1]` に `a[1:] += a[:-1]` を適用すると、前者は各要素が元の隣の値と足されて `[1, 2, 2, 2, 2]` に、後者は書き換えたばかりの `a[1]=2` を `a[2]` に足し、`a[2]=3` を `a[3]` に足し、と進んで `[1, 2, 3, 4, 5]` になります。証明の 1 は、読む番地と書く番地が重ならなければこの差が出ないこと、2 はこの例の一般化です。 **1.** 逐次実行の $t$ 番目の操作の直前のメモリを $\mu_{t-1}$（$\mu_0=\mu$）とします。$t$ についての帰納法で、「$\mu_{t-1}$ は、番地 $\omega(j_1),\ldots,\omega(j_{t-1})$ 以外では $\mu$ に一致する」ことがわかります（書き込むのは $\omega(j_s)$ だけです）。$t$ 番目の操作が読む番地 $\rho_l(j_t)$ は、$s\lt t$ について $j_s\ne j_t$ なので、仮定から $\omega(j_s)$ のどれとも異なります。したがって読む値は $\mu(\rho_l(j_t))$ で、書く値は $f(\mu(\rho_1(j_t)),\ldots,\mu(\rho_m(j_t)))=\tau_{j_t}$ です。$\omega$ は単射なので、番地 $\omega(j_t)$ はほかの操作では書き換えられず、最終的な値は $\tau_{j_t}$ です。これは一時配列を使う実行の結果で、順序 $j_1,j_2,\ldots$ によりません。 出力と入力の番地の集合が交わらなければ、$\omega(j)\ne\rho_l(j')$ はすべての $j,j'$ で成り立ちます。同じ配置の場合は、$\omega(j)=\rho_1(j)$ で、$j\ne j'$ なら $\omega$ の単射性から $\omega(j)\ne\omega(j')=\rho_1(j')$、ほかの入力とは重ならないので、やはり条件を満たします。 **2.** $\omega(1)=1=\rho_2(2)$ なので、$j=1$, $j'=2$ が条件に反します。一時配列を使う実行は、元の値から $\tau_j=a_j+a_{j-1}$ を作って書くので $a'_j=a_j+a_{j-1}$ です。添字の小さい順の逐次実行では、$j$ についての帰納法で $a'_j=a_0+\cdots+a_j$ です。$j=1$ は $a_1+a_0$。$j-1$ まで成り立てば、$j$ 番目の操作は、すでに書き換えられた番地 $j-1$ の値 $a_0+\cdots+a_{j-1}$ と、まだ書き換えられていない番地 $j$ の値 $a_j$ を読んで足すので、$a'_j=a_0+\cdots+a_j$ です。$j=2$ での差は $(a_0+a_1+a_2)-(a_2+a_1)=a_0$ です。 NumPy の ufunc は（1.13 以降）、入力と出力の重なりを検出すると、一時配列を使う実行と同じ結果になるように計算します。したがって `a = np.array([1, 1, 1, 1, 1]); a[1:] += a[:-1]` の結果は `[1 2 2 2 2]` で、累積和 `[1 2 3 4 5]` にはなりません（手元の NumPy 2.0.2 で確認）。同じ式を Python の `for` ループで添字の小さい順に書くと累積和になります。ループをベクトル化するときは、ループが「書き換えたばかりの値を読む」形になっていないかを確かめてください。累積和そのものが欲しいなら `np.cumsum` を使います。 もう 1 つ、出力の番地が重複する場合があります。`idx = np.array([0, 0, 2])` に対する `d[idx] += 1` は、`d[idx]` を取り出し（コピー）、$1$ を足し、`d[idx]` に書き戻す、という一時配列を使う実行です。番地 $0$ には同じ値が 2 回書かれるだけなので、`d` が零ベクトルなら結果は `[1 0 1]` で、`[2 0 1]` にはなりません（ の仮定「$\omega$ は単射」が破れている場合です）。重複する添字ごとに足し込みたいときは `np.add.at(d, idx, 1)` を使います。こちらは `[2 0 1]` を返します。 ## 5. 縮約の順序 ### 5.1. Einstein の縮約記法 配列 $T^{(1)},\ldots,T^{(m)}$ の各軸に**添字の名前**（文字）を付け、同じ名前の軸は同じ長さをもつとする。出力に残す名前の列を指定したとき、**縮約**の結果は、出力の名前の値を固定し、残りのすべての名前について和をとった $$ R[\text{出力の添字}]=\sum_{\text{出力にない添字}}\ \prod_{r=1}^mT^{(r)}[\text{その配列の添字}] $$ である。`np.einsum("ij,jk->ik", A, B)` は行列の積 $R_{ik}=\sum_jA_{ij}B_{jk}$ を表す。 $m\ge3$ のとき、配列を 2 つずつ縮約して中間結果に置き換えていく手順を**縮約順序**という。2 つの配列の縮約の**費用**を、定義の式どおりに計算したときのスカラーの乗算の回数、すなわち 2 つの配列に現れるすべての名前の長さの積とする。$p\times q$ 行列と $q\times r$ 行列の積の費用は $pqr$ である。 添字に名前を付ける書き方は、軸の順序を覚えなくてよい点でも有用です。`np.einsum("bij,bjk->bik", A, B)` は「先頭の軸 `b` ごとの行列の積」であることが式から読めます。 ### 5.2. 行列の鎖の最適な括弧づけ 行列の積は結合的なので、$A_1A_2\cdots A_n$ はどう括弧を付けても同じ行列ですが、費用は括弧づけで変わります。 $A_i$ を $p_{i-1}\times p_i$ 行列（$i=1,\ldots,n$。$p_0,\ldots,p_n$ は正の整数）とする。$A_i\cdots A_j$（$i\le j$）の**括弧づけ**を、$i=j$ のときは $A_i$ そのもの、$i\lt j$ のときはある $k$（$i\le k\lt j$）と $A_i\cdots A_k$ の括弧づけ $\mathcal{L}$、$A_{k+1}\cdots A_j$ の括弧づけ $\mathcal{R}$ の組 $(\mathcal{L}\,\mathcal{R})$ と再帰的に定め、その費用を $$ \mathrm{cost}(A_i)=0,\qquad\mathrm{cost}\bigl((\mathcal{L}\,\mathcal{R})\bigr)=\mathrm{cost}(\mathcal{L})+\mathrm{cost}(\mathcal{R})+p_{i-1}p_kp_j $$ と定める。$A_i\cdots A_j$ の括弧づけの費用の最小値を $m(i,j)$ とする。 1. $m(i,i)=0$ であり、$i\lt j$ のとき $$ m(i,j)=\min_{i\le k\lt j}\bigl\{m(i,k)+m(k+1,j)+p_{i-1}p_kp_j\bigr\}. $$ 2. すべての $m(i,j)$ と、費用 $m(1,n)$ を達成する括弧づけは、$O(n^3)$ 時間・$O(n^2)$ 領域で求まる。 証明の方針は、最後の積を作る直前の分割点 $k$ で場合分けすることです。1 は、左右それぞれが最適な括弧づけであるほど全体の費用が小さいことから漸化式を導きます。2 は、区間の短いものから順に表を埋めれば、各 $m(i,j)$ の計算に必要な値がすでに表にあることから、時間と領域を見積もります。 費用の定義の第 3 項について。$A_i\cdots A_k$ は括弧づけによらず $p_{i-1}\times p_k$ 行列、$A_{k+1}\cdots A_j$ は $p_k\times p_j$ 行列なので、最後の積の費用は $p_{i-1}p_kp_j$ で、$\mathcal{L}$, $\mathcal{R}$ の選び方によりません。 **1.** $m(i,i)=0$ は定義です。$i\lt j$ とし、右辺を $M$ とおきます。 （$m(i,j)\ge M$）費用 $m(i,j)$ を達成する括弧づけを $(\mathcal{L}\,\mathcal{R})$、その分割点を $k$ とします。$\mathcal{L}$ は $A_i\cdots A_k$ の括弧づけの 1 つなので $\mathrm{cost}(\mathcal{L})\ge m(i,k)$、同様に $\mathrm{cost}(\mathcal{R})\ge m(k+1,j)$ です。よって $m(i,j)\ge m(i,k)+m(k+1,j)+p_{i-1}p_kp_j\ge M$ です。 （$m(i,j)\le M$）$i\le k\lt j$ を任意にとり、費用 $m(i,k)$ の括弧づけ $\mathcal{L}$ と費用 $m(k+1,j)$ の括弧づけ $\mathcal{R}$ をとります。$(\mathcal{L}\,\mathcal{R})$ は $A_i\cdots A_j$ の括弧づけで、費用は $m(i,k)+m(k+1,j)+p_{i-1}p_kp_j$ です。$m(i,j)$ は最小値なのでこれ以下で、$k$ は任意だったので $m(i,j)\le M$ です。 この式は、区間 $(i,j)$ の最適値が、より短い区間の最適値だけから決まることを示しています。問題の族 $\{(i,j)\}$ は最適部分構造（）をもちます。 **2.** $S=\{(i,j):1\le i\le j\le n\}$、$V(i,j)=m(i,j)$ とし、$(i,j)$（$i\lt j$）の依存先を $(i,k)$, $(k+1,j)$（$i\le k\lt j$）とします。依存先の区間の長さ $k-i$, $j-k-1$ はどちらも $j-i$ より小さいので、依存グラフの辺は区間の長さが真に増える向きだけに伸び、閉路はありません。よって $(S,D,g)$ は部分問題族（）です。$|S|=n(n+1)/2$ で、依存先の値が表にあるとき、1 の右辺の評価は $j-i\le n$ 個の候補それぞれに定数回の演算なので $c(i,j)\le Cn$（$C$ は定数）です。動的計画法の計算量の定理（）から、全体は $$ O\Bigl(|S|+\sum_{(i,j)\in S}c(i,j)\Bigr)=O\Bigl(\frac{n(n+1)}{2}+\frac{n(n+1)}{2}\cdot Cn\Bigr)=O(n^3) $$ 時間（）、$O(|S|)=O(n^2)$ 領域です。 括弧づけを得るには、1 の最小値を与える $k$ を $\kappa(i,j)$ として同じ表に記録しておきます（領域は $O(n^2)$ のまま）。$\mathrm{paren}(i,i)=A_i$、$\mathrm{paren}(i,j)=(\mathrm{paren}(i,\kappa)\,\mathrm{paren}(\kappa+1,j))$（$\kappa=\kappa(i,j)$）と再帰的に組み立てると、1 の（$\le$）の証明のとおり費用 $m(i,j)$ の括弧づけになります（区間の長さについての帰納法）。再帰の呼び出しは区間ごとに 1 回で、現れる区間は $2n-1$ 個（葉 $n$ 個と分割 $n-1$ 個）なので、$O(n)$ 時間で済みます。 6 個の行列（大きさ $40\times8$, $8\times30$, $30\times5$, $5\times60$, $60\times12$, $12\times25$）の鎖について、漸化式で最小の乗算回数を求め、`np.einsum_path` が選ぶ縮約順序と比べます。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64） import numpy as np def matrix_chain(p): # p[0..n]: A_i は p[i-1] × p[i] n = len(p) - 1 m = [[0] * (n + 1) for _ in range(n + 1)] cut = [[0] * (n + 1) for _ in range(n + 1)] for length in range(2, n + 1): # 区間の長さの小さい順 for i in range(1, n - length + 2): j = i + length - 1 m[i][j], cut[i][j] = min( (m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j], k) for k in range(i, j)) def paren(i, j): if i == j: return f"A{i}" k = cut[i][j] return f"({paren(i, k)}{paren(k + 1, j)})" return m[1][n], paren(1, n) def left_to_right(p): # ((A1 A2) A3) ... の乗算回数 return sum(p[0] * p[k] * p[k + 1] for k in range(1, len(p) - 1)) p = [40, 8, 30, 5, 60, 12, 25] best, expr = matrix_chain(p) print("最小の乗算回数", best, expr) print("左から順に掛ける", left_to_right(p)) ops = [np.empty((p[i], p[i + 1])) for i in range(len(p) - 1)] path, text = np.einsum_path("ab,bc,cd,de,ef,fg->ag", *ops, optimize="optimal") print(path) for line in text.splitlines(): if "FLOP" in line or "intermediate" in line: print(line.strip()) ``` ```text 最小の乗算回数 12900 ((A1(A2A3))((A4A5)A6)) 左から順に掛ける 68400 ['einsum_path', (1, 2), (0, 4), (0, 1), (0, 2), (0, 1)] Naive FLOP count: 5.184e+09 Optimized FLOP count: 2.580e+04 Largest intermediate: 1.000e+03 elements ``` **検算.** (a) 手計算: 得られた括弧づけの費用は、$A_2A_3$ が $8\cdot30\cdot5=1200$、$A_1(A_2A_3)$ が $40\cdot8\cdot5=1600$、$A_4A_5$ が $5\cdot60\cdot12=3600$、$(A_4A_5)A_6$ が $5\cdot12\cdot25=1500$、最後の積が $40\cdot5\cdot25=5000$ で、合計 $12900$ です。(b) 別手法: 6 個の行列の括弧づけは 42 通りで、全部を列挙して費用を計算すると最小は $12900$、最大は $168600$ でした（別のスクリプトで確認）。(c) NumPy との照合: `einsum_path` の順序は、被演算子の番号を読み解くと「$A_2A_3$、それと $A_1$、$A_4A_5$、それと $A_6$、最後に 2 つ」で、漸化式の括弧づけと同じです。NumPy の FLOP 数 $2.580\times10^4$ は乗算と加算を 1 回ずつ数えるので、乗算回数 $12900$ のちょうど 2 倍です。定義どおりの 7 重の和（Naive）は $5.184\times10^9$ で、縮約順序を選ぶだけで 5 桁減ります。 ### 5.3. 順序で指数が変わる 行列の鎖では、どの括弧づけでも中間結果は行列でした。一般のテンソルの縮約では、順序によって中間結果の軸の本数が増え、演算量の**指数**が変わります。 $n\times n$ 行列 $A,B,C$ について $R_{il}=\sum_{j,k}A_{ij}B_{jk}C_{kl}$ を、2 つずつの縮約で計算します（`"ij,jk,kl->il"`）。 - **$A$ と $B$ を先に縮約する。** $X_{ik}=\sum_jA_{ij}B_{jk}$ の費用は、現れる名前 $i,j,k$ の長さの積で $n^3$。次に $R_{il}=\sum_kX_{ik}C_{kl}$ が $n^3$。合計 $2n^3$ 回の乗算で、中間結果の要素数は $n^2$ です。 - **$A$ と $C$ を先に縮約する。** $A$ と $C$ は共通の名前をもたないので、縮約は外積 $Y_{ijkl}=A_{ij}C_{kl}$ になり、費用は $n^4$。次に $R_{il}=\sum_{j,k}Y_{ijkl}B_{jk}$ が $n^4$。合計 $2n^4$ 回の乗算で、中間結果の要素数は $n^4$ です。 $n=40$ なら、前者は $1.28\times10^5$ 回・中間結果 $1600$ 要素、後者は $5.12\times10^6$ 回・$2.56\times10^6$ 要素（binary64 で約 20 MB）です。手元の NumPy 2.0.2 の `np.einsum_path` に 2 つの順序を指定すると、中間結果の最大の要素数はそれぞれ $1.600\times10^3$, $2.560\times10^6$ と表示されます。`np.einsum` は既定では縮約順序の最適化をしない（`optimize=False`）ので、3 個以上の配列を渡すときは `optimize=True` を付けるか、自分で 2 つずつに分けて書きます。
ABX ACY → → X = AB：n² 要素（n³ 回の乗算） Y = A⊗C：n⁴ 要素（n⁴ 回の乗算）
行列の鎖では最適な順序が $O(n^3)$ 時間で求まりましたが、一般のテンソルネットワークについて費用が最小の縮約順序を求める問題は NP 困難であることが知られています（Lam, Sadayappan, Wenger の 1997 年の論文。参考文献）。実用上は、配列が少なければ全探索、多ければ貪欲法などの発見的方法が使われます（Pfeifer ほかの論文に探索法の整理があります）。この記事の命題と証明は、この事実を使いません。 ## 6. ベクトル化の費用 「ベクトル化すると速くなる」のは、要素ごとの演算の**回数**が減るからではありません。 $X$ を $n\times d$ の配列（$n$ 個の点の座標）とし、全対の距離の 2 乗 $D_{ij}=\sum_{k=1}^d(X_{ik}-X_{jk})^2$ を次の 3 通りで計算する。 - (a) **二重ループ**: $i,j$ についてのループの中で、$k=1,\ldots,d$ について差・2 乗・足し込みをスカラーで行う。 - (b) **ブロードキャスト**: `((X[:, None, :] - X[None, :, :]) ** 2).sum(axis=2)`。 - (b′) **ブロック化**: 行を $b$ 個ずつ（$1\le b\le n$）に分け、各ブロック `X[s:s+b]` について (b) と同じ式を計算する。 このとき、 1. 3 通りとも、引き算を $n^2d$ 回、乗算を $n^2d$ 回、足し算を $n^2(d-1)$ 回行う。RAM モデルでの計算時間はどれも $\Theta(n^2d)$ である。 2. 入力と出力 $D$ のほかに使う作業領域は、(a) が $O(1)$、(b) が $n^2d$ 要素の一時配列、(b′) が高々 $bnd$ 要素の一時配列である。 **1.** (a) は $(i,j,k)$ の $n^2d$ 通りについて引き算と 2 乗（乗算）を 1 回ずつ行い、各 $(i,j)$ で $d$ 個の値を足し合わせるのに $d-1$ 回足します。(b) では、`X[:, None, :]` の shape は $(n,1,d)$、`X[None, :, :]` は $(1,n,d)$ で、結合は $(n\vee1,1\vee n,d\vee d)=(n,n,d)$ です（）。引き算の ufunc はこの shape の各要素について 1 回、2 乗もそれぞれ 1 回で $n^2d$ 回ずつ、第 3 軸に沿った和の縮約は $n^2$ 個の出力それぞれに $d-1$ 回の足し算です。(b′) は、行数 $b_1,b_2,\ldots$（$\sum_rb_r=n$）のブロックごとに shape $(b_r,n,d)$ で同じ計算をするので、引き算と乗算は $\sum_rb_rnd=n^2d$ 回ずつ、足し算は $\sum_rb_rn(d-1)=n^2(d-1)$ 回です。 どの方法でも、基本演算（）の回数は、上の算術演算と、それに比例する回数の読み書き・添字の計算の和です。これは $n^2d$ の定数倍以上・定数倍以下なので $\Theta(n^2d)$ です（）。 **2.** (a) が使うのは、ループ変数と足し込み用のスカラーだけです。(b) の引き算の結果は shape $(n,n,d)$ の新しい配列で、$n^2d$ 要素です。(b′) の一時配列は shape $(b_r,n,d)$ で、$b_r\le b$ から高々 $bnd$ 要素です。ブロックの計算が終われば一時配列は不要になるので、同時に必要な量はこれで足ります。 ベクトル化で速くなるのは、演算 1 回あたりの定数が小さくなるからです。Python のループは 1 回の演算ごとにインタプリタの処理（型の判定、オブジェクトの生成）をはさみますが、ufunc はその部分をコンパイル済みのループで行います。次数は同じなので、$n$ を 2 倍にしたときの時間の伸び方は変わりません。 一方、作業領域は次数が変わります。$n=10^5$, $d=3$ の全対の差は $3\times10^{10}$ 要素、binary64 で $240$ GB で、普通の計算機には載りません。(b′) で $b=100$ とすれば $3\times10^7$ 要素（$240$ MB）で済み、演算の回数は 1 のとおり同じです。「1 行で書ける」ことと「メモリに載る」ことは別の問題です。ブロックの大きさは、キャッシュの大きさ（[計算機アーキテクチャ](/computer-science/cs-basics/computer-architecture)）にも関係しますが、ここでは立ち入りません。 全対の距離には、もう 1 つよく使われる式があります。$\|x-y\|^2=\|x\|^2+\|y\|^2-2\langle x,y\rangle$ を使って、Gram 行列 $XX^{\mathsf{T}}$ から計算する方法 (c) です。行列の積は高度に最適化されたライブラリが担うので速く、一時配列も $n^2$ 要素で済みます。ところがこの式は、近い 2 点について大きな数どうしの引き算になります。原点から遠く（座標が $1000$ 前後）、互いに近い 300 個の点で比べます。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64）／乱数種 20261008 import numpy as np rng = np.random.default_rng(20261008) n, d = 300, 3 X = 1000.0 + rng.random((n, d)) # 原点から遠く、互いに近い点 X[1] = X[0] + 1e-7 # 0 番と 1 番はとくに近い def loops(X): # (a) 二重ループ D = np.empty((n, n)) for i in range(n): for j in range(n): D[i, j] = ((X[i] - X[j]) ** 2).sum() return D def broadcast(X): # (b) ブロードキャスト（一時配列 n×n×d） return ((X[:, None, :] - X[None, :, :]) ** 2).sum(axis=2) def blocked(X, b=64): # (b') 行を b 個ずつ（一時配列 b×n×d） D = np.empty((n, n)) for s in range(0, n, b): D[s:s + b] = ((X[s:s + b, None, :] - X[None, :, :]) ** 2).sum(axis=2) return D def gram(X): # (c) |x|^2 + |y|^2 - 2 x·y sq = (X ** 2).sum(axis=1) return sq[:, None] + sq[None, :] - 2 * (X @ X.T) Da, Db, Dbb, Dc = loops(X), broadcast(X), blocked(X), gram(X) print("(a)=(b):", np.array_equal(Da, Db), " (b)=(b'):", np.array_equal(Db, Dbb)) print("(c) の負の成分の数", int((Dc < 0).sum()), " 最小値", f"{Dc.min():.3e}") print("(c) の対角の最大絶対値", f"{np.abs(np.diag(Dc)).max():.3e}") print("D[0,1]: (b)", f"{Db[0, 1]:.6e}", " (c)", f"{Dc[0, 1]:.6e}") u = 2.0 ** -53 print("|(c)-(b)| の最大", f"{np.abs(Dc - Db).max():.3e}", " u·(|x|^2+|y|^2+2|x·y|) の目安", f"{u * 4 * (X**2).sum(axis=1).max():.3e}") print("一時配列の要素数: (b)", n * n * d, " (b')", 64 * n * d) ``` ```text (a)=(b): True (b)=(b'): True (c) の負の成分の数 30 最小値 -9.313e-10 (c) の対角の最大絶対値 9.313e-10 D[0,1]: (b) 2.999998e-14 (c) 0.000000e+00 |(c)-(b)| の最大 1.737e-09 u·(|x|^2+|y|^2+2|x·y|) の目安 1.335e-09 一時配列の要素数: (b) 270000 (b') 57600 ``` **検算.** (a) 別手法: 二重ループ・ブロードキャスト・ブロック化の 3 つは、同じ演算を同じ順序で行うのでビット単位で一致しました。(b) 既知値: 0 番と 1 番の点の距離の 2 乗は $3\times(10^{-7})^2=3\times10^{-14}$ で、(b) の値はこれと 6 桁一致します（残りの違いは、$1000$ の近くで $10^{-7}$ を足したときの丸めです）。(c) の値は $0$ で、1 桁も合っていません。(c) 誤差の大きさ: (c) では $a=\|x\|^2+\|y\|^2\approx6\times10^6$ と $b=2\langle x,y\rangle\approx6\times10^6$ の差をとります。それぞれが相対 $u$ 程度の丸め誤差を含むので、桁落ちの評価（）から、差の絶対誤差は $u(|a|+|b|)\approx1.3\times10^{-9}$ 程度になりえます。実際の (b) との差の最大 $1.7\times10^{-9}$ はこの大きさで、距離の 2 乗が $10^{-9}$ より小さい対では符号すら信用できず、負の値が 30 個出ています。このまま平方根をとると `nan` になります。 対策は、引き算の前に誤差の入らない形にすることです。差を直接とる (b) か (b′) を使うか、(c) を使うなら先に重心を引いて $\|x\|$ を小さくします。負の値を $0$ に切り上げるだけでは、符号が直るだけで桁は戻りません。 ## 7. 固定幅の整数 NumPy の整数は Python の `int`（桁数に制限がない）と違い、決まったビット数で表されます。 $w\ge1$ を整数、$I_w=\{-2^{w-1},\ldots,2^{w-1}-1\}$ とする。整数 $v$ に対し、$v\equiv r\pmod{2^w}$（）を満たす $r\in I_w$ を $\mathrm{wrap}_w(v)$ と書く。**符号つき $w$ ビット整数**の加算・減算・乗算を、$a,b\in I_w$ に対し $$ a\oplus b=\mathrm{wrap}_w(a+b),\qquad a\ominus b=\mathrm{wrap}_w(a-b),\qquad a\otimes b=\mathrm{wrap}_w(ab) $$ と定める。`int32` は $w=32$、`int64` は $w=64$ である。 $\mathrm{wrap}_w(v)$ がただ 1 つ定まることは、次の命題の 1 で示します。$I_w$ は剰余環 $\mathbb{Z}/2^w\mathbb{Z}$（）の各剰余類から代表を 1 つずつ選んだもので、$\oplus,\otimes$ はこの環の演算を代表元で書いたものです。 1. 任意の整数 $v$ に対し、$v\equiv r\pmod{2^w}$ となる $r\in I_w$ がただ 1 つ存在する。$\mathrm{wrap}_w(v)=v$ となるのは $v\in I_w$ のとき、そのときに限る。 2. $a,b\in I_w$ について、$a\oplus b=a+b$ となるのは真の和 $a+b$ が $I_w$ に属するとき、そのときに限る。$\ominus$ と $a-b$、$\otimes$ と $ab$ についても同じである。 3. $I_w$ の元を変数に代入し、$+,-,\times$ だけからなる式を $\oplus,\ominus,\otimes$ で計算した結果は、同じ式を整数として計算した真の値 $v$ に対する $\mathrm{wrap}_w(v)$ に等しい。とくに、最終的な真の値が $I_w$ に属していれば、途中の値が $I_w$ を出ていても、結果は正しい。 **1.** 存在: $q=\lfloor(v+2^{w-1})/2^w\rfloor$ とおくと $q\le(v+2^{w-1})/2^w\lt q+1$ なので、$r=v-2^wq$ は $-2^{w-1}\le r\lt2^{w-1}$ を満たし、$r\in I_w$、$2^w\mid(v-r)$ です。一意性: $r,r'\in I_w$ がともに $v$ と合同なら、$2^w\mid(r-r')$ です。一方 $|r-r'|\le(2^{w-1}-1)-(-2^{w-1})=2^w-1\lt2^w$ なので $r-r'=0$ です。後半: $v\in I_w$ なら $v$ 自身が条件を満たす $I_w$ の元なので、一意性から $\mathrm{wrap}_w(v)=v$ です。逆に $\mathrm{wrap}_w(v)=v$ なら、$\mathrm{wrap}_w(v)\in I_w$ から $v\in I_w$ です。 **2.** $a\oplus b=\mathrm{wrap}_w(a+b)$ に 1 の後半を $v=a+b$ として適用します。$\ominus$, $\otimes$ も $v=a-b$, $v=ab$ として同じです。 **3.** 式の構造についての帰納法で、「計算結果 $\hat v$ は真の値 $v$ と $2^w$ を法として合同で、$\hat v\in I_w$」を示します。変数では $\hat v=v\in I_w$ です。式が $E_1\circ E_2$（$\circ\in\{+,-,\times\}$）で、部分式の計算結果と真の値が $\hat v_k\equiv v_k$（$k=1,2$）を満たすとします。合同式の四則（）から $\hat v_1\circ\hat v_2\equiv v_1\circ v_2\pmod{2^w}$ です。計算結果は $\hat v=\mathrm{wrap}_w(\hat v_1\circ\hat v_2)$ で、定義から $\hat v\in I_w$ かつ $\hat v\equiv\hat v_1\circ\hat v_2\equiv v_1\circ v_2=v$ です。式全体について $\hat v\in I_w$, $\hat v\equiv v$ なので、1 の一意性から $\hat v=\mathrm{wrap}_w(v)$ です。$v\in I_w$ なら 1 の後半から $\hat v=v$ です。 2 が言うのは、あふれは**真の値が範囲を出たときにだけ**起き、そのとき結果は誤差を含む近似値ではなく、$2^w$ の整数倍だけずれたまったく別の数になる、ということです。浮動小数点数のあふれは `inf` になって目に見えますが、整数のあふれは見た目がもっともらしい数になるので気づきにくくなります。 $20!=2432902008176640000\lt2^{63}=9223372036854775808\lt21!=51090942171709440000$ なので、`int64` の階乗は $21!$ で初めて $I_{64}$ を出ます。 ```python # Python 3.9.6 / NumPy 2.0.2（arm64） import math import numpy as np def wrap(v, w): # 2^w を法とする代表元を [-2^(w-1), 2^(w-1)) にとる return (v + 2 ** (w - 1)) % 2 ** w - 2 ** (w - 1) f = np.cumprod(np.arange(1, 23, dtype=np.int64)) # 1!, 2!, ..., 22! for k in (20, 21, 22): true = math.factorial(k) print(k, int(f[k - 1]), int(f[k - 1]) == true, int(f[k - 1]) == wrap(true, 64)) a = np.full(3, 2 * 10**9, dtype=np.int32) # 真の和は 6·10^9 > 2^31 - 1 print(a.sum(dtype=np.int32), wrap(6 * 10**9, 32)) # int32 のまま足す print(a.sum(), a.sum().dtype) # 既定では int64 に上げて足す print(a + a, (a + a).dtype) # 要素ごとの和は int32 のまま c = np.full(3, 50000, dtype=np.int32) # 真の内積は 7.5·10^9 print(c @ c, wrap(75 * 10**8, 32)) # 内積も int32 のまま print(c.astype(np.int64) @ c.astype(np.int64)) # 先に int64 にする b = np.array([2**31 - 1, 1, -5], dtype=np.int32) # 途中であふれるが真の和は範囲内 print(b.sum(dtype=np.int32), 2**31 - 1 + 1 - 5) ``` ```text 20 2432902008176640000 True True 21 -4249290049419214848 False True 22 -1250660718674968576 False True 1705032704 1705032704 6000000000 int64 [-294967296 -294967296 -294967296] int32 -1089934592 -1089934592 7500000000 2147483643 2147483643 ``` **検算.** (a) 既知値: `math.factorial` は Python の整数で計算した真の値です。$20!$ は一致し、$21!$, $22!$ は一致しません。(b) 命題との照合: あふれた値は、真の値を $2^{64}$ または $2^{32}$ を法として $I_w$ に戻した `wrap` の値と、階乗・和・内積のすべてで一致します（ の 3）。最後の行は、途中の和 $2^{31}$ が $I_{32}$ を出ても、真の和 $2147483643$ が範囲内なので正しい結果になる例です。(c) 既定の型: この環境の `np.sum` は、`int32` の配列を既定で `int64` に上げてから足すので $6\times10^9$ を正しく返します。しかし要素ごとの和 `a + a` と内積 `c @ c` は `int32` のまま計算し、警告なしに負の数を返します。和・積・内積の結果が $2^{31}\approx2.1\times10^9$ を超えうるときは、先に `astype(np.int64)` で広げるか、`dtype=np.int64` を指定します。`int64` でも $2^{63}\approx9.2\times10^{18}$ を超えれば同じことが起きるので、階乗や二項係数のように急に大きくなる量は Python の `int` か浮動小数点数で計算します。 ## 8. 実務の注意 ここまでの命題にもとづく注意を、型の昇格とあわせてまとめます。 | 落とし穴 | 何が起きるか | 対処 | |---|---|---| | `float32` の混入 | `float32` の配列だけからなる演算の結果は `float32` のままで、単位丸めは $2^{-24}\approx6\times10^{-8}$。ファイルから読んだ配列が `float32` だと、以後の計算の精度が 7 桁になる | 読み込んだ直後に `dtype` を確かめ、必要なら `astype(np.float64)` | | 型の昇格の規則 | 手元の NumPy 2.0.2 では、`float32` の配列に Python の `2.0` を掛けると `float32`、`np.float64(2.0)` を掛けると `float64`、`int32` や `int64` の配列と `float32` の配列の和は `float64` になる（NumPy 2 の規則。1.x とは一部異なる） | 結果の `dtype` を仮定せず確かめる。版を記録する | | ビューを通じた書き換え | スライス・転置はビューなので、書き込むと元の配列が変わる（） | 独立な配列が要るなら `.copy()` | | マスク・整数配列による添字 | コピーを返す（）。`a[mask][0] = 1` は一時的なコピーに書くだけで `a` は変わらない。`a[mask] = 1` は `a` に書く | 代入は 1 回の添字で書く | | 重複する添字への足し込み | `d[idx] += 1` は重複を 1 回しか数えない（4 節） | `np.add.at` か `np.bincount` | | 書き換えたばかりの値を読むループのベクトル化 | 結果が変わる（） | `np.cumsum` などの専用の関数を使う | | 一時配列 | ブロードキャストで $n^2d$ 要素の配列ができる（） | ブロック化する | | メモリ配置と走査の向き | C 順の配列を列方向に走査すると、番地が $s_1$ ずつ飛び、キャッシュに載りにくい | 最後の軸を内側のループにする。必要なら `np.ascontiguousarray` | 配列の型・shape・連続性は、コードの中で仮定している箇所に `assert a.dtype == np.float64` のように書いておくと、仮定が破れたときにその場で止まります。数値の精度をどう選ぶかは、次の[精度の選び方](/computer-science/research-computing/precision-and-error-budget)で扱います。 ## 9. 演習 (1) shape が $(5,1,4)$, $(3,1)$, $()$（0 次元）の 3 つの配列の和の shape を、 に従って求めよ。(2) shape $(2,3)$ と $(3,2)$ の配列は足せるか。足せない場合、片方をどう変形すれば足せるか。 (1) 3 次元にそろえると $(5,1,4)$, $(1,3,1)$, $(1,1,1)$ です。成分ごとの結合は $5\vee1\vee1=5$、$1\vee3\vee1=3$、$4\vee1\vee1=4$ で、$(5,3,4)$ です。 (2) 第 1 成分が $2\vee3=\top$ なので両立しません。$(3,2)$ の配列を転置すれば $(2,3)$ になり、足せます（要素の対応が意図どおりかは別に確かめます）。あるいは、すべての組合せの和（shape $(2,3,3,2)$）が欲しいのなら、`a[:, :, None, None] + b[None, None, :, :]` のように長さ 1 の軸を足します。$(2,3,1,1)$ と $(1,1,3,2)$ の結合は $(2,3,3,2)$ です。 `a = np.arange(24, dtype=np.int32).reshape(2, 3, 4)`（C 順）とする。`b = a[:, ::-1, 1::2].transpose(2, 0, 1)` の shape・strides・offset を から求め、`b[1, 0, 2]` の値を番地の式から計算せよ。 $w=4$ で、`a` の strides は $(4\cdot3\cdot4,\ 4\cdot4,\ 4)=(48,16,4)$、offset は $0$ です。 第 2 軸の逆順（$\alpha=2$, $t=-1$, $m=3$）で stride が $-16$、offset が $2\cdot16=32$ になります。第 3 軸の `1::2`（$\alpha=1$, $t=2$, $m=2$）で stride が $8$、offset が $32+1\cdot4=36$ になります。ここまでで shape $(2,3,2)$、strides $(48,-16,8)$、offset $36$ です。`transpose(2, 0, 1)` は、新しい第 1, 2, 3 軸を元の第 3, 1, 2 軸にするので、shape $(2,2,3)$、strides $(8,48,-16)$、offset $36$ です。 `b[1, 0, 2]` の番地は $36+1\cdot8+0\cdot48+2\cdot(-16)=12$ です。`a` は `int32` の C 順で、番地 $12$ は先頭から $12/4=3$ 番目の要素、値は $3$ です。確かめると、`b[1, 0, 2]` はスライス後の配列の添字 $(0,2,1)$、すなわち `a[0, 2-2, 1+2·1] = a[0, 0, 3] = 3` です。NumPy でも同じ shape・strides・値が得られます。 $A_1,A_2,A_3,A_4$ の大きさを $20\times2$, $2\times30$, $30\times4$, $4\times10$ とする。 の漸化式で $m(1,4)$ と最適な括弧づけを求め、左から順に掛けた場合の乗算回数と比べよ。 $p=(20,2,30,4,10)$ です。長さ 2 の区間: $m(1,2)=20\cdot2\cdot30=1200$、$m(2,3)=2\cdot30\cdot4=240$、$m(3,4)=30\cdot4\cdot10=1200$。 長さ 3 の区間: $m(1,3)=\min\{0+240+20\cdot2\cdot4,\ 1200+0+20\cdot30\cdot4\}=\min\{400,3600\}=400$（$k=1$）。$m(2,4)=\min\{0+1200+2\cdot30\cdot10,\ 240+0+2\cdot4\cdot10\}=\min\{1800,320\}=320$（$k=3$）。 長さ 4: $m(1,4)=\min\{0+320+20\cdot2\cdot10,\ 1200+1200+20\cdot30\cdot10,\ 400+0+20\cdot4\cdot10\}=\min\{720,8400,1200\}=720$（$k=1$）。 最適な括弧づけは $A_1((A_2A_3)A_4)$ で $720$ 回です。左から順に掛けると $1200+20\cdot30\cdot4+20\cdot4\cdot10=1200+2400+800=4400$ 回で、約 6 倍です。細い行列 $A_1$（$20\times2$）を最後に掛けると、中間結果が小さいまま進むからです。 `int8`（$w=8$, $I_8=\{-128,\ldots,127\}$）で次を計算した結果を から求めよ。(1) $100\oplus100$。(2) $(100\oplus100)\oplus(-120)$。(3) 絶対値が $2^{31}$ 未満の整数 $N$ 個の和を `int64` で計算するとき、あふれが起きないことを保証できる $N$ の範囲。 (1) 真の和 $200$ は $I_8$ の外です。$200-256=-56\in I_8$ なので $100\oplus100=-56$ です。 (2) 命題の 3 から、結果は真の値 $100+100-120=80$ に対する $\mathrm{wrap}_8(80)$ で、$80\in I_8$ なので $80$ です。途中の $-56$ は誤っていますが、最終結果は正しくなります。実際 $-56\oplus(-120)=\mathrm{wrap}_8(-176)=-176+256=80$ です。 (3) 和の絶対値は $N\cdot2^{31}$ 未満です。$N\cdot2^{31}\le2^{63}$、すなわち $N\le2^{32}\approx4.3\times10^9$ なら、和は $-2^{63}$ より大きく $2^{63}$ より小さいので $I_{64}$ に属し、命題の 3 から（途中の部分和も同じ評価で範囲内ですが、仮に出たとしても）結果は正しいです。 ## 参考文献 - C. R. Harris ほか, "Array programming with NumPy", *Nature* 585 (2020), 357–362. [doi:10.1038/s41586-020-2649-2](https://doi.org/10.1038/s41586-020-2649-2) — 配列の構造（データ・shape・strides）、ビュー、ブロードキャストの概説。 - S. van der Walt, S. C. Colbert, G. Varoquaux, "The NumPy array: a structure for efficient numerical computation", *Computing in Science & Engineering* 13 (2011), 22–30. [doi:10.1109/MCSE.2011.37](https://doi.org/10.1109/MCSE.2011.37) - NumPy 公式ドキュメント, [Broadcasting](https://numpy.org/doc/stable/user/basics.broadcasting.html)、[Copies and views](https://numpy.org/doc/stable/user/basics.copies.html)、[NEP 50 — Promotion rules for Python scalars](https://numpy.org/neps/nep-0050-scalar-promotion.html)。 - T. H. Cormen, C. E. Leiserson, R. L. Rivest, C. Stein, *Introduction to Algorithms*, 3rd ed., MIT Press, 2009 — 動的計画法の章（行列の鎖の積）。 - C.-C. Lam, P. Sadayappan, R. Wenger, "On optimizing a class of multi-dimensional loops with reductions for parallel execution", *Parallel Processing Letters* 7 (1997) — 縮約順序の最適化の計算困難性。 - R. N. C. Pfeifer, J. Haegeman, F. Verstraete, "Faster identification of optimal contraction sequences for tensor networks", *Physical Review E* 90 (2014), 033315. [doi:10.1103/PhysRevE.90.033315](https://doi.org/10.1103/PhysRevE.90.033315)