私はNumPyが嫌いです

• 筆者は NumPy への不満をテーマに、複数の例とともに問題点を説明している
• 簡単な配列演算は NumPy で容易に書けるが、次元が増えると複雑さと混乱が急激に増大する
• ブロードキャスト や高度なインデックス指定など、NumPy の設計は明確さと抽象化の面で不十分である
• 明示的に軸を指定する代わりに、推測と試行錯誤に頼るコードを書くことが不可欠になる
• 改善された配列言語に関するアイデアを提示し、具体的な代替案は次回の記事で紹介する予定である

---

序論: NumPyへの愛憎

• 筆者は長年 NumPy を使ってきたが、その限界に大きく失望してきたと述べている
• NumPy は Python における配列演算のための必須かつ影響力の大きいライブラリである
• PyTorch などの現代的な 機械学習ライブラリ にも、NumPy と似た問題が存在する

NumPyの簡単な点と難しい点

• 基本的な線形方程式の求解のような単純な演算は、明快でエレガントな文法 で記述できる
• しかし、配列の次元が高くなったり演算が複雑になったりすると、for ループ なしで一括処理する必要が出てくる
• ループを使えない環境（GPU 演算など）では、特殊なベクトル化文法 や特別な関数呼び出しの方法が必要になる
• しかし、こうした関数の 正確な使い方 は曖昧で、ドキュメントだけでは明確に理解しにくい
• 実際、NumPy の linalg.solve 関数は、高次元配列の場合にどう使うのが正しいのか、誰も確信を持ちにくい

NumPyの問題点

• NumPy には、多次元配列 の一部や特定の軸に演算を適用するための一貫した理論が不足している
• 配列の次元が 2 以下のときは明快だが、3 次元以上では各配列ごとにどの軸が演算対象なのかが不明確になる
• 明示的に次元を合わせるために None の使用、ブロードキャスト、np.tensordot など複雑な方法を強いられる
• こうしたやり方は、ミスを誘発し、コードの可読性を下げ、バグの可能性を高める

反復処理と明確さ

• 実際には、反復処理を許せば、より 簡潔で明確なコード を書くことができる
• 反復処理のコードは洗練されていないように見えるかもしれないが、明確さの面では大きな利点がある
• 一方で、配列の次元が変わると transpose や軸の順序を逐一考えなければならず、複雑さが増してしまう

np.einsum: 例外的に優れた関数

• np.einsum は、軸の名前を指定できる柔軟なドメイン特化言語を提供しており、非常に強力である
• einsum は演算の意図が明確で一般化もしやすく、複雑な軸演算を明示的に実装できる
• しかし、einsum に似た方式の演算サポートは 一部の演算 に限られており、たとえば linalg.solve では使えない

ブロードキャストの問題点

• NumPy の中核的なトリックである ブロードキャスト は、次元が合わないときに自動的に合わせてくれる機能である
• 単純なケースでは便利だが、実際には次元を明確に把握しにくくし、誤りの原因になることが多い
• ブロードキャストは暗黙的であるため、コードを読むたびに演算がどう働くかを毎回確認しなければならない

インデックス指定の不明確さ

• NumPy の 高度なインデックス指定 は、配列の shape を予測するのが非常に難しく、不明確である
• インデックス指定の組み合わせによって結果配列の shape が変わるため、実際に扱った経験がなければ予測が困難である
• インデックス規則の説明文書も長く複雑で、習得に大きな時間を要する
• 単純なインデックス指定だけを使いたくても、特定の演算ではどうしても高度なインデックス指定を使わざるを得ない

NumPy関数設計の限界

• 多くの NumPy 関数は、特定の配列 shape にだけ最適化されている
• 高次元配列では追加の axes 引数、別の関数名、慣習などを使う必要があり、関数ごとに一貫性がない
• 抽象化と再利用を基本とするプログラミング原則に逆行する構造である
• 特定の問題を解く関数を使っても、さまざまな配列や軸に再適用するには、まったく別のコードとして書き直さなければならない

実例: self-attention の実装

• self-attention の実装を NumPy で書くとき、反復処理を使えば明快だが、ベクトル化を強制するとコードが複雑になる
• マルチヘッド attention のように高次元演算が必要な場合、einsum と軸変換を組み合わせて使う必要があり、コードが難解になる

