計算科学のための大規模数値解法

大規模数値計算では、同じ「方程式を解く」でも、支配要因が小規模計算とは質的に変化する。計算量のオーダーだけでなく、メモリ階層やデータ移動、通信コストが性能を決定し、さらに丸め誤差や悪条件性が解の信頼性を左右する。本書は、学部3年生が研究・卒業研究に入る前段階で、数値計算が「なぜ速くならないのか」「なぜ壊れるのか」を説明できる土台を与えることを目的とする。

全体は、線形方程式を軸に、疎行列を「データ構造＋グラフ」として捉える視点から始め、直接法と反復法の役割分担、Krylov部分空間法（CG/GMRES/BiCGStab）と前処理の設計思想へ進む。さらに、固有値問題、非線形方程式（Newton-Krylov）、大規模最適化（勾配・準Newton・近接・確率的手法）、PDE離散化（FDM/FEM）を一つの流れに統合し、「多くの科学技術計算は疎な線形系を繰り返し解く問題である」という共通構造を繰り返し確認する。

後半では、時間発展問題における安定性・剛性・陰解法を扱い、「陰解法＝毎ステップ線形（非線形）ソルバが必要」という実務上の核心を押さえる。最後に、並列化・スケーリング・プロファイル・再現性・テスト・ログ設計といった研究コードの作法を整理し、学生が“動くプログラム”を“使える研究ソフトウェア”へ引き上げる道筋を示す。各章には、手計算での見積もりと、最小限の実装で体験できる演習を配置し、理解を概念で終わらせない構成とする。

第1章　大規模数値計算とは何か：問題分類と計算資源

狙い：計算規模が増えると何が支配的になるか（計算量・メモリ・通信）を掴む

1.0 この章で学ぶこと

「大規模数値計算」とは、単に N（データ数・未知数の数）が大きいだけでなく、

計算時間が現実的でなくなる
メモリに乗らない（そもそも保持できない）
並列化しても通信がボトルネックになって速くならない

といった計算資源の制約が支配的になる状況を指す。

本章の到達目標は次の3点である。

数値問題を「線形・非線形・最適化・PDE」などに分類し、典型的に何が重くなるかを説明できる
計算量（時間）とメモリ量を、オーダー記法（O( )）で見積もれる
大規模化で支配的になる「メモリ階層」と「通信（データ移動）」を、簡単なモデル式で説明できる

1.1 「大規模」とは何か：3つのボトルネック

大規模計算で問題になるのは、主に次の3つである。

(A) 計算量（Compute）

「演算回数が多すぎて時間がかかる」という問題である。
例：密行列のLU分解は概ね O(N^3) の演算が必要になる。

(B) メモリ（Memory）

「必要な配列や行列をメモリに保持できない」という問題である。
例：N×N の密行列を倍精度（8 byte）で持つだけで 8N^2 byte 必要になる。

(C) 通信・データ移動（Communication / Data movement）

並列計算やGPU計算では「演算そのもの」よりも、データを動かすコストが支配的になりやすい。
特に、多数ノードの並列では、ネットワーク通信（MPIなど）が効いてくる。

重要：近年の高性能計算（HPC）やGPU計算では、演算性能（FLOPS）が伸びる一方で、
メモリ帯域や通信は相対的に制約になりやすく、
「どれだけ速く計算できるか」は「どれだけ効率よくデータを運べるか」に強く依存する。

1.2 数値計算で扱う「問題」の代表的な型

本講義では、数値問題を次の4つに大別する。

1.2.1 線形問題（連立一次方程式）

最も基本で、かつ大規模化の中心になる問題が

A x = b

である。

$A$ ：行列（サイズ $n \times n$ ）
$x$ ：未知ベクトル（サイズ $n$ ）
$b$ ：既知ベクトル（サイズ $n$ ）

多くの問題（PDEの離散化、最小二乗、非線形のNewton法など）は、最終的に「何度も $A x = b$ を解く」形に帰着する。

1.2.2 非線形問題（方程式 $f (x) = 0$ ）

f (x) = 0

の形で、未知 $x$ を求める。代表例はNewton法であり、各反復で

J (x_{k}) Δ x = - f (x_{k}), x_{k + 1} = x_{k} + Δ x

を解く（ここで $J$ はヤコビアン）。
つまり非線形でも結局「大規模線形問題を繰り返し解く」ことが多い。

1.2.3 最適化問題（目的関数の最小化）

min_{x} ϕ (x)

の形。材料設計、逆問題、機械学習などで頻出する。
勾配法などでは、反復ごとに勾配 $\nabla ϕ (x)$ を評価し、場合によってはヘッセ行列（2階微分）に関わる線形系も出てくる。

1.2.4 PDE（偏微分方程式）の数値解

例としてPoisson方程式：

- \nabla^{2} u = f

を格子やメッシュで離散化すると、巨大な連立一次方程式

A u = b

が得られる。PDE → 離散化 → 疎行列の線形系、という流れは大規模計算の王道である。

1.3 オーダー記法（O( )）で何を見積もるのか

大規模化の見積もりでは、係数の細部よりまず「増え方（スケーリング）」を見る。

1.3.1 時間計算量（Time complexity）

入力サイズを $n$ とするとき、実行時間が

$O (n)$ ：線形に増える（理想的）
$O (n^{2})$ ：二乗で増える
$O (n^{3})$ ：三乗で増える（大規模では破綻しやすい）

などで表される。

例：

ベクトルの内積： $O (n)$
密行列×ベクトル： $O (n^{2})$
密行列のLU分解：おおむね $O (n^{3})$

1.3.2 メモリ計算量（Memory complexity）

保存すべきデータ量の増え方を見る。

例：

密行列 $n \times n$ （倍精度）：必要メモリは $8 n^{2} byte$
ベクトル $n$ ： $8 n byte$

大規模計算では「計算時間」以前に「そもそもメモリに乗らない」ことが頻繁に起きる。

1.4 密行列と疎行列：大規模計算の分岐点

1.4.1 疎行列（sparse matrix）とは

行列 $A$ の成分の多くがゼロで、ゼロでない要素（nonzero）の数が少ない行列を疎行列という。
疎行列では、ゼロ要素を全部保存するのは無駄なので、非ゼロ要素だけを保存する形式（CSR/CSCなど）を使う。

$n$ ：行列の次元
$nnz$ ：非ゼロ要素数（number of non-zeros）

密行列： $nnz \approx n^{2}$
疎行列： $nnz ≪ n^{2}$

1.4.2 保存量の比較（密 vs 疎）

密行列（倍精度）

memory \approx 8 n^{2} byte

CSR形式（概算）

CSRは「値配列」「列番号配列」「行ポインタ配列」を持つ。概算として（列番号などを4 byte整数と仮定すると）

memory \approx (8 + 4) nnz + 4 (n + 1) byte

となる（実装によって8 byte indexの場合もあり、ここは環境依存である）。

重要：疎行列は「保存」だけでなく、行列×ベクトル（SpMV）の計算量も
$O (nnz)$
に落ちることが多い。ここが大規模で決定的に効く。

1.4.3 具体例：もし $n = 10^{6}$ なら？

密行列：
$8 n^{2} = 8 \times (10^{6})^{2} = 8 \times 10^{12} byte \approx 8 TB$
これは普通の計算機どころか、多くの計算ノードでも単体では厳しい。
疎行列：仮に1行あたり平均7個だけ非ゼロ（ $nnz \approx 7 n$ ）なら、
$nnz \approx 7 \times 10^{6}$
CSR概算で
$(12) nnz \approx 12 \times 7 \times 10^{6} \approx 84 \times 10^{6} byte \approx 84 MB$
となり、現実的なサイズになる。

このように、「疎性があるかどうか」で問題の難易度が桁違いに変わる。

1.5 メモリ階層：なぜ“データを動かす”のが遅いのか

計算機は一様なメモリではなく、一般に次のような階層を持つ（概念図）。

レジスタ：最速だが極小
キャッシュ（L1/L2/L3）：速いが容量は小さい
主記憶（DRAM）：容量は大きいが遅い
ストレージ（SSD/HDD）：さらに遅い（桁違いに遅い）

この階層があるため、同じ演算でも

データがキャッシュにうまく乗る → 速い
毎回DRAMから読み込む → 遅いとなる。

したがって大規模計算では「演算回数」だけでなく
「メモリ帯域（どれだけ速くデータを供給できるか）」が支配的になる。

1.6 “計算が遅い”の正体：Rooflineの考え方（直観）

大規模計算では、しばしば「CPU/GPUは速いはずなのに、実効性能が出ない」ことが起きる。
その直観を与えるのが Rooflineモデルである。

1.6.1 算術強度（Arithmetic Intensity）

「読み書きするデータ量に対して、どれだけ演算するか」を

I = \frac{FLOPs}{Bytes}

で定義する（算術強度）。

$I$ が小さい：少し計算してすぐ次のデータを取りに行く → メモリ帯域律速（memory-bound）
$I$ が大きい：データを一度取ったらたくさん計算する → 演算律速（compute-bound）

1.6.2 性能の上限（概念式）

Rooflineの基本は

P \leq min (P_{peak}, B_{mem} \times I)

である。

$P$ ：実効性能（FLOPs/s）
$P_{peak}$ ：ピーク演算性能
$B_{mem}$ ：メモリ帯域（Bytes/s）
$I$ ：算術強度（FLOPs/Bytes）

つまり、算術強度が低い計算（例：疎行列×ベクトル）はしばしばメモリ帯域で頭打ちになる。

1.7 並列計算と通信：強スケーリング・弱スケーリング

大規模計算では並列化が必須だが、並列化しても必ず速くなるわけではない。
その理由を、まず概念として整理する。

1.7.1 通信時間の単純モデル

通信時間はしばしば

T_{comm} \approx α + \frac{m}{β}

と書ける（概念モデル）。

$α$ ：レイテンシ（通信開始の固定コスト）
$β$ ：帯域（Bytes/s）
$m$ ：送るデータ量（Bytes）

小さいデータを何回も送ると $α$ が効いて遅い。
大きいデータを送ると $m / β$ が効いて遅い。
つまり「通信はタダではない」。

1.7.2 強スケーリング（Strong scaling）とAmdahlの法則

問題サイズ固定でプロセッサ数を増やし、どこまで速くなるかを見るのが強スケーリングである。
Amdahlの法則（概念）：

S (P) = \frac{1}{s + \frac{1 - s}{P}}

$S (P)$ ：P並列のスピードアップ
$s$ ：並列化できない割合（逐次割合）

逐次割合 $s$ が少しでも残ると、Pを増やしても頭打ちになる。

1.7.3 弱スケーリング（Weak scaling）とGustafsonの法則

プロセッサ数に比例して問題サイズも増やす見方が弱スケーリングである。
Gustafsonの法則（概念）：

S (P) = P - (P - 1) s

実務では「計算資源が増えたら、より大きい問題を解く」ことが多く、弱スケーリングの視点が重要になる。

1.8 本章のまとめ：大規模計算の“設計図”

大規模数値計算を考えるときは、次の順序で整理するとギャップが少ない。

問題の型を決める（線形／非線形／最適化／PDE）
支配的データ構造を見極める（密か疎か、 $nnz$ はどの程度か）
時間とメモリをオーダーで見積もる（O( )）
“データ移動”が支配的かを疑う（メモリ帯域・通信）
強/弱スケーリングのどちらを目標にするかを決める

この見取り図ができると、第2章以降で学ぶ「数値解法の選択」が合理的にできるようになる。

演習

演習1：密行列 vs 疎行列のメモリ見積もり

次を計算して、結果をMB/GB/TBで書け。

倍精度の密行列 $n \times n$ の保存量： $8 n^{2}$ byte
CSRの概算保存量： $(8 + 4) nnz + 4 (n + 1)$ byte（列indexが4 byteの場合）

設問

(a) $n = 10^{4}, 10^{5}, 10^{6}$ の密行列は何byteか？
(b) 1行あたり平均 $k$ 個の非ゼロ（ $nnz \approx k n$ ）として、
$k = 5, 10, 50$ のときCSRは何byteか？
(c) 「この規模なら密は無理だが疎なら可能」と言える境目を文章で説明せよ。

演習2：計算時間のスケーリング測定（体感）

ベクトルサイズ $n$ を変えながら、次の処理の時間を測り、log-logでプロットせよ（Python等でよい）。

ベクトル内積（期待： $O (n)$ ）
密行列×ベクトル（期待： $O (n^{2})$ ）※小さめのnで
疎行列（CSR）×ベクトル（期待： $O (nnz) \approx O (n)$ ）※1行あたりk固定なら

レポートの観点

実測が理想のO( )から外れる理由を、メモリ階層の観点で説明せよ
疎行列×ベクトルが「計算回数は少ないのに遅い」場合、Rooflineの観点で原因を推測せよ

第2章　浮動小数点と誤差：丸め誤差・桁落ち・悪条件

狙い：数値計算が「なぜ壊れるか」を説明できるようにする

2.0 この章で学ぶこと

この章のテーマは、「数学的には正しいはずの計算が、コンピュータ上では壊れる」理由を、順序立てて理解することである。結論から言うと、壊れる理由は主に次の3つに集約される。

表現できる数が有限（浮動小数点の制約）
演算のたびに丸めが入る（丸め誤差・桁落ち）
問題自体が不安定（悪条件：条件数が大きい）

本章を読み終えると、次の問いに答えられるようになることを目指す。

「なぜ 0.1 + 0.2 == 0.3 が必ずしも真にならないのか？」
「なぜ sqrt(1+x) - 1 のような式は危険なのか？」
「同じアルゴリズムでも、行列の条件数が大きいと解が信用できなくなるのはなぜか？」
「“安定なアルゴリズム”とは何を意味するのか？」

2.1 実数と浮動小数点数は別物：まず“数の世界観”を揃える

数学では、実数は無限に細かく表現できる。しかし計算機では、数はビット列で表されるため、表現できる数は有限である。この有限集合の中で実数を近似する代表方式が IEEE 754 浮動小数点である（[S1]）。

2.1.1 浮動小数点数の一般形

2進の浮動小数点数は概念的に

x = (- 1)^{s} \times m \times 2^{e}

と表される（実際には規格化やビット配置のルールがある）。

$s$ ：符号（0なら正、1なら負）
$m$ ：仮数（significand, mantissa に相当）…「有効桁」を担う
$e$ ：指数（exponent）…「桁の位置（スケール）」を担う

このため、浮動小数点には次の特徴がある。

一定の有効桁数しか持てない（細かい差が表現できない）
非常に大きい数・非常に小さい数は表現限界がある（オーバーフロー／アンダーフロー）
数直線上で等間隔ではない（数が大きいほど刻みが粗い）

2.2 IEEE 754 の重要ポイント：NaN・∞・±0・サブノーマル

IEEE 754 では、通常の有限値だけでなく、例外的な値も定義されている（[S1]）。

2.2.1 特殊値の意味

$+ \infty, - \infty$ ：オーバーフローやゼロ除算などで現れる
NaN (Not a Number)：無効演算（例： $\sqrt{- 1}$ ）などの結果
$+ 0, - 0$ ：ゼロにも符号がある（符号付きゼロ）
サブノーマル（subnormal）：極めて小さい数を“徐々に0へ”落とすための表現（急に0になりにくい）

これらは、「計算を止める」のではなく「値として持ち回って、後で検出できる」ことを意図している。現場では NaN が混入すると、その後の計算が全滅することが多いので、途中で NaN/inf をチェックする習慣が重要である。

2.3 丸め（rounding）とは何か：計算は毎回“近似”される

浮動小数点演算では、実数としての正確な演算結果を、そのまま表現できない場合が多い。そのため、結果は最も近い表現可能な数へ丸められる（[S1]）。

2.3.1 丸めモード（代表例）

IEEE 754 には複数の丸め規則があるが、最も一般的なのは

最近接丸め（ties to even）：最も近い値へ。ちょうど中間なら“偶数側”へ（統計的偏りを減らす）

この「偶数側」は、一見変だが、長い計算で偏りが蓄積しにくいという利点がある。

2.4 機械イプシロンと unit roundoff：誤差を式で扱う準備

数値誤差を議論するには、「丸め誤差の大きさ」を記号で表しておくと便利である。

2.4.1 機械イプシロン（machine epsilon）

機械イプシロン $ε$ は、ざっくり言えば

1 の次に大きい浮動小数点数との差

で定義されることが多い（定義流儀に注意）（[S2]）。倍精度（binary64）では概ね

ε \approx 2^{- 52} \approx 2.22 \times 10^{- 16}

程度である（“1付近の相対刻み”の代表値）。

2.4.2 unit roundoff $u$

「最近接丸め」のもとでは、丸め誤差を扱う基本定数として

u = \frac{ε}{2}

を用いる流儀も多い（[S2][S3]）。この $u$ を使うと、非常に重要な“計算モデル”が書ける。

2.5 浮動小数点の基本モデル：fl( ) と $δ$

「計算機が返す演算結果」を $fl (\cdot)$ と書く（floating の略）。代表的なモデルは

fl (a \circ b) = (a \circ b) (1 + δ), | δ | \leq u

（ $\circ$ は $+, -, \times, /$ など）

である。これは次を意味する。

真の結果 $(a \circ b)$ に対して
相対誤差が高々 $u$ 程度の範囲で
近い数へ丸められる

このモデルは「誤差を定量化する」ための出発点であり、以後の章でも何度も登場する。

注意：これは“だいたいこうなる”ではなく、「正しく丸められる（correct rounding）」という設計思想のもとで成り立つ近似モデルである（[S1]）。

2.6 桁落ち（catastrophic cancellation）：危険な差分

2.6.1 何が起きるのか

桁落ちとは、

近い2つの数を引き算すると、有効桁が大量に失われる

現象である。例えば $a$ と $b$ が 15桁くらい一致していて、差が 1桁しかないとき、

a - b

の結果は元の有効桁をほとんど失う。

