雑に読む「Numerical Recipes in C 日本語版」
雑に読む「問題解決力を鍛える！アルゴリズムとデータ構造」にインスパイアされた
雑に読むシリーズ増えてきて草
1993
動機
2022後期に数値計算のレポートがあったのだが、数値計算になれてなさすぎたせいか、何もアイデアが浮かばなかった
そもそも数値計算の解像度の粗い地図がなかった
数値計算の手法をざっくり把握しておきたい
これだと数値計算の常識も定番の予感
目標
継続
何日かに1回読む
井戸端を開いたときに必ず書き込む
ちょっと強すぎるかも？
分厚い本なので、机に座らないと読めない
1章までは写真scanしたので、スマホでも読める
読んだ日
| 2023-03-20 |
| 2023-03-22 |
コードはJSで書く
Python版も書きたい
Python苦手～><
通読だけなら、7章までざっと読んである
解像度の粗い地図はもう手に入った感はあるが……
具体的に数式やコードを書いて動かしたわけではない
少なくとも、2章は数式やコードまで全部追いたい
完了
一通りコードを書いて動かしたいが、分量が多すぎて終わりそうにない
数式の理解は飛ばしていい
分量が多すぎて終わらない
第2章あたりの計算はさっとコードを書けるようになる
数値計算コードを書くのに慣れたい
関数の計算あたりは考え方がわかればいい
いつかやる
Rustで書き直す
第１章 準備
1.0 はじめに
1.1 プログラムの構成と制御構造
知ってるので略
(説明できるとは言っていない)
1.2 科学計算のためのC言語プログラム作法
C言語の地雷のことしか書かれてないので略
1.3 誤差，正確さ，安定性
2023-02-16
ここから書いてく
int や long は固定小数点としている整数\subseteq固定小数点ということだろうか
まあ間違ってはいない
浮動小数点表現の構成要素
\mathsf{value}=s\cdot M\cdot B^{e-E}
s：符号ビット
M：仮数部(仮数)
e：指数部(指数)
B：基数(基底)
>ここでBは基数(基底)で,通常B=2であるが,B=16のこともある
B=16の処理系なんてあるのか
2023-02-26
例：float
s=\pm1
0\le M<2^{23}=8388608
32bit符号付き整数の最大値が2^{31}だから、正確に表せる数値の範囲が整数より狭いだいぶ狭いことがわかる
0\le e-E<2^8=256
10進数に換算すると何桁くらい？
2^m={10}^n
\iff m\ln2=n\ln10
\iff n=m\frac{\ln2}{\ln10}=0.30102999566\cdots\times m
\therefore B^{e-E}<2^{256}\simeq{10}^{77.06}
77桁くらいか
あれ？256のうちどのくらいが負の指数に使われるんだ？
指数の下駄Eがそれを決めているのかな？
floatの例
| s | 指数 | 仮数 | |
| 1/4 | 0 | 01111111 | 1000000000000 |
| 1/2 | 0 | 10000000 | 1000000000000 |
| 3 | 0 | 10000000 | 1100000000000 |
はい？なんで1/2=1\times2^{-1}と3=3\times2^0の指数部の値が同じなの？？
1/4=1\times2^{-2}の指数部が1/2の指数部を反転した値になっているのも謎
>浮動小数点表現では, ビットの並びが異なっても同じ数値を表すことがありうる。
>例えば B = 2 で, 仮数部 (浮動小数点表示の) の左側 (高位の桁) 0のビットがある場合は、仮数部全体を左にシフトして(つまり2の累乗を乗算して)、その分だけ指数部 (浮動小数点表示の)を減らしても,同じ数値を表す。
>この操作により仮数部の最も左側に0のビットが来ないようにすることを正規化 (normalization) という｡
このあたり、数値計算のアルゴリズムとはあまり関係ないから、飛ばして2章を先にやったほうがいい気がしてきた
2023-03-17 13:11:50 飛ばします
『精度保証付き数値計算 (現代非線形科学シリーズ)』の最初の方に、数式でわかりやすく書かれた解説があった
いずれ読む
第２章 連立１次方程式の解法
>ここから開始
2.0 はじめに
特異な連立方程式、非特異な連立方程式
連立方程式でうまく解が求まらない場合がある
理論的に求まらない場合
特異な方程式
>退化している方程式は特異(singular)であるという。
すべて、解が不定になる？
連立方程式の場合は、以下に分けられる
行退化した方程式
1つ以上の方程式が、他の方程式の線型結合で表されるもの
例
\begin{pmatrix}1&2&3\\7&-1&4\\8&1&7\end{pmatrix}\begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}-4\\5\\1\end{pmatrix}
1行目と2行目を足すと3行目になる
\begin{pmatrix}1&2&3\\7&-1&4\end{pmatrix}\begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}-4\\5\end{pmatrix}と同値
列退化した方程式
> また、どの方程式にも、いくつかの変数が全く同じ線形結合となって現れる(この条件を列退化 column degeneracyという)ならば、やはり解は一意には定まらない
言っていることがわからない
>いくつかの変数が全く同じ線形結合となって現れる
「なんの」線形結合？？
解が不能になる方程式は特異ではない
例\begin{pmatrix}1&2\\7&-1\\8&2\end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix}=\begin{pmatrix}-4\\5\\1\end{pmatrix}
後ろの解説から推測するに、これは優決定方程式と呼ぶっぽい？
計算誤差で求まらない場合
いくつかの方程式が非常に線型従属に近い場合
丸め誤差で、特異な方程式になってしまう場合がある
誤差が累積して真の解から離れる場合
未知数が非常に多いときに特に注意
連立方程式を解くpackageは、これらの病理の検出・防止のせいで複雑になってしまう
大雑把に言って、以下の場合なら単純に計算して問題ないらしい
20~50個の未知数をfloatで解く
数百個の未知数をdoubleで解く
数千個の未知数をもつ疎行列を解く
行列
記法の解説
計算線形代数の仕事
2章で扱う問題の紹介
これは解が不能にならない場合もそうなる？
例：M-N個の方程式が、それ以外の方程式の線型結合で表せるとき
この際厳密解は存在しないが、できるだけ正確に満たす値を最小二乗法で求めることはできる
右辺と左辺の差の平方和の最小化に帰着させた問題を「線型最小2乗問題」と呼ぶ
正規方程式とよぶ以下の式を解く問題に帰着する
\pmb{A}^\top\cdot\pmb{A}\cdot\pmb{x}=\pmb{A}^\top\cdot\pmb{b}
標準的なサブルーチンパッケージ
C言語の場合
linpack以外は化石ライブラリと化していそう
pythonならnumpyか
Rustは前調べてたことがある。まだデファクトスタンダードが決まっていなかったはず
JSはあるのかな
スクリプト言語で行列計算する意味を問うてはならない
2.1 Gauss-Jordan法
この辺からprogramの出番が出る
逆行列を求める最も効率のいい方法の一つらしい
メリット
\pmb{A}\cdot\pmb{C}=\pmb{B}から逆行列\pmb{A}^{-1}と方程式の解\pmb{C}を同時に求められる
algorithmが単純で理解しやすい
デメリット
>すべての右辺を同時に格納・操作しなければならない
意味がわからなかった
何と同時？
どこに格納？
\pmb{A}\cdot\pmb{C}=\pmb{B}から解\pmb{C}を求めるだけなら、これより3倍速い方法があるらしい
連立方程式を解くのにも逆行列を求めるにも、これよりLU分解のほうがいい
ただ、何か異常が発生して連立方程式を解くルーチンが怪しくなったとき、Gauss-Jordan法が予備の手法として心強い
algorithmが単純で理解しやすいから
列拡大行列の消去法
任意個の連立方程式と逆行列を一つの連立方程式に合体し、一度に解く
任意個の連立方程式\pmb{A}\cdot\pmb{x}_i=\pmb{b}_iと\pmb{A}\cdot\pmb{Y}=\pmb{I}を組み合わせて、
\pmb{A}\cdot(\bigsqcup_i\pmb{x}_i)\sqcup\pmb{Y}=(\bigsqcup_i\pmb{b}_i)\sqcup\pmb{I}
\sqcupは列拡大演算子
e.g. \begin{pmatrix}0\\3\\4\end{pmatrix}\sqcup\begin{pmatrix}5\\7\\-1\end{pmatrix}=\begin{pmatrix}0&5\\3&7\\4&-1\end{pmatrix}
を作る
これを解くと
(\bigsqcup_i\pmb{x}_i)\sqcup\pmb{Y}=\pmb{A}^{-1}\cdot(\bigsqcup_i\pmb{b}_i)\sqcup\pmb{I}
=(\bigsqcup_i\pmb{A}^{-1}\cdot\pmb{b}_i)\sqcup\pmb{A}^{-1}
となる。
ピボット選択
pivot選択なる単語が登場
しらない言葉だ
>ダンスやバスケットボールにおける片足旋回のこと
関係なさそう
まだ理解してない
CコードをTSに書き直す
うわ要素番号1始まりかよ
まあ直すの大して難しくないけど
2次元配列と行列との対応メモ
bはこの関係が逆転しているので注意
mod.tsexport type Matrix = number[][];
export type Vector = number[];
/** 完全pivot選択のGauss-Jordan法で複数の連立方程式の解と逆行列を求める
*
* @param a 連立方程式の係数
* @param bs 求めたい連立方程式の右辺
* @return 連立方程式の解とaの逆行列
export const pivotGaussJordan = (a: Matrix, ...bs: Vector[]): [Matrix, ...Vector[]] => {
// 次元チェック
const n = a.length;
if (!a.every((row) => row.length === n)) throw Error(`a must be a square matrix`);
const m = bs.length;
if (!b.every((column) => column.length === n)) throw Error(`b must be the same column length as a`);
/** pivot選択記録用 */
const indxc: number[] = [];
/** pivot選択記録用 */
const indxr: number[] = [];
/** pivot選択記録用 */
const ipiv: number[] = new Array(n).fill(0);
let big = 0;
for (let i = 0; i < n; i++) {
big = 0;
for (let j = 0; j < n; j++) {
if (ipiv[j] === 1) continue;
for (let k = 0; k < n; k++) {
}>ここまで読んでる
疑問点
掃き出し法との関係は？
2.3 LU分解
2.4 逆行列
LU分解から秒で求められる
2.5 行列式
LU分解から秒で求められる
2.6 3重対角な連立1次方程式
2.7 連立方程式の解の反復改良
2.9 特異値分解
非正方行列でもできる分解
これ使えばジョルダン標準形なんて作らなくてもいいのかな？
2.10 疎な連立1次方程式
2.11 逆行列の計算はN^3乗の過程か？
おまけ節
O(N^3)より計算量を減らすことができるが、計算が複雑になるしそこまでする必要ある？というお話
第３章 補間と補外
3.0 はじめに
3.2 有理関数による補間と補外
3.3 3次スプライン補間
3.4 整列した表の探索法
3.5 補間多項式の係数
3.6 ２次元以上の補間
第４章 関数の積分
4.0 はじめに
4.1 等間隔分点の古典公式
4.2 初等的なアルゴリズム
4.3 Romberg積分
4.4 変格積分
4.5 Gaussの求積法
4.6 多次元積分
第５章 関数の計算
5.0 はじめに
5.1 級数と収束性
5.2 連分数の計算
5.3 多項式と有理関数
5.4 漸化式と Clenshawの漸化公式
5.5 ２次方程式，３次方程式
5.6 Chebyshev近似
5.7 Chebyshev 近似した関数の微分と積分
5.8 Chebyshev係数からの多項式近似
第６章 特殊関数
6.0 はじめに
6.4 整数次の Bessel関数
6.5 整数次の変形Bessel関数
6.6 球面調和関数
6.7 楕円積分と Jacobiの楕円関数
第７章 乱数
まだメルセンヌ・ツイスターすら登場していなかったころの話
7.0 はじめに
7.1 一様乱数
7.4 ランダムなビットの生成
7.5 データ暗号化に基づく乱数列
7.6 モンテカルロ積分
第８章 ソーティング
8.0 はじめに
8.2 ヒープソート
8.3 索引づけと順位づけ
8.4 クイックソート
8.5 同値類の決定
第９章 非線形方程式と非線形連立方程式の解法
9.0 はじめに
9.4 微分を利用した Newton-Raphson 法
9.5 多項式の根
9.6 非線形連立方程式のNewton-Raphson 法
第１０章 関数の最大・最小
10.0 はじめに
10.1 １次元の黄金分割法
10.3 １階導関数を使う１次元探索
10.4 多次元の滑降シンプレックス法
10.5 多次元の方向集合
10.6 多次元の共役勾配法
10.7 多次元の可変計量法
10.8 線形計画法とシンプレックス(単体)法
10.9 アニーリング法
第１１章 固有値問題の数値計算法
11.0 はじめに
11.1 対称行列のJacobi変換
11.3 ３重対角行列の固有値，固有ベクトル
11.4 Hermite行列
11.5 一般行列のHessenberg形への変換
11.6 実Hessenberg行列のQR法
11.7 逆反復法による固有値の改良と固有ベクトルの計算
第１２章 フーリエ変換
ここまでやると、来年度の講義の予習になる
12.0 はじめに
12.1 離散データのFourier変換
12.2 高速Fourier変換
12.3 実関数の FFT，サイン・コサイン変換
12.5 FFTによる相関と自己相関
12.6 FFTによる最適(Wiener)フィルタ
12.7 FFTによるパワースペクトルの推定
12.8 最大エントロピー(全極)法によるパワースペクトル推定
12.9 時間領域でのディジタルフィルタ
12.11 ２次元以上の FFT
第１３章 データの統計的記述
13.0 はじめに
13.2 メディアンの効率的な探索
13.3 連続データの最頻値の推定
13.4 二つの分布が同じ平均値・分散を持つかどうかの検定
13.5 二つの分布は異なるのか？
13.6 2分布の分割表による解析
13.7 線形相関
13.8 ノンパラメトリック(順位)相関
13.9 データのスムージング
第１４章 データのモデル化
14.0 はじめに
14.1 最尤推定量としての最小2乗法
14.2 データの直線への当てはめ
14.3 一般の線形最小2乗法
14.4 非線形モデル
14.5 モデルパラメータの推定値の信頼限界
14.6 ロバスト推定
第１５章 常微分方程式の数値解法
15.0 はじめに
15.1 Runge-Kutta
15.3 修正中点法
15.6 硬い連立方程式
第１６章 2点境界値問題
16.0 はじめに
16.1 ねらい撃ち法
16.2 適合点へのねらい撃ち
16.3 緩和法
16.4 実例：回転楕円体調和関数
16.5 メッシュ点の自動割当て
16.6 内点での境界条件，特異点の処理
第１７章 偏微分方程式
17.0 はじめに
17.1 流束保存の初期値問題
17.2 拡散初期値問題
17.3 多次元の初期値問題
17.5 境界値問題の緩和法
付録A 参考文献
付録B プログラムの依存関係
付録C プロトタイプ宣言一覧
付録D ユーティリティルーチン
付録E 複素数演算パッケージ
索引
掲載プログラム販売のお知らせ
訳者あとがき
観客席
辞書を通読するのは自分には難しいのでどう読むのか気になる
これ辞書だったんですね
そういえばタイトルにレシピと書いてある