結論と代替案

• 筆者は、NumPy は「他の配列言語より悪い点が多いにもかかわらず、市場で重要になった唯一の選択肢」だと述べている
• NumPy のさまざまな問題点（ブロードキャスト、インデックス指定の不明確さ、関数の非一貫性など）を克服するため、改良された配列言語 のプロトタイプを作ったことを予告している
• 具体的な改善案（新しい配列言語 API）は、今後別の記事で紹介する予定である

[youn17]

Juliaがなぜ生まれたのかという話のようですね。ライブラリ群を学ぶ必要はありますが、NumPyの多くの問題を解決してくれるという点で、本当に魅力的な選択肢だと思います。

[ahwjdekf]

NumPy はベクトル化をうまく使えないと、性能は悲惨になりますよね。そういうことを考慮して書くのはストレスですし、難しいですよね。

[domino]

少し古い Python ライブラリは、どれも似たような問題を抱えている気がします

[GN⁺]

Hacker News のコメント

• 最初の例では b の型だけを見てドキュメントを読むと分かりづらいが、返される shape の説明があるので、実際に b ベクトルが行列形式なのか（特に K=1 の場合）確認する必要がある
• 配列の次元が 2 つを超えるなら、NumPy 配列に次元名を追加してくれる Xarray を使うのがおすすめ。次元合わせや transpose 作業なしでブロードキャストや整列が自動化されるため、この種の問題の多くが解消される。Xarray は線形代数の面では NumPy より弱いが、簡単に NumPy に戻れるし、補助関数を作るだけでよい。Xarray を使うと 3 次元以上のデータを扱うときの生産性が大きく上がる
  • Xarray は Pandas と NumPy の長所を合わせた感じ。da.sel(x=some_x).isel(t=-1).mean(["y", "z"]) のようなインデックス指定が簡単で、次元名が尊重されるのでブロードキャストも明確。複数の CRS を持つ地理空間データの処理に強く、Arviz との相性も抜群なのでベイズ分析で追加次元を扱うのも簡単。複数の配列を 1 つの dataset にまとめて共通の座標も共有でき、ds.isel(t=-1) のように時間軸を持つすべての配列に対して簡単に操作できる
  • Xarray のおかげで初歩的な NumPy の使い方をかなり減らせて、生産性がずっと上がった
  • Tensorflow、Keras、Pytorch のようなフレームワークにも似たものがあるのか気になる。以前、ここで触れられていた内容をデバッグするのに苦労した記憶がある
  • 紹介ありがとう、ぜひ使ってみるつもり。array[:, :, None] みたいな文法が不便だと感じていたのが自分だけではなかったと分かってうれしい
  • biosignal 分野では NeuroPype が NumPy 上で n 次元テンソル向けの名前付き軸をサポートし、各軸ごとに per-element データ（チャネル名、位置など）も保存できる
  • NumPy が昔の Numeric と Numarray ライブラリから派生していた頃を思い出す。Numarray 派が 20 年間ずっと主張を続け、資金提供を受けて Xarray に名前を変え、ついに NumPy に勝ったと想像してみる（もちろん大半は作り話）
• Julia を使い始めた理由の 1 つは、NumPy の文法があまりに難しかったから。MATLAB から NumPy に移ると、プログラミングが上達するどころか下手になったように感じ、数学よりも性能トリックを覚えることに時間を使っていた。Julia ではベクトル化もループもうまく動くので、コードの可読性だけ気にしていればいい。そうした経験や感情をこの記事からそのまま感じた。np.linalg.solve のようなものを最速だと思い込んで無条件にそれに合わせろという「ブラックボックス」的な発想は正しくないと思う。問題特化のカーネルを自分で書いた方がよい理由もいくつもある
  • 原因は、Julia が科学計算向けに設計された言語である一方、NumPy は科学計算向けではない言語の上に無理やり載せられたライブラリだからだ。いつか Julia が勝って、ネットワーク効果のせいで Python を使っている人たちが解放されることを願う
  • MATLAB もベクトル化せずにループを回すと Python と同じくらい遅い。Python の遅さこそが最大の問題で、Julia に明確な利点があるのは確かだが、実際には用途がかなり限定される。Python には JIT hack のようなものも出てきたが、まだ不完全。Python の代替が切実に必要だ
  • MATLAB は本当に違うのだろうか。ループが遅いのは変わらず、最速なのは \\ 演算子のような完全に最適化されたブラックボックスだ
  • Fortran の最新バージョンも Julia のようにベクトル化とループの両方が高速に動作するので、可読性だけに集中できる