2.6.2 典型例： $\sqrt{1 + x} - 1$ （ $x$ が小さい）

$x$ が小さいとき、 $\sqrt{1 + x} \approx 1 + \frac{x}{2}$ なので

\sqrt{1 + x} - 1 \approx \frac{x}{2}

となり、非常に小さい値になる。しかし計算機上では $\sqrt{1 + x}$ が 1 に非常に近い数として丸められ、そこから 1 を引くため、桁落ちが起きやすい。

回避策（有理化）

\sqrt{1 + x} - 1 = \frac{(\sqrt{1 + x} - 1) (\sqrt{1 + x} + 1)}{\sqrt{1 + x} + 1} = \frac{x}{\sqrt{1 + x} + 1}

右辺は引き算が消えており、 $x$ が小さくても安定に計算できる。

このように、同じ数学的量でも“計算の形”で精度が変わる。
大規模計算では誤差が何百万回も積み重なるので、こうした「形の工夫」が効いてくる。

2.7 誤差の種類：絶対誤差・相対誤差・ULP

数値結果 $\hat{x}$ （x-hat）が真値 $x$ を近似するとき、

絶対誤差： $| \hat{x} - x |$
相対誤差： $\frac{| \hat{x} - x |}{| x |}$ （ただし $x \neq 0$ ）

を用いる。

また浮動小数点では ULP（unit in the last place）という概念もよく使う。 ULPは「その値付近での最小刻み幅」に相当し、誤差を“何ULPか”で測ることがある（実装や規格の議論で有用）。

2.8 誤差伝播：足し算でも誤差は積み上がる

2.8.1 逐次和の誤差（なぜ和が危ないのか）

s = \sum_{i = 1}^{n} x_{i}

を逐次に足すと、各ステップで丸めが入るため、誤差が蓄積する。しかも足し算は結合法則が厳密には成り立たない：

(x_{1} + x_{2}) + x_{3} \neq x_{1} + (x_{2} + x_{3})

これは“数学が壊れた”のではなく、丸めが入るからである。

2.8.2 最近の計算で重要な点：並列和と再現性

並列計算では、和を「木構造」で足す（reduction）ため、足し順が変わる。足し順が変われば丸めも変わるので、同じコードでも実行のたびに最終値が微妙に変わることが起こり得る。この問題は近年のGPU/大規模並列で特に顕在化しやすい（再現性の管理が重要になる）。

2.8.3 対策例：Kahan和（補償和）

精度改善の有名手法として Kahan summation がある（[S1]）。「失われた下位桁」を補償変数に蓄えておき、次のステップで回収する考え方である。

2.9 悪条件（ill-conditioned）と条件数：問題そのものが不安定

計算が壊れるとき、原因は「アルゴリズムが悪い」だけとは限らない。問題自体が不安定で、入力が少し変わるだけで解が大きく変わる場合がある。この不安定さを定量化するのが条件数（condition number）である。

2.9.1 関数の条件数（1変数の例）

y = f (x)

で $x$ が少し変わったときの相対変化を考えると、

\frac{| Δ y |}{| y |} \approx κ_{f} (x) \frac{| Δ x |}{| x |}

の形になることが多い。ここで

