$100 \times 100$ の連立1次方程式を手計算で解くことは現実的ではありません。しかしコンピュータなら一瞬です。では $10^6 \times 10^6$ の場合はどうでしょうか。実は、行列のサイズが大きくなると解法の選択が計算時間を桁違いに左右します。$O(n^3)$ のアルゴリズムで $n = 10^6$ を処理しようとすると、$10^{18}$ 回の演算が必要になり、現在のスーパーコンピュータでも数年以上かかります。しかし行列の構造を活かした適切なアルゴリズムを使えば、同じ問題が数秒で解ける場合もあります。
数値線形代数（numerical linear algebra）は、連立1次方程式 $A\bm{x} = \bm{b}$ を数値的に効率よく解く手法を研究する分野です。ここで「効率よく」とは、単に速いだけでなく、数値的に安定（丸め誤差が蓄積しにくい）であることも含みます。浮動小数点演算では避けられない丸め誤差が、不適切なアルゴリズムでは致命的に拡大する場合があるためです。
有限要素法の離散化、画像処理、機械学習のパラメータ推定など、科学技術計算のほぼすべての場面で連立方程式の求解が必要になります。「数値計算の時間の80%以上は線形代数に費やされている」とも言われるほど、数値線形代数は計算科学の根幹をなす分野です。
数値線形代数を理解すると、以下のような場面で活用できます。
- 構造解析: 有限要素法で生じる大規模疎行列の連立方程式。航空機の翼の応力解析では数百万自由度の系を解く必要があります
- 機械学習: 正規方程式 $X^TX\bm{w} = X^T\bm{y}$ の数値的に安定な解法。直接 $X^TX$ を計算するのは悪条件になりやすく、QR分解やSVDを使うのが推奨されます
- 流体力学: CFDの圧力ポアソン方程式の反復求解。3次元流体シミュレーションでは $10^7$ 以上の未知数を持つ疎行列が生じます
- グラフ解析: ページランクの固有値問題。Googleの初期のウェブページランキングは、数十億次元の疎行列の主固有ベクトルを反復法で計算していました
本記事の内容
- ガウス消去法とピボット選択
- LU分解の理論と計算量
- 条件数と数値的安定性
- ヤコビ法・ガウス=ザイデル法（反復法）
- 共役勾配法（CG法）の理論
- Pythonでの実装と比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
ガウス消去法
直感的な理解
ガウス消去法の考え方は、中学校で学んだ連立方程式の「加減法」の一般化です。2元連立方程式では、一方の式を定数倍してもう一方から引くことで変数を1つ消去しました。ガウス消去法はこの操作を $n$ 変数の連立方程式に体系的に適用し、変数を1つずつ消去していくことで、最終的に1変数の方程式に帰着させます。
たとえば3元連立方程式なら、まず第1式を使って第2式・第3式から $x_1$ を消去し、次に（$x_1$ が消えた）第2式を使って第3式から $x_2$ を消去します。すると第3式は $x_3$ だけの方程式になるので、$x_3$ の値が直ちに求まります。あとは逆順に $x_2$, $x_1$ を求めればよいのです。
アルゴリズム
ガウス消去法は連立方程式 $A\bm{x} = \bm{b}$ を解く最も基本的な直接法です。前進消去（forward elimination）で行列を上三角行列に変換し、後退代入（back substitution）で解を求めます。
前進消去: 第 $k$ ステップで、$k$ 行目より下の行から $k$ 列の成分を消去します。
$$ a_{ij}^{(k+1)} = a_{ij}^{(k)} – \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}} a_{kj}^{(k)}, \quad i > k $$
乗数 $l_{ik} = a_{ik}^{(k)} / a_{kk}^{(k)}$ はLU分解の $L$ 行列の要素になります。
後退代入: 上三角系 $U\bm{x} = \bm{c}$ を下から順に解きます。
$$ x_i = \frac{1}{u_{ii}}\left(c_i – \sum_{j=i+1}^n u_{ij}x_j\right), \quad i = n, n-1, \ldots, 1 $$
計算量
前進消去の演算回数（乗算・除算）は
$$ \sum_{k=1}^{n-1} \sum_{i=k+1}^n (n – k + 1) \approx \frac{n^3}{3} $$
後退代入は $n^2/2$ 程度なので、全体の計算量は $O(n^3/3)$ です。この $n^3/3$ という値は、$n$ が大きくなると急速に増大します。たとえば $n = 1000$ では約 $3.3 \times 10^8$ 回、$n = 10000$ では約 $3.3 \times 10^{11}$ 回の演算が必要です。現代のCPUは毎秒 $10^{10}$ 程度の浮動小数点演算が可能なので、$n = 10000$ でも数十秒で解けますが、$n = 10^6$ になると現実的な時間では解けなくなります。
ピボット選択
$a_{kk}^{(k)} = 0$（または非常に小さい値）の場合、消去は失敗するか数値的に不安定になります。部分ピボット選択（partial pivoting）は、$k$ 列の $k$ 行以下で最大の絶対値を持つ行を探し、$k$ 行と交換します。
$$ |a_{pk}^{(k)}| = \max_{i \geq k} |a_{ik}^{(k)}| $$
部分ピボット選択により、乗数 $|l_{ik}| \leq 1$ が保証され、数値的安定性が大幅に向上します。なぜ乗数の大きさが安定性に関係するかを直感的に理解しておきましょう。乗数 $l_{ik}$ が非常に大きいと、$l_{ik}$ を掛けた行を引く操作で、元の値に比べてはるかに大きな値が現れます。大きな値から小さな値を引く際に桁落ちが発生し、有効桁数が失われるのです。ピボット選択により乗数を1以下に抑えることで、この桁落ちのリスクを最小化できます。
なお、完全ピボット選択（行と列の両方を探索して最大要素を選ぶ）はさらに安定ですが、探索コストが $O(n^2)$ になるため、部分ピボット選択が実用上ほとんどの場合に十分です。
ガウス消去法を体系的に整理すると、LU分解が自然に導かれます。
LU分解
直感的な理解
ガウス消去法は「計算してその場で使い捨てる」手法でしたが、消去の過程を記録しておけば再利用できます。この「消去過程の記録」がLU分解です。
料理に例えると、ガウス消去法は「レシピを見ながら1回だけ料理を作る」ようなものですが、LU分解は「レシピ自体を整理して保存しておく」ようなものです。レシピ（分解結果）が手元にあれば、材料（右辺ベクトル $\bm{b}$）を変えて何度でも素早く料理（求解）ができます。
理論
LU分解は、行列 $A$ を下三角行列 $L$（対角成分が1）と上三角行列 $U$ の積に分解します。
$$ A = LU $$
ガウス消去法の前進消去で得られる上三角行列が $U$、乗数を集めた下三角行列が $L$ です。
$$ L = \begin{pmatrix} 1 & & & \\ l_{21} & 1 & & \\ l_{31} & l_{32} & 1 & \\ \vdots & & & \ddots \end{pmatrix}, \quad U = \begin{pmatrix} u_{11} & u_{12} & u_{13} & \cdots \\ & u_{22} & u_{23} & \cdots \\ & & u_{33} & \cdots \\ & & & \ddots \end{pmatrix} $$
LU分解の利点
$A\bm{x} = \bm{b}$ を解くには
- $L\bm{y} = \bm{b}$（前進代入, $O(n^2)$）
- $U\bm{x} = \bm{y}$（後退代入, $O(n^2)$）
LU分解自体は $O(n^3/3)$ ですが、一度分解すれば異なる $\bm{b}$ に対して $O(n^2)$ で解けます。これは $\bm{b}$ が複数ある場合に非常に効率的です。たとえば時間発展問題で各時刻の右辺ベクトルが異なる場合や、行列の逆行列を列ごとに計算する場合（$A^{-1}$ の第 $j$ 列は $A\bm{x} = \bm{e}_j$ の解）などで、この利点が活きます。
ピボット付きLU分解
部分ピボット選択を含むLU分解は $PA = LU$（$P$ は置換行列）と表されます。これが実用上の標準的なLU分解です。
コレスキー分解
$A$ が対称正定値行列の場合、$A = LL^T$（$L$ は下三角行列）と分解できます。コレスキー分解はLU分解の特殊ケースで、計算量は $n^3/6$（LU分解の約半分）です。計算量が半分になる理由は、対称性 $A = A^T$ から $L$ と $U = L^T$ が独立でなくなり、保存すべき情報量が半分で済むためです。さらにピボット選択も不要で（対称正定値行列のすべての主座小行列式が正であるため）、対称正定値行列に対しては常に最良の選択です。
コレスキー分解が適用できる条件は「対称正定値」ですが、この条件は実用上非常によく満たされます。有限要素法の剛性行列、正規方程式 $X^TX$、共分散行列など、多くの重要な行列が対称正定値です。コレスキー分解の存在自体が、行列が正定値であることの数値的なテストとしても使えます（分解中に負の平方根が現れたら正定値でない）。
直接法は正確な解を有限ステップで得られますが、大規模疎行列に対しては反復法の方が効率的な場合があります。しかし直接法と反復法のどちらが適切かを判断するためには、まず行列の「数値的な扱いやすさ」を定量的に評価する指標が必要です。それが条件数です。
条件数と数値的安定性
直感的な理解
条件数は、連立方程式が「小さな揺れに対してどれだけ敏感か」を数値化した指標です。日常的な例でいえば、2本のほぼ平行な直線の交点を求める問題を考えてください。直線の傾きをほんの少し変えるだけで、交点は大きく移動してしまいます。このような問題は「悪条件」であり、条件数が大きいことに対応します。一方、2本の直線がほぼ直交していれば、傾きの小さな変化に対して交点はほとんど動きません。これが「良条件」です。
条件数の定義
$A\bm{x} = \bm{b}$ において、$\bm{b}$ に小さな摂動 $\delta\bm{b}$ が加わったとき、解の相対誤差は
$$ \frac{\|\delta\bm{x}\|}{\|\bm{x}\|} \leq \kappa(A) \frac{\|\delta\bm{b}\|}{\|\bm{b}\|} $$
ここで $\kappa(A) = \|A\| \cdot \|A^{-1}\|$ が条件数（condition number）です。
2ノルムの場合、$\kappa_2(A) = \sigma_{\max} / \sigma_{\min}$（最大特異値と最小特異値の比）です。
条件数の意味
- $\kappa(A) \approx 1$: 良条件（well-conditioned）。入力の誤差が拡大されない
- $\kappa(A) \gg 1$: 悪条件（ill-conditioned）。入力の小さな誤差が大きく拡大される
- $\kappa(A) \approx 10^k$: おおよそ $k$ 桁の精度が失われる
倍精度浮動小数点（約16桁の精度）で条件数 $10^{10}$ の行列を解くと、有効数字は約6桁しか残りません。これは「入力に16桁の精度があるのに、出力は6桁しか信頼できない」ということです。
ヒルベルト行列の例
$H_{ij} = 1/(i + j – 1)$ で定義されるヒルベルト行列は悪条件行列の典型例です。$n = 10$ で $\kappa \approx 10^{13}$、$n = 20$ で $\kappa \approx 10^{28}$ になり、実質的に求解不可能です。
ヒルベルト行列がこれほど悪条件になる理由は直感的に理解できます。$H$ の各行は $1/(i+j-1)$ という「少しずつずれた」なだらかな関数であり、行間の差が非常に小さいのです。ほぼ同じ方向を向いたベクトルの集まりから独立な情報を引き出すのは困難であり、これが大きな条件数として現れます。
悪条件への対処法
条件数が大きい場合、以下のアプローチが有効です。
- スケーリング: 行列の行や列を適切にスケーリングして、条件数を改善する。物理単位の不統一が悪条件の原因になっていることがあります
- 正則化: チホノフ正則化 $(A^TA + \lambda I)\bm{x} = A^T\bm{b}$ のように、正則化項を追加して問題を安定化する
- 高精度演算: 多倍長精度演算を用いて、丸め誤差自体を減らす。ただし計算コストが大幅に増加します
- 問題の再定式化: SVDを用いた擬似逆行列や、QR分解による最小二乗解法に切り替える
条件数が大きい場合、直接法でも反復法でも精度に限界があることを認識しておく必要があります。次に、直接法とは異なるアプローチである反復法を見ていきましょう。反復法は条件数が適度な問題で、行列の疎構造を活かして効率的に解を求める手法です。
反復法
動機
直接法の計算量 $O(n^3)$ は、$n$ が非常に大きい（$10^5$ 以上）場合に問題になります。たとえば $n = 10^6$ の密行列にLU分解を適用すると $10^{18}$ 回の演算が必要ですが、同じ行列が疎行列で各行に平均5個の非零要素しかなければ、反復法では行列-ベクトル積が $O(5n) = O(n)$ で計算でき、数百回の反復で十分な精度が得られる場合があります。
有限要素法や差分法で生じる行列は多くの場合疎行列（sparse matrix, ほとんどの成分が0）です。たとえば3次元ラプラシアンの差分近似では、$n \times n \times n$ の格子で $N = n^3$ 個の未知数がありますが、各方程式に関与する変数は高々7個（自分自身と6つの隣接格子点）です。このような疎構造を活かせるのが反復法の強みです。
反復法の基本的な考え方は、適当な初期推測 $\bm{x}^{(0)}$ から出発して、ある漸化式に従って $\bm{x}^{(1)}, \bm{x}^{(2)}, \ldots$ を生成し、真の解 $\bm{x}^*$ に収束させるというものです。各反復で必要な主要な演算は行列-ベクトル積 $A\bm{v}$ であり、$A$ が疎であればこの演算は $O(\text{nnz})$（非零要素数に比例）で済みます。
ヤコビ法
$A = D + L + U$（$D$: 対角成分、$L$: 狭義下三角、$U$: 狭義上三角）として
$$ \bm{x}^{(k+1)} = D^{-1}(\bm{b} – (L + U)\bm{x}^{(k)}) $$
成分で書くと
$$ x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i – \sum_{j \neq i} a_{ij}x_j^{(k)}\right) $$
各成分の更新に古い値のみを使うため、並列計算に向いています。直感的には、「各変数について、他の変数を固定して自分だけを解く」操作を繰り返しているのがヤコビ法です。
ガウス=ザイデル法
ヤコビ法では全変数の更新に「前の反復の値」を使いましたが、すでに更新した変数の最新の値を使ったほうが、より正確な情報に基づく更新になるはずです。この自然な改良がガウス=ザイデル法です。
$$ x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i – \sum_{j < i} a_{ij}x_j^{(k+1)} - \sum_{j > i} a_{ij}x_j^{(k)}\right) $$
一般にヤコビ法より収束が速いですが、逐次的な更新のため並列化が困難です。
収束条件
反復法 $\bm{x}^{(k+1)} = B\bm{x}^{(k)} + \bm{c}$ が収束するための必要十分条件は
$$ \rho(B) < 1 $$
ここで $\rho(B) = \max |\lambda_i|$ はスペクトル半径です。スペクトル半径が小さいほど収束が速くなります。
対角優位条件 $|a_{ii}| > \sum_{j \neq i} |a_{ij}|$ が満たされれば、ヤコビ法もガウス=ザイデル法も収束します。
SOR法（逐次過緩和法）
ガウス=ザイデル法に緩和パラメータ $\omega$ を導入します。
$$ x_i^{(k+1)} = (1 – \omega)x_i^{(k)} + \frac{\omega}{a_{ii}}\left(b_i – \sum_{j < i} a_{ij}x_j^{(k+1)} - \sum_{j > i} a_{ij}x_j^{(k)}\right) $$
$\omega = 1$ でガウス=ザイデル法に一致し、$1 < \omega < 2$ で過緩和（加速）、$0 < \omega < 1$ で不足緩和になります。最適な $\omega$ の選択で収束を大幅に加速できますが、最適値の決定自体が難しい問題です。
特殊な場合として、1次元ポアソン方程式の差分近似行列では最適緩和パラメータが $\omega_{\text{opt}} = 2/(1 + \sin(\pi h))$ と厳密に求まることが知られています。この場合、SOR法の収束率はガウス=ザイデル法より大幅に改善されます。しかし一般の行列に対して最適な $\omega$ を求めるのは困難であり、そもそもスペクトル半径の情報が必要になるという循環的な問題に直面します。
ヤコビ法、ガウス=ザイデル法、SOR法はいずれも定常反復法（各ステップで同じ操作を繰り返す）に分類されます。これらの手法は収束が遅い場合がありますが、探索方向をより賢く選ぶクリロフ部分空間法はこれを劇的に改善します。
共役勾配法（CG法）
動機と基本思想
定常反復法が「固定的な更新ルールを繰り返す」のに対し、共役勾配法は「これまでの計算結果をフルに活用して最適な方向と距離を選ぶ」という、より洗練された戦略を取ります。
$A$ が対称正定値のとき、$A\bm{x} = \bm{b}$ の解は2次関数
$$ \phi(\bm{x}) = \frac{1}{2}\bm{x}^T A\bm{x} – \bm{b}^T\bm{x} $$
の最小化問題と等価です。この等価性は、$\phi$ の勾配を計算すると $\nabla\phi = A\bm{x} – \bm{b}$ となり、$\nabla\phi = \bm{0}$ の解がまさに $A\bm{x} = \bm{b}$ の解であることから確認できます。$A$ が正定値であるという条件は、$\phi$ が下に凸（唯一の極小値を持つ）であることを保証しています。
幾何学的には、$\phi(\bm{x})$ は $n$ 次元空間での楕円体状の等高面を持つ「お椀」のような形をしており、その最も低い点が解 $\bm{x}^*$ です。
最急降下法（steepest descent）は $-\nabla\phi$ の方向（残差 $\bm{r} = \bm{b} – A\bm{x}$ の方向）に進みますが、楕円体が扁平な場合にジグザグ軌道になるため収束が遅くなりがちです。条件数 $\kappa$ が大きい（楕円体が扁平な）場合、最急降下法は同じ方向を何度も行ったり来たりする非効率な動きをします。
共役勾配法（CG法）は、$A$ に関して共役な方向を使うことでこの問題を解決します。
共役方向とは
2つのベクトル $\bm{d}_i, \bm{d}_j$ が $A$-共役であるとは
$$ \bm{d}_i^T A \bm{d}_j = 0 \quad (i \neq j) $$
を満たすことです。通常の直交条件 $\bm{d}_i^T \bm{d}_j = 0$ の内積を $A$-内積 $\langle \bm{d}_i, \bm{d}_j \rangle_A = \bm{d}_i^T A \bm{d}_j$ に置き換えたものです。
$n$ 個の $A$-共役ベクトルが存在すれば、各方向の1次元最小化を $n$ 回行うだけで厳密解が得られます（丸め誤差を無視すれば）。
CG法のアルゴリズム
初期値 $\bm{x}_0$, 残差 $\bm{r}_0 = \bm{b} – A\bm{x}_0$, 探索方向 $\bm{d}_0 = \bm{r}_0$ として
$$ \alpha_k = \frac{\bm{r}_k^T \bm{r}_k}{\bm{d}_k^T A \bm{d}_k} $$
$$ \bm{x}_{k+1} = \bm{x}_k + \alpha_k \bm{d}_k $$
$$ \bm{r}_{k+1} = \bm{r}_k – \alpha_k A\bm{d}_k $$
$$ \beta_k = \frac{\bm{r}_{k+1}^T \bm{r}_{k+1}}{\bm{r}_k^T \bm{r}_k} $$
$$ \bm{d}_{k+1} = \bm{r}_{k+1} + \beta_k \bm{d}_k $$
各ステップで必要なのは行列-ベクトル積 $A\bm{d}_k$ 1回と内積数回のみです。
収束性
CG法は理論上 $n$ ステップ以内に厳密解に収束します。実用的には、条件数 $\kappa$ が大きくない場合は $n$ よりはるかに少ないステップで十分な精度が得られます。
$$ \frac{\|\bm{x}_k – \bm{x}^*\|_A}{\|\bm{x}_0 – \bm{x}^*\|_A} \leq 2\left(\frac{\sqrt{\kappa} – 1}{\sqrt{\kappa} + 1}\right)^k $$
条件数が大きい場合、前処理（preconditioning）$M^{-1}A\bm{x} = M^{-1}\bm{b}$ により実効的な条件数を下げて収束を加速します。前処理行列 $M$ は $A$ の近似逆行列のようなもので、$M^{-1}A$ の条件数が $A$ の条件数より小さくなることを目指します。理想的には $M = A$ とすれば条件数は1になりますが、そのためには $A^{-1}$ を計算する必要があり、元の問題と同じ困難さに戻ってしまいます。実用的な前処理としては、不完全LU分解（ILU）、不完全コレスキー分解（ICC）、代数的マルチグリッド法などがよく使われます。
ここまでの理論を実際にPythonで確認してみましょう。直接法と反復法の収束特性、条件数の影響、計算時間のスケーリングを可視化します。
Pythonによる実装
以下のコードでは、1次元離散ラプラシアン行列（有限差分法でポアソン方程式を離散化した際に現れる典型的な行列）をテスト問題として使い、各手法の特性を比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import lu, cho_factor, cho_solve, hilbert
import time
# 実装
def jacobi(A, b, x0=None, max_iter=1000, tol=1e-10):
n = len(b)
x = np.zeros(n) if x0 is None else x0.copy()
D_inv = 1.0 / np.diag(A)
residuals = []
for k in range(max_iter):
r = b - A @ x
residuals.append(np.linalg.norm(r))
if residuals[-1] < tol:
break
x_new = D_inv * (b - (A - np.diag(np.diag(A))) @ x)
x = x_new
return x, residuals
def gauss_seidel(A, b, x0=None, max_iter=1000, tol=1e-10):
n = len(b)
x = np.zeros(n) if x0 is None else x0.copy()
residuals = []
for k in range(max_iter):
r = b - A @ x
residuals.append(np.linalg.norm(r))
if residuals[-1] < tol:
break
for i in range(n):
x[i] = (b[i] - A[i, :i] @ x[:i] - A[i, i+1:] @ x[i+1:]) / A[i, i]
return x, residuals
def conjugate_gradient(A, b, x0=None, max_iter=1000, tol=1e-10):
n = len(b)
x = np.zeros(n) if x0 is None else x0.copy()
r = b - A @ x
d = r.copy()
residuals = [np.linalg.norm(r)]
for k in range(max_iter):
Ad = A @ d
alpha = (r @ r) / (d @ Ad)
x = x + alpha * d
r_new = r - alpha * Ad
residuals.append(np.linalg.norm(r_new))
if residuals[-1] < tol:
break
beta = (r_new @ r_new) / (r @ r)
d = r_new + beta * d
r = r_new
return x, residuals
# テスト行列: 対称正定値（離散ラプラシアン）
def make_laplacian_1d(n):
A = np.zeros((n, n))
for i in range(n):
A[i, i] = 2.0
if i > 0:
A[i, i-1] = -1.0
if i < n - 1:
A[i, i+1] = -1.0
return A
fig, axes = plt.subplots(2, 2, figsize=(14, 11))
# --- (a) 反復法の収束比較 ---
ax = axes[0, 0]
n = 50
A = make_laplacian_1d(n)
b = np.ones(n)
_, res_jac = jacobi(A, b)
_, res_gs = gauss_seidel(A, b)
_, res_cg = conjugate_gradient(A, b)
ax.semilogy(res_jac, "b-", linewidth=2, label="Jacobi")
ax.semilogy(res_gs, "r-", linewidth=2, label="Gauss-Seidel")
ax.semilogy(res_cg, "g-", linewidth=2, label="CG")
ax.set_xlabel("Iteration", fontsize=12)
ax.set_ylabel("Residual $\\|\\mathbf{r}\\|$", fontsize=12)
ax.set_title(f"(a) Convergence Comparison ($n = {n}$, 1D Laplacian)", fontsize=12)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, which="both")
# --- (b) 条件数と収束速度 ---
ax = axes[0, 1]
sizes = [10, 20, 50, 100]
for n_size in sizes:
A_test = make_laplacian_1d(n_size)
b_test = np.ones(n_size)
_, res = conjugate_gradient(A_test, b_test, max_iter=200)
kappa = np.linalg.cond(A_test)
ax.semilogy(res, linewidth=1.5, label=f"$n={n_size}$, $\\kappa={kappa:.0f}$")
ax.set_xlabel("CG Iteration", fontsize=12)
ax.set_ylabel("Residual $\\|\\mathbf{r}\\|$", fontsize=12)
ax.set_title("(b) CG Convergence vs Condition Number", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")
# --- (c) 条件数の成長（ヒルベルト行列） ---
ax = axes[1, 0]
ns_hilbert = range(2, 16)
kappas = [np.linalg.cond(hilbert(n_h)) for n_h in ns_hilbert]
ax.semilogy(list(ns_hilbert), kappas, "ro-", linewidth=2, markersize=8)
ax.axhline(1e16, color="gray", linestyle="--", linewidth=1,
label="Double precision limit ($\\approx 10^{16}$)")
ax.set_xlabel("Matrix size $n$", fontsize=12)
ax.set_ylabel("Condition number $\\kappa_2$", fontsize=12)
ax.set_title("(c) Hilbert Matrix Condition Number", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which="both")
# --- (d) 計算時間の比較 ---
ax = axes[1, 1]
sizes_time = [50, 100, 200, 500, 1000]
times_lu = []
times_chol = []
times_cg = []
for n_t in sizes_time:
A_t = make_laplacian_1d(n_t)
b_t = np.ones(n_t)
# LU分解
t0 = time.perf_counter()
for _ in range(3):
np.linalg.solve(A_t, b_t)
times_lu.append((time.perf_counter() - t0) / 3)
# コレスキー分解
t0 = time.perf_counter()
for _ in range(3):
c, low = cho_factor(A_t)
cho_solve((c, low), b_t)
times_chol.append((time.perf_counter() - t0) / 3)
# CG法
t0 = time.perf_counter()
for _ in range(3):
conjugate_gradient(A_t, b_t, tol=1e-10)
times_cg.append((time.perf_counter() - t0) / 3)
ax.loglog(sizes_time, times_lu, "bo-", linewidth=2, markersize=6, label="LU (numpy)")
ax.loglog(sizes_time, times_chol, "rs-", linewidth=2, markersize=6, label="Cholesky")
ax.loglog(sizes_time, times_cg, "g^-", linewidth=2, markersize=6, label="CG")
# 参考線
ns_ref = np.array(sizes_time, dtype=float)
ax.loglog(ns_ref, 1e-7 * ns_ref**3, "k:", linewidth=1, alpha=0.5, label="$O(n^3)$")
ax.set_xlabel("Matrix size $n$", fontsize=12)
ax.set_ylabel("Time (sec)", fontsize=12)
ax.set_title("(d) Computation Time", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which="both")
plt.tight_layout()
plt.savefig("numerical_linear_algebra.png", dpi=150, bbox_inches="tight")
plt.show()
この4つの図から、数値線形代数の手法の特性が明確に読み取れます。
-
左上（反復法の収束比較）: CG法が圧倒的に速く収束しています。ヤコビ法は最も遅く、ガウス=ザイデル法はその約2倍速い収束を示していますが、どちらもCG法には遠く及びません。50次元の1次元ラプラシアンに対して、CG法は約50ステップ以内に収束しています
-
右上（条件数とCG収束）: 行列サイズ $n$ が大きくなると条件数 $\kappa$ も増加し、CG法の収束に必要なステップ数が増えています。条件数と収束速度の関係 $(\sqrt{\kappa} – 1)/(\sqrt{\kappa} + 1)$ が反映されています
-
左下（ヒルベルト行列の条件数）: ヒルベルト行列の条件数は $n$ に対して指数的に増大しています。$n = 12$ 程度で倍精度の限界 $10^{16}$ を超え、$n = 15$ では $10^{17}$ 以上に達しています。この行列に対する数値解は信頼できないことがわかります
-
右下（計算時間）: LU分解とコレスキー分解は $O(n^3)$ の参考線に沿っています。CG法は密行列でも行列-ベクトル積が $O(n^2)$ のため、直接法より速い場合があることが確認できます。疎行列の場合、この差はさらに顕著になります。コレスキー分解がLU分解より一貫して速いのは、対称性を活用して演算量が約半分に削減されているためです
手法選択のガイドライン
実際の問題で「どの手法を使うべきか」は、行列の性質と問題の要件によって異なります。以下に実践的な指針をまとめます。
| 状況 | 推奨手法 | 理由 |
|---|---|---|
| 小〜中規模の密行列 ($n \leq 10^4$) | LU分解 or コレスキー分解 | $O(n^3)$ でも十分高速。厳密解が得られる |
| 大規模疎行列、対称正定値 | 前処理付きCG法 | 疎構造を活用。1反復 $O(\text{nnz})$ |
| 大規模疎行列、非対称 | GMRES, BiCGSTAB | CG法の非対称行列への拡張 |
| 複数の右辺ベクトル | LU分解を保存して再利用 | 分解は1回、求解は $O(n^2)$ |
| 悪条件行列 | SVDまたは正則化 | 条件数に応じた打ち切りが可能 |
まとめ
本記事では、連立1次方程式を解く直接法と反復法の理論を体系的に解説しました。
- ガウス消去法は $O(n^3/3)$ の計算量で厳密解を得る直接法で、部分ピボット選択が数値的安定性に重要。乗数を1以下に抑えることで桁落ちのリスクを最小化する
- LU分解はガウス消去法を体系化したもので、$A = LU$ の分解を保存すれば複数の右辺に $O(n^2)$ で対応できる。対称正定値行列にはコレスキー分解が半分の計算量で適用できる
- 条件数 $\kappa(A)$ は解の精度の上限を決め、悪条件行列では直接法でも精度に限界がある。ヒルベルト行列のように条件数が指数的に増大する行列には特に注意が必要
- ヤコビ法・ガウス=ザイデル法は定常反復法で実装は容易だが収束が遅い。SOR法による加速も可能だが、最適パラメータの決定が課題
- CG法は対称正定値行列に対する最適な反復法で、$n$ ステップ以内に厳密解に収束し、疎行列に特に有効。前処理により収束を大幅に加速できる
数値線形代数は「単に方程式を解く」だけではなく、計算量・数値安定性・メモリ効率のバランスを考慮した手法選択が求められる奥深い分野です。現代の科学技術計算の多くは、最終的にはこの分野の知見に支えられています。
次のステップとして、以下の記事も参考にしてください。
- QR分解の理論 — 最小二乗問題の数値解法
- コレスキー分解の理論 — 対称正定値行列の分解
- 常微分方程式の数値解法 — 陰解法での線形代数