• Matlab、Julia と比べた numpy への不満をまとめると、関数ごとに軸関連の引数や命名、ベクトル化の提供方法がばらばらで、ある軸に関数を適用したいだけでもコードを完全に書き直さなければならない点。プログラミングの基本は抽象化なのに、NumPy はそれを難しくしている。Matlab ではベクトル化コードはほぼそのまま動くか、修正点が明確だが、NumPy では毎回ドキュメントを掘り返し、transpose や reshape などで型合わせをするやり方に一貫性がなく曖昧
  • Matlab の 3 次元以上の配列サポートがあまりに弱いため、むしろ記事で触れられたような問題は起きにくい
  • 2 つ目の問題には jax の vmap を試してみる価値がある
  • 2x2 配列向けの特定関数を書いた後、それを 3x2x2 配列の一部に適用したいという話は slice や squeeze などで可能。この問題自体があまりに曖昧で理解しにくい
  • reshape で処理できる
• numpy で最も混乱する点は、どの演算がベクトル化されて動くのかが明確でなく、Julia のように dot 構文で明示することもできないこと。返り値の型に関しても落とし穴が多い。たとえば poly1d オブジェクト P を右から z0 で掛けると poly1d が返るが、左から z0*P の形で掛けると配列しか返らず、型変換が静かに起きる。quadratic の leading coefficient も P.coef[0] と P[2] の 2 通りで取得でき、混乱しやすい。公式には poly1d は「古い」API で、新しいコードには Polynomial クラスが推奨されているが、実際には deprecated 警告すら出ない。こうした型変換やデータ型の不一致のような地雷がライブラリの至る所にあり、デバッグが悪夢になる
• 著者の指摘に共感する。Matlab から Numpy に移るとき不便なことが多く、データのスライシングも Numpy の方が Matlab や Julia より不便だと感じた。ただ、Matlab の toolbox ライセンス費用を考えれば、Numpy の欠点は相殺される。記事で挙げられた問題は主に 2 次元超のテンソルで発生し、Numpy はもともと行列（2D）ベースなのでその限界は当然とも言える。Torch のような専用ライブラリの方がよいが、それも簡単ではない。結局のところ「NumPy は他のどの array 言語より少しだけひどいが、だからといって使えるものが他にあまりない」が正解という感じだ
  • Numpy は当初から N 次元配列を目標とし、numarray の延長線上にあったので、2D にとどまっていたわけではない
• Python データサイエンス生態系の最大の問題は、すべてが非標準であること。10 個あまりのライブラリが 4 つの言語分くらいばらばらに動き、せいぜい to_numpy() くらいしか統一されていない。結局、問題を解く時間よりデータ形式の変換に時間を取られる。Julia にも長所しかないわけではないが、単位や不確かさなどさまざまなライブラリ間の連携がうまく、Python は常にボイラープレートコードを大量に必要とする
  • array-api プロジェクトが Python 生態系全体で配列操作 API の標準化を進めようとしている
  • R はむしろ 4 つのクラスシステムがあるせいでさらに複雑だ
• どうしてみんな sage ではなく numpy を使うのか不思議だ
• いくつかの問題は numpysane と gnuplotlib を使えば解決する。この組み合わせができてから numpy をあらゆる作業に積極的に使うようになった。これがなければ到底使えないレベルだった
  • numpysane は結局 Python ループであって、本当のベクトル化ではない
  • 紹介ありがとう。こういう問題でたびたび愚痴っていたのに、そんな単純な上位ライブラリがあるとは思いつかなかった
• vectorized multi-head attention のために、すべての行列積を einsum に入れ、optimize="optimal" で matrix chain 積アルゴリズムを使って性能向上を試した。実際、一般的なベクトル化実装と比べて 2 倍程度速くはなったが、驚いたことにループ方式の素朴な実装の方がさらに速かった。理由が気になる人はコードを参照してほしい。einsum 内部の cache coherency にはまだ改善の余地があるのではないかと推測している