κ_{f} (x) = | \frac{x f^{'} (x)}{f (x)} |

を条件数と呼ぶことがある。 $κ_{f} (x)$ が大きいほど、入力誤差が出力に大きく増幅される。

2.9.2 行列の条件数（線形方程式）

線形系

A x = b

では、（適切なノルム $∥ \cdot ∥$ に対して）

κ (A) = ∥ A ∥ ∥ A^{- 1} ∥

を条件数と呼ぶ。 $κ (A)$ が大きいほど、 $b$ のわずかな誤差が $x$ に大きく増幅される。

重要：条件数が大きいとき、どんなに良いアルゴリズムでも「解が不安定」になり得る。
つまり “解けた”と“信用できる”は別である。

2.10 前進誤差と後退誤差：安定性を測る2つの物差し

誤差の議論には2つの視点がある（[S4]）。

2.10.1 前進誤差（forward error）

∥ \hat{x} - x ∥

のように、「真の解 $x$ と計算結果 $\hat{x}$ の差」を直接見る。

2.10.2 後退誤差（backward error）

「計算結果 $\hat{x}$ が、少しだけ違う問題の厳密解になっている」と見なす考え方。たとえば線形系なら

(A + Δ A) \hat{x} = b + Δ b

を満たすような $Δ A, Δ b$ が小さいかどうかを見る。

残差による直観（線形系）

残差（residual）

r = b - A \hat{x}

が小さいほど、「 $\hat{x}$ は元の式に近い」という意味で後退誤差が小さいと言える。実務では「残差が小さい＝良い解」と短絡しがちだが、条件数が大きい場合は注意が必要である。

2.11 数値的安定性：後退安定なら（だいたい）安心

アルゴリズムが後退安定（backward stable）であるとは、「計算結果が、元の問題の“わずかな摂動”に対する厳密解になっている」と言える性質である（[S4]）。

そして大まかに

前進誤差 ≲ κ (問題) \times 後退誤差

という関係が成り立つ。つまり

後退誤差が小さく（アルゴリズムが良く）
条件数が小さければ（問題が良く）

前進誤差も小さくなり、結果を信用しやすい。

2.12 最近の大規模計算での重要トピック：混合精度（mixed precision）

近年、GPUやAI向けハードウェアの普及により、FP64（倍精度）より低精度の FP32/FP16/bfloat16/TF32 などを混ぜて高速化する混合精度が重要になっている（[S5]）。

低精度：速い・メモリ転送も軽いが、丸め誤差が大きい
高精度：遅い・重いが、誤差に強い

実務では「全部低精度」ではなく、

主要な反復は低精度で走らせ
必要なところだけ高精度で補正する
（例：反復改良 iterative refinement の考え方）のように、精度と速度の折衷を行う。

この視点は、数値解析と計算資源の両方を理解して初めて扱える“現代的テーマ”である。

2.13 演習

演習2-1：桁落ちの体験（式変形の威力）

$x = 10^{- 8}, 10^{- 12}, 10^{- 16}$ について次を倍精度で評価し、結果を比較せよ。

$\sqrt{1 + x} - 1$
$\frac{x}{\sqrt{1 + x} + 1}$

どの $x$ で差が顕著になるか？理由を「丸め」と「差分」の観点で説明せよ。

演習2-2：条件数の大きい行列で解が不安定になることの確認

ヒルベルト行列など、条件数が大きい行列 $A$ を用意して

A x = b

を解く。次に $b$ に微小摂動 $Δ b$ を加えたとき、解 $\hat{x}$ がどれだけ変わるかを調べよ。「残差は小さいのに解が大きく変わる」ことが起きる理由を述べよ。

演習2-3：和の順序で結果が変わることの確認

同じ配列 ${x_{i}}$ について

前から順に足す
小さい順に並べて足す
Kahan和（可能なら）で結果を比較し、どれが最も安定かを観察せよ。

第3章　行列とグラフ：疎構造・帯行列・再順序化

狙い：疎行列を「データ構造＋グラフ」として理解し、解法設計につなげる

3.0 この章で学ぶこと（なぜグラフが出てくるのか）

大規模数値計算では、行列はしばしば疎行列であり、しかもその疎構造が解法の性能を支配する。ここで重要になる発想は次の2つである。

疎行列は「数値の表」ではなく、非ゼロ位置のパターン（構造）が第一級の情報である
その構造は、グラフとして理解すると直観的で、再順序化や分割が設計できる

本章では、疎行列を扱うための “道具立て” を揃える。

3.1 疎行列の基本： $nnz$ とスパース演算

疎行列 $A$ のサイズが $n \times n$ のとき、非ゼロ要素数を

nnz = # {(i, j) | A_{i j} \neq 0}

と書く。

密行列では $nnz \approx n^{2}$ だが、PDE離散化やネットワーク問題では

nnz \approx O (n)

となることも多い（1行あたり定数個の結合）。

疎行列計算の代表操作は疎行列×ベクトル（SpMV）

y = A x

であり、理想的には計算量が

O (nnz)

となる。

ただし注意：SpMV は演算回数が少なくても、メモリアクセスが散らばるため、実効性能は“メモリ帯域律速”になりやすい（第1章のRooflineの話と接続）。

3.2 疎行列のデータ構造：COO / CSR / CSC

疎行列は、ゼロ要素を保存しないために専用フォーマットを使う。代表が CSR（Compressed Sparse Row）と CSC（Compressed Sparse Column）である（[S6][S7]）。

3.2.1 CSR（行圧縮）：3つの配列で表す

CSRは次の3配列で表す（[S6]）。

data：非ゼロ要素の値を並べた配列（長さ $nnz$ ）
indices：各 data が属する列番号（長さ $nnz$ ）
indptr：各行の開始位置を指すポインタ（長さ $n + 1$ ）

行 $i$ の非ゼロは

indices[indptr[i] : indptr[i+1]]

に列番号が、

data[indptr[i] : indptr[i+1]]

に値が入る。

直観：CSRは「行ごとに必要な非ゼロだけ連続で持つ」ので、SpMV（行方向アクセス）に相性が良い。

3.2.2 CSC（列圧縮）：列ごとにまとめる

CSCはCSRの「列版」であり、列 $j$ の非ゼロが連続に格納される（[S7]）。直観：CSCは列方向の操作（例：列スライス、特定列のアクセス）に向く。

3.2.3 最近の動向：sparse “matrix” から sparse “array” へ

近年のSciPyでは csr_matrix に加えて csr_array のような sparse array も整備されている（[S8]）。ユーザー視点では「NumPy配列に近いふるまい」を目指しており、将来的には

より自然なブロードキャスト
型・index配列の統一など、扱いやすさが増している。

3.3 疎行列をグラフとして見る：隣接関係＝非ゼロ

疎行列の非ゼロパターンはグラフで表現できる。

3.3.1 対称パターンの場合（無向グラフ）

$A$ の疎パターンが対称（ $(i, j)$ が非ゼロなら $(j, i)$ も非ゼロ）なとき、

頂点：変数 $1, \dots, n$
辺： $A_{i j} \neq 0$ （ $i \neq j$ ）なら $i$ と $j$ を結ぶ

で無向グラフになる。 PDE離散化では「隣接格子点が結合」するので、格子グラフに近い構造が出る。

3.3.2 非対称の場合（二部グラフなど）

非対称行列では、行と列を別集合として持つ二部グラフ表現などが便利になる（ここでは詳細は第4章以降で扱う）。

3.4 帯行列・帯幅（bandwidth）と“近くに寄せる”発想

疎行列でも、非ゼロが対角付近に集中していると計算に有利なことがある。その“近さ”を測る代表が帯幅である。

3.4.1 帯幅の定義（代表例）

行列 $A$ の（上側）帯幅を

β = max {| i - j | : A_{i j} \neq 0}

で定義することがある（流儀あり）。 $β$ が小さいほど、非ゼロが対角の近くに集まっている。

重要：帯幅は「行列を並べ替える（再順序化する）」ことで変えられる。
つまり同じ問題でも、変数の番号の付け方で計算が速くなる。

3.5 fill-in（フィルイン）：疎直接法を重くする“新しい非ゼロ”

疎行列を直接法（LU/Cholesky）で分解すると、元々ゼロだった場所に非ゼロが現れることがある。これを fill-in という。

3.5.1 なぜfill-inが問題か

fill-inが増えると

追加のメモリが必要
演算回数が増える（疎なのに密っぽくなる）
並列化の効率が落ちる場合がある

つまり「疎だから楽なはず」が崩れる。

3.5.2 fill-inは順序に依存する

同じ $A$ でも、変数の並び（行・列の順序）を変えると fill-in は大きく変わる。そのため、疎直接法の前処理として ordering（再順序化）が極めて重要になる（[S9][S10]）。

3.6 再順序化（ordering）の代表：RCM と AMD

再順序化は大きく2系統ある。

帯幅を減らす（bandwidth reduction）
fill-in を減らす（fill-reducing ordering）

3.6.1 RCM（Reverse Cuthill–McKee）：帯幅を減らす

RCMは、グラフの幅優先探索（BFS）に基づき、頂点番号を付け替えて帯幅を小さくする手法である（[S9]）。直観的には

近い頂点同士が近い番号になるように並べ替えることで、非ゼロを対角近くへ寄せる。

最近は「開始点の選び方」で品質が変わる問題に対し、改良提案も出ている（[S11]）。

3.6.2 AMD（Approximate Minimum Degree）：fill-inを減らす

AMDは、対称疎行列の因子分解に先立って、fill-in を減らすための順序を与える手法である（[S10][S12]）。考え方の核は「（近似的に）次数が小さい頂点から消去する」ことで、消去のたびに増える辺（＝fill-in）を抑えることにある。

実装では、 $A + A^{T}$ のパターンを作って順序を決めるなどの工夫がある（[S10]）。

3.6.3 最近の動向：ordering自体の並列化

大規模化が進むと、「分解」だけでなく「orderingの計算」も無視できないコストになる。 AMDの並列化や高速化を狙った研究・実装も活発である（[S13][S14]）。

3.7 グラフ分割と nested dissection：さらに大規模へ

fill-in削減の有力戦略として nested dissection（入れ子分割）がある。これは「グラフをセパレータ（separator）で分けて、階層的に順序を作る」発想である（[S15][S16]）。

3.7.1 直観：セパレータで“切ってから解く”

大きいグラフを、少数の頂点（セパレータ）を取り除くことで2つに分け、

左側を先に
右側を次に
最後にセパレータの順に処理することで、fill-in を抑えつつ並列性も確保しやすい。

3.7.2 METIS：実務でよく使われるツール

METISはグラフ分割や疎行列orderingを提供するソフトウェアとして広く使われている（[S17]）。実務では「AMD/COLAMD/Metis のどれが良いか」は問題構造に依存し、実際に試して選ぶことも多い。

3.8 実務的な設計指針：どの順序を選ぶか

初学者が最初に迷いやすい点なので、最低限の指針を置く。

目的が帯幅削減（視覚化・局所性改善） → RCM系をまず試す
疎直接法（Cholesky/LU）のfill-in削減 → AMD / nested dissection（METIS等）を試す
反復法（CG/GMRES）中心でも、前処理（ILU/IC）に ordering が効くことがある（第7章で再登場）

大規模計算では「解法」だけでなく、「順序（ordering）」が“アルゴリズムの一部”になる。

3.9 演習

演習3-1：疎行列の可視化（spyプロット）

PDE由来の疎行列（Poissonの離散化など）を作り、非ゼロパターンを可視化せよ。

元の順序のパターン
RCM後のパターンを並べて、非ゼロが対角へ寄ることを確認せよ。

演習3-2：再順序化で fill-in が減ることの体験

小さめの疎対称行列 $A$ を用意し、次を比較せよ。

orderingなしで Cholesky（または疎LU）
AMD ordering を入れてから分解
nested dissection（可能ならMETIS）で ordering してから分解

「因子の非ゼロ数（nnz）」や「分解時間」がどう変わるかを計測し、理由を説明せよ。

演習3-3：CSR/CSC の中身を手で追う

小さな疎行列を1つ作り、CSRの data/indices/indptr を自分で書き出して説明せよ。行 $i$ の非ゼロがどのスライスに入っているかを、式と対応づけて述べよ。

参考文献

[S1] IEEE 754 の概要（特殊値・丸め・FMA・例外など）
[S2] machine epsilon / unit roundoff の定義の注意点
[S3] 浮動小数点の入門（単位丸め・丸め誤差モデル）
[S4] backward error の直観的説明（前進誤差との違い）
[S5] TF32/混合精度が高速化に効く理由（帯域・Tensor Core）
[S6] SciPy CSR の定義（data/indices/indptr の関係）
[S7] CSC の定義（columns の indptr、rows の indices）
[S8] SciPy sparse arrays（csr_array 等）とdtype/index配列の注意
[S9] RCM（帯幅削減、BFS、reverseの意図）
[S10] AMD の目的と基本説明（fill-in削減、A+A^T パターン）
[S11] RCMの改良（開始点問題など、近年の研究例）
[S12] AMD 原論文（Approximate Minimum Degree ordering）
[S13] ParAMD 等：ordering自体の高速化・並列化の動き
[S14] AMD並列化に関する近年の報告（orderingコストが無視できない背景）
[S15] nested dissection と fill-reducing ordering の概念
[S16] nested dissection の実務的解説（再訪・実装観点）
[S17] METIS マニュアル（分割・orderingの提供）

第4章　直接法（基礎と限界）：LU/Cholesky と疎直接法の考え方

狙い：直接法が強い場面と、大規模で厳しくなる理由を把握する

4.0 この章で学ぶこと

直接法（direct method）は、線形方程式

A x = b

を「有限回の手順」で（丸め誤差はあっても）原理的には一気に解く方法である。代表は LU分解と Cholesky分解である。

直接法の強み：
1. 収束判定が不要で、一度分解できれば安定に解が出る
2. 同じ $A$ に対して右辺 $b$ が多数あるとき（多点計算・時間発展の各ステップなど）、因子を再利用できて強い
直接法の限界（大規模で厳しくなる理由）：
1. 密行列では計算量が概ね $O (n^{3})$ 、メモリが $O (n^{2})$ と重い
2. 疎行列でも fill-in（フィルイン）により、因子が疎でなくなりメモリが爆発し得る
3. 並列化すると、演算より通信・同期・メモリ帯域が支配的になりやすい（第1章の話と接続）

本章の到達目標は、次を自分の言葉で説明できることである。

LU/Cholesky が「何をしているか」（ガウス消去との対応）
なぜピボット（pivot）が必要で、どのように安定性に関係するか
疎直接法が「(1) 解析（ordering）→ (2) 分解 → (3) 三角解法」から成り、特に fill-in が本質である理由
直接法と反復法（次章）を「どの場面で選ぶか」の判断基準

4.1 ガウス消去：直接法の原点（式変形として理解する）

まず密行列の基本として、ガウス消去を「何をしているか」から押さえる。

4.1.1 目的：連立一次方程式を“上三角”にしてから解く

連立一次方程式

A x = b, A \in R^{n \times n}, x, b \in R^{n}

は、 $A$ が上三角行列 $U$ であれば後退代入で解ける。したがって狙いは

A \Rightarrow U

に変形することにある。

4.1.2 消去の1ステップ（2×2の例で直観）

例として

(\begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) = (\begin{matrix} b_{1} \\ b_{2} \end{matrix})

を考える。下の式から $a_{21}$ を消したいので、上の行を係数

m_{21} = \frac{a_{21}}{a_{11}}

倍して引く（行基本変形）：

{row}_{2} \leftarrow {row}_{2} - m_{21} {row}_{1}

これにより $(2, 1)$ 成分が0になり、上三角化が進む。

4.1.3 一般のアルゴリズム（密行列の典型形）

ステップ $k = 1, 2, \dots, n - 1$ について、

ピボット（割る数）： $a_{k k}$
乗数（multiplier）： $m_{i k} = \frac{a_{i k}}{a_{k k}} (i = k + 1, \dots, n)$
更新： $a_{i j} \leftarrow a_{i j} - m_{i k} a_{k j} (i = k + 1, \dots, n, j = k, \dots, n)$ $b_{i} \leftarrow b_{i} - m_{i k} b_{k}$

この操作は「行列の左から下三角行列を掛けていく」ことと同値であり、それが LU 分解へつながる。

4.2 LU分解：ガウス消去を“分解”として保存する

ガウス消去で現れる乗数 $m_{i k}$ を下三角に保存すると、行列 $A$ を

A = L U

と因子分解できる。ここで

$L$ ：単位下三角行列（対角が1）
$U$ ：上三角行列

である。

4.2.1 なぜ分解がうれしいのか（解く手順が2段）

A x = b \Rightarrow L U x = b

と書けるので、まず

L y = b (前進代入)

を解いてから

U x = y (後退代入)

を解けばよい。ここで $y$ は中間ベクトルである。

分解（factorization）は高価
代入（triangular solve）は比較的安い

という性質があるため、「同じ $A$ に対して多数の右辺を解く」用途で直接法は特に強い。

4.2.2 計算量（密行列の目安）

密行列のLU分解は概ね

演算量 \sim O (n^{3}), メモリ \sim O (n^{2})

であり、 $n$ が大きいと急速に厳しくなる。

4.3 ピボット（pivot）と安定性：なぜ入れ替えが必要か

ガウス消去では $a_{k k}$ で割るため、もし $a_{k k} = 0$ なら計算不能である。さらに $a_{k k}$ が非常に小さいと

m_{i k} = \frac{a_{i k}}{a_{k k}}

が巨大になり、丸め誤差が増幅されやすい（第2章：誤差・悪条件と直結する）。

4.3.1 部分ピボット（partial pivoting）

実務で標準的なのは、列 $k$ の $k$ 行以降で絶対値最大の要素をピボットに選び、行交換する方法である。これを繰り返すと

P A = L U

という形になる。

$P$ ：置換行列（行の入れ替えを表す）
$L, U$ ：分解結果

4.3.2 ピボット成長（pivot growth）という危険信号

分解の途中で $U$ の成分が非常に大きくなると、丸め誤差が増幅されやすい。実務では「ピボット成長」や関連指標を監視し、解が信用できるかを判断することがある。LAPACK には、成長に関係する情報（例：reciprocal pivot growth factor）を返すルーチン群がある（[D7]）。

初学者向けの要点：
「分解できた＝常に信用できる」ではない。
ただし“適切なピボット＋後述の反復改良”により、実務上は安定に解ける場面が多い。

4.4 Cholesky分解：対称正定値（SPD）では最強に近い

もし $A$ が対称正定値（symmetric positive definite, SPD）なら、

A = L L^{T}

（または $A = R^{T} R$ ）と分解できる。これが Cholesky 分解である。

4.4.1 SPDの定義

A = A^{T}, x^{T} A x > 0 (\forall x \neq 0)

を満たすとき $A$ はSPDである。

4.4.2 Choleskyが強い理由

LUより演算量が少ない（概ね半分程度）
ピボットが不要（SPDなら理論的に安定）
因子が三角1つ（ $L$ だけ）で済むため、メモリも有利

実務では、PDE（拡散・ポアソン）や最小二乗の正規方程式などで SPD が現れ、Cholesky が第一候補になる。

4.5 疎直接法：3段階（解析→分解→解法）で理解する

大規模計算で重要なのは「疎行列に対する直接法」である。疎直接法は、概念的に次の3段階で設計される。

解析（symbolic analysis）：非ゼロ構造だけ見て ordering と fill-in を見積もる
数値分解（numeric factorization）：実際に因子（L,U あるいは L）を計算する
解（solve）：前進・後退代入（多右辺なら繰り返す）

ここで「1.解析」が重要である理由は、疎行列では fill-in が順序に強く依存し、メモリと計算量を支配するからである（第3章で扱った再順序化がここで主役になる）

4.6 fill-in（フィルイン）：疎直接法の“真のコスト”

4.6.1 fill-inとは

疎行列 $A$ を分解すると、元はゼロだった位置が非ゼロになることがある。これを fill-in と呼ぶ。

元の非ゼロ数： $nnz (A)$
分解後の非ゼロ数： $nnz (L) + nnz (U)$ （または $nnz (L)$ ）

が大きくなるほど、メモリも時間も増える。

4.6.2 なぜfill-inが起きるのか（グラフの直観）

対称パターンの場合、疎行列はグラフに対応し、消去操作は「隣接頂点同士を結ぶ辺を追加する」操作に対応する。消去順序が悪いと、辺が大量に増え（＝fill-inが増え）、最終的に“ほぼ密”になってしまう。

4.6.3 直接法が厳しくなる典型ケース

3次元PDEの離散化（近傍結合が多く、分解で結合が広がりやすい）
グラフのセパレータが大きい構造（切り分けにくい）

このため大規模では、直接法単独ではメモリが持たず、反復法＋前処理へ移行する場面が多い（次章の主題）。

4.7 疎直接法の実装の中身：multifrontal / supernodal という考え方

疎直接法は「ゼロを保存しない」だけでなく、キャッシュや並列性を活かすために、内部で“まとまり”を作って計算する。

4.7.1 multifrontal（多前面法）の直観

消去木（elimination tree）に沿って、小さな密行列（フロント、front）を作り、

フロント内は密演算（BLAS-3）で高速に処理し
部分結果を親へ“組み立てる”

という方法である。MUMPS は並列多前面法として広く知られる（[D3]）。

4.7.2 supernodal（スーパーノード）Cholesky の直観

Choleskyでは、同じ非ゼロパターンを持つ連続列を束ねて「スーパーノード」を作り、密ブロック演算としてまとめて処理する。CHOLMOD などで採用されている（[D1]）。

初学者向けの要点：
疎直接法は「疎をそのまま扱う」だけでなく、
“疎の中に現れる密の塊”を見つけて、そこは密として高速化する設計である。

4.8 近年の実務的トレンド：GPU・低ランク・混合精度

「最近の情報」を、初学者が押さえるべき範囲で整理する。

4.8.1 GPU支援の疎直接法

疎直接法でも、フロント（密塊）や三角解法の一部をGPUへ載せて高速化する流れがある。SuiteSparse の CHOLMOD/SPQR には CUDA 加速に関する記述がある（[D1]）。また SuperLU も MPI/OpenMP/CUDA を使う実装系を持つ（[D2]）。

4.8.2 メモリ削減：BLR（Block Low-Rank）など

MUMPS のユーザーガイドには、ある種の問題で因子のメモリ削減や計算削減を狙う BLR 機能が言及されている（[D3]）。これは「フロント内部が（近似的に）低ランク構造を持つ」ことを利用する発想である。

直接法で得た近似解を、残差を使って高精度化する反復改良が現代的に重要である（[D8]）。特に、低精度（高速）で初期解を作り、高精度で残差を評価して補正する「混合精度反復改良」が普及している（[D8][D9]）。

4.9 直接法はいつ使うべきか

初学者が実務で迷わないための、最低限の判断基準を置く。

4.9.1 直接法が第一候補になりやすい

$n$ が中規模まで（例：メモリが十分で、因子が保持できる）
SPDで Cholesky が使える
$A$ が固定で、多数の右辺 $b$ を解く（再利用が効く）
高精度が必要で、反復法の収束が不安／前処理が難しい

4.9.2 直接法が厳しくなりやすい

fill-in が大きく、 $nnz (L) + nnz (U)$ が巨大になる
3D問題などで因子が“準密”化してメモリ- (a) 毎回 LU して解く
(b) 一度だけ LU して、三角解法を100回するの時間を比較せよ。なぜ (b) が強いかを「分解 vs 代入」の観点で説明せよ。

演習4-2：疎化すると逆に遅くなる例を観察する

サイズは中程度だが、fill-in が大きくなる疎行列（順序が悪いもの）を用意し、

ordering なし
AMD/ND 等の ordering ありで分解時間・因子の nnz を比較せよ。「疎行列なのに遅い」原因を fill-in とメモリ階層（第1章）で説明せよ。

演習4-3：ピボットの役割を体験する

ピボットなしのガウス消去が破綻する行列（小さいピボットが現れる例）を作り、

ピボットなし
部分ピボットありで解の誤差や残差を比較せよ。なぜピボットが必要かを、丸め誤差（第2章）と関連づけて述べよ。

参考文献

[D1] SuiteSparse（UMFPACK/CHOLMOD/SPQR、CUDA加速の記述を含む）
[D2] SuperLU（MPI/OpenMP/CUDA を含む疎直接ソルバ）
[D3] MUMPS Users’ guide（並列多前面法、BLR などの機能）
[D7] LAPACK（ピボット成長などの診断情報）
[D8] 混合精度反復改良（iterative refinement）に関する近年の代表論文
[D9] cuSOLVER（反復改良と高精度パッケージの記述）

第5章　反復法の基本：固定点反復とKrylov部分空間の直観

狙い：大規模の主役「反復法」の全体像を掴む

5.0 この章で学ぶこと

大規模線形問題

A x = b, A \in R^{n \times n}

で $n$ が非常に大きいとき、直接法は

因子のメモリが持たない
fill-in が増えすぎる
並列化しても通信が重い

などの理由で厳しくなることがある。そこで主役になるのが反復法（iterative method）である。

反復法の核心は次の一文に尽きる。

「行列 $A$ そのものを分解して保存するのではなく、
行列-ベクトル積 $y = A x$ を何度も計算しながら、解へ近づく」

この発想により、

必要メモリは主に $nnz (A)$ に比例（疎なら軽い）
1ステップの計算は概ね $O (nnz)$
“行列を明示しない” matrix-free 設計が可能

となり、超大規模へ伸びる。

本章の到達目標は以下である。

固定点反復と反復行列（iteration matrix）の意味を説明できる
収束条件が「スペクトル半径 $ρ (T) < 1$ 」で表される理由を説明できる
Krylov部分空間法が「多項式近似」と「最小化」の問題だと直観できる固定点反復は

x^{(k + 1)} = g (x^{(k)})

で表される。解 $x^{\*}$ は固定点

x^{\*} = g (x^{\*})

として定義される。

5.1.2 1次元での収束条件

1次元で $x^{\*}$ の近くにいるとき、テイラー展開より

x^{(k + 1)} - x^{\*} \approx g^{'} (x^{\*}) (x^{(k)} - x^{\*})

となる。したがって

| g^{'} (x^{\*}) | < 1

なら誤差が縮み、収束する。これが「局所収束条件」の基本である。

5.2 線形方程式を反復に落とす：行列分割（splitting）

線形方程式

A x = b

を固定点反復にする標準ルートは、行列分割である。

5.2.1 分割の定義

A = M - N

と分割し、 $M$ が解きやすい（例えば対角・下三角など）ように選ぶ。すると

(M - N) x = b \Rightarrow M x = N x + b

x = M^{- 1} N x + M^{- 1} b

となるので、反復は

x^{(k + 1)} = T x^{(k)} + c

T = M^{- 1} N, c = M^{- 1} b

で与えられる。ここで $T$ を反復行列（iteration matrix）と呼ぶ。

5.3 残差と誤差：反復法の“計測器”

反復法では、真の誤差 $e^{(k)} = x - x^{(k)}$ は通常わからない（真の解 $x$ を知らないため）。そこで観測可能な量として残差を使う。

5.3.1 残差（residual）

r^{(k)} = b - A x^{(k)}

で定義する。 $r^{(k)} = 0$ なら厳密解である。

5.3.2 残差と誤差の関係

真の誤差

e^{(k)} = x - x^{(k)}

について、

A e^{(k)} = A (x - x^{(k)}) = b - A x^{(k)} = r^{(k)}

すなわち

r^{(k)} = A e^{(k)}

が成り立つ。よって、もし $A^{- 1}$ が“穏やか”なら

∥ e^{(k)} ∥ \leq ∥ A^{- 1} ∥ ∥ r^{(k)} ∥

となり、残差が小さいほど誤差も小さいと期待できる。ただし $∥ A^{- 1} ∥$ が大きい（悪条件）と、残差が小さくても誤差が小さいとは限らない（第2章の条件数と接続）。

5.4 収束条件： $ρ (T) < 1$ がなぜ出るのか

反復

x^{(k + 1)} = T x^{(k)} + c

の真の解 $x$ は

x = T x + c

を満たす。両者の差を取ると

e^{(k + 1)} = x - x^{(k + 1)} = T (x - x^{(k)}) = T e^{(k)}

つまり誤差は

e^{(k)} = T^{k} e^{(0)}

で伝播する。よって $T^{k} \to 0$ なら収束する。

線形代数の結果として、十分な条件として

ρ (T) < 1

（ $ρ (T)$ は $T$ の固有値の絶対値最大＝スペクトル半径）が重要になる。直観的には「誤差のモード（固有ベクトル成分）が、反復を重ねるごとに縮む」ことを意味する。

5.5 Jacobi法とGauss–Seidel法：最初に触るべき2つの反復

5.5.1 Jacobi法

$A$ を

A = D + (L + U)

（ $D$ ：対角、 $L$ ：下三角、 $U$ ：上三角）と分け、

M = D, N = - (L + U)

とする。すると

x^{(k + 1)} = D^{- 1} (b - (L + U) x^{(k)})

となる。

1ステップが軽い
並列化しやすい（各成分が同時更新できる）
収束は遅いことが多い（ $ρ (T)$ に依存）

5.5.2 Gauss–Seidel法

M = D + L, N = - U

とすると

x^{(k + 1)} = (D + L)^{- 1} (b - U x^{(k)})

となる。各成分が逐次更新されるので、Jacobiより収束が速いことが多い一方で、並列化は難しくなる（依存関係があるため）。

5.6 反復法の“本命”：Krylov部分空間法

Jacobi/GS は理解の入口として重要だが、超大規模で強いのは Krylov 部分空間法である。現代の実務コードの中心には「Krylov法＋前処理」がある（[I1]）。

5.6.1 Krylov部分空間の定義

初期値 $x^{(0)}$ を取り、初期残差

r^{(0)} = b - A x^{(0)}

を定義する。Krylov部分空間は

K_{k} (A, r^{(0)}) = span {r^{(0)}, A r^{(0)}, A^{2} r^{(0)}, \dots, A^{k - 1} r^{(0)}}

で定義される。

5.6.2 何がうれしいのか：多項式近似

Krylov空間にいるベクトルは $r^{(0)}$ に $A$ を繰り返し作用させた線形結合であるから、残差は

r^{(k)} = p_{k} (A) r^{(0)}

のような形（ $p_{k}$ は多項式）で表される。
よって「うまい多項式 $p_{k}$ を選んで、 $p_{k} (A)$ が誤差成分を消す」ことができれば、少ない反復で収束する。

5.7 代表的Krylov法の考え方

5.7.1 CG（共役勾配法）：SPDの場合

$A$ がSPDなら CG が第一候補になる。CGは

あるエネルギー（2次形式）を最小化する
3項漸化式で進む（メモリが少ない）などの特徴があり、巨大SPD問題に強い。

5.7.2 GMRES：非対称・一般の場合

一般の非対称行列に対して代表的なのが GMRES である。
GMRESは「残差ノルム $∥ r^{(k)} ∥_{2}$ を最小化する」反復であり、そのために Arnoldi 法で直交基底を作る（詳細は後章で扱う）。

ただし GMRES は反復とともに記憶が増えるため、実務では

再スタート GMRES(m)
柔軟GMRES（前処理が変化する場合）などが用いられる。

5.8 なぜ「行列-ベクトル積が中心」になるのか：計算資源

Krylov法の主なコストは

1回の行列-ベクトル積 $y = A x$ （疎なら $O (nnz)$ ）
内積・ノルム計算（並列では同期が必要）である。

5.8.1 反復法が大規模に向く理由

行列因子を保存しないのでメモリが軽い
行列-ベクトル積は局所計算に分解しやすく、並列化しやすい
PDEのような局所結合では、通信は主に境界データ交換に限定されやすい

5.9 近年の実務的トレンド：通信削減（pipelined Krylov）と混合精度

「最近の情報」を、初学者が理解できる形で要点だけ整理する。

5.9.1 大規模並列で効く：pipelined Krylov

Krylov法は内積を多用し、そのたびに全体同期（global reduction）が発生する。ノード数が大きいと同期が支配的になりうる。この問題に対し、同期回数を減らす pipelined Krylov が実装されている（例：PETSc の PIPEFGMRES）（[I2]）。

要点：
反復法は「計算量」だけでなく、通信（同期）の回数が性能を支配し得る。
したがって“速い反復法”は“同期の少ない反復法”でもある。

5.9.2 混合精度＋反復改良

反復法でも、低精度演算で高速に進めつつ、高精度で誤差を補正する設計が広がっている（[D8][D9]）。 GPUライブラリでも、必要に応じて高精度反復改良を使う旨が説明されている（[D9]）。

5.10 第5章のまとめ：反復法を設計する最初の一歩

この章の核心は以下である。

反復法は「残差を観測しながら解へ近づく」方法である
収束性は反復行列のスペクトル（ $ρ (T)$ ）と強く関係する
大規模で主役になるのは Krylov法であり、行列-ベクトル積が中心になる
近年は、通信削減（pipelining）と混合精度が性能の鍵になりやすい

次章以降では、Krylov法の具体アルゴリズム（CG/GMRES等）、前処理、停止判定を体系的に扱う。

5.11 演習

演習5-1：Jacobi/Gauss–Seidel の収束性をスペクトルと結びつける

小さな行列 $A$ をいくつか用意し、Jacobi と Gauss–Seidel の反復行列 $T$ を作り、

$ρ (T)$ を計算（もしくは近似）
実際の収束の速さ（反復回数）と比較せよ。なぜ $ρ (T)$ が小さいと速いのかを説明せよ。

演習5-2：残差は減るのに誤差が減らない例を探す

条件数が大きい行列で、

残差 $∥ r^{(k)} ∥$ は小さくなるが
解の誤差 $∥ e^{(k)} ∥$ が思ったほど小さくならない例を観察せよ。その理由を「 $r = A e$ 」「 $∥ A^{- 1} ∥$ 」の観点で述べよ。

演習5-3：matrix-free の発想を体験する

明示的に $A$ を保存せず、関数として $y = A x$ を返す実装を用意し、Krylov法（利用可能な範囲で）を適用せよ。「保存するのは何で、何が不要になるか」を言語化せよ。

参考文献

[I1] PETSc KSP：Krylov法＋前処理が現代コードの中心であるという整理
[I2] PETSc の PIPEFGMRES：pipelined Krylov の実装例
[D8] 混合精度反復改良（iterative refinement）の代表論文
[D9] cuSOLVER：必要に応じて高精度反復改良を使うという説明

第6章　代表的Krylov法：CG・GMRES・BiCGStab

6.0 この章の狙い

この章では、大規模線形方程式

A x = b

を解くときに実戦で頻出する Krylov部分空間法の代表例として、CG, GMRES, BiCGStab を学ぶ。特に次の3点を説明できるようになることを狙う。

行列 $A$ の性質（対称性・正定値性・非対称性）に応じて解法を選べる
反復法の停止条件（「どこで止めるのが正しいか」）を数式で説明できる
大規模計算で支配的になるコスト（行列–ベクトル積・内積の通信・メモリ）を意識して実装・運用できる

6.1 Krylov部分空間とは何か（反復法の“作業空間”）

反復法は「少しずつ近づける」方法であり、その“近づけ方”を体系化したのがKrylov法である。

6.1.1 残差と誤差

ある近似解 $x_{k}$ に対して

残差（residual） $r_{k} = b - A x_{k}$
誤差（error） $e_{k} = x^{*} - x_{k} (x^{*} : 真の解)$

を定義する。真の解は未知なので、実際の計算では残差 $r_{k}$ を監視して収束判定するのが基本である。

6.1.2 Krylov部分空間の定義

初期残差 $r_{0} = b - A x_{0}$ を起点に、次の空間を考える：

K_{k} (A, r_{0}) = span {r_{0}, A r_{0}, A^{2} r_{0}, \dots, A^{k - 1} r_{0}}

この $K_{k}$ を Krylov部分空間と呼ぶ。
直観的には「 $A$ を何回も掛けて得られる方向の組を集めた空間」であり、反復法は

x_{k} \in x_{0} + K_{k} (A, r_{0})

という制約のもとで「最もよい $x_{k}$ を選ぶ」戦略を取る。

6.2 停止条件：相対残差・絶対残差・最大反復

大規模計算では「無限に回せば良い」わけではない。止め方はアルゴリズムの一部であり、目的精度と計算コストを結びつける。

6.2.1 典型的な停止条件

残差ノルム $∥ r_{k} ∥_{2}$ を用いて、例えば

∥ r_{k} ∥_{2} \leq max (rtol \cdot ∥ r_{0} ∥_{2}, atol)

のように判定することが多い。

rtol（relative tolerance）: 初期残差に対する相対精度
atol（absolute tolerance）: 絶対精度
さらに安全のため最大反復回数 $k_{max}$ を設定する

6.2.2 大規模並列で重要になる注意：内積＝通信

Krylov法の多くは $⟨ u, v ⟩$ のような内積を使う。分散メモリ並列（MPI）では、内積は全プロセス集約（Allreduce）を伴い、計算より通信待ちが支配的になることがある。
この背景が、後述する pipelined（通信隠蔽）Krylov法が近年重要になっている理由である。

6.3 CG法（Conjugate Gradient）：SPD行列の“最適解”

6.3.1 いつ使うか：SPD条件

CGは次の条件で理論的に保証が強い。

$A$ が対称： $A = A^{⊤}$
$A$ が正定値： $\forall z \neq 0, z^{⊤} A z > 0$

このとき $A$ は SPD (Symmetric Positive Definite) と呼ばれ、CGは非常に強力な選択になる。

6.3.2 何を最適化しているか：二次関数最小化

SPDのとき、線形方程式 $A x = b$ は二次関数

ϕ (x) = \frac{1}{2} x^{⊤} A x - b^{⊤} x

の最小化と等価である（ $\nabla ϕ (x) = A x - b$ ）。
CGはこの $ϕ (x)$ を Krylov空間上で効率的に最小化する方法だと理解できる。

6.3.3 アルゴリズム（基本形）

初期値 $x_{0}$ から開始し

r_{0} = b - A x_{0}, p_{0} = r_{0}

として、反復ごとに

ステップ長 $α_{k} = \frac{r_{k}^{⊤} r_{k}}{p_{k}^{⊤} A p_{k}}$
解の更新 $x_{k + 1} = x_{k} + α_{k} p_{k}$
残差の更新 $r_{k + 1} = r_{k} - α_{k} A p_{k}$
方向係数 $β_{k} = \frac{r_{k + 1}^{⊤} r_{k + 1}}{r_{k}^{⊤} r_{k}}$
探索方向の更新 $p_{k + 1} = r_{k + 1} + β_{k} p_{k}$

という短い再帰（short recurrence）で進む。
保存するベクトルが少ない（数本）ため、メモリ面でも有利。

6.3.4 収束の直観：条件数が効く

SPDでは条件数

κ (A) = \frac{λ_{max} (A)}{λ_{min} (A)}

が収束に強く影響する（ $λ$ は固有値）。
ざっくり言えば $κ (A)$ が大きいほど収束は遅い。
この「 $κ$ を下げる」役割が次章の前処理である。

6.3.5 大規模並列の話題：pipelined CG

CGは反復ごとに内積（通信）を複数回行う。近年は通信待ちが支配的になりやすいため、内積を減らし計算と通信を重ねる pipelined CG が使われることがある。

6.4 GMRES：非対称でも“残差最小”を狙う王道

6.4.1 いつ使うか：非対称・非正定値

$A$ が非対称（ $A \neq A^{⊤}$ ）な場合、CGは原理的に使えない（保証が壊れる）。
そこで登場するのがGMRESである。

6.4.2 GMRESの基本思想：残差ノルム最小化

GMRESは

x_{k} \in x_{0} + K_{k} (A, r_{0})

の中で

∥ b - A x_{k} ∥_{2} = ∥ r_{k} ∥_{2}

を最小にする $x_{k}$ を選ぶ方法である。
つまり「Krylov空間の中で最も残差が小さい解」を毎回選ぶので、直観的にも分かりやすい。

6.4.3 Arnoldi過程と直交化（なぜ重いか）

GMRESは $K_{k}$ の正規直交基底 ${v_{1}, \dots, v_{k}}$ を作るために Arnoldi過程を使う。
このとき新しいベクトルを過去の全ベクトルに対して直交化する必要があり、反復が進むほど

メモリ：基底ベクトルを全部保存（ $k$ 本）
計算：直交化のコストが増大
通信：内積が増える（Allreduceが増える）

という理由で重くなる。

6.4.4 Restart（GMRES(m)）とトレードオフ

この問題を避けるため、実務では Restart GMRES が標準的：

$m$ 回反復したらいったん止めて、得た $x_{m}$ を新しい初期値にして再スタート
これを GMRES( $m$ ) と書く

利点：メモリと計算が上限 $m$ で抑えられる
欠点：再スタートで情報を捨てるため、収束が鈍る（停滞する）ことがある

6.4.5 直交化の実務：MGS vs CGS

直交化は数値安定性に直結する。一般に

MGS（修正グラム・シュミット）：安定だがやや遅い
CGS（古典グラム・シュミット）：速いが不安定になりやすい（必要なら反復改良で補う）

という性質があり、実装やライブラリのオプションで選べることが多い。

6.5 BiCGStab：短い再帰で非対称を解く“軽量”ルート

6.5.1 いつ使うか：非対称で、メモリを抑えたい

GMRESは堅牢だが、反復が増えるとメモリが重い。
BiCGStabは短い再帰で進むため、保存ベクトルが少なく、メモリが軽い。

6.5.2 直観：BiCGの不安定さを“滑らか”にした

BiCGStabは BiCG（双共役勾配法）を基にし、残差の振る舞いを“安定化（stabilize）”することで、より滑らかな収束を狙う。

注意点として

収束が速いこともあるが、問題によっては停滞・破綻もあり得る
前処理の設計次第で大きく性能が変わる

という「当たり外れ」があり、実戦ではGMRESと比較して使い分ける。

6.5.3 pipelined BiCGStab（通信隠蔽）

BiCGStabも内積（通信）を含むため、大規模並列では通信隠蔽（pipelined）版が研究・実装されている。
ただし通信隠蔽版は、数値安定性の工夫（残差置換など）が重要になることがある。

6.6 どう選ぶか：CG / GMRES / BiCGStab の意思決定

同じ $A x = b$ でも、行列の性質・メモリ制約・並列環境で最適解が変わる。

$A$ がSPD：まず CG（+前処理）
非対称・頑健性重視：GMRES(m)（+前処理）
非対称・メモリ節約・軽量志向：BiCGStab（+前処理）
前処理が反復ごとに変わる／非線形：FGMRES など柔軟型を検討

6.7 演習

演習6A：同一問題をCG/GMRESで比較

2種類の行列を用意する
- SPD（例：Poisson離散化）
- 非対称（例：移流拡散の離散化）
各々に対して
- CG（可能なら）
- GMRES(m)
- BiCGStab
  を適用し、反復回数・時間・最終残差を比較せよ。
「内積回数（＝通信回数）」が多いほど並列で不利になる、という観点で結果を考察せよ。

演習6B：停止条件の影響

rtol を $10^{- 2}, 10^{- 4}, 10^{- 6}, 10^{- 8}$ に変えたとき、反復回数・時間がどう変わるか観察せよ。
目的（可視化用途／設計最適化の内側ループ／高精度検証）に応じて「適切な停止条件」が変わる理由を文章で説明せよ。

第7章　前処理：なぜ効くのか、どう作るのか

7.0 この章の狙い

前処理（preconditioning）は「反復法の性能を決める主役」である。
この章では次の3点を身につける。

前処理が何を変えているのか（方程式の同値変形）を数式で説明できる
「効く前処理」の直観：条件数・固有値分布・クラスタリングを言葉で説明できる
Jacobi/SSOR/ILU/IC/近似逆行列/分割型などの代表例を、コストと並列性まで含めて選べる

7.1 前処理の基本：同じ解を保ったまま“解きやすい形”にする

線形方程式

A x = b

に対して「解 $x$ は同じだが、反復法が速く進む形」に変形するのが前処理。

7.1.1 左前処理（Left preconditioning）

$M \approx A$ を前処理行列（または作用素）として

M^{- 1} A x = M^{- 1} b

を解く。反復法は実際には $M^{- 1}$ を明示的に作らず、「 $y = M^{- 1} v$ を解く操作」を繰り返し使う。

7.1.2 右前処理（Right preconditioning）

A M^{- 1} y = b, x = M^{- 1} y

と変形して $y$ を解く。
右前処理は「残差の定義や監視の仕方」に注意が必要だが、柔軟型Krylov法（FGMRESなど）と相性がよい。

7.1.3 対称（左右）前処理（Symmetric preconditioning）

対称性を保ちたい場合に

\tilde{A} = D^{- 1} A D^{- T}

のような形を考えることがある（ただし、適用可能な手法や実装に制約が出る）。

7.2 なぜ効くのか：条件数とスペクトルの整理

7.2.1 条件数を下げる＝反復回数を減らす

特にSPDの場合、CGの収束が $κ (A)$ に依存するため

κ (M^{- 1} A) ≪ κ (A)

となるように $M$ を作れれば、反復回数が大きく減る。

7.2.2 固有値の“まとまり”が重要

GMRES等では単純に条件数だけでは語れず、固有値が

ある点の周りにクラスタしている
0付近に危険な固有値が少ない
スペクトルが“良い形”に整う

などが効いてくる（ここは数値線形代数の発展的内容）。
ただし初学者としては、次の経験則をまず押さえる：

前処理の狙いは「反復法が苦手な成分（遅い誤差モード）を先に潰す」こと
そのためには、問題の物理（PDEのタイプ、係数の変化、メッシュ品質）に合う前処理を選ぶ必要がある

7.3 代表的前処理（1）：Jacobi（対角スケーリング）

7.3.1 定義

最も簡単な前処理は

M = diag (A)

つまり対角成分のみを使う方法（Jacobi前処理、スケーリング）。

7.3.2 何が嬉しいか

実装が簡単
並列化が容易（各成分が独立）
メモリが軽い

一方で、強く結合した問題（PDEの強い結合、異方性、悪条件）では効果が小さいことが多い。

7.4 代表的前処理（2）：(S)SOR / SSOR（Gauss–Seidel系）

7.4.1 基本イメージ

SOR/SSORは「近傍の情報を使って順次更新する」タイプで、Jacobiよりは局所的に強い。

7.4.2 並列での注意

通常のGauss–Seidelは逐次性が強く、分散メモリ並列では工夫が必要。
実務では「局所（各サブドメイン内）でSSOR」など、分割型（domain decomposition）の部品として使われることが多い。

7.5 代表的前処理（3）：ILU / IC（不完全分解）

7.5.1 発想：直接法の“近似”を前処理として使う

直接法（LUやCholesky）は強力だが、疎行列では fill-in（非ゼロ増殖）が増えてメモリが厳しくなる。
そこで「fill-inを制限した分解」を作り、前処理として用いるのが

ILU（Incomplete LU）
IC / ICC（Incomplete Cholesky）

である。

7.5.2 レベル・オブ・フィル（fill level）

ILU(k), ICC(k) のように、許すfillのレベルを増やすと一般に

前処理は強くなる（反復回数が減る）
しかし構築コスト・適用コスト・メモリが増える

というトレードオフが生じる。

7.5.3 並列での注意：ILUは“そのまま”では難しいことが多い

ILU/ICはデータ依存の逐次性を含むため、素朴な並列化は難しい場合がある。
実務では

ブロックJacobi（領域分割）＋各ブロックでILU
あるいは並列ILU（外部ライブラリ）などにすることが多い。

7.6 代表的前処理（4）：近似逆行列（Sparse Approximate Inverse, SPAI）

7.6.1 ねらい： $M \approx A^{- 1}$ を直接“疎”に近似する

前処理の理想は $M = A^{- 1}$ だが、逆行列は一般に密になる。
そこで「疎な形を保ったまま逆行列を近似する」考え方がSPAI。

7.6.2 強みと弱み

強み：適用（ $M v$ ）が疎行列–ベクトル積に近く、並列化しやすい場合がある
弱み：構築が重くなりやすく、パラメータ調整が要る

7.7 代表的前処理（5）：分割型（Domain Decomposition / Block preconditioning）

7.7.1 大規模計算で頻出する理由

大規模並列では「全体を一気に解く」よりも「分割して解いて協調させる」方が現実的なことが多い。
分割型前処理は

物理領域を分割（サブドメイン）
サブドメイン内の解法（ILU/IC/直接法/局所多重格子など）を組み合わせる

ことで、並列効率と収束を両立させる設計思想である。

7.8 現代の実務で重要な補足：多重格子（MG/AMG）とGPU対応

教科書の基本としてはILU/ICを学べばよいが、実務（特にPDE由来の大規模問題）では

幾何学的多重格子（GMG）
代数的多重格子（AMG）

が「前処理の王道」になっていることが多い。
さらに近年はGPU利用が一般化し、AMGもGPU対応が進んでいる。
初学者としてはまず「多重格子は“遅い誤差モードを粗い格子で一気に潰す”」という直観を持てば十分である。

7.9 前処理の設計手順

初学者が現場で迷わないために、まずは次の順で考えるとよい。

スケーリング（Jacobi）：係数の桁差が大きいなら必須
SPDなら IC + CG を試す（可能なら）
非対称なら ILU + GMRES(m) / ILU + BiCGStab を比較
PDE由来で巨大なら（A）MG/AMG系や（B）分割型を検討
重要：比較は「反復回数」ではなく $総時間 = 前処理構築 + 反復（適用×回数）$ で判断する

7.10 演習

演習7A：前処理あり/なしで収束が激変する例

Poisson型SPD問題を作り、CGで解く
前処理なし／Jacobi／IC(0)／IC(k) を比較し、
- 反復回数
- 総計算時間
- 最終残差
  をまとめよ。

演習7B：前処理コストも含めた最適点探索

ILUのレベル（k）や、SPAIの許容誤差などパラメータを振って
「前処理構築時間」と「反復時間」を分解し
総時間が最小となる設定を探し、その理由を説明せよ。

演習7C（発展）：大規模並列で“通信が支配的”になる状況を想像せよ

内積（Allreduce）が律速になる状況を文章で説明し、
pipelined Krylov法がなぜ効くのかを、自分の言葉でまとめよ。

第8章　固有値問題（大規模）：Lanczos・Arnoldi・反復固有値解法

8.0 この章の狙い

材料・物理・工学の数値計算では、しばしば「行列のすべての固有値」ではなく、端（最大/最小）や、あるエネルギー近傍の少数の固有値（固有対）だけが欲しくなる。例えば、安定性判定（最小固有値の符号）、振動モード（低周波固有値）、電子状態計算でのスペクトル近傍などが典型である。

この章では、大規模（疎）行列に対して「少数固有値」を求めるための基本を身につける。

固有値問題の種類（対称/非対称、標準/一般化）を整理できる
パワー法・逆反復・Rayleigh商を、数式の意味まで説明できる
Lanczos/Arnoldi（Krylov型固有値解法）が「なぜ大規模向き」か説明できる
シフト＆インバートの発想が「内部固有値」を取り出す鍵であることを理解する
固有値計算が本質的に「線形方程式を解く能力（前処理・反復法）」と結びつくことを理解する

8.1 固有値問題の定義：標準形と一般化形

8.1.1 標準固有値問題

最も基本は次の形である。

A x = λ x

$A \in R^{n \times n}$ （または $C^{n \times n}$ ）は行列
$x \neq 0$ は固有ベクトル
$λ$ は固有値

「 $A$ を作用させても方向が変わらず、スケールだけ変わるベクトル」が固有ベクトルである。

8.1.2 一般化固有値問題（Generalized EVP）

有限要素法などでは次が現れる。

A x = λ B x

ここで $B$ は質量行列のような役割を持つことが多い。

$B$ が対称正定値（SPD）なら、適切な変換で標準形に近い形へ写せる
実務では「 $B$ を明示的に逆行列にしない」形で処理することが多い

8.2 大規模固有値計算で「全部求めない」理由

8.2.1 密行列の全固有値は重い

密行列の全固有値・全固有ベクトル計算（例えばQR法）は、計算量が概ね

O (n^{3})

となり、 $n$ が大きいと現実的でない。

8.2.2 疎行列では「行列–ベクトル積」が主役

疎行列 $A$ は非ゼロ要素数 $nnz$ が $n^{2}$ より遥かに小さい。疎行列–ベクトル積 $y = A x$ は概ね

O (nnz)

で計算できる。したがって大規模では、

行列を密にしない
行列–ベクトル積だけで進む
欲しい固有値だけ取り出すという設計が自然になる。

8.3 Rayleigh商：固有値の「自然な推定量」

ベクトル $x \neq 0$ に対し、Rayleigh商を

ρ (x) = \frac{x^{⊤} A x}{x^{⊤} x}

（複素なら $ρ (x) = \frac{x^{*} A x}{x^{*} x}$ ）と定義する。

8.3.1 直観

もし $x$ が固有ベクトル $x = v$ なら、

ρ (v) = \frac{v^{⊤} (λ v)}{v^{⊤} v} = λ

となる。つまり「固有ベクトルに近いベクトルを入れると、固有値に近い値が返る」。

8.3.2 対称行列での重要性

$A$ が対称なら、Rayleigh商の極値が固有値に対応し、最大固有値・最小固有値と深く結びつく。初学者としては次の感覚を持つと良い：

$ρ (x)$ を大きくする方向へ $x$ を“育てる”と最大固有値に近づく
$ρ (x)$ を小さくする方向へ“育てる”と最小固有値に近づく

8.4 パワー法：最大固有値を当てに行く最初の一歩

8.4.1 基本アイデア

最大固有値（絶対値最大）に対応する固有ベクトル方向は、行列を繰り返し掛けると支配的になることがある。

初期ベクトル $x_{0}$ を与え、反復

x_{k + 1} = \frac{A x_{k}}{∥ A x_{k} ∥}

を行う。

8.4.2 なぜ効くのか

固有ベクトルが基底を作る状況をイメージし、

x_{0} = \sum_{i} c_{i} v_{i}, A v_{i} = λ_{i} v_{i}

と書けるとすると、

A^{k} x_{0} = \sum_{i} c_{i} λ_{i}^{k} v_{i}

となる。もし $| λ_{1} | > | λ_{2} | \geq \dots$ で、かつ $c_{1} \neq 0$ なら、 $λ_{1}^{k}$ の項が支配的になって方向が $v_{1}$ に寄っていく。

8.4.3 限界：固有値の分離が悪いと遅い

$| λ_{2} | / | λ_{1} |$ が 1 に近い（最大固有値が“他と近い”）ほど収束は遅くなる。この「収束を速める」一般手段が、後で出てくるKrylov法（Lanczos/Arnoldi）である。

8.5 逆反復とシフト：欲しい固有値

8.5.1 逆反復（Inverse Iteration）

逆反復は

x_{k + 1} = \frac{A^{- 1} x_{k}}{∥ A^{- 1} x_{k} ∥}

である。これは「 $A^{- 1}$ の最大固有値」を取りに行くので、「 $A$ の最小固有値」を取りやすい。

ただし、実装上は $A^{- 1}$ を作らない。毎回

A y = x_{k}

という線形方程式を解いて $y = A^{- 1} x_{k}$ を得る。

ここで重要：
固有値計算は、結局「線形方程式を何回も解く問題」になり得る。
したがって、前章までのKrylov法・前処理がそのまま効いてくる。

8.5.2 シフト付き逆反復（Shifted Inverse Iteration）

ある $σ$ 近傍の固有値（内部固有値）を狙うには、

(A - σ I)^{- 1}

を使うのが自然である。反復は

(A - σ I) y = x_{k}, x_{k + 1} = \frac{y}{∥ y ∥}

となる。これは「 $λ$ が $σ$ に近いものほど $(λ - σ)^{- 1}$ が大きい」ことを利用している。

8.6 Lanczos法（対称/Hermitian向け）：短い再帰で強い

8.6.1 Krylov部分空間での固有値近似

第6章で学んだKrylov部分空間

K_{m} (A, v_{1}) = span {v_{1}, A v_{1}, \dots, A^{m - 1} v_{1}}

に固有ベクトルの“よい近似”が含まれることが多い。

Lanczos法は、 $A$ が対称（Hermitian）であるとき、この空間の正規直交基底を作りつつ、小さな三重対角行列の固有値問題へ落とし込む。

8.6.2 三項間漸化式（short recurrence）

Lanczosでは、直交化が“うまく簡約”され、概ね

β_{k + 1} v_{k + 1} = A v_{k} - α_{k} v_{k} - β_{k} v_{k - 1}

のような三項間漸化式になる（ $α_{k}, β_{k}$ はスカラー）。この「短い再帰」が、メモリと計算を抑える鍵である。

8.6.3 実務上の注意：直交性の崩れと再直交化

有限精度では $v_{k}$ の直交性が少しずつ崩れ、“ゴースト固有値”のような挙動が出ることがある。実装では

部分再直交化
あるいは暗黙再スタート（implicit restart）などが使われる。

8.7 Arnoldi法（一般行列向け）：GMRESの親戚

8.7.1 何が違うか

非対称（非Hermitian）の場合、Lanczosの簡約は成立しない。そこでArnoldi法を使い、直交基底 ${v_{1}, \dots, v_{m}}$ を作って

A V_{m} = V_{m} H_{m} + h_{m + 1, m} v_{m + 1} e_{m}^{⊤}

の形（ $H_{m}$ は上ヘッセンベルグ行列）を得る。

8.7.2 なぜ重いか

Arnoldiは新しいベクトルを過去すべてに直交化するので、反復が進むほど

内積が増える（並列では通信が増える）
保存するベクトルが増える（メモリ増）という課題が出る。

8.7.3 暗黙再スタート（Implicit Restart）

実務でよく使われるARPACK系の方法は、Arnoldiをそのまま伸ばし続けず、

必要な情報だけ残して
再スタートすることでメモリと計算を制御する（“Implicitly Restarted Arnoldi Method”の考え方）。

8.8 シフト＆インバート（Shift-and-Invert）：内部固有値の定番ルート

端の固有値（最大・最小）はLanczos/Arnoldiで比較的取りやすいが、内部固有値は難しい。そこでスペクトル変換を使う。

8.8.1 変換の定義

T = (A - σ I)^{- 1}

の固有値問題を考えると、

A x = λ x \Rightarrow T x = \frac{1}{λ - σ} x

となる。

8.8.2 何が嬉しいか

$λ$ が $σ$ に近いほど $| (λ - σ)^{- 1} |$ が大きくなり、“端の固有値”のように扱える。つまり、内部固有値問題を端の固有値問題へ変換できる。

8.8.3 代償：毎反復で線形方程式を解く必要

ただし $T$ を作用させるには

(A - σ I) y = x

を解かないといけない。ここで

直接法（疎直接法）
反復法＋前処理の選択が支配的になる。

結論：
「固有値を解く」＝「線形方程式を（何度も）解く」
したがって、線形ソルバと前処理が性能の中心になる。

8.9 大規模固有値計算の現代的な実務感

現場では自作より、信頼できる反復固有値ソルバを使うことが多い。例として、並列疎行列の固有値問題に特化した枠組み（例：SLEPcのEPS）は、

Lanczos/Arnoldi系
シフト＆インバートなどのスペクトル変換
反復法の再スタート
収束判定・誤差評価を統一インターフェースで扱える。

また、内部固有値を大量に求めたい場合は、スペクトルを区間に分けて並列に探索する「スペクトルスライシング」の系統も重要になる（発展）。

8.10 演習

演習8A：最大固有値の推定（パワー法）

疎行列 $A$ を用意し、パワー法で $λ_{max}$ を推定せよ。
反復ごとに Rayleigh商 $ρ (x_{k})$ を計算し、推定値がどう収束するかプロットせよ。
$| λ_{2} | / | λ_{1} |$ が 1 に近い例（固有値が近い例）を作り、収束が遅くなることを確認せよ。

演習8B：疎行列の端の固有値を少数だけ求める

Lanczos/Arnoldi系の関数（例：Pythonの疎固有値ソルバ）で、最大/最小の固有値を数個だけ求めよ。
反復回数・時間・精度を比較し、「全部求める」方法が不利な理由を文章で説明せよ。

演習8C（発展）：内部固有値とシフト＆インバート

$σ$ を選び、 $(A - σ I)^{- 1}$ を使う発想を、式変形で説明せよ。
そのとき「線形方程式を解くコスト」が支配的になる理由を、前処理の観点も含めてまとめよ。

第9章　非線形方程式：Newton法と大規模化（Newton-Krylov）

9.0 この章の狙い

大規模計算では、最終的に解きたい問題が非線形であることが多い（非線形PDE、相平衡、最適化のKKT条件など）。しかし重要な事実は次である：

非線形問題を解く“中心”は、各ステップで現れる線形方程式を速く安定に解く能力である。

この章では

Newton法の導出（線形化）
ヤコビアン（Jacobian）の意味
収束を保証する工夫（線形探索・信頼領域）
大規模向けの Newton-Krylov / Jacobian-free Newton-Krylov（JFNK）を、数式のギャップが出ないように丁寧に学ぶ。

9.1 非線形方程式の定義：残差形式 $F (u) = 0$

非線形方程式を、未知ベクトル $u \in R^{n}$ に対して

F (u) = 0

という形で書く。ここで $F : R^{n} \to R^{n}$ は残差ベクトル（各式の“ずれ”）である。

例：離散化した非線形Poisson型

連続系： $- Δ u + u^{3} = f$
離散化すると： $K u + u^{\circ 3} - f = 0$
- $K$ ：ラプラシアン由来の疎行列
- $u^{\circ 3}$ ：成分ごとの3乗（ベクトル関数）

このとき

F (u) = K u + u^{\circ 3} - f

が残差である。

9.2 Newton法の導出：Taylor展開による線形化

9.2.1 ヤコビアンの定義

ヤコビアン $J (u)$ は

J (u) = \frac{\partial F (u)}{\partial u}

で定義される $n \times n$ 行列で、成分表示では

J_{i j} (u) = \frac{\partial F_{i} (u)}{\partial u_{j}}

である。

9.2.2 1次Taylor展開

現在の近似解 $u_{k}$ の近傍で

F (u_{k} + δ) \approx F (u_{k}) + J (u_{k}) δ

と近似する。

解きたいのは $F (u_{k} + δ) = 0$ なので、近似式に代入して

F (u_{k}) + J (u_{k}) δ \approx 0

すなわち

J (u_{k}) δ_{k} = - F (u_{k})

という線形方程式を解いて更新する：

u_{k + 1} = u_{k} + δ_{k}

これがNewton法である。

9.2.3 重要な直観：Newton法は「線形方程式を逐次解く」方法

Newton法の1反復は

残差 $F (u_{k})$ を評価する
ヤコビアン $J (u_{k})$ を用意する（または作用を計算できるようにする）
線形方程式 $J (u_{k}) δ = - F (u_{k})$ を解く
更新 $u_{k + 1} = u_{k} + δ$

で構成される。大規模化すると、コストは 3（線形ソルバ）が支配することが多い。

9.3 収束が「うまくいかない」理由：局所法の限界

Newton法は、初期値が十分近いときには非常に速い（しばしば2次収束）一方で、初期値が悪いと

発散する
変な方向へ飛ぶ
途中で停滞することがある。

そこで必要になるのが「大域化（globalization）」である。

9.4 大域化（Globalization）：線形探索と信頼領域

9.4.1 線形探索（Line Search）

Newton方向 $δ_{k}$ をそのまま足すのではなく、係数 $α_{k} \in (0, 1]$ を入れて

u_{k + 1} = u_{k} + α_{k} δ_{k}

とする。

$α_{k}$ は、残差ノルム $∥ F (u) ∥$ が十分減るように選ぶ（例：バックトラッキング）。直観は「大きく踏み出して失敗するくらいなら、少しずつ進む」。

9.4.2 信頼領域（Trust Region）

「線形化が信頼できる範囲（半径）」を設け、その中で最適な更新を選ぶ考え方。線形化が当たっていれば大きく進め、当たっていなければ半径を縮める。

9.5 Inexact Newton：線形方程式を“厳密に”解かない

大規模では、毎回 $J δ = - F$ を高精度で解くとコストが重い。そこで「ある程度の精度で線形ソルバを止める」Inexact Newton を用いる。

9.5.1 典型的な停止条件（フォーシングターム）

線形ソルバで得られた近似解 $δ_{k}$ が

∥ J (u_{k}) δ_{k} + F (u_{k}) ∥ \leq η_{k} ∥ F (u_{k}) ∥

を満たすところで止める。ここで $η_{k} \in (0, 1)$ がフォーシングタームであり、

解から遠いとき： $η_{k}$ を大きめ（粗く解く）
解に近いとき： $η_{k}$ を小さめ（精密に解く）とすることで総計算量を減らす。

大事な直観：
まだ非線形残差が大きい段階で、線形方程式を過剰に高精度で解くのは無駄になりやすい。

9.6 Newton-Krylov：Newtonの中でKrylov法を回す

Newtonステップの線形方程式

J (u_{k}) δ = - F (u_{k})

を、直接法ではなく Krylov反復法（GMRESなど）で解く構成を Newton-Krylov と呼ぶ。

9.6.1 どのKrylov法を使うか

ヤコビアン $J$ は一般に非対称になり得るため、実務では

GMRES（またはFGMRES）
BiCGStab などがよく使われる。

9.6.2 前処理はどうするか（ここが実戦の核心）

Krylov法の収束は前処理に強く依存する。Newton-Krylovでも同様で、前処理 $M^{- 1}$ を用いて

M^{- 1} J (u_{k}) δ = - M^{- 1} F (u_{k})

の形で解くのが典型である。

前処理の設計方針（初学者向けの指針）：

まずはヤコビアンの近似（ブロック構造、拡散項だけ、線形化の一部）で $M$ を作る
$M$ は“厳密でなくてよい”が、“適用が安い”必要がある
非線形反復で $J$ が変化するので、前処理の更新頻度もチューニング対象

9.7 Jacobian-free Newton-Krylov（JFNK）：行列を作らずに $J v$ を得る

大規模ではヤコビアン行列 $J$ を組み立て・保存すること自体が重い場合がある。そこで、Krylov法が必要とするのは本質的に「行列そのもの」ではなく「行列–ベクトル積 $J v$ 」である点を利用する。

9.7.1 有限差分近似による $J v$

J (u) v \approx \frac{F (u + ε v) - F (u)}{ε}

ここで $ε$ は十分小さいスカラー。この式は、 $F (u)$ の方向微分の定義に一致し、Taylor展開から自然に出てくる。

9.7.2 何が嬉しいか

ヤコビアンの行列を明示的に作らない（メモリ節約）
実装が「残差評価 $F (u)$ さえ書ければ」成立しやすい

9.7.3 注意：前処理は別途必要になりやすい

JFNKは $J v$ を行列なしで計算できるが、前処理 $M^{- 1}$ まで自動で良くなるわけではない。むしろ実務では

$J$ は行列フリー
前処理 $M$ は“行列ベースの近似”で持つという混合構成がよく使われる。

9.8 非線形Poissonのミニ例：Newton + Krylov の流れ

例として

- Δ u + u^{3} = f

を離散化して

F (u) = K u + u^{\circ 3} - f = 0

とする。

9.8.1 ヤコビアン

成分ごとに微分すると

J (u) = \frac{\partial (K u)}{\partial u} + \frac{\partial (u^{\circ 3})}{\partial u} = K + diag (3 u^{\circ 2})

ここで $diag (3 u^{\circ 2})$ は、対角に $3 u_{i}^{2}$ を持つ対角行列。

9.8.2 Newtonステップ

各反復で

(K + diag (3 u_{k}^{\circ 2})) δ_{k} = - (K u_{k} + u_{k}^{\circ 3} - f)

を解き、 $u_{k + 1} = u_{k} + α_{k} δ_{k}$ （線形探索あり）で更新する。

9.8.3 大規模では何が支配的か

行列–ベクトル積（疎行列 $K$ と対角行列）
Krylov法の内積（並列なら通信）
前処理の構築・適用が支配的になりやすい。

9.9 現代的な実務感

非線形ソルバの枠組み（例：SNESなど）は、

Newton（線形探索・信頼領域）
Inexact Newton（フォーシングターム）
Jacobian-free（行列フリー）
Newton-Krylov（内部でKrylov法）を統一的に扱う設計になっていることが多い。

「非線形を解く」実務は、しばしば

残差評価 $F (u)$ の実装
物理に即した前処理 $M$ の設計
停止条件（非線形/線形）の整合の3点セットで考えるのがコツである。

9.10 演習

演習9A：非線形PoissonでNewtonの流れを実装

1次元でもよいので格子を切り、 $- u^{″} + u^{3} = f$ を差分で離散化せよ。
$F (u)$ と $J (u)$ を実装し、Newton法で解け。
初期値を変えて、収束/発散の違いを観察し、線形探索の効果を確認せよ。

演習9B：Newton + Krylov（Newton-Krylov）で線形ソルバを差し替え比較

Newtonステップの線形方程式を、直接法ではなくGMRESで解くようにせよ。
前処理なし/あり（Jacobiなど簡単で良い）を比較し、反復回数と総時間がどう変わるかまとめよ。
Inexact Newton の停止条件（ $η_{k}$ ）を変えて、総反復コストがどう変わるか観察せよ。

演習9C（発展）：JFNK（行列フリー）を試す

$J v \approx (F (u + ε v) - F (u)) / ε$ を使ってGMRESに渡せ（行列を組み立てない）。
$ε$ の値を変えると挙動がどう変わるか観察し、丸め誤差（第2章）との関係を考察せよ。
行列フリーでも前処理が必要になる理由を文章で説明せよ。

第10章　最適化（大規模）：勾配法・共役勾配・準Newton・確率的手法

10.0 この章の狙い

逆問題・機械学習・材料設計（パラメータ同定、設計変数探索、PDE制約付き最適化など）では、最終的に

min_{x} f (x)

という最適化問題に帰着することが多い。本章では、大規模（高次元）で現実に動く最適化を理解するために、次を到達目標とする。

目的関数・勾配・ヘッセ行列（2階微分）の役割を、式の意味まで説明できる
大規模では「行列を作る」より「作用（ベクトルとの積）」が重要になる理由を説明できる
勾配法・共役勾配（最適化版）・準Newton（L-BFGS）・近接勾配（ISTA/FISTA）・SGD系を使い分ける判断軸を持つ
停止条件（いつ止めるか）とスケーリング（単位・正規化）が性能を左右することを理解する

10.1 最適化問題の基本形：何を最小化しているのか

10.1.1 変数・目的関数・制約

最も基本は無制約最適化：

min_{x \in R^{n}} f (x)

$x$ ：設計変数／未知パラメータ（材料組成、モデル係数、メッシュ上の自由度など）
$f (x)$ ：目的関数（損失、誤差、エネルギー、コスト）

制約がある場合は

min_{x} f (x) s.t. c (x) = 0, g (x) \leq 0

のようになるが、本章ではまず無制約〜単純制約（境界付き）を中心に扱う。

10.1.2 勾配とヘッセ行列（微分の意味）

勾配（gradient）：

\nabla f (x) = {[\frac{\partial f}{\partial x_{1}}, \dots, \frac{\partial f}{\partial x_{n}}]}^{⊤}

勾配は「その点で最も増える方向」。したがって $- \nabla f (x)$ は「最も減る方向」。

ヘッセ行列（Hessian）：

\nabla^{2} f (x) = {[\frac{\partial^{2} f}{\partial x_{i} \partial x_{j}}]}_{i, j}

ヘッセ行列は「曲率（どれだけ曲がっているか）」を表し、Newton法・準Newton法の要になる。

10.1.3 例：最小二乗（データ同化・回帰の基本）

観測 $y \in R^{m}$ 、モデル $A x$ （ $A \in R^{m \times n}$ ）に対し

min_{x} \frac{1}{2} ∥ A x - y ∥_{2}^{2}

これは線形最小二乗。目的関数を

f (x) = \frac{1}{2} ∥ A x - y ∥_{2}^{2} = \frac{1}{2} (A x - y)^{⊤} (A x - y)

とおくと、

\nabla f (x) = A^{⊤} (A x - y), \nabla^{2} f (x) = A^{⊤} A

となる。ここから

勾配が「残差 $A x - y$ を $A^{⊤}$ で引き戻したもの」
ヘッセ行列が $A^{⊤} A$ （しばしば大規模・疎・条件が悪い）であることが分かる。

10.2 勾配降下法（Gradient Descent）：最も基本で大規模向き

10.2.1 更新式（アルゴリズムの核）

勾配降下法は

x_{k + 1} = x_{k} - α_{k} \nabla f (x_{k})

$x_{k}$ ：反復 $k$ 回目の解の近似
$α_{k} > 0$ ：ステップサイズ（学習率）
$\nabla f (x_{k})$ ：勾配

ポイントは、必要なのが基本的に「勾配ベクトル」だけで、ヘッセ行列を明示的に作らない点にある。大規模ではここが非常に重要になる。

10.2.2 ステップサイズがなぜ重要か（壊れる典型原因）

直観的には

$α_{k}$ が大きすぎる：最小値を飛び越えて発散しやすい
$α_{k}$ が小さすぎる：収束が遅すぎる

よく使う設計として

固定ステップ（簡単だが調整が必要）
直線探索（line search）： $α_{k}$ をその場で決めるがある。直線探索は「今の方向にどれだけ進むのが安全か」を見積もる仕組みで、非凸でも暴走を抑える。

10.2.3 「条件数」と収束の関係（大規模で効く基本概念）

二次関数

f (x) = \frac{1}{2} x^{⊤} Q x - b^{⊤} x, Q ≻ 0

（ $Q$ は対称正定値）では、最適解は $Q x = b$ の解である。勾配降下法の収束は概ね

κ (Q) = \frac{λ_{max} (Q)}{λ_{min} (Q)}

（条件数）が大きいほど遅くなる。つまり「曲率が方向によって極端に違う（細長い谷）」ほど難しい。

この事実は次の大きな教訓に繋がる：

大規模最適化の性能は、ステップサイズの工夫だけでなく、前処理（スケーリング）で条件数を下げることが決定的に重要。

10.3 共役勾配法（最適化版）：二次問題なら“最短級”に強い

10.3.1 線形方程式のCGと最適化のCGは同じ根っこ

二次関数最小化は線形方程式 $Q x = b$ に等価で、ここで $Q$ がSPDなら CG（共役勾配法）が使える。

CGは「勾配法より賢い方向」を作っていく。直観的には

勾配法：毎回「一番下りやすい方向」に進むが、ジグザグしやすい
CG：過去の情報を使ってジグザグを抑え、より少ない反復で進む

10.3.2 何が“大規模向き”か

CGの主要計算は

行列–ベクトル積 $Q p$
内積（スカラー）
ベクトル更新（axpy）であり、疎行列なら $Q p$ が $O (nnz)$ で済む。

10.3.3 前処理付きCG（PCG）

条件数が大きいとCGでも遅くなるため、前処理 $M^{- 1}$ を用いて

M^{- 1} Q x = M^{- 1} b

の形で解く（前処理付きCG）。ここで $M \approx Q$ を「解きやすい形」で用意できると収束が劇的に改善する。 ※前処理の設計思想は第7章と接続する。

10.4 準Newton法：曲率を“賢く近似”して速くする

10.4.1 Newton法の理想と現実

Newton法は

x_{k + 1} = x_{k} - {(\nabla^{2} f (x_{k}))}^{- 1} \nabla f (x_{k})

で、局所的に非常に速い（うまくいけば二次収束）。しかし大規模では

$\nabla^{2} f$ の構築が重い
逆行列は作らずとも、線形方程式を解くのが重いという問題がある。

10.4.2 準Newton（BFGS）とL-BFGS（大規模の定番）

準Newtonはヘッセ行列の逆 ${(\nabla^{2} f)}^{- 1}$ を逐次更新で近似する。特に L-BFGS は「limited-memory（少量メモリ）」で、過去の更新情報を $m$ ステップ分だけ保持する。

更新に使う差分は

s_{k} = x_{k + 1} - x_{k}, y_{k} = \nabla f (x_{k + 1}) - \nabla f (x_{k})

で、これらのペアを少数保持し、二重ループ（two-loop recursion）で近似的に $(H_{k}) \nabla f$ を計算する。

行列 $H_{k}$ を明示的に保持しない
主要計算はベクトルの内積・線形結合という点で大規模向き。

10.4.3 実務での使いどころ

中規模〜大規模で高精度に詰めたい：L-BFGSが強い
目的関数評価が重い（PDEを解く必要がある等）：反復回数が少ない手法が有利になりやすい
ただし非凸では直線探索や初期値依存が重要（暴走防止が必要）

10.5 近接勾配（Proximal Gradient）：正則化・非滑らか最適化の入口

10.5.1 合成最適化（smooth + nonsmooth）

次の形が非常によく出る：

min_{x} F (x) = f (x) + g (x)

$f$ ：滑らか（微分可能）で勾配が計算できる（例：最小二乗）
$g$ ：非滑らか（例： $ℓ_{1}$ 正則化、制約の指示関数）

例：LASSO（ $ℓ_{1}$ 正則化）

min_{x} \frac{1}{2} ∥ A x - y ∥_{2}^{2} + λ ∥ x ∥_{1}

ここで $∥ x ∥_{1} = \sum_{i} | x_{i} |$ は微分不可能（0で尖っている）なので、普通の勾配降下がそのまま使いづらい。

10.5.2 近接作用素（prox）の定義

近接作用素を

{prox}_{λ g} (v) = \arg min_{x} (g (x) + \frac{1}{2 λ} ∥ x - v ∥_{2}^{2})

と定義する。意味は

$v$ の近くで
$g (x)$ を小さくするように
“ほどよく”引き寄せた点を返す操作である。

10.5.3 ISTA（基本）とFISTA（加速）

近接勾配法（ISTA）は

x_{k + 1} = {prox}_{λ g} (x_{k} - λ \nabla f (x_{k}))

と書ける。

まず滑らかな部分 $f$ に対して勾配ステップ
次に非滑らかな部分 $g$ をproxで処理という分離（splitting）になっている。

FISTAはこれを加速し、（凸・滑らか条件のもとで）目的関数値の収束が概ね

ISTA： $O (1 / k)$
FISTA： $O (1 / k^{2})$ と速くなることが知られている。

初学者の直観：
“普通の勾配法 + 正則化の処理”を上手に組み合わせると、大規模でも回る強いアルゴリズムになる。

10.6 確率的手法（SGD系）：データが巨大なときの標準

10.6.1 有限和（finite-sum）最適化

機械学習・統計では

min_{x} \frac{1}{N} \sum_{i = 1}^{N} ℓ_{i} (x)

の形が多い。全データで勾配を計算するのは重いので、ミニバッチ $B$ を使い

g_{k} \approx \nabla f (x_{k}) = \frac{1}{| B |} \sum_{i \in B} \nabla ℓ_{i} (x_{k})

という“雑音入り勾配”で更新するのがSGDである。

10.6.2 代表的な改良：Momentum・Adam・AdamW

Momentum：過去の更新方向を蓄積して谷のジグザグを抑える
Adam：勾配の一次・二次統計（移動平均）を使って座標ごとに学習率を調整
AdamW：重み減衰（weight decay）を勾配項から分離（decoupled）して扱う

ここで重要なのは、実装上「L2正則化」と「weight decay」が同じだと思い込むと、適応的手法（Adam系）では振る舞いがズレることがある、という点である。AdamWはこの点を整理した代表例で、現代の大規模学習で標準的に使われることが多い。

10.6.3 大規模での停止条件

確率的手法では勾配がノイズを含むため、単純に $∥ \nabla f (x) ∥$ が0に近いかだけで止めるのは難しい。典型的には

目的関数（訓練損失）の移動平均が改善しない
検証データの指標が改善しない（早期終了）
学習率を段階的に下げた後も改善が止まるなどの基準を併用する。

10.7 スケーリングと正則化：最適化を“壊さない”設計

10.7.1 スケーリング（単位・正規化）の重要性

変数のスケールがバラバラだと、等方的なステップ（勾配法）が極端に非効率になる。例：ある成分は $10^{- 6}$ の変化が重要、別の成分は $10^{2}$ の変化が重要、など。

対策は

変数の無次元化（物理なら特に重要）
入力・特徴量の標準化（平均0・分散1など）
前処理（対角スケーリング、近似ヘッセの利用）である。これは「条件数を下げる」実務そのもの。

10.7.2 正則化（Regularization）

最小二乗にL2正則化を足すと

min_{x} \frac{1}{2} ∥ A x - y ∥_{2}^{2} + \frac{λ}{2} ∥ x ∥_{2}^{2}

となり、これはRidge回帰に相当する。勾配は

\nabla f (x) = A^{⊤} (A x - y) + λ x

で、ヘッセは

\nabla^{2} f (x) = A^{⊤} A + λ I

となる。ここで $λ I$ を足すと最小固有値が上がりやすく、条件が改善することがある（ただし過大な $λ$ は過平滑化＝精度低下）。

10.8 演習

演習10A：L2正則化付き最小二乗（Ridge）を解く

ランダムな疎行列 $A$ とデータ $y$ を作り、

min_{x} \frac{1}{2} ∥ A x - y ∥_{2}^{2} + \frac{λ}{2} ∥ x ∥_{2}^{2}

を（i）勾配降下、（ii）CG（正規方程式 $ (A^\top A+\lambda I)x=A^\top y$ を解く）、（iii）L-BFGS のいずれかで解け。
2. 反復回数・時間・最終誤差を比較し、「どこが支配的か」を説明せよ。

演習10B：条件数と収束の関係を観察する

二次関数 $f (x) = \frac{1}{2} x^{⊤} Q x$ を用意し、 $Q$ の固有値の比（条件数）を変える。
勾配法・CGでの収束（残差や目的関数値）がどう変わるかをプロットし、理由を言葉で説明せよ。

演習10C：ミニバッチ（SGD系）の挙動比較

有限和目的関数（回帰・分類の簡単なもの）で、バッチサイズを変える。
目的関数の減り方、振動（ノイズ）、最終精度を比較し、停止条件の難しさをまとめよ。
可能なら Adam と AdamW を比較し、weight decay の扱いの違いが何に効くか考察せよ。

第11章　PDE離散化：FDM/FEMの入口と線形系への帰着

11.0 この章の狙い

偏微分方程式（PDE）は、物理現象（拡散・弾性・電磁場・流体など）を記述する標準言語である。大規模計算では、PDEを

PDE（連続） \Rightarrow 離散化（格子/メッシュ） \Rightarrow 疎行列の連立方程式 \Rightarrow ソルバ

という一本の流れで理解することが極めて重要になる。本章では次を目標とする。

FDM（有限差分）とFEM（有限要素）の違いを「何を近似しているか」で説明できる
Poisson方程式を例に、PDEが疎行列 $A x = b$ に落ちる導出を自分で追える
境界条件、弱形式、メッシュ、離散化誤差（精度）と安定性の基本用語を理解する
実装で何がボトルネックになり、なぜ現代は“行列を作らない（matrix-free）”発想が重要なのかを掴む

11.1 例題としてのPoisson方程式

領域 $Ω$ で次を考える：

- Δ u = f in Ω

ここで

$u (x)$ ：未知関数（温度、電位、変位の成分など）
$f (x)$ ：ソース項
$Δ = \nabla \cdot \nabla$ ：ラプラシアン

境界条件として典型的に

Dirichlet： $u = g$ on $\partial Ω$ （値を固定）
Neumann： $\frac{\partial u}{\partial n} = h$ on $\partial Ω$ （法線微分＝流束を固定）を扱う。

11.2 離散化で必ず出てくる3語：整合性・安定性・収束

離散化した解を $u_{h}$ （格子幅 $h$ の解）とするとき、基本は次の三点セットで考える。

整合性（consistency）：連続方程式を離散近似したとき、 $h \to 0$ で元のPDEに近づくか
安定性（stability）：誤差や擾乱が計算で暴走しないか
収束（convergence）： $u_{h} \to u$ （真の解）に近づくか

初学者はまず「整合でも安定でないと壊れる」「安定でも整合でないと間違う」という直観を持つとよい。

11.3 有限差分法（FDM）：微分をテイラー展開で置き換える

11.3.1 1次元Poissonの導出

区間 $[0, 1]$ で

- u^{″} (x) = f (x), u (0) = u (1) = 0

を考える。格子点 $x_{i} = i h$ （ $h = 1 / (n + 1)$ ）を置き、未知を $u_{i} \approx u (x_{i})$ とする。

テイラー展開から

u^{″} (x_{i}) \approx \frac{u_{i - 1} - 2 u_{i} + u_{i + 1}}{h^{2}}

よって

- \frac{u_{i - 1} - 2 u_{i} + u_{i + 1}}{h^{2}} = f_{i}

両辺に $h^{2}$ を掛けて

- u_{i - 1} + 2 u_{i} - u_{i + 1} = h^{2} f_{i}

これは $i = 1, \dots, n$ について成り立つので、行列でまとめると

A u = b

$u = [u_{1}, \dots, u_{n}]^{⊤}$
$b = h^{2} [f_{1}, \dots, f_{n}]^{⊤}$
$A$ は三重対角の疎行列（対称正定値）

11.3.2 2次元Poisson（5点差分）への拡張

2次元格子では

- Δ u \approx - \frac{u_{i + 1, j} + u_{i - 1, j} + u_{i, j + 1} + u_{i, j - 1} - 4 u_{i, j}}{h^{2}}

となり、非ゼロが各行あたり最大5個程度の疎行列が出る（境界を除く）。

11.3.3 FDMの特徴（向き・不向き）

長所：導出が直感的、実装が簡単、構造格子に強い
短所：複雑形状・局所精密化に弱い（格子生成が難しい）
大規模化：結局は疎行列 $A$ を解く問題なので、第6〜7章のKrylov＋前処理が本体になる

11.4 有限要素法（FEM）：弱形式で“積分方程式”にしてから離散化

FEMは「微分方程式を、そのまま差分で置換する」のではなく、

弱形式（weak form）に変換して
関数空間の中で近似する
という流れを取る。これが複雑形状に強い理由の中心である。

11.4.1 弱形式の導出（Poisson）

- Δ u = f in Ω, u = g on \partial Ω

を考える。テスト関数 $v$ （任意の“揺らぎ”）を掛けて積分：

\int_{Ω} (- Δ u) v d x = \int_{Ω} f v d x

左辺を部分積分（グリーンの公式）すると

\int_{Ω} \nabla u \cdot \nabla v d x - \int_{\partial Ω} \frac{\partial u}{\partial n} v d s = \int_{Ω} f v d x

Dirichlet境界が強制される設定では、境界項の扱いを整理して最終的に典型形：

（弱形式） \int_{Ω} \nabla u \cdot \nabla v d x = \int_{Ω} f v d x (\forall v \in V_{0})

の形になる（ $V_{0}$ は境界で0になるテスト関数空間、という意味だと思えばよい）。

初学者の直観：
微分（局所の情報）を、積分（全体の平均的な情報）に変換してから近似するので、複雑形状でも扱いやすい。

11.4.2 基底関数による離散化：疎行列が出る理由

近似解を

u_{h} (x) = \sum_{j = 1}^{n} U_{j} ϕ_{j} (x)

と展開する。ここで

$ϕ_{j}$ ：メッシュ上の基底関数（要素ごとに局所的に非ゼロ）
$U_{j}$ ：未知係数（自由度）

テスト関数も $v = ϕ_{i}$ と取る（Galerkin法）と、

\sum_{j = 1}^{n} U_{j} \int_{Ω} \nabla ϕ_{j} \cdot \nabla ϕ_{i} d x = \int_{Ω} f ϕ_{i} d x

よって

K U = F

剛性行列（stiffness matrix）：

K_{i j} = \int_{Ω} \nabla ϕ_{j} \cdot \nabla ϕ_{i} d x

右辺：

F_{i} = \int_{Ω} f ϕ_{i} d x

ここで $ϕ_{i}$ は局所的にしか非ゼロでないため、 $K_{i j} \neq 0$ となる $i, j$ の組はごく一部に限られ、結果として疎行列が自然に生まれる。

11.4.3 境界条件の扱い（入口レベルの整理）

Dirichlet（値固定）：自由度を固定して未知から外す（強制法）、あるいは弱く課す方法もある
Neumann（流束固定）：弱形式の境界項として自然に入ることが多い

11.5 メッシュと要素：h法・p法・適応の考え方

11.5.1 要素の種類

2D：三角形（triangle）、四角形（quad）
3D：四面体（tet）、六面体（hex）

幾何が複雑なら三角形/四面体が扱いやすく、構造的な格子なら四角形/六面体が強いことが多い。

h法：メッシュを細かくする（ $h ↓$ ）
p法：要素内の多項式次数を上げる（ $p ↑$ ）
hp法：両方使う

高次（大きな $p$ ）では、行列組み立てよりも「作用を高速に計算する（matrix-free）」が効いてくる場面が増える。

11.6 現代の大規模FEMで重要な発想：行列を“作らない”・GPUに載せる

大規模化すると、単に疎行列を作って解く以外に、

行列の組み立て（assembly）のコスト
メモリ帯域（メモリから読む量）の制約が支配的になることがある。

そこで近年の重要トレンドとして

partial assembly
matrix-free operator evaluation（行列を明示せずに $y = A x$ を計算）
GPUを意識したデータ配置が広く使われている。

初学者はここで「第5章で出た“行列–ベクトル積が中心”が、PDE離散化でもそのまま効いてくる」ことを掴むとよい。

11.7 PDE → 線形系（/非線形系）→ ソルバ：実装の全体像

離散化して得るのは典型的に

定常線形： $A U = b$
定常非線形： $F (U) = 0$ （第9章のNewton-Krylovへ）
時間発展： $M \dot{U} + K U = F$ （質量行列 $M$ 、剛性 $K$ ）

したがって、PDE計算の実務フローは次の形で整理できる。

PDEと境界条件を決める
メッシュ（格子）を作る
離散化（FDM or FEM）して疎行列・残差を作る（または matrix-free で作用だけ作る）
連立方程式を解く（Krylov＋前処理、または疎直接法）
誤差評価（メッシュ収束、残差、物理量の収束）を行う

11.8 演習

演習11A：Poisson方程式をFDMで離散化して解く

1Dまたは2DでPoissonを差分化し、疎行列 $A$ と右辺 $b$ を構成せよ。
CG（必要なら前処理付き）で解き、反復回数・時間を記録せよ。
格子を細かくし、解がどう変化するか（収束）を確認せよ。

演習11B：Poisson方程式をFEM（弱形式）で組み立てて解く（可能なら）

弱形式から $K_{i j}$ と $F_{i}$ を書き下し、一次要素で実装せよ。
メッシュを変えて非ゼロ構造（疎パターン）を可視化せよ。
境界条件の実装方法（強制 or 弱い課し方）の違いを文章で説明せよ。

演習11C（発展）：matrix-freeの発想を言葉で説明する

「行列を組まずに $y = A x$ を計算できる」理由を、要素積分と局所計算の観点からまとめよ。
行列組み立てと比較して、どのコスト（計算量/メモリ/通信）が支配的になり得るかを整理せよ。

第13章　並列計算と実務：性能評価・再現性・ソフトウェア設計

この章の狙い

この章は「動くコード」から「研究で使えるコード」へ進むための実務章である。
大規模計算では、数値解法そのものの良さに加えて、

並列化して速くなるか（スケーリング）
どこが遅いのか（プロファイル）
結果が再現できるか（再現性）
将来の自分や共同研究者が読めるか（設計・テスト）

が研究成果の信頼性を左右する。ここでは“定石”を体系化する。

13.1 並列化の基本パターン：何を分けるか

13.1.1 並列計算の2つの世界：共有メモリと分散メモリ

共有メモリ（shared memory）：同じメモリ空間をスレッドが共有（例：OpenMP）
分散メモリ（distributed memory）：プロセスごとにメモリが分かれ、通信が必要（例：MPI）

大規模HPCでは「ノード内は共有（OpenMP）、ノード間は分散（MPI）」のハイブリッドが典型である。

13.1.2 データ並列・領域分割・タスク並列

データ並列：同じ処理をデータの分割に対して適用（最頻出）
領域分割（domain decomposition）：空間領域を分割し、境界で“のりしろ（halo）通信”
タスク並列：異なる仕事（例：I/O、前処理構築、可視化）を同時に進める

PDEや疎行列計算は多くが領域分割で自然に並列化できるが、通信が性能の鍵になる。

13.2 通信コスト：なぜ“速くならない”のかを数式で掴む

13.2.1 $α$ - $β$ モデル

通信時間を

T_{comm} (n) \approx α + β n

と近似することが多い。

$α$ ：レイテンシ（送受信開始の固定コスト）
$β$ ：1バイトあたりの時間（帯域の逆数）
$n$ ：送るデータ量（バイト）

小さいメッセージを頻繁に送ると $α$ が支配的になり、帯域が活かせない。
“まとめて送る”“回数を減らす”“計算と重ねる（overlap）”が重要になる。

13.2.2 境界通信（halo exchange）が支配的になる理由

領域分割では、各プロセスは自分の領域内部の計算（計算量）に加え、境界の更新（通信）を行う。
プロセス数を増やすほど、1プロセスあたりの領域は小さくなるので

計算量（体積） $\propto$ 小さくなる
通信量（表面） $\propto$ 相対的に減りにくいとなり、通信が相対的に支配的になる（“表面／体積問題”）。

13.3 強スケーリングと弱スケーリング：測り方を間違えない

13.3.1 強スケーリング（Strong Scaling）

問題サイズ固定でプロセス数 $p$ を増やす。
スピードアップを

S (p) = \frac{T (1)}{T (p)}

効率を

E (p) = \frac{S (p)}{p} = \frac{T (1)}{p T (p)}

で定義する。理想は $S (p) = p$ （完全並列）だが、通信・同期・負荷不均衡で崩れる。

13.3.2 弱スケーリング（Weak Scaling）

1プロセスあたりの問題サイズを一定にし、 $p$ に応じて全体問題を大きくする。
理想は $T (p)$ がほぼ一定。実務では「どの程度までサイズを伸ばせるか」の評価に向く。

13.3.3 AmdahlとGustafson（何が限界を決めるか）

Amdahl：固定問題の強スケーリング限界（逐次割合が効く）
Gustafson：問題を大きくする弱スケーリングの見方（“大きい問題なら並列が効く”）

どちらも「何を固定して議論しているか」を意識する。

13.4 ルーフライン（Roofline）：計算かメモリかを一瞬で分類する

13.4.1 算術強度（Arithmetic Intensity）

I = \frac{FLOPs}{Bytes}

を算術強度という。

$I$ が小さい：メモリからデータを運ぶコストが支配（メモリ律速）
$I$ が大きい：演算器が支配（計算律速）

疎行列-ベクトル積（SpMV）は典型的に $I$ が小さく、メモリ律速になりやすい。
この分類ができると、「SIMD最適化より、データ局所性や圧縮形式の方が効く」など、方針が明確になる。

13.5 プロファイルの作法：最短で“遅い理由”に到達する

13.5.1 プロファイルは3段階でやる

計測（measurement）：どこに時間が使われているか（関数別、区間別）
原因仮説（hypothesis）：計算律速か、メモリ律速か、通信律速か
検証（validation）：変更→再測定（必ずベースラインと比較）

13.5.2 計測の種類：プロファイルとトレース

プロファイル（統計・集計）：軽量で全体像が掴みやすい
トレース（タイムライン）：詳細だがデータ量が大きい、オーバーヘッドも増えやすい

“まずプロファイル、必要な部分だけトレース”が実務の定石。

13.5.3 代表的ツール（現代の標準セット）

CPU側：ホットスポット、スレッド効率、メモリ帯域
MPI/OpenMP混在：通信時間・同期・ロードバランス
GPU：カーネル時間、転送、占有率、メモリ効率

研究室レベルでも、少なくとも「区間タイマ＋ログ」は必須で、外部ツールは段階的に導入するとよい。

13.6 再現性（Reproducibility）：結果を“出し直せる”設計にする

13.6.1 何が再現性を壊すか

依存ライブラリのバージョン差（コンパイラ・BLAS・MPI・Python環境）
並列還元（reduction）の順序違いによる丸め誤差の差（浮動小数は結合律が成り立たない）
乱数（seed未固定、並列乱数のストリーム設計不備）
GPUの非決定性（アルゴリズムや実装により結果が微小に揺れることがある）

再現性は「完全一致」だけでなく、「許容誤差内で一致」も含めて設計する（回帰テストで扱う）。

13.6.2 コンテナと環境固定

研究コードの再現性には、実行環境を“丸ごと固定する”手法が強力である。
HPCでは、管理者権限なしで使えるコンテナが広く使われる。

実験に必要なソフト一式をコンテナ化
実行ログに「コンテナID・git commit・入力ファイルhash」を残す
結果と一緒に環境情報を保存

13.7 ソフトウェア設計

13.7.1 構成管理（Configuration）

“コードを変えずに条件を変える”ために、設定をファイル（例：YAML/TOML/JSON）に分離する。

物理パラメータ
ソルバ（収束判定、前処理種）
並列設定
出力（頻度、フォーマット）

13.7.2 テスト（最低限の3種）

単体テスト：小さな関数（例：行列組立、境界条件適用）
回帰テスト：同じ入力で同等の出力（許容誤差つき）
収束テスト：格子・刻みを細かくすると誤差が理論通り下がるか

13.7.3 ログとメタデータ

最低限ログに残すべきもの：

実行日時、計算機情報、MPIプロセス数、スレッド数
入力ファイル（またはハッシュ）
反復回数、残差履歴、停止理由
主要な性能指標（総時間、主要カーネル時間、通信時間など）

13.8 演習（性能と再現性を“レポート化”する）

演習13-A：強/弱スケーリングを測る（簡易ベンチ）

疎行列-ベクトル積（SpMV）を用意（例：Poisson行列）
プロセス数 $p$ を変えて強スケーリングを測定
1プロセスあたりの格子数を一定にして弱スケーリングを測定
$S (p), E (p)$ を計算し、どこから頭打ちになるかを議論

演習13-B：プロファイル→仮説→改善→再測定

タイマで主要区間（組立、SpMV、通信、前処理）を計測
ルーフライン的に「メモリ律速/通信律速」を仮説化
1つだけ改善（例：通信回数削減、データ構造変更、前処理更新頻度の変更）
再測定して改善が本物か確認（“速くなった気がする”は禁止）

演習13-C：再現性パッケージを作る

入力一式、実行スクリプト、ログ出力、環境情報をまとめる
別マシン（または別環境）で再実行し、同等結果が出るか確認
「一致しないときの切り分け手順」を文章化する（研究実務ではここが価値になる）

まとめ

本書を通して得られる最も重要な見取り図は、「大規模計算の主戦場は疎な線形代数である」という事実である。PDEを離散化すれば疎行列が現れ、非線形問題は線形化により疎線形系を反復的に解く形に落ち、最適化も勾配計算や線形サブ問題を通じて疎線形代数と結びつく。したがって、疎構造の理解（データ構造・グラフ・再順序化）と、Krylov法＋前処理という設計技術が、分野横断で最も再利用性の高い基盤になる。

第二に、数値計算の信頼性は「正しい式」と同じくらい「正しい計算過程」に依存する。浮動小数点の丸め誤差、桁落ち、条件数、安定性は、単なる注意事項ではなく、計算結果の解釈を左右する“理論の一部”である。時間発展ではCFLや剛性が刻み幅を支配し、陰解法では線形（非線形）ソルバが計算の中核になる。どの章でも、誤差と安定性の視点を一貫して持つことで、「計算できた」から「正しく計算できた」へと視点が進む。

第三に、研究コードはアルゴリズムだけでは完成しない。並列化での通信コスト、強/弱スケーリング、プロファイルによるボトルネック同定、再現性（環境固定・乱数・ログ）、テストと設定管理は、研究成果を“検証可能な形”で残すための必須要件である。本書では最終章でそれらを体系化し、学生が将来、共同研究・論文・学会発表の場で説得力ある計算結果を提示できるよう、数値解法と実務を同じ地平で接続する。

展望

今後の大規模数値計算は、単に計算機が速くなるだけでなく、計算機アーキテクチャの多様化（CPUの多コア化、GPU/アクセラレータ、異種混在）とともに発展する。これにより、「浮動小数点演算を減らす」よりも「データ移動を減らす」「行列を作らない（matrix-free）」「通信を隠す（overlap）」といった発想の重要性が増す。疎行列演算や前処理は、メモリ帯域・通信・並列還元の影響を強く受けるため、アルゴリズム設計と性能モデル（例：算術強度や通信モデル）を往復する能力が、ますます価値を持つ。

また、PDEや最適化の世界では、ブラックボックス化されたソルバを“使う”だけでなく、“作る・組み合わせる”ことが競争力になる。Newton-Krylov、IMEX、近接勾配、分割法などの枠組みは、方程式の構造（対称性、保存則、スケール分離、疎性）を利用して性能と安定性を引き出す設計思想である。今後は、物理モデルの複雑化（マルチフィジックス、マルチスケール）に伴い、単一手法ではなく、複数ソルバを整合的に接続する“ソルバ工学”が中心課題になっていく。

さらに、データ駆動・機械学習の拡大は、最適化と数値解法の境界を曖昧にしつつある。大規模最適化は逆問題や材料設計に直結し、PDE制約付き最適化や学習とシミュレーションの融合（ハイブリッドモデリング）も一般化していく。こうした流れの中で、本書が提供する「疎線形代数」「誤差・安定性」「前処理」「時間積分」「並列実務」という基盤は、分野が変わっても長く効く“共通言語”になる。本書の次の発展としては、GPU前提のmatrix-free実装、領域分割型前処理の実装演習、PDE制約付き最適化のミニプロジェクトなどを追加し、研究現場の標準形へさらに近づけることができる。

付録

用語集：本書で頻出する専門用語（定義・直観・どこで使うか）

目的：数式の「記号」と「計算上の役割」を結びつけ、読み飛ばしが起きないようにする。

1. 勾配（gradient） $\nabla f (x)$

定義

スカラー値関数 $f : R^{n} \to R$ に対して

\nabla f (x) = [\begin{matrix} \frac{\partial f}{\partial x_{1}} \\ ⋮ \\ \frac{\partial f}{\partial x_{n}} \end{matrix}] .

直観

$\nabla f (x)$ は「その点で $f$ が最も増える方向」を向くベクトル。
したがって $- \nabla f (x)$ は「最も減る方向」なので、勾配降下法は $x_{k + 1} = x_{k} - α_{k} \nabla f (x_{k})$ と進む。

どこで使うか

最適化（第10章）、Newton法（第9章）、PDEの変分（弱形式）をエネルギー最小化として見るときなど。

2. ヤコビアン（Jacobian） $J (x)$

定義（ベクトル値関数の微分）

ベクトル値関数 $F : R^{n} \to R^{m}$ に対して

J (x) = \frac{\partial F}{\partial x} = [\begin{matrix} \frac{\partial F_{1}}{\partial x_{1}} & \dots & \frac{\partial F_{1}}{\partial x_{n}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial F_{m}}{\partial x_{1}} & \dots & \frac{\partial F_{m}}{\partial x_{n}} \end{matrix}] .

行数 $m$ ：出力の次元
列数 $n$ ：入力の次元

直観

$F (x)$ を $x$ の周りで一次近似すると $F (x + δ) \approx F (x) + J (x) δ .$
「入力の微小変化 $δ$ が、出力をどの方向にどれだけ変えるか」を線形写像として表す。

どこで使うか

非線形方程式 $F (x) = 0$ をNewton法で解くとき（第9章）：

J (x_{k}) δ_{k} = - F (x_{k}), x_{k + 1} = x_{k} + δ_{k} .

陰解法で出る非線形方程式（第12章）でも同様に現れる。

3. ヘッセ行列（Hessian） $\nabla^{2} f (x)$

定義（スカラー値関数の2階微分）

スカラー値関数 $f : R^{n} \to R$ に対して

\nabla^{2} f (x) = [\begin{matrix} \frac{\partial^{2} f}{\partial x_{1} \partial x_{1}} & \dots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{n}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial^{2} f}{\partial x_{n} \partial x_{1}} & \dots & \frac{\partial^{2} f}{\partial x_{n} \partial x_{n}} \end{matrix}] .

直観

ヘッセ行列は「曲率（2次の変化）」を表す。
1次元では $f^{″} (x)$ が曲率。多次元ではそれが行列になる。
2次近似（テイラー展開）： $f (x + δ) \approx f (x) + \nabla f (x)^{⊤} δ + \frac{1}{2} δ^{⊤} \nabla^{2} f (x) δ .$

どこで使うか

Newton法（最適化版）：

x_{k + 1} = x_{k} - {(\nabla^{2} f (x_{k}))}^{- 1} \nabla f (x_{k}) .

大規模では $\nabla^{2} f$ を作る・逆を取るのが重く、準Newton（L-BFGSなど）に繋がる（第10章）。

4. 残差（residual） $r$

定義

線形方程式 $A x = b$ に対して、近似解 $\tilde{x}$ のズレを

r = b - A \tilde{x}

で定義する。

直観

「いまの解で、方程式をどれだけ満たしていないか」の指標。
Krylov法ではこの残差のノルム $∥ r ∥$ が減るように反復を進める。

どこで使うか

反復法全般（第5〜7章）、Newton法の線形サブ問題（第9章）、陰解法での線形ソルバ（第12章）。

5. 条件数（condition number） $κ (A)$

定義（2ノルム条件数）

正則行列 $A$ に対して

κ_{2} (A) = ∥ A ∥_{2} ∥ A^{- 1} ∥_{2} .

対称正定値（SPD）なら固有値で

κ_{2} (A) = \frac{λ_{max} (A)}{λ_{min} (A)} .

直観

$κ (A)$ が大きいほど「入力（右辺）の小さな誤差が解に大きく増幅される」。
反復法では、条件数が大きいほど収束が遅くなりやすい（特にSPDのCGで顕著）。

どこで使うか

誤差と安定性（第2章）、前処理の意義（第7章）、最適化の収束（第10章）。

6. 安定性（stability）と絶対安定領域

定義（テスト方程式）

y^{'} = λ y

に対して、数値法が

y_{n + 1} = R (z) y_{n}, z = h λ

と書けるとき、 $| R (z) | \leq 1$ を満たす $z$ の集合が絶対安定領域。

直観

真の解が減衰する系（ $Re (λ) < 0$ ）でも、刻み幅 $h$ が大きいと数値解が発散し得る。
陽解法は安定領域が小さく、陰解法は広い（剛性に強い）。

どこで使うか

時間発展問題（第12章）。

7. CFL条件（Courant–Friedrichs–Lewy condition）

定義（代表例）

移流方程式で

CFL = \frac{| c | Δ t}{Δ x} ≲ 1

のような条件として現れる（スキーム依存）。

拡散ではしばしば

Δ t \leq C \frac{Δ x^{2}}{D}

（ $C$ は定数）という形になる。

直観

1ステップで情報が格子を“飛び越えない”ようにする制約。
空間分解能を上げる（ $Δ x ↓$ ）ほど $Δ t$ が厳しくなる。

どこで使うか

陽解法の破綻、剛性の理解（第12章）。

8. 剛性（stiffness）

定義（厳密定義は文脈依存）

一般に「解がゆっくり変化して見えるのに、安定性のために極端に小さい $Δ t$ が必要」な状況。

線形系 $y^{'} = A y$ では、固有値の実部のスケール差が大きいと剛性が出やすい。

直観

速い減衰モードが“すぐ消える”のに、陽解法はその速いモードに合わせて刻みを小さくしないと発散する。
だから陰解法が必要になる。

どこで使うか

時間発展問題、拡散優勢、反応拡散など（第12章）。

9. 前処理（preconditioning）と前処理行列 $M$

定義

線形系 $A x = b$ を、解きやすい形に変換する。典型は

左前処理： $M^{- 1} A x = M^{- 1} b$
右前処理： $A (M^{- 1} y) = b, x = M^{- 1} y$

直観

目的は「条件数を下げる（収束を速くする）」こと。
$M \approx A$ だが、 $M$ は解くのが簡単（ILU/IC、Jacobiなど）。

どこで使うか

Krylov法（第6章）を実戦投入する鍵（第7章）。陰解法の各ステップ線形解法にも直結（第12章）。

10. Krylov部分空間（Krylov subspace） $K_{k} (A, r_{0})$

定義

初期残差 $r_{0} = b - A x_{0}$ に対し

K_{k} (A, r_{0}) = span {r_{0}, A r_{0}, A^{2} r_{0}, \dots, A^{k - 1} r_{0}} .

直観

「行列 $A$ を何回か当てたベクトルで張られる空間」の中で最適な近似解を探す。
大規模で行列を因数分解しない代わりに、行列-ベクトル積を繰り返して解を改善する。

どこで使うか

CG/GMRES/BiCGStabなど（第6章）、Newton-Krylov（第9章）、陰解法の線形ソルバ（第12章）。

11. 弱形式（weak form）とテスト関数

定義（Poissonの例）

- Δ u = f

に対し、任意のテスト関数 $v$ を掛けて積分し、部分積分で

\int_{Ω} \nabla u \cdot \nabla v d x = \int_{Ω} f v d x

のような形に変換したもの。

直観

微分方程式を、積分（平均的関係）として書き直してから近似する。
FEMが複雑形状に強い理由の中心。

どこで使うか

FEM（第11章）、変分法、エネルギー最小化との接続。

12. ルーフライン（Roofline）と算術強度 $I$

定義

算術強度：

I = \frac{FLOPs}{Bytes}

性能は概ね

計算律速：演算性能の上限
メモリ律速：メモリ帯域の上限のどちらで頭打ちになる。

直観

$I$ が小さい（例：SpMV）と、演算器を増やしても速くならず、メモリ帯域が支配する。
大規模計算では「データ移動を減らす」が重要になる。

どこで使うか

性能評価（第13章）、疎行列演算の最適化方針の判断。

13. スケーリング：強スケーリング・弱スケーリング

定義

強スケーリング：問題サイズ固定でプロセス数 $p$ を増やす $S (p) = \frac{T (1)}{T (p)}, E (p) = \frac{S (p)}{p}$
弱スケーリング：1プロセスあたりの仕事量一定で $p$ を増やす（理想は $T (p)$ 一定）

直観

強スケーリングは「同じ問題をどこまで速くできるか」
弱スケーリングは「問題をどこまで大きくできるか」

どこで使うか

並列計算の性能評価（第13章）。

14. 再現性（reproducibility）と決定性（determinism）

定義（実務的）

再現性：同じ入力・同等環境で、同等の結果（許容誤差内）を得られる性質
決定性：完全に同じビット列の結果が毎回得られる性質（並列・GPUでは難しい場合がある）

直観

浮動小数点は「足し算の順序」で結果が微小に変わり得る（結合律が成り立たない）。
したがって研究では「許容誤差つき再現性」を設計し、ログ・設定・環境を残す。

どこで使うか

研究コード設計・論文の信頼性（第13章）。

15. 代表的なノルム（norm） $∥ \cdot ∥$

定義

ベクトル $x \in R^{n}$ に対して

2ノルム： $∥ x ∥_{2} = \sqrt{\sum_{i = 1}^{n} x_{i}^{2}}$
1ノルム： $∥ x ∥_{1} = \sum_{i = 1}^{n} | x_{i} |$
$\infty$ ノルム： $∥ x ∥_{\infty} = max_{i} | x_{i} |$

直観

残差の大きさ $∥ r ∥$ を測って「収束したか」を判断するのに使う。
どのノルムを使うかで停止条件の感度が変わる。

どこで使うか

反復法の停止条件、誤差評価、最適化の収束判定。

計算科学のための大規模数値解法 ​

第1章 大規模数値計算とは何か：問題分類と計算資源 ​

1.0 この章で学ぶこと ​

1.1 「大規模」とは何か：3つのボトルネック ​

(A) 計算量（Compute） ​

(B) メモリ（Memory） ​

(C) 通信・データ移動（Communication / Data movement） ​

1.2 数値計算で扱う「問題」の代表的な型 ​

1.2.1 線形問題（連立一次方程式） ​

1.2.2 非線形問題（方程式 f(x)=0） ​

1.2.3 最適化問題（目的関数の最小化） ​

1.2.4 PDE（偏微分方程式）の数値解 ​

1.3 オーダー記法（O( )）で何を見積もるのか ​

1.3.1 時間計算量（Time complexity） ​

1.3.2 メモリ計算量（Memory complexity） ​

1.4 密行列と疎行列：大規模計算の分岐点 ​

1.4.1 疎行列（sparse matrix）とは ​

1.4.2 保存量の比較（密 vs 疎） ​

密行列（倍精度） ​

CSR形式（概算） ​

1.4.3 具体例：もし n=106 なら？ ​

1.5 メモリ階層：なぜ“データを動かす”のが遅いのか ​

1.6 “計算が遅い”の正体：Rooflineの考え方（直観） ​

1.6.1 算術強度（Arithmetic Intensity） ​

1.6.2 性能の上限（概念式） ​

1.7 並列計算と通信：強スケーリング・弱スケーリング ​

1.7.1 通信時間の単純モデル ​

1.7.2 強スケーリング（Strong scaling）とAmdahlの法則 ​

1.7.3 弱スケーリング（Weak scaling）とGustafsonの法則 ​

1.8 本章のまとめ：大規模計算の“設計図” ​

演習 ​

演習1：密行列 vs 疎行列のメモリ見積もり ​

設問 ​

演習2：計算時間のスケーリング測定（体感） ​

レポートの観点 ​

第2章 浮動小数点と誤差：丸め誤差・桁落ち・悪条件 ​

2.0 この章で学ぶこと ​

2.1 実数と浮動小数点数は別物：まず“数の世界観”を揃える ​

2.1.1 浮動小数点数の一般形 ​

2.2 IEEE 754 の重要ポイント：NaN・∞・±0・サブノーマル ​

2.2.1 特殊値の意味 ​

2.3 丸め（rounding）とは何か：計算は毎回“近似”される ​

2.3.1 丸めモード（代表例） ​

2.4 機械イプシロンと unit roundoff：誤差を式で扱う準備 ​

2.4.1 機械イプシロン（machine epsilon） ​

2.4.2 unit roundoff u ​

2.5 浮動小数点の基本モデル：fl( ) と δ ​

2.6 桁落ち（catastrophic cancellation）：危険な差分 ​

2.6.1 何が起きるのか ​

2.6.2 典型例：1+x−1（x が小さい） ​

回避策（有理化） ​

2.7 誤差の種類：絶対誤差・相対誤差・ULP ​

2.8 誤差伝播：足し算でも誤差は積み上がる ​

2.8.1 逐次和の誤差（なぜ和が危ないのか） ​

2.8.2 最近の計算で重要な点：並列和と再現性 ​

2.8.3 対策例：Kahan和（補償和） ​

2.9 悪条件（ill-conditioned）と条件数：問題そのものが不安定 ​

2.9.1 関数の条件数（1変数の例） ​

2.9.2 行列の条件数（線形方程式） ​

2.10 前進誤差と後退誤差：安定性を測る2つの物差し ​

2.10.1 前進誤差（forward error） ​

2.10.2 後退誤差（backward error） ​

残差による直観（線形系） ​

2.11 数値的安定性：後退安定なら（だいたい）安心 ​

2.12 最近の大規模計算での重要トピック：混合精度（mixed precision） ​

2.13 演習 ​

演習2-1：桁落ちの体験（式変形の威力） ​

演習2-2：条件数の大きい行列で解が不安定になることの確認 ​

演習2-3：和の順序で結果が変わることの確認 ​

第3章 行列とグラフ：疎構造・帯行列・再順序化 ​

3.0 この章で学ぶこと（なぜグラフが出てくるのか） ​

3.1 疎行列の基本：nnz とスパース演算 ​

3.2 疎行列のデータ構造：COO / CSR / CSC ​

3.2.1 CSR（行圧縮）：3つの配列で表す ​

3.2.2 CSC（列圧縮）：列ごとにまとめる ​

3.2.3 最近の動向：sparse “matrix” から sparse “array” へ ​

3.3 疎行列をグラフとして見る：隣接関係＝非ゼロ ​

3.3.1 対称パターンの場合（無向グラフ） ​

3.3.2 非対称の場合（二部グラフなど） ​

3.4 帯行列・帯幅（bandwidth）と“近くに寄せる”発想 ​

計算科学のための大規模数値解法

第1章　大規模数値計算とは何か：問題分類と計算資源

1.0 この章で学ぶこと

1.1 「大規模」とは何か：3つのボトルネック

(A) 計算量（Compute）

(B) メモリ（Memory）

(C) 通信・データ移動（Communication / Data movement）

1.2 数値計算で扱う「問題」の代表的な型

1.2.1 線形問題（連立一次方程式）

1.2.2 非線形問題（方程式 $f (x) = 0$ ）

1.2.3 最適化問題（目的関数の最小化）

1.2.4 PDE（偏微分方程式）の数値解

1.3 オーダー記法（O( )）で何を見積もるのか

1.3.1 時間計算量（Time complexity）

1.3.2 メモリ計算量（Memory complexity）

1.4 密行列と疎行列：大規模計算の分岐点

1.4.1 疎行列（sparse matrix）とは

1.4.2 保存量の比較（密 vs 疎）

密行列（倍精度）

CSR形式（概算）

1.4.3 具体例：もし $n = 10^{6}$ なら？

1.5 メモリ階層：なぜ“データを動かす”のが遅いのか

1.6 “計算が遅い”の正体：Rooflineの考え方（直観）

1.6.1 算術強度（Arithmetic Intensity）

1.6.2 性能の上限（概念式）

1.7 並列計算と通信：強スケーリング・弱スケーリング

1.7.1 通信時間の単純モデル

1.7.2 強スケーリング（Strong scaling）とAmdahlの法則

1.7.3 弱スケーリング（Weak scaling）とGustafsonの法則

1.8 本章のまとめ：大規模計算の“設計図”

演習

演習1：密行列 vs 疎行列のメモリ見積もり

設問

演習2：計算時間のスケーリング測定（体感）

レポートの観点

第2章　浮動小数点と誤差：丸め誤差・桁落ち・悪条件

2.0 この章で学ぶこと

2.1 実数と浮動小数点数は別物：まず“数の世界観”を揃える

2.1.1 浮動小数点数の一般形

2.2 IEEE 754 の重要ポイント：NaN・∞・±0・サブノーマル

2.2.1 特殊値の意味

2.3 丸め（rounding）とは何か：計算は毎回“近似”される

2.3.1 丸めモード（代表例）

2.4 機械イプシロンと unit roundoff：誤差を式で扱う準備

2.4.1 機械イプシロン（machine epsilon）

2.4.2 unit roundoff $u$

2.5 浮動小数点の基本モデル：fl( ) と $δ$

2.6 桁落ち（catastrophic cancellation）：危険な差分

2.6.1 何が起きるのか

2.6.2 典型例： $\sqrt{1 + x} - 1$ （ $x$ が小さい）

回避策（有理化）

2.7 誤差の種類：絶対誤差・相対誤差・ULP

2.8 誤差伝播：足し算でも誤差は積み上がる

2.8.1 逐次和の誤差（なぜ和が危ないのか）

2.8.2 最近の計算で重要な点：並列和と再現性

2.8.3 対策例：Kahan和（補償和）

2.9 悪条件（ill-conditioned）と条件数：問題そのものが不安定

2.9.1 関数の条件数（1変数の例）

2.9.2 行列の条件数（線形方程式）

2.10 前進誤差と後退誤差：安定性を測る2つの物差し

2.10.1 前進誤差（forward error）

2.10.2 後退誤差（backward error）

残差による直観（線形系）

2.11 数値的安定性：後退安定なら（だいたい）安心

2.12 最近の大規模計算での重要トピック：混合精度（mixed precision）

2.13 演習

演習2-1：桁落ちの体験（式変形の威力）

演習2-2：条件数の大きい行列で解が不安定になることの確認

演習2-3：和の順序で結果が変わることの確認

第3章　行列とグラフ：疎構造・帯行列・再順序化

3.0 この章で学ぶこと（なぜグラフが出てくるのか）

3.1 疎行列の基本： $nnz$ とスパース演算

3.2 疎行列のデータ構造：COO / CSR / CSC

3.2.1 CSR（行圧縮）：3つの配列で表す

3.2.2 CSC（列圧縮）：列ごとにまとめる

3.2.3 最近の動向：sparse “matrix” から sparse “array” へ

3.3 疎行列をグラフとして見る：隣接関係＝非ゼロ

3.3.1 対称パターンの場合（無向グラフ）

3.3.2 非対称の場合（二部グラフなど）

3.4 帯行列・帯幅（bandwidth）と“近くに寄せる”発想