数値線形代数
概要
数値線形代数（numerical linear algebra）とは、連立 1 次方程式・最小二乗問題・固有値問題・特異値分解を計算機の上で解く方法を、演算の回数・問題の条件・方法の安定性の 3 つの観点から研究する分野である。行列 $A$ の条件数 $\kappa(A)=\|A\|\,\|A^{-1}\|$ は、データの相対誤差が解の相対誤差にどれだけ拡大されうるかを測り、残差が小さくても誤差は $\kappa(A)$ 倍まで大きくなりうる。大きな疎行列に使う Jacobi 法・Gauss–Seidel 法は反復行列のスペクトル半径が $1$ より小さいときに限って収束し、行について狭義の対角優位な行列では必ず収束する。冪乗法は、絶対値が最大の固有値がただ 1 つのとき、その固有ベクトルの方向へ収束する。
数値線形代数 は、連立 1 次方程式 $Ax=b$、最小二乗問題、固有値問題、特異値分解といった線形代数の問題を、計算機の上で速く正確に解く方法を研究する分野である。手で解ける大きさの問題なら 掃き出し法 で十分だが、$n$ が数千から数百万になると、演算の回数・記憶の量・丸め誤差の影響が問題の中心になる。
数値線形代数の議論は、次の 3 つの問いに分けると見通しがよい。
- 手間：演算の回数がどれだけか。行列の構造（対称、帯、疎）を使うと減らせるか。
- 問題の条件：データの小さな誤差が、答えにどれだけ拡大されうるか。これは問題そのものの性質で、解き方によらない。行列の 条件数 で測る。
- 方法の安定性：方法が計算の途中で入る丸め誤差をどれだけ拡大するか。
この記事では、行列のノルムと条件数を定義し、連立 1 次方程式の解の誤差が条件数で抑えられること（摂動定理、thm-nla-perturbation）を証明する。さらに、大きな疎行列に使う反復法が収束する条件（thm-nla-stationary・thm-nla-diagdom）と、最大の固有値を求める冪乗法の収束（thm-nla-power）を証明する。個々の分解（LU 分解・Cholesky 分解・QR 分解・特異値分解）は、それぞれの記事へ案内する（節「分野の地図」）。
定義と使える条件
以下、$\mathbb{K}$ は $\mathbb R$ または $\mathbb C$ とし、行列とベクトルの成分は $\mathbb{K}$ の元とする。
$x\in\mathbb{K}^n$ について
$$
\|x\|_1:=\sum_{i=1}^n\lvert x_i\rvert,\qquad \|x\|_2:=\Bigl(\sum_{i=1}^n\lvert x_i\rvert^2\Bigr)^{1/2},\qquad \|x\|_\infty:=\max_i\lvert x_i\rvert
$$
と定める。$\mathbb{K}^n$ のノルム $\|\cdot\|$ を 1 つ固定したとき、$A\in M_n(\mathbb{K})$ の 作用素ノルム（誘導されたノルム）を
$$
\|A\|:=\max_{x\ne0}\frac{\|Ax\|}{\|x\|}=\max_{\|x\|=1}\|Ax\|
$$
と定める。
単位球面 $\{\|x\|=1\}$ は有界閉集合で $x\mapsto\|Ax\|$ は連続なので、最大値は存在する。定義から $\|Ax\|\le\|A\|\,\|x\|$、$\|AB\|\le\|A\|\,\|B\|$、$\|I\|=1$ である。
$A=(a_{ij})\in M_n(\mathbb{K})$ について次が成り立つ。
- $\|A\|_\infty=\max_i\sum_j\lvert a_{ij}\rvert$（行の絶対値の和の最大）。
- $\|A\|_1=\max_j\sum_i\lvert a_{ij}\rvert$（列の絶対値の和の最大）。
- $\|A\|_2$ は $A$ の最大の特異値 $\sigma_1$ に等しく、$A$ が正則なら $\|A^{-1}\|_2=1/\sigma_n$ である（特異値分解 の定理「特異値とノルム・行列式・固有値」）。
要点：1 は各行の和で上から抑え、最大の行に合わせて符号（複素数では偏角）をそろえたベクトルで等号を作る。2 は単位ベクトルで等号を作る。
詳しい証明を開く
1：$\|x\|_\infty=1$ なら $\lvert(Ax)_i\rvert\le\sum_j\lvert a_{ij}\rvert\lvert x_j\rvert\le\sum_j\lvert a_{ij}\rvert$ なので、$\|A\|_\infty$ は右辺以下である。逆に、右辺の最大を与える行を $i$ とし、$x_j:=\overline{a_{ij}}/\lvert a_{ij}\rvert$（$a_{ij}=0$ なら $x_j:=1$）とおくと、$\|x\|_\infty=1$ で $(Ax)_i=\sum_j\lvert a_{ij}\rvert$ である。
2：$\|x\|_1=1$ なら $\|Ax\|_1\le\sum_j\lvert x_j\rvert\sum_i\lvert a_{ij}\rvert\le\max_j\sum_i\lvert a_{ij}\rvert$ である。最大を与える列を $j$ とし、$x$ を第 $j$ 単位ベクトルにとると等号が成り立つ。3 は引用した定理そのものである。$\square$
正則な $A\in M_n(\mathbb{K})$ の、ノルム $\|\cdot\|$ に関する 条件数 を
$$
\kappa(A):=\|A\|\,\|A^{-1}\|
$$
と定める。ノルムを明示するときは $\kappa_\infty(A)$、$\kappa_2(A)$ などと書く。正則でない $A$ では $\kappa(A):=\infty$ とする。
$1=\|I\|\le\|A\|\,\|A^{-1}\|$ なので $\kappa(A)\ge1$ であり、$c\ne0$ について $\kappa(cA)=\kappa(A)$ である。prop-nla-norm-formulas の 3 により $\kappa_2(A)=\sigma_1/\sigma_n$ で、これは 最小二乗法 の記事が正方でない行列にも定めている条件数 $\sigma_{\max}/\sigma_{\min}$ と、正方の場合に一致する。$\mathbb{K}^n$ のどの 2 つのノルムも定数倍で比べられるので、ノルムを替えても条件数は $n$ だけによる定数倍しか変わらない（たとえば $\|x\|_\infty\le\|x\|_2\le\sqrt n\,\|x\|_\infty$ から $\kappa_2(A)\le n\,\kappa_\infty(A)$）。
条件数は、スカラーの問題での相対条件数（誤差(数値解析) の定義「スカラー問題の相対条件数」）を、写像 $b\mapsto A^{-1}b$ に対して最悪の方向で測ったものにあたる。使える条件は次のとおりである。
- 以下の誤差の評価は、ノルムを 1 つ固定し、ベクトルと行列で同じノルム（作用素ノルム）を使ったときに成り立つ。
- 評価は 相対誤差 についての上界であり、$x=0$ や $b=0$ では使えない。
- 条件数は問題（行列 $A$）の性質であり、どの方法で解いたかによらない。方法の良し悪しは後退誤差（誤差(数値解析) の定義「前進誤差と後退誤差」）で別に測る。
分野の地図
数値線形代数の主な問題と、標準的な方法を表にする。演算の回数は、加減乗除の 1 回を 1 と数えた主要項である（BV04 Appendix C）。
| 問題 | 行列の条件 | 標準的な方法 | 演算の回数 | 根拠・詳しい記事 |
|---|---|---|---|---|
| $Ax=b$ | 正則（一般） | 部分ピボット選択つき LU 分解 $PA=LU$、前進・後退代入 | $\frac23n^3$（分解）$+2n^2$（代入） | LU分解、BV04 C.3.1（p. 668） |
| $Ax=b$ | 正定値対称（Hermite） | Cholesky 分解 $A=LL^{*}$ | $\frac13n^3$ | 正定値行列 の命題「Cholesky 分解の一意性」、BV04 C.3.2（p. 669） |
| $Ax=b$ | 帯幅 $k\ll n$ | 帯行列の LU 分解 | 約 $4nk^2$ | BV04 C.3.1（p. 669） |
| $Ax=b$ | 大きく疎、対角優位など | Jacobi 法・Gauss–Seidel 法などの反復法 | 1 回の反復で非零成分の数の定数倍 | thm-nla-stationary、thm-nla-diagdom |
| $\min\Vert Ax-b\Vert _2$ | 列が 1 次独立 | QR 分解 $A=QR$ と $R\hat x=Q^{*}b$ | — | 最小二乗法 の定理「QR 分解による最小二乗解」、Gram–Schmidtの直交化 |
| 最大の固有値 | 絶対値が最大の固有値が 1 つ | 冪乗法、Rayleigh 商 | 1 回の反復で行列とベクトルの積 1 回 | thm-nla-power、Rayleigh商 |
| 低階数近似・階数 | 任意 | 特異値分解 | — | 特異値分解 の定理「作用素ノルムでの最良の低階数近似」 |
ほかに、行列とベクトルの積は $2n^2$ 回、行列どうしの積は $2n^3$ 回（BV04 C.1.2、pp. 663–664）、三角行列の連立方程式は $n^2$ 回（同 C.2.1、pp. 664–665）である。逆行列を求めてから掛けるより、分解して代入で解くほうが手間が少ない（逆行列の計算は $\frac83n^3$ 回程度、同 p. 668）。
条件数と摂動定理
まず、単位行列に近い行列が正則であることを確かめる。
$M\in M_n(\mathbb{K})$ が $\|M\|<1$ を満たすならば、$I+M$ は正則であり
$$
\|(I+M)^{-1}\|\le\frac1{1-\|M\|}
$$
である。
任意の $y$ について $\|y\|\le\|(I+M)y\|+\|My\|\le\|(I+M)y\|+\|M\|\,\|y\|$ なので
$$
(1-\|M\|)\,\|y\|\le\|(I+M)y\|
$$
である。$1-\|M\|>0$ なので、$(I+M)y=0$ なら $y=0$ であり、正方行列 $I+M$ は正則である。$y=(I+M)^{-1}z$ を代入すると $\|(I+M)^{-1}z\|\le\|z\|/(1-\|M\|)$ である。$\square$
$A\in M_n(\mathbb{K})$ を正則、$b\ne0$、$Ax=b$ とする。$E\in M_n(\mathbb{K})$、$f\in\mathbb{K}^n$ が
$$
\|A^{-1}\|\,\|E\|<1
$$
を満たすとき、$A+E$ は正則であり、$(A+E)\hat x=b+f$ の解 $\hat x$ は
$$
\frac{\|\hat x-x\|}{\|x\|}\le\frac{\kappa(A)}{1-\kappa(A)\,\varepsilon_A}\,\bigl(\varepsilon_A+\varepsilon_b\bigr),\qquad \varepsilon_A:=\frac{\|E\|}{\|A\|},\quad\varepsilon_b:=\frac{\|f\|}{\|b\|}
$$
を満たす。とくに $E=0$ なら $\dfrac{\|\hat x-x\|}{\|x\|}\le\kappa(A)\,\dfrac{\|f\|}{\|b\|}$ である。
$A+E=A(I+A^{-1}E)$ であり、$\|A^{-1}E\|\le\|A^{-1}\|\,\|E\|<1$ なので、lem-nla-identity により $I+A^{-1}E$ は正則で、$A+E$ も正則である。
$(A+E)\hat x=b+f$ から $Ax=b$ を引くと $(A+E)(\hat x-x)=f-Ex$ である。よって
$$
\hat x-x=(I+A^{-1}E)^{-1}A^{-1}(f-Ex)
$$
であり、lem-nla-identity により
$$
\|\hat x-x\|\le\frac{\|A^{-1}\|\,(\|f\|+\|E\|\,\|x\|)}{1-\|A^{-1}\|\,\|E\|}
$$
である。$b\ne0$ なので $x\ne0$ で、$\|b\|=\|Ax\|\le\|A\|\,\|x\|$ から $\|f\|=\varepsilon_b\|b\|\le\varepsilon_b\|A\|\,\|x\|$ である。また $\|E\|=\varepsilon_A\|A\|$、$\|A^{-1}\|\,\|E\|=\kappa(A)\,\varepsilon_A$ である。これらを代入して $\|x\|$ で割ると
$$
\frac{\|\hat x-x\|}{\|x\|}\le\frac{\|A^{-1}\|\,\|A\|\,(\varepsilon_b+\varepsilon_A)}{1-\kappa(A)\,\varepsilon_A}=\frac{\kappa(A)}{1-\kappa(A)\,\varepsilon_A}\,(\varepsilon_A+\varepsilon_b)
$$
を得る。$\square$
この定理は、データ $A,b$ の相対誤差が $\varepsilon$ 程度なら、解の相対誤差はおよそ $\kappa(A)\,\varepsilon$ 程度まで拡大されうることを言っている。$\kappa(A)\approx10^k$ なら、10 進でおよそ $k$ 桁を失う覚悟が要る。
$A$ を正則、$Ax=b\ne0$ とし、$\hat x$ を任意のベクトル、$r:=b-A\hat x$ を残差とする。このとき
$$
\frac{\|\hat x-x\|}{\|x\|}\le\kappa(A)\,\frac{\|r\|}{\|b\|}
$$
である。
$A\hat x=b-r$ なので、thm-nla-perturbation を $E=0$、$f=-r$ として使う。$\square$
計算した $\hat x$ の残差は手元で計算できる量なので、この系は「残差が小さいこと」から「誤差が小さいこと」を言うための道具である。ただし倍率 $\kappa(A)$ が掛かるので、条件数が大きいと残差が小さくても誤差は大きくなりうる（ex-nla-nearly-singular）。
後退安定性との組み合わせ：方法の計算した $\hat x$ が、ある $E$ について $(A+E)\hat x=b$ を正確に満たし、$\|E\|/\|A\|$ が丸めの単位 $u$ の定数倍以下であるとき、その方法を 後退安定 という。このとき thm-nla-perturbation により、相対誤差はおよそ $\kappa(A)$ と $u$ の定数倍の積で抑えられる。「前進誤差は条件数と後退誤差の積」という一般の原理（誤差(数値解析) の命題「前進誤差は条件数と後退誤差の積で決まる」）の、連立 1 次方程式での形である。部分ピボット選択つきの LU 分解では、この定数が途中の成分の増大に依存する（LU分解 の定理「部分ピボット選択での増大率の上界」）。
反復法の収束
大きな疎行列では、LU 分解の途中で零だった成分が非零になり（fill-in）、記憶も手間も増える。そこで、行列とベクトルの積だけを繰り返す 反復法 を使う。$A=M-N$（$M$ は正則で、$My=c$ が解きやすい行列）と分け、
$$
x_{k+1}:=M^{-1}(Nx_k+b)\qquad(k=0,1,2,\dots)
$$
とする方法を 定常反復法 といい、$T:=M^{-1}N$ を 反復行列 という。$A=D+L+U$（$D$ は対角部分、$L$ は狭義下三角部分、$U$ は狭義上三角部分）と書き、$D$ の対角成分はすべて $0$ でないとする。
| 方法 | $M$ | 反復行列 $T$ | 成分ごとの式 |
|---|---|---|---|
| Jacobi 法 | $D$ | $T_J=-D^{-1}(L+U)$ | $x^{(k+1)}_i=\bigl(b_i-\sum_{j\ne i}a_{ij}x^{(k)}_j\bigr)/a_{ii}$ |
| Gauss–Seidel 法 | $D+L$ | $T_{GS}=-(D+L)^{-1}U$ | $x^{(k+1)}_i=\bigl(b_i-\sum_{j< i}a_{ij}x^{(k+1)}_j-\sum_{j>i}a_{ij}x^{(k)}_j\bigr)/a_{ii}$ |
Gauss–Seidel 法は、新しく求めた成分をその場で使う。$x_*$ を $Ax=b$ の解とすると $Mx_*=Nx_*+b$ なので、誤差 $e_k:=x_k-x_*$ は $e_{k+1}=Te_k$、すなわち $e_k=T^ke_0$ を満たす。
$A=M-N$ を正則、$M$ を正則とし、$T:=M^{-1}N$ とする。次は同値である。
- どの $b$ とどの初期値 $x_0$ についても、$x_k$ は $Ax=b$ の解に収束する。
- $T$ のスペクトル半径 $\rho(T):=\max\{\lvert\lambda\rvert:\lambda\text{ は }T\text{ の複素固有値}\}$ が $1$ より小さい。
要点：誤差は $e_k=T^ke_0$ を満たす。$\rho(T)<1$ なら Gelfand の公式で $\|T^k\|\to0$、$\rho(T)\ge1$ なら絶対値 $1$ 以上の固有値の固有ベクトル（実の場合はその実部か虚部）を $e_0$ にとる。
詳しい証明を開く
$e_k=T^ke_0$ であり、$e_0=x_0-x_*$ は $x_0$ を選べば任意のベクトルにできる。よって 1 は「すべての $e_0$ で $T^ke_0\to0$」と同値である。
2 ⇒ 1：$T$ を $\mathbb C^n$ の作用素とみると、スペクトル半径 の定理「Gelfand のスペクトル半径公式」により $\lim_k\|T^k\|^{1/k}=\rho(T)<1$ である（行列のノルムはどの 2 つも定数倍で比べられ、定数の $k$ 乗根は $1$ に近づくので、この極限は行列のノルムの選び方によらない。$\mathbb{K}=\mathbb R$ で $\mathbb R^n$ の作用素ノルムを使うときも同じである）。$\rho(T)< q< 1$ となる $q$ をとると、十分大きい $k$ で $\|T^k\|\le q^k$ となり、$\|e_k\|\le q^k\|e_0\|\to0$ である（Neumann級数 の命題「線形反復の収束」と同じ結論）。
1 ⇒ 2：$\rho(T)\ge1$ とし、$Tv=\lambda v$、$\lvert\lambda\rvert\ge1$、$v\in\mathbb C^n$、$v\ne0$ とする。$T^kv=\lambda^kv$ は $0$ に収束しない。$\mathbb{K}=\mathbb C$ ならこれで 1 が崩れる。$\mathbb{K}=\mathbb R$ なら $v=v'+iv''$（$v',v''\in\mathbb R^n$）と分けると、$T$ が実行列なので $T^kv=T^kv'+iT^kv''$ であり、$T^kv'$ と $T^kv''$ の少なくとも一方は $0$ に収束しない。その実ベクトルを $e_0$ にとれば 1 が崩れる。$\square$
実際に条件 2 を確かめるには固有値を求める必要があるが、行列の形から直接分かる十分条件がある。$A$ が 行について狭義の対角優位 であるとは、すべての $i$ で $\lvert a_{ii}\rvert>\sum_{j\ne i}\lvert a_{ij}\rvert$ であることをいう。
$A$ が行について狭義の対角優位であるとし、
$$
q:=\max_i\frac{\sum_{j\ne i}\lvert a_{ij}\rvert}{\lvert a_{ii}\rvert}<1
$$
とおく。このとき $A$ は正則で、Jacobi 法と Gauss–Seidel 法はどの初期値からも $Ax=b$ の解に収束し、誤差は
$$
\|e_k\|_\infty\le q^k\,\|e_0\|_\infty
$$
を満たす。
Jacobi 法：$T_J$ の $(i,j)$ 成分は $j\ne i$ で $-a_{ij}/a_{ii}$、$j=i$ で $0$ なので、prop-nla-norm-formulas の 1 により $\|T_J\|_\infty=q$ である。よって $\|e_k\|_\infty=\|T_J^ke_0\|_\infty\le q^k\|e_0\|_\infty$ である。
Gauss–Seidel 法：$e':=T_{GS}e$ とおくと、成分ごとの式から $a_{ii}e'_i=-\sum_{j< i}a_{ij}e'_j-\sum_{j>i}a_{ij}e_j$ である。$\alpha_i:=\sum_{j< i}\lvert a_{ij}\rvert/\lvert a_{ii}\rvert$、$\beta_i:=\sum_{j>i}\lvert a_{ij}\rvert/\lvert a_{ii}\rvert$ とおくと $\alpha_i+\beta_i\le q$ である。$\lvert e'_i\rvert$ が最大となる $i$ をとると
$$
\|e'\|_\infty=\lvert e'_i\rvert\le\alpha_i\|e'\|_\infty+\beta_i\|e\|_\infty
$$
であり、$\alpha_i\le q<1$ なので $\|e'\|_\infty\le\frac{\beta_i}{1-\alpha_i}\|e\|_\infty$ である。$\beta_i\le q-\alpha_i\le q-q\alpha_i=q(1-\alpha_i)$ なので、$\|e'\|_\infty\le q\|e\|_\infty$ である。これを $k$ 回使えばよい。
正則性：$Ay=0$ なら、$b=0$ の方程式の解 $y$ を初期値とする Jacobi 法の列はつねに $y$ だが、上の評価から $0$ に収束するので $y=0$ である。$\square$
$\|T\|<1$ は $\rho(T)<1$ の十分条件にすぎない。Jacobi 法と Gauss–Seidel 法のどちらか一方だけが収束する行列もあり、一般の行列についてどちらが速いとは言えない。
冪乗法
固有値問題では、大きな行列の固有値のうち絶対値が最大のものだけが欲しいことが多い。冪乗法 は、ベクトルに $A$ を掛けて正規化することを繰り返す方法である：
$$
x_{k+1}:=\frac{Ax_k}{\|Ax_k\|_2},\qquad \mu_k:=x_k^{*}Ax_k\qquad(k=0,1,2,\dots).
$$
$\mu_k$ は $x_k$ での Rayleigh商 である（$\|x_k\|_2=1$ なので分母は $1$）。
$A\in M_n(\mathbb C)$ が固有ベクトルの基底 $v_1,\dots,v_n$（$Av_j=\lambda_jv_j$）をもち、
$$
\lvert\lambda_1\rvert>\lvert\lambda_2\rvert\ge\cdots\ge\lvert\lambda_n\rvert
$$
とする。$\|x_0\|_2=1$ を $x_0=\sum_jc_jv_j$ と書いたとき $c_1\ne0$ とし、$r:=\lvert\lambda_2\rvert/\lvert\lambda_1\rvert<1$、$u:=v_1/\|v_1\|_2$ とおく。このとき $Ax_k\ne0$ がすべての $k$ で成り立ち、絶対値 $1$ の複素数 $\theta_k$ と、$k$ によらない定数 $C$ があって
$$
\|\theta_kx_k-u\|_2\le C\,r^k,\qquad\lvert\mu_k-\lambda_1\rvert\le2\|A\|_2\,C\,r^k
$$
が成り立つ。とくに $x_k$ の張る直線は固有ベクトル $v_1$ の直線に近づき、$\mu_k\to\lambda_1$ である。
帰納法で $x_k=A^kx_0/\|A^kx_0\|_2$ である（$A^kx_0\ne0$ は下で示す）。
$$
A^kx_0=\sum_jc_j\lambda_j^kv_j=c_1\lambda_1^k\,w_k,\qquad w_k:=v_1+\sum_{j\ge2}\frac{c_j}{c_1}\Bigl(\frac{\lambda_j}{\lambda_1}\Bigr)^kv_j
$$
と書ける（$\lambda_1\ne0$ に注意）。$C_0:=\sum_{j\ge2}\lvert c_j/c_1\rvert\,\|v_j\|_2$ とおくと $\|w_k-v_1\|_2\le C_0r^k$ である。$v_1,\dots,v_n$ は 1 次独立で $c_1\ne0$ なので $A^kx_0\ne0$、すなわち $w_k\ne0$ である。
$\phi_k:=c_1\lambda_1^k/\lvert c_1\lambda_1^k\rvert$ とおくと $x_k=\phi_k\,w_k/\|w_k\|_2$ なので、$\theta_k:=\overline{\phi_k}$ とすれば $\theta_kx_k=w_k/\|w_k\|_2$ である。$0$ でないベクトル $a,b$ について
$$
\Bigl\|\frac a{\|a\|}-\frac b{\|b\|}\Bigr\|\le\frac{\|a-b\|}{\|b\|}+\Bigl\lvert\frac1{\|a\|}-\frac1{\|b\|}\Bigr\rvert\,\|a\|=\frac{\|a-b\|+\bigl\lvert\|b\|-\|a\|\bigr\rvert}{\|b\|}\le\frac{2\|a-b\|}{\|b\|}
$$
なので、$a=w_k$、$b=v_1$ として $\|\theta_kx_k-u\|_2\le2C_0r^k/\|v_1\|_2$ である。$C:=2C_0/\|v_1\|_2$ とおく。
$\lvert\theta_k\rvert=1$ なので $\mu_k=(\theta_kx_k)^{*}A(\theta_kx_k)$ であり、$u^{*}Au=\lambda_1u^{*}u=\lambda_1$ である。$y:=\theta_kx_k$ とおくと $\|y\|_2=\|u\|_2=1$ で
$$
\lvert\mu_k-\lambda_1\rvert=\lvert y^{*}Ay-u^{*}Au\rvert=\lvert y^{*}A(y-u)+(y-u)^{*}Au\rvert\le2\|A\|_2\,\|y-u\|_2\le2\|A\|_2\,C\,r^k
$$
である。$\square$
収束の速さは比 $r=\lvert\lambda_2\rvert/\lvert\lambda_1\rvert$ で決まり、$r$ が $1$ に近いと遅い。$A$ が Hermite 行列なら、Rayleigh 商の誤差は固有ベクトルの誤差の 2 乗の程度になり（Rayleigh商 の定理「2 乗の精度」）、$\mu_k$ は $r^{2k}$ の速さで $\lambda_1$ に近づく（ex-nla-power）。
計算の手順
連立 1 次方程式 $Ax=b$ を解くときの方法の選び方と、結果の確かめ方は次のとおりである。
- 行列の構造を調べる。正定値対称（Hermite）なら Cholesky 分解、帯行列なら帯の LU 分解、一般の密な行列なら部分ピボット選択つきの LU 分解を使う（節「分野の地図」の表）。逆行列は求めない。
- 同じ $A$ で右辺が何本もあるときは、分解を 1 回だけ行い、右辺ごとに前進・後退代入をする。
- 行列が大きく疎で、分解すると成分が埋まって記憶が足りないときは反復法を使う。対角優位なら Jacobi 法・Gauss–Seidel 法が収束する（thm-nla-diagdom）。一般には反復行列のスペクトル半径が $1$ より小さいことが収束の条件である（thm-nla-stationary）。
- 反復を止める判定：相対残差 $\|b-Ax_k\|/\|b\|$ が許容値 $\tau$ 以下になったら止める。cor-nla-residual により、そのときの相対誤差は $\kappa(A)\,\tau$ 以下である。反復の上限回数も決めておく。
- 条件数を見積もる。$\kappa(A)$ と丸めの単位 $u$ の積が $1$ に近いなら、どの方法でも答えの桁はほとんど信用できない。問題の定式化（変数の尺度、基底の選び方）を見直す。
固有値の最大のものだけが欲しいときは冪乗法を使い、$\mu_k$ の変化が許容値以下になったら止める。比 $\lvert\lambda_2\rvert/\lvert\lambda_1\rvert$ が $1$ に近いと遅いので、その場合はほかの方法を考える。
例題
$H_n:=\bigl(\frac1{i+j-1}\bigr)_{i,j=1}^n$ とする。
使う公式：def-nla-condition、prop-nla-norm-formulas の 1、thm-nla-perturbation。
答え：$\kappa_\infty(H_n)$ は $n=2,3,4,5,6$ で $27$、$748$、$28375$、$943656$、$29070279$ であり、$n=8$ で約 $3.39\times10^{10}$ である。
$n=3$ での摂動の計算を開く
$H_3^{-1}=\begin{pmatrix}9&-36&30\\-36&192&-180\\30&-180&180\end{pmatrix}$ なので $\|H_3\|_\infty=1+\frac12+\frac13=\frac{11}6$、$\|H_3^{-1}\|_\infty=36+192+180=408$、$\kappa_\infty(H_3)=748$ である。$x=(1,1,1)^{\mathsf T}$ のとき $b=H_3x=(\frac{11}6,\frac{13}{12},\frac{47}{60})^{\mathsf T}$ である。右辺の第 3 成分だけを $10^{-3}$ 増やすと、解の変化は $H_3^{-1}(0,0,10^{-3})^{\mathsf T}=(0.03,-0.18,0.18)^{\mathsf T}$ である。右辺の相対変化 $\frac{10^{-3}}{11/6}\approx5.4545\times10^{-4}$ に対して、解の相対変化は $0.18$ で、拡大率は $330$ である。これは thm-nla-perturbation の上界 $\kappa_\infty(H_3)=748$ 以下である。
$A=\begin{pmatrix}1&1\\1&1+10^{-4}\end{pmatrix}$、$x=(1,1)^{\mathsf T}$、$b=Ax=(2,\ 2.0001)^{\mathsf T}$ とする。
使う公式：cor-nla-residual。
答え：$A^{-1}=\begin{pmatrix}10001&-10000\\-10000&10000\end{pmatrix}$、$\kappa_\infty(A)=2.0001\times20001=40004.0001$ である。近似解 $\hat x=(2,0)^{\mathsf T}$ の残差は $r=b-A\hat x=(0,\ 10^{-4})^{\mathsf T}$ で、相対残差は $\|r\|_\infty/\|b\|_\infty\approx5.0\times10^{-5}$ と小さいが、相対誤差は $\|\hat x-x\|_\infty/\|x\|_\infty=1$ である。cor-nla-residual の上界は $\kappa_\infty(A)\times\frac{10^{-4}}{2.0001}=2.0001$ で、誤差 $1$ はこの上界以下である。
$A=\begin{pmatrix}4&-1&0\\-1&4&-1\\0&-1&4\end{pmatrix}$、$b=(3,2,3)^{\mathsf T}$（解は $x_*=(1,1,1)^{\mathsf T}$）、$x_0=0$ とする。
使う公式：thm-nla-diagdom。$A$ は行について狭義の対角優位で、$q=\max(\frac14,\frac24,\frac14)=\frac12$ である。
答え：Jacobi 法の反復は $x_1=(\frac34,\frac12,\frac34)$、$x_2=(\frac78,\frac78,\frac78)$、$x_3=(\frac{31}{32},\frac{15}{16},\frac{31}{32})$ と進み、$\|e_{10}\|_\infty=2^{-15}\approx3.0518\times10^{-5}$ である。定理の上界は $q^{10}\|e_0\|_\infty=2^{-10}\approx9.7656\times10^{-4}$ である。Gauss–Seidel 法では $x_1=(0.75,\ 0.6875,\ 0.921875)$ で、$\|e_{10}\|_\infty\approx4.6566\times10^{-9}$ とずっと速い。
スペクトル半径との比較を開く
$T_J=\frac14\begin{pmatrix}0&1&0\\1&0&1\\0&1&0\end{pmatrix}$ の固有値は $0,\pm\frac{\sqrt2}4$ なので $\rho(T_J)=\frac{\sqrt2}4\approx0.35355$ である。$T_{GS}=-(D+L)^{-1}U$ の固有値は $0,0,\frac18$ で $\rho(T_{GS})=\frac18$ である。どちらも $q=\frac12$ より小さく、thm-nla-diagdom の評価は上界にすぎない。Jacobi 法の誤差は 2 回の反復でちょうど $\frac18=\rho(T_J)^2$ 倍になっている（$\|e_2\|_\infty=\frac18$、$\|e_4\|_\infty=\frac1{64}$）。
$A=\begin{pmatrix}2&1\\1&3\end{pmatrix}$、$x_0=(1,0)^{\mathsf T}$ とする。固有値は $\lambda_1=\frac{5+\sqrt5}2\approx3.618034$、$\lambda_2=\frac{5-\sqrt5}2\approx1.381966$ で、$r\approx0.381966$ である。
使う公式：thm-nla-power。
答え：$A^kx_0$ は $(2,1)$、$(5,5)$、$(15,20)$ と進み、$\mu_1=3$、$\mu_2=\frac72$、$\mu_3=\frac{18}5$ である。$\mu_5\approx3.6176471$（誤差 $3.87\times10^{-4}$）、$\mu_8\approx3.6180328$（誤差 $1.20\times10^{-6}$）である。$A$ は対称なので、誤差はおよそ反復ごとに $r^2\approx0.1459$ 倍になる。
使える条件と反例
| 外す条件 | 反例 | 成り立たなくなること |
|---|---|---|
| $\Vert A^{-1}\Vert \,\Vert E\Vert <1$（thm-nla-perturbation） | $A=I$、$E=-I$ | $A+E$ が正則であること |
| 行について狭義の対角優位（thm-nla-diagdom） | $A=\begin{pmatrix}1&2\\2&1\end{pmatrix}$ の Jacobi 法 | 反復が収束すること |
| $\lvert\lambda_1\rvert>\lvert\lambda_2\rvert$（thm-nla-power） | $A=\operatorname{diag}(1,-1)$、$x_0=\frac1{\sqrt2}(1,1)^{\mathsf T}$ | $x_k$ の直線が固有ベクトルの直線に近づくこと |
| $c_1\ne0$（thm-nla-power） | $A=\operatorname{diag}(2,1)$、$x_0=(0,1)^{\mathsf T}$ | $\mu_k\to\lambda_1$ |
4 つの確かめを開く
1 行目：$\|A^{-1}\|\,\|E\|=1$ で、$A+E=0$ は正則でない。2 行目：$q=2>1$ である。$T_J=\begin{pmatrix}0&-2\\-2&0\end{pmatrix}$ の固有値は $\pm2$ で $\rho(T_J)=2$ なので、thm-nla-stationary により収束しない初期値がある。実際 $b=(3,3)^{\mathsf T}$（解は $(1,1)^{\mathsf T}$）、$x_0=0$ とすると、反復は $(3,3)$、$(-3,-3)$、$(9,9)$、$(-15,-15)$、$(33,33)$ と発散する。3 行目：固有ベクトルの基底 $e_1,e_2$ があり $c_1=\frac1{\sqrt2}\ne0$ だが、$\lvert\lambda_1\rvert=\lvert\lambda_2\rvert=1$ である。$x_k$ は $\frac1{\sqrt2}(1,-1)^{\mathsf T}$ と $\frac1{\sqrt2}(1,1)^{\mathsf T}$ を交互にとり、どちらの直線も固有ベクトルの直線ではない。$\mu_k=0$ はどちらの固有値でもない。4 行目：$x_k=(0,1)^{\mathsf T}$ がずっと続き、$\mu_k=1\ne\lambda_1=2$ である。丸めのある計算では、$c_1$ の方向の成分が丸めで入って、やがて $\lambda_1$ に向かうことが多い。
よくある誤り
| 誤り | 正しくは | なぜ |
|---|---|---|
| $A^{-1}$ を求めてから $x=A^{-1}b$ とする | LU 分解などで分解し、代入で解く | 手間が約 4 倍になる（節「分野の地図」） |
| 残差が小さいので解は正確だとする | 相対誤差は $\kappa(A)$ 倍まで大きくなりうる | cor-nla-residual、ex-nla-nearly-singular |
| $\det A$ が小さいので条件が悪いとする | 条件数で判断する。$\det(10^{-1}I_n)=10^{-n}$ だが $\kappa(10^{-1}I_n)=1$ | $\kappa(cA)=\kappa(A)$ |
| 最小二乗問題を正規方程式 $A^{*}Ax=A^{*}b$ で解く | QR 分解を使う | 正規方程式の条件数は $\kappa_2(A)^2$（最小二乗法 の命題「正規方程式の条件数」） |
| Jacobi 法はどの行列でも収束するとする | $\rho(T)<1$ が必要十分 | thm-nla-stationary、ex-nla-counterexamples |
補足
この記事で扱った方法のほかに、数値線形代数の中心には次のものがある（この記事では扱わない）。正定値対称な大きな疎行列には共役勾配法、一般の行列には GMRES などの Krylov 部分空間法が使われ、反復法を速くするための前処理が重要である。すべての固有値を求めるには、直交変換で行列を簡単な形に直してから反復する QR 法などが使われ、特異値分解も反復で計算される。精度の保証が必要なときは、区間演算 で解を含む区間を計算する方法がある。
丸め誤差の解析の基本（標準浮動小数点モデル、前進誤差と後退誤差）は 誤差(数値解析)、行列のスペクトル半径とノルムの関係は スペクトル半径、反復の収束の一般論は Neumann級数 にある。
関連項目
参考文献
Mathpediaは寄付と、参考文献の書籍リンク（Amazonアソシエイト）の紹介料で運営されています。 支援について ／ 寄付する
