Python の numpy の裏では FORTRAN のライブラリが動いているらしい

平鍋さんのツイートで、Python の numpy の裏では FORTRAN で書かれたライブラリが動いているよ、という資料のリンクがあったのでリンクしておく。
参考になる。

【参考】
(1) Kenji HiranabeさんはTwitterを使っています 「Python の numpy の裏では FORTRAN で書かれた BLAS, LAPACK が現役で動いていますよ！ 行列数値計算は自分で書いてはダメ．これだけの歴史の蓄積がある．これはいい資料． https://t.co/G8UaXisxn8」 / Twitter

線形代数演算ライブラリBLASとLAPACKの基礎と実践 (I)BLAS, LAPACK入門編

興味深い点をメモしておく。

行列の使い道　→　連立一次方程式
量子コンピュータは無限次元の線形代数≒Hilbert空間論

コンピュータでの数値計算に再現性はあるか?
違う結果が出る場合がある。最近のコンピュータはマルチスレッドで、足し算の順序がコンピュータの都合で変わることがある。従って、結合法則が成り立たないことがあることより、違う結果が出る場合がある。

自作して失敗する例:連立一次方程式をクラメールの公式で解く
→理論的にはクラメールの公式を使ってプログラムを作ると、コンピュータで解ける。
→しかしながら、行列のサイズ n が大きくなるとすぐ解けなくなる
　あっというまに宇宙時間より長い時間が必要になる。

立場によって観点が違う。

・数学者の意識 : 原理的に可能, 解の存在のみ興味ある場合が多い
・情報系の数学より : アルゴリズムが多項式程度なものを考えたがる

・自然科学系研究者 : ともかく答えが求まる方法なら何でも良い。問題が出るまで放置。よく指数関数的なアルゴリズムを意識せずにゴリ押しする。

・HPC or 数値解析系研究者 : 1 clockでも速い方法 1bitでも転送量が少ない方法、1桁でも精度の良い方法などから選択する。アルゴリズム重視？

⇒数学者は解の存在だけ分かればいい。
⇒情報系学徒は、計算量オーダーが多項式程度しか興味がない。プログラムが実行できないから。
⇒自然科学者は、とにかく実際に計算できればいい。計算量オーダーの問題が出てくるまで計算ゴリ押し。

連立一次方程式をどう解くか？

・数学者の意識 : 行列式が0でないなら解ける。終わり
・情報系の数学より : LU分解も逆行列求める方法もオーダーは同じO(n^3)なので違いはない、クラメールはO(n!)なんで論外。
・自然科学系研究者 : とりあえず逆行列求めとこう。LU分解って何? 解きたい問題は別だし他に考えるべきことが多いし後まわしにしよう。

・HPC or 数値計算系研究者 : LU分解経由一択。基本であるdgemmだからキャッシュヒットも高いし、逆行列を使うより誤差も少ない。よってLU分解経由以外あり得ない。クラメールは誤差も大きい。

⇒数学者は理論1本だけ。
⇒科学者や情報系学徒はQCDのトレードオフを考える。情報系学徒は計算量オーダーが大事。

コンピュータで線形代数演算するならBLAS+LAPACKを使いましょう。
品質、信頼性がとても高く、無料で入手できる。

BLASはBasic Linear Algebra Subprogramsの略。
基礎的な線形代数の「サブ」プログラム Level 1 ～ 3まである。
FORTRAN77でさまざまなルーチンの仕様を提供している。
BLASのルーチンを「ブロック」にしてより高度なことをする。

LAPACKとはBLASをビルディングブロックとして使いつつ、より高度な問題である連立一次方程式、 最小二乗法固有値問題、固有値問題、特異値問題を解くことができる.
下請けルーチン群も提供する: 行列の分解(LU分解, コレスキー分解, QR分解, 特異値分解, Schur分解, 一般化Schur分解)、条件数の推定ルーチン, 逆行列計算など。

BLAS, LAPACKを利用したソフトウェアには、Ruby, Python (numpy), Perl, Java, C, Mathematica, Maple, Matlab, R, octave, SciLabがある。

⇒Numpyにも使われているのは面白い。