# Gauss 消去と LU 分解：ピボット選択・成長因子・後退誤差解析 > LU 分解の存在と一意性、部分ピボット選択つき分解 PA = LU の存在と演算量を示し、三角方程式・LU 分解・Gauss 消去の後退誤差 |ΔA| ≤ γ|L||U| を証明する。成長因子の上界 2^(n−1)、列対角優位な行列、混合精度の反復改良の収束も扱う。 > https://rikai.mugen-giken.com/mathematics/numerical-linear-algebra/lu-and-pivoting ## 0. この記事の要点 連立一次方程式をコンピュータで解くとき、答えは丸め誤差のせいで少しずれます。その「少し」がどれくらいなのか、問題が難しいせいなのか、解き方が悪いせいなのかを、きちんと見分けるのがこの記事の話です。Gauss の消去法を「行列を 2 つの三角形の積に分ける」手順として読み直し、行の入れ替え（ピボット選択）がなぜ必要か、入れ替えても数値が極端に膨らむ例があるか、を順に調べます。読み終えると、`numpy.linalg.solve` が返す解をどこまで信頼してよいかを、式を使って説明できるようになります。 - 先頭の主小行列式がすべて $0$ でないことと、$A=LU$（$L$ は対角が $1$ の下三角、$U$ は正則な上三角）と一意に書けることは同値です。**部分ピボット選択**を使えば、どんな正方行列も $PA=LU$, $|l_{ij}|\le1$ と分解できます。演算量は $\tfrac23n^3+O(n^2)$ です。 - 浮動小数点数で計算した因子は $\hat L\hat U=A+\Delta A$, $|\Delta A|\le\gamma_n|\hat L||\hat U|$ を満たし、前進代入・後退代入と合わせた解は $(A+\Delta A)\hat{\boldsymbol{x}}=\boldsymbol{b}$, $|\Delta A|\le\gamma_{3n}|\hat L||\hat U|$ を満たします。丸めの基本は既存の「1 ステップの丸め誤差」の補題に帰着します。 - 後退誤差を $|A|$ の尺度で小さくできるかは $|\hat L||\hat U|$ の大きさ、つまり**成長因子** $\rho_n$ で決まります。部分ピボット選択では $\rho_n\le2^{n-1}$ で、この上界に達する行列（Wilkinson の例）があります。ただし乱数行列では $\rho_n$ はずっと小さく、実用上は後退安定とみなされます。 - 列対角優位な行列はピボット選択なしで消去でき、$\rho_n\le2$ です。Schur 補行列が列対角優位を保つからです。 - **反復改良**: 補正方程式を誤差 $E_k$ のある分解で解いても、$\|A^{-1}E_k\|\le\theta\lt1/2$ なら誤差は 1 回ごとに $\theta/(1-\theta)$ 倍以下に縮みます。binary32 の LU と binary64 の残差を組み合わせた混合精度の解法の根拠です。 ## 1. 動機 連立一次方程式 $A\boldsymbol{x}=\boldsymbol{b}$ を解く標準的な方法は、高校以来の掃き出し法（Gauss の消去法）です。数値計算では、これを「$A$ を下三角 $L$ と上三角 $U$ の積に分解し、$L\boldsymbol{y}=\boldsymbol{b}$, $U\boldsymbol{x}=\boldsymbol{y}$ を順に解く」と読み直します。分解は $O(n^3)$、三角方程式は $O(n^2)$ なので、同じ $A$ で右辺を何本も解くときに分解を使い回せます。SciPy の `scipy.linalg.solve`、NumPy の `numpy.linalg.solve` の中身は LAPACK の `getrf`（部分ピボット選択つき LU 分解）と `getrs`（三角方程式）です。 この章の問いは、浮動小数点数で計算した解がどれだけ信頼できるかです。[行列のノルムと線形方程式の条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning)の章で、前進誤差は条件数 × 後退誤差で抑えられることを見ました（）。したがって示すべきは、Gauss 消去の**後退誤差**が小さいことです。答えは「ピボット選択をすれば、ほとんどの場合に小さい」で、その「ほとんど」の中身が成長因子です。 ## 2. LU 分解の存在と一意性 ### 2.1. ピボット選択なしの LU 分解 $A\in M_n(\mathbb{K})$ の左上の $k$ 次の正方部分行列を $A_k$ と書き、**先頭主小行列**と呼びます。 $A\in M_n(\mathbb{K})$ について、対角成分がすべて $1$ の下三角行列 $L$（**単位下三角行列**）と上三角行列 $U$ で $A=LU$ と書くことを、$A$ の **LU 分解**という。置換行列 $P$ を使って $PA=LU$ と書くことを**ピボット選択つき LU 分解**という。 三角行列について、次の 2 つの事実を何度も使います。 - 三角行列の行列式は対角成分の積です。$n$ についての帰納法で、ブロック三角行列の行列式の公式（）を $1+(n-1)$ の分割に使えば出ます。 - 単位下三角行列 $L$ は正則で、$L^{-1}$ も単位下三角です。$LX=I$ の第 $j$ 列を前進代入で解くと、第 $1$〜$j-1$ 成分は $0$、第 $j$ 成分は $1$ になるからです。単位下三角行列どうしの積も単位下三角です（積の $(i,j)$ 成分 $\sum_kl_{ik}m_{kj}$ は $i\lt j$ なら $k\le i\lt j$ の項しか残らず $0$）。上三角についても同様です。 $A\in M_n(\mathbb{K})$ について、次は同値である。 1. 先頭主小行列式がすべて $0$ でない: $\det A_k\ne0$（$k=1,\ldots,n$）。 2. 単位下三角行列 $L$ と、対角成分がすべて $0$ でない上三角行列 $U$ で $A=LU$ と書ける。 さらにこのとき、2 の $L$, $U$ は一意で、$u_{11}\cdots u_{kk}=\det A_k$ である。 **2 ⇒ 1.** $A=LU$ をブロックに分けると、$L$ が下三角、$U$ が上三角なので左上のブロックについて $A_k=L_kU_k$ です（$(LU)_{ij}=\sum_ml_{im}u_{mj}$ で、$i,j\le k$ なら $l_{im}\ne0$ には $m\le i\le k$ が必要）。積の定理（）と三角行列の行列式から $\det A_k=\det L_k\det U_k=u_{11}\cdots u_{kk}\ne0$ です。 **1 ⇒ 2.** $n$ についての帰納法です。$n=1$ は $L=(1)$, $U=(a_{11})$ です。$n-1$ で正しいとし、 $$ A=\begin{pmatrix}A_{n-1}&\boldsymbol{b}\\\boldsymbol{c}^{\mathsf{T}}&a_{nn}\end{pmatrix} $$ と分けます。$A_{n-1}$ の先頭主小行列は $A_1,\ldots,A_{n-1}$ なので、帰納法の仮定から $A_{n-1}=L'U'$（$L'$ は単位下三角、$U'$ は対角が $0$ でない上三角）です。$U'$, $L'$ は正則なので $$ L=\begin{pmatrix}L'&\boldsymbol{0}\\\boldsymbol{c}^{\mathsf{T}}U'^{-1}&1\end{pmatrix},\qquad U=\begin{pmatrix}U'&L'^{-1}\boldsymbol{b}\\\boldsymbol{0}^{\mathsf{T}}&\delta\end{pmatrix},\qquad\delta=a_{nn}-\boldsymbol{c}^{\mathsf{T}}A_{n-1}^{-1}\boldsymbol{b} $$ とおくと、ブロックの積を計算して $LU=A$ です（右下は $\boldsymbol{c}^{\mathsf{T}}U'^{-1}L'^{-1}\boldsymbol{b}+\delta=\boldsymbol{c}^{\mathsf{T}}A_{n-1}^{-1}\boldsymbol{b}+\delta=a_{nn}$）。積の定理とブロック三角行列の行列式（）から $\det A=\det L\det U=\det U'\cdot\delta$ で、$\det A=\det A_n\ne0$ なので $\delta\ne0$ です。 **一意性.** $L_1U_1=L_2U_2$ なら $L_2^{-1}L_1=U_2U_1^{-1}$ です。左辺は単位下三角、右辺は上三角なので、両辺は対角行列で、しかも対角が $1$、すなわち $I$ です。よって $L_1=L_2$, $U_1=U_2$。最後の等式は 2 ⇒ 1 の計算です。 構成から、$U$ の対角成分 $u_{kk}=\det A_k/\det A_{k-1}$ は消去の $k$ 段目の**ピボット**です。$\delta$ は $A$ の $A_{n-1}$ に関する Schur 補行列（）で、消去の各段で残る右下のブロックは、すべて Schur 補行列です（左上のブロックを消去したあとに右下に残る、「$B-\boldsymbol{l}\boldsymbol{w}^{\mathsf{T}}$」の形の小さい行列のことです）。 ### 2.2. 部分ピボット選択 $\det A_1=a_{11}=0$ なら、正則な $A=\begin{pmatrix}0&1\\1&0\end{pmatrix}$ でも LU 分解はありません。行を入れ替えれば解決します。数値計算では、さらに**ピボットの絶対値が最大になるように**行を入れ替えます。 任意の $A\in M_n(\mathbb{K})$ に対し、置換行列 $P$、単位下三角行列 $L$、上三角行列 $U$ があって $$ PA=LU,\qquad|l_{ij}|\le1\quad(i\gt j) $$ が成り立つ。$A$ は特異でもよい（そのとき $U$ の対角に $0$ が現れる）。 $n$ についての帰納法です。$n=1$ は $P=L=(1)$, $U=A$ です。$n\ge2$ とし、第 1 列で絶対値が最大の成分の行 $r$ を選び（**部分ピボット選択**）、第 1 行と第 $r$ 行を入れ替える置換行列を $P_1$ とします。 $$ P_1A=\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{v}&B\end{pmatrix},\qquad|v_i|\le|\alpha| $$ です。$\alpha\ne0$ なら $\boldsymbol{l}=\boldsymbol{v}/\alpha$、$\alpha=0$ なら（第 1 列は $0$ なので $\boldsymbol{v}=\boldsymbol{0}$）$\boldsymbol{l}=\boldsymbol{0}$ とおくと、$|l_i|\le1$ で $$ P_1A=\begin{pmatrix}1&\boldsymbol{0}^{\mathsf{T}}\\\boldsymbol{l}&I\end{pmatrix}\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{0}&S\end{pmatrix},\qquad S=B-\boldsymbol{l}\boldsymbol{w}^{\mathsf{T}} $$ です（右辺の左下は $\alpha\boldsymbol{l}=\boldsymbol{v}$、右下は $\boldsymbol{l}\boldsymbol{w}^{\mathsf{T}}+S=B$）。帰納法の仮定から $P'S=L'U'$, $|l'_{ij}|\le1$ です。$P=\operatorname{diag}(1,P')P_1$ とすると $$ PA=\begin{pmatrix}1&\boldsymbol{0}^{\mathsf{T}}\\P'\boldsymbol{l}&P'\end{pmatrix}\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{0}&S\end{pmatrix}=\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\alpha P'\boldsymbol{l}&P'\boldsymbol{l}\boldsymbol{w}^{\mathsf{T}}+P'S\end{pmatrix}=\begin{pmatrix}1&\boldsymbol{0}^{\mathsf{T}}\\P'\boldsymbol{l}&L'\end{pmatrix}\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{0}&U'\end{pmatrix} $$ です（最後は $P'S=L'U'$ を使いました）。右辺の左の因子は単位下三角で、成分は $P'\boldsymbol{l}$（$\boldsymbol{l}$ の並べ替え）と $L'$ の成分なので絶対値は $1$ 以下です。 証明の手順がそのまま算法です。$k$ 段目では、残っている右下のブロック $A^{(k)}$（$A^{(1)}=A$）の第 1 列から絶対値最大の行を選んで入れ替え、乗数 $l_{ik}=a^{(k)}_{ik}/a^{(k)}_{kk}$ を作り、 $$ a^{(k+1)}_{ij}=a^{(k)}_{ij}-l_{ik}a^{(k)}_{kj}\qquad(i,j\gt k) $$ で次のブロックを作ります。$U$ の第 $k$ 行は $a^{(k)}_{kj}$（$j\ge k$）です。行の入れ替えは浮動小数点数でも厳密なので、丸め誤差の解析では、**最終的な $P$ を先に $A$ に施した行列をピボット選択なしで消去した**と考えてかまいません。 上の算法（ピボット選択の比較を除く）の浮動小数点演算の回数は $\tfrac23n^3-\tfrac12n^2-\tfrac16n$ である。分解が済めば、$L\boldsymbol{y}=P\boldsymbol{b}$ と $U\boldsymbol{x}=\boldsymbol{y}$ は合わせて $2n^2-n$ 回で解ける。 $k$ 段目（$m=n-k$ とおく）は、乗数の割り算が $m$ 回、$a^{(k+1)}_{ij}$ の更新が $m^2$ 個でそれぞれ乗算と減算の 2 回です。$m=1,\ldots,n-1$ について足すと $$ \sum_{m=1}^{n-1}(2m^2+m)=\frac{(n-1)n(2n-1)}{3}+\frac{n(n-1)}{2}=\frac23n^3-\frac12n^2-\frac16n $$ です。前進代入（対角は $1$ なので割り算なし）は第 $i$ 行で乗算と減算が $i-1$ 回ずつで計 $n^2-n$ 回、後退代入は第 $i$ 行でさらに割り算 1 回が加わって計 $n^2$ 回です。 $A=\begin{pmatrix}\varepsilon&1\\1&1\end{pmatrix}$, $\boldsymbol{b}=(1,2)^{\mathsf{T}}$ の解は $x_1=1/(1-\varepsilon)$, $x_2=(1-2\varepsilon)/(1-\varepsilon)$ で、$\varepsilon\to0$ ではどちらも $1$ に近く、$\kappa_\infty(A)\approx4$ と条件のよい問題です。binary64 で $\varepsilon=10^{-20}$ としてピボット選択なしで消去すると、$l_{21}=10^{20}$, $u_{22}=\mathrm{fl}(1-10^{20})=-10^{20}$ で、$1$ の情報が丸めで消えます。後退代入の結果は $\hat x_2=1$, $\hat x_1=\mathrm{fl}((1-1)/\varepsilon)=0$ で、$x_1$ がまったく違います。部分ピボット選択なら 2 行を入れ替えて $l_{21}=\varepsilon$ となり、`numpy.linalg.solve` は $(1,1)$ を返します。$\hat L\hat U$ の大きさを見ると、ピボット選択なしでは $|\hat L||\hat U|$ の成分が $10^{20}$ に達し、次節の後退誤差の評価 $\gamma|\hat L||\hat U|$ が $|A|$ の尺度ではまったく小さくないことがわかります。 $A$ が三重対角（$|i-j|\ge2$ で $a_{ij}=0$）で、ピボット選択なしで消去できるとします。$k$ 段目で第 $k$ 列の非零は $a^{(k)}_{k+1,k}$ だけで、第 $k$ 行の非零は $a^{(k)}_{k,k+1}$ だけなので、更新されるのは $a^{(k+1)}_{k+1,k+1}$ の 1 成分だけです。したがって $L$ は対角と 1 本下の副対角、$U$ は対角と 1 本上の副対角だけをもち、分解は割り算・乗算・減算 1 回ずつの $3(n-1)$ 回で済みます。$\tfrac23n^3$ と比べると、比はおよそ $\tfrac29n^2$ で、$n=10^6$ では約 $2\times10^{11}$ 倍の差です。一般の帯行列で帯幅が保たれることは、Cholesky 分解と帯行列の章で定理として扱います。 ## 3. 後退誤差解析 ### 3.1. 算法と記号 以下、行列は浮動小数点数を成分とし、途中でオーバーフローもアンダーフローも起きず、$3nu\lt1$ とします。行列やベクトルの絶対値 $|A|$ と不等号 $|E|\le F$ は成分ごとです。丸めの記号 $\mathrm{fl}$、$\gamma_k=ku/(1-ku)$ は[浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)の章のとおりです。 ピボット選択なしの消去を、浮動小数点数では次の順序で行うとします（行の入れ替えは厳密なので、部分ピボット選択の場合はあらかじめ入れ替えた $PA$ に適用したと考えます）。$\hat a^{(1)}_{ij}=a_{ij}$ とし、$k=1,\ldots,n-1$ について $$ \hat l_{ik}=\mathrm{fl}\bigl(\hat a^{(k)}_{ik}/\hat a^{(k)}_{kk}\bigr),\qquad\hat a^{(k+1)}_{ij}=\mathrm{fl}\bigl(\hat a^{(k)}_{ij}-\mathrm{fl}(\hat l_{ik}\hat a^{(k)}_{kj})\bigr)\qquad(i,j\gt k) $$ を計算し、$\hat u_{kj}=\hat a^{(k)}_{kj}$（$j\ge k$）とします。成分 $(i,j)$ に注目すると、$s_0=a_{ij}$ から始めて $s_m=\mathrm{fl}(s_{m-1}-\mathrm{fl}(\hat l_{im}\hat u_{mj}))$ を $m=1,2,\ldots$ と更新しており、これは既存の「1 ステップの丸め誤差」の補題（）の形そのものです。つまり、消去の 1 成分の計算は「引き算と掛け算を何回か繰り返し、最後に割り算をする」だけの計算で、その丸め誤差は 1 つの補題でまとめて扱えます。下の引用は、その補題の主張だけを再掲したものです。 > **1 ステップの丸め誤差**（ の 1）: $c,a_k,b_k,d\in\mathbb{F}$、$s_0=c$, $s_k=\mathrm{fl}(s_{k-1}-\mathrm{fl}(a_kb_k))$（$k=1,\ldots,m$）、$y=\mathrm{fl}(s_m/d)$ なら、$c=dy(1+\theta^{(0)})+\sum_{k=1}^ma_kb_k(1+\theta^{(k)})$, $|\theta^{(k)}|\le\gamma_{m+1}$ と書ける。 ### 3.2. 三角方程式 $T\in\mathbb{F}^{n\times n}$ を対角成分が $0$ でない下三角行列、$\boldsymbol{b}\in\mathbb{F}^n$ とし、前進代入 $$ \hat x_i=\mathrm{fl}\Bigl(\bigl(b_i-\textstyle\sum_{k=1}^{i-1}t_{ik}\hat x_k\bigr)/t_{ii}\Bigr) $$ （和は $s\leftarrow\mathrm{fl}(s-\mathrm{fl}(t_{ik}\hat x_k))$ で $k=1,\ldots,i-1$ の順に計算する）で $\hat{\boldsymbol{x}}$ を求める。このとき $$ (T+\Delta T)\hat{\boldsymbol{x}}=\boldsymbol{b},\qquad|\Delta T|\le\gamma_n|T| $$ となる下三角行列 $\Delta T$ がある。上三角行列の後退代入（$k$ を任意の順に計算してよい）についても同じことが成り立つ。 第 $i$ 行に「1 ステップの丸め誤差」の補題（ の 1）を $c=b_i$, $m=i-1$, $a_kb_k=t_{ik}\hat x_k$, $d=t_{ii}$, $y=\hat x_i$ で使うと $$ b_i=t_{ii}\hat x_i(1+\theta^{(0)}_i)+\sum_{k=1}^{i-1}t_{ik}\hat x_k(1+\theta^{(k)}_i),\qquad|\theta^{(k)}_i|\le\gamma_i\le\gamma_n $$ です（$\gamma_k$ は $k$ について増加。 の 1）。$\Delta t_{ii}=t_{ii}\theta^{(0)}_i$, $\Delta t_{ik}=t_{ik}\theta^{(k)}_i$（$k\lt i$）とおけば主張です。後退代入では和の項の番号を付け替えれば同じ議論です。 ### 3.3. LU 分解と Gauss 消去 この節の流れは次のとおりです。まず分解 $\hat L\hat U$ が元の $A$ から成分ごとにどれだけずれるかを示し、次に前進代入・後退代入の誤差（3.2 節）を重ねて、最終的な解が「$A$ を少し変えた方程式の正確な解」になっていることを示します。どちらも、各成分に入る丸めの回数を数えて $\gamma_k$ で上から抑えるだけです。 $A\in\mathbb{F}^{n\times n}$ の 3.1 節の消去がピボット $\hat a^{(k)}_{kk}\ne0$ のまま最後まで進み、$\hat L$（単位下三角、$(i,k)$ 成分 $\hat l_{ik}$）と $\hat U$ を返したとする。このとき $$ \hat L\hat U=A+\Delta A,\qquad|\Delta A|\le\gamma_n|\hat L||\hat U|. $$ $(\hat L\hat U)_{ij}=\sum_{k\le\min(i,j)}\hat l_{ik}\hat u_{kj}$（$\hat l_{ii}=1$）です。 $i\le j$ のとき、$\hat u_{ij}=\hat a^{(i)}_{ij}$ は $s_0=a_{ij}$ から $s_m=\mathrm{fl}(s_{m-1}-\mathrm{fl}(\hat l_{im}\hat u_{mj}))$（$m=1,\ldots,i-1$）で得た $s_{i-1}$ です。これは $d=1$, $y=\mathrm{fl}(s_{i-1}/1)=s_{i-1}$ とした補題（ の 1）の形なので $$ a_{ij}=\hat u_{ij}(1+\theta^{(0)})+\sum_{k=1}^{i-1}\hat l_{ik}\hat u_{kj}(1+\theta^{(k)}),\qquad|\theta^{(k)}|\le\gamma_i\le\gamma_n. $$ $i\gt j$ のとき、$\hat l_{ij}=\mathrm{fl}(s_{j-1}/\hat u_{jj})$ で、$s_{j-1}$ は同じ更新を $m=1,\ldots,j-1$ まで行った値です。補題を $d=\hat u_{jj}$, $y=\hat l_{ij}$ で使うと $$ a_{ij}=\hat l_{ij}\hat u_{jj}(1+\theta^{(0)})+\sum_{k=1}^{j-1}\hat l_{ik}\hat u_{kj}(1+\theta^{(k)}),\qquad|\theta^{(k)}|\le\gamma_j\le\gamma_n. $$ どちらの場合も $a_{ij}=\sum_{k\le\min(i,j)}\hat l_{ik}\hat u_{kj}(1+\theta^{(k)})$ なので、$\Delta a_{ij}=(\hat L\hat U)_{ij}-a_{ij}=-\sum_k\hat l_{ik}\hat u_{kj}\theta^{(k)}$ で、$|\Delta a_{ij}|\le\gamma_n\sum_k|\hat l_{ik}||\hat u_{kj}|=\gamma_n(|\hat L||\hat U|)_{ij}$ です。 の $\hat L$, $\hat U$ を使い、前進代入 $\hat L\hat{\boldsymbol{y}}=\boldsymbol{b}$ と後退代入 $\hat U\hat{\boldsymbol{x}}=\hat{\boldsymbol{y}}$ を の算法で解いたとする。このとき $$ (A+\Delta A)\hat{\boldsymbol{x}}=\boldsymbol{b},\qquad|\Delta A|\le\gamma_{3n}|\hat L||\hat U|. $$ 部分ピボット選択の場合は、$\hat L\hat U$ が $PA$ の分解で $\hat L\hat{\boldsymbol{y}}=P\boldsymbol{b}$ を解くとして、$(A+\Delta A)\hat{\boldsymbol{x}}=\boldsymbol{b}$, $|\Delta A|\le\gamma_{3n}P^{\mathsf{T}}|\hat L||\hat U|$ である。 から $(\hat L+\Delta L)\hat{\boldsymbol{y}}=\boldsymbol{b}$, $|\Delta L|\le\gamma_n|\hat L|$、$(\hat U+\Delta U)\hat{\boldsymbol{x}}=\hat{\boldsymbol{y}}$, $|\Delta U|\le\gamma_n|\hat U|$ です。代入して の $\hat L\hat U=A+\Delta A_1$ を使うと $$ \boldsymbol{b}=(\hat L+\Delta L)(\hat U+\Delta U)\hat{\boldsymbol{x}}=(A+\Delta A_1+\hat L\Delta U+\Delta L\hat U+\Delta L\Delta U)\hat{\boldsymbol{x}} $$ なので、$\Delta A=\Delta A_1+\hat L\Delta U+\Delta L\hat U+\Delta L\Delta U$ とおきます。成分ごとの三角不等式から $$ |\Delta A|\le(\gamma_n+2\gamma_n+\gamma_n^2)|\hat L||\hat U| $$ です。 の 2 から $2\gamma_n+\gamma_n^2\le\gamma_{2n}$、同じく $\gamma_{2n}+\gamma_n\le\gamma_{3n}$ なので、係数は $\gamma_{3n}$ 以下です。部分ピボット選択の場合は $PA$, $P\boldsymbol{b}$ に適用して $(PA+\Delta A')\hat{\boldsymbol{x}}=P\boldsymbol{b}$ を得、左から $P^{\mathsf{T}}$ を掛けて $\Delta A=P^{\mathsf{T}}\Delta A'$ とします（置換は成分の絶対値を並べ替えるだけです）。 この定理は**成分ごと**の後退誤差の評価ですが、尺度は $|A|$ ではなく $|\hat L||\hat U|$ です。$|\hat L||\hat U|$ が $|A|$ と同じ程度の大きさなら、Gauss 消去は（成分ごと・ノルムともに）後退安定です。しかし のように $|\hat L||\hat U|\gg|A|$ となることがあり、そのときは評価が役に立ちません。ピボット選択の目的は、$|\hat L||\hat U|$ を $|A|$ の尺度に抑えることです。 なお LAPACK の `getrf` はキャッシュを活かすためにブロック化した順序で計算し、内積の順序も上の算法と異なりますが、各成分が受ける丸めの回数は同じ程度で、同じ形の評価（定数が少し変わる）が成り立つことが知られています（Higham『Accuracy and Stability of Numerical Algorithms』第 9 章・第 13 章）。以下のコードの数値は LAPACK の結果なので、定理の評価との比較は「同じ大きさか」を見るものです。 ## 4. 成長因子 ### 4.1. 部分ピボット選択の成長因子 3.1 節の消去で現れるすべての中間値の絶対値の最大と、$A$ の成分の絶対値の最大の比 $$ \hat\rho_n(A)=\frac{\max_{i,j,k}|\hat a^{(k)}_{ij}|}{\max_{i,j}|a_{ij}|} $$ を（計算した消去の）**成長因子**という。厳密な算術で行った消去の中間値 $a^{(k)}_{ij}$ で同じ比をとったものを $\rho_n(A)$ と書く。 $\hat u_{kj}=\hat a^{(k)}_{kj}$ なので $\max|\hat u_{kj}|\le\hat\rho_n\max|a_{ij}|$ です。部分ピボット選択では $|\hat l_{ik}|\le1$ なので、$|\hat L||\hat U|$ の大きさは成長因子で決まります。次の定理は、その成長因子が最悪でもどこまで大きくなりうるかを与えます。証明の方針は、各段で成分の絶対値が高々 2 倍にしかならないことを示し、それを段数だけ繰り返すことです。 部分ピボット選択を行う消去について、 1. 厳密な算術では $\rho_n(A)\le2^{n-1}$。 2. 浮動小数点数の消去（3.1 節の算法を、各段で $|\hat a^{(k)}_{ik}|$ が最大の行を選んで行う）では $|\hat l_{ik}|\le1$ であり、$\hat\rho_n(A)\le\bigl(2(1+u)^2\bigr)^{n-1}$。 1. $k$ 段目のピボットは第 $k$ 列で絶対値最大なので $|l_{ik}|=|a^{(k)}_{ik}/a^{(k)}_{kk}|\le1$ です（ピボットが $0$ なら列は $0$ で、乗数を $0$ とします）。よって $|a^{(k+1)}_{ij}|\le|a^{(k)}_{ij}|+|l_{ik}||a^{(k)}_{kj}|\le2\max_{p,q}|a^{(k)}_{pq}|$ で、$k$ についての帰納法から $\max|a^{(k)}_{pq}|\le2^{k-1}\max|a_{pq}|$ です。$k\le n$ なので主張が出ます。 2. $t=\hat a^{(k)}_{ik}/\hat a^{(k)}_{kk}$ は $|t|\le1$ です。$\hat l_{ik}=\mathrm{fl}(t)$ が $1$ より大きいとすると、$1\in\mathbb{F}$ は $t$ との距離 $1-t$ が $\mathrm{fl}(t)-t$ より小さく、最も近い浮動小数点数を選ぶ丸めの定義に反します。$-1$ についても同様で、$|\hat l_{ik}|\le1$ です。標準モデル（）から $$ |\hat a^{(k+1)}_{ij}|\le(1+u)\bigl(|\hat a^{(k)}_{ij}|+(1+u)|\hat l_{ik}||\hat a^{(k)}_{kj}|\bigr)\le2(1+u)^2\max_{p,q}|\hat a^{(k)}_{pq}| $$ で、1 と同じ帰納法です。 部分ピボット選択つき Gauss 消去で計算した解 $\hat{\boldsymbol{x}}$ は $$ (A+\Delta A)\hat{\boldsymbol{x}}=\boldsymbol{b},\qquad\|\Delta A\|_\infty\le n^2\gamma_{3n}\,\hat\rho_n(A)\,\|A\|_\infty $$ を満たす。したがってノルムに関する後退誤差（、$\infty$ ノルム）は $n^2\gamma_{3n}\hat\rho_n(A)$ 以下である。 から $|\Delta A|\le\gamma_{3n}P^{\mathsf{T}}|\hat L||\hat U|$ です。$\|\cdot\|_\infty$ は行和の最大（ の 2）なので、成分ごとの不等式 $|E|\le F$ から $\|E\|_\infty\le\|F\|_\infty$ が従い、行の並べ替え $P^{\mathsf{T}}$ はノルムを変えません。また誘導ノルムは $\|XY\boldsymbol{z}\|\le\|X\|\|Y\boldsymbol{z}\|\le\|X\|\|Y\|\|\boldsymbol{z}\|$ から劣乗法的です。よって $$ \|\Delta A\|_\infty\le\gamma_{3n}\,\||\hat L|\|_\infty\,\||\hat U|\|_\infty $$ です。$|\hat l_{ik}|\le1$（ の 2）から $\||\hat L|\|_\infty\le n$、$|\hat u_{kj}|\le\hat\rho_n\max|a_{ij}|$ から $\||\hat U|\|_\infty\le n\hat\rho_n\max|a_{ij}|$、そして $\max|a_{ij}|\le\|A\|_\infty$ です。後半は $\Delta\boldsymbol{b}=\boldsymbol{0}$ として定義に当てはめたものです。 したがって条件数との積で、前進誤差は $\|\hat{\boldsymbol{x}}-\boldsymbol{x}\|_\infty/\|\boldsymbol{x}\|_\infty\lesssim2\kappa_\infty(A)\,n^2\gamma_{3n}\hat\rho_n$ です（）。$n^2$ の係数は最悪の場合の数え方で、実際の誤差はずっと小さいのが普通です。問題は $\hat\rho_n$ です。 対角成分が $1$、対角より下がすべて $-1$、最後の列がすべて $1$、他が $0$ の行列 $W_n$ を考えます。第 1 列は $1,-1,\ldots,-1$ で絶対値がすべて等しく、部分ピボット選択は行を入れ替えません（同点なら上の行を選ぶ約束）。1 段目の消去で最後の列は $1,2,\ldots,2$ になり、残りのブロックは同じ形で最後の列が 2 倍になったものです。これを繰り返すと $u_{nn}=2^{n-1}$ となり、$\rho_n(W_n)=2^{n-1}$ で の上界に達します。一方 $\kappa_\infty(W_n)=n$ で、問題の条件はよいのです。 次のコードは $n=10,20,30,50$ について、`scipy.linalg.lu` の $U$ から測った成長因子（$=2^{n-1}$ かの判定とその値）、$P|\hat L||\hat U||\hat{\boldsymbol{x}}|$ を尺度とした成分ごとの後退誤差 $\omega=\max_i|r_i|/(P|\hat L||\hat U||\hat{\boldsymbol{x}}|)_i$ と評価 $\gamma_{3n}$、Rigal–Gaches の後退誤差 $\eta=\|\boldsymbol{r}\|_\infty/(\|A\|_\infty\|\hat{\boldsymbol{x}}\|_\infty+\|\boldsymbol{b}\|_\infty)$（）、$\kappa_\infty$、mpmath の 50 桁の解に対する前進誤差を並べます（SciPy の `lu` は $A=PLU$ の形で返すので、尺度の $P$ は本文の $P^{\mathsf{T}}$ にあたります）。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 / mpmath 1.3.0 # （x86-64, OpenBLAS）／乱数種 20261006 import numpy as np import scipy.linalg as sla import mpmath as mp def wilkinson(n): # 対角 1, 対角より下 -1, 最後の列 1 A = np.eye(n) - np.tril(np.ones((n, n)), -1) A[:, -1] = 1.0 return A mp.mp.dps = 50 rng = np.random.default_rng(20261006) g = lambda k: k * 2.0**-53 / (1 - k * 2.0**-53) for n in (10, 20, 30, 50): A = wilkinson(n) P, L, U = sla.lu(A) # A = P L U rho = np.abs(U).max() / np.abs(A).max() b = rng.standard_normal(n) x = sla.lu_solve(sla.lu_factor(A), b) r = b - A @ x om = np.max(np.abs(r) / (P @ np.abs(L) @ np.abs(U) @ np.abs(x))) eta = np.abs(r).max() / (n * np.abs(x).max() + np.abs(b).max()) kap = np.linalg.cond(A, np.inf) xm = mp.lu_solve(mp.matrix(A.tolist()), mp.matrix(b.tolist())) fe = max(abs(xm[i] - x[i]) for i in range(n)) / mp.norm(xm, mp.inf) print(f"{n:>2} {rho == 2.0**(n-1)} {rho:.2e} {om:.1e} {g(3*n):.1e}" f" {eta:.2e} {kap:.0f} {float(fe):.2e}") ``` ```text 10 True 5.12e+02 6.9e-17 3.3e-15 6.04e-16 10 2.47e-15 20 True 5.24e+05 1.2e-16 6.7e-15 3.38e-14 20 4.91e-13 30 True 5.37e+08 8.4e-17 1.0e-14 8.78e-11 30 9.45e-10 50 True 5.63e+14 6.2e-17 1.7e-14 9.65e-05 50 2.59e-03 ``` **検算.** (a) 既知値: 成長因子は 4 つの $n$ すべてで $2^{n-1}$ に一致し（`True`）、$\kappa_\infty(W_n)=n$ も閉じた形のとおりです。(b) 別精度: 前進誤差は mpmath の 50 桁の解との比較で、 の評価 $2\kappa_\infty\eta/(1-\kappa_\infty\eta)$ の内側にあります（$n=50$ で $2.6\times10^{-3}\le9.7\times10^{-3}$）。(c) 次数: $|\hat L||\hat U|$ を尺度とした後退誤差 $\omega$ はどの $n$ でも $u$ 程度で の $\gamma_{3n}$ の内側ですが、$|A|$ の尺度の後退誤差 $\eta$ は、$n=10$ から $20$ では約 $56$ 倍ですが、その先は $n$ が 10 増えるごとにおよそ $2^{10}\approx10^3$ 倍に増え、成長因子の増え方 $2^{n-1}$ と同じ割合で悪化します。条件のよい問題なのに答えが 3 桁しか合わないのは、アルゴリズムの後退誤差が大きいからです。 残りのブロック全体から絶対値最大の成分をピボットに選ぶ（行と列を入れ替える）**完全ピボット選択**では、Wilkinson が Hadamard の不等式を使って $\rho_n\le n^{1/2}(2\cdot3^{1/2}\cdots n^{1/(n-1)})^{1/2}$ を示しました（Wilkinson 1961、Higham 第 2 版 §9.4）。列だけでなく行の最大でもあるピボットを探す**ルークピボット選択**も、成長因子の上界が完全ピボット選択と同様に指数関数よりずっと遅く増えることが知られています（Foster 1997、Higham §9.1）。これらは比較の手間が増えるので、LAPACK の標準の解法は部分ピボット選択です。部分ピボット選択で成長因子が大きくなる行列は実際にはまれで、乱数行列の成長因子の分布は Trefethen–Schreiber（1990）が数値実験で調べています。この記事の証明はこれらの評価を使いません。 乱数種を固定した標準正規乱数の $n\times n$ 行列を各 $n$ で 30 本作り、`scipy.linalg.lu` の $U$ から測った $\max|\hat u_{kj}|/\max|a_{ij}|$（成長因子の下界）の平均と最大を出します。最後の列は最後の 1 本の $\max|\hat l_{ik}|$ です。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64, OpenBLAS） # 乱数種 20261006 import numpy as np import scipy.linalg as sla rng = np.random.default_rng(20261006) for n in (100, 200, 400, 800): rho = [] for _ in range(30): A = rng.standard_normal((n, n)) P, L, U = sla.lu(A) rho.append(np.abs(U).max() / np.abs(A).max()) rho = np.array(rho) print(f"{n:>4} 平均 {rho.mean():6.2f} 最大 {rho.max():6.2f}" f" 平均/√n {rho.mean() / n**0.5:.3f}" f" |L|最大 {np.abs(L).max():.2f}") ``` ```text 100 平均 5.51 最大 10.11 平均/√n 0.551 |L|最大 1.00 200 平均 7.57 最大 12.48 平均/√n 0.535 |L|最大 1.00 400 平均 11.21 最大 14.44 平均/√n 0.560 |L|最大 1.00 800 平均 15.65 最大 24.16 平均/√n 0.553 |L|最大 1.00 ``` **検算.** (a) 既知値: $|\hat l_{ik}|\le1$（ の 2）が成り立っています。(b) 別の物差し: 成長因子は上界 $2^{n-1}$（$n=100$ で $6\times10^{29}$）より 28 桁以上小さく、 とは対照的です。(c) 増え方: この範囲では平均は $\sqrt n$ にほぼ比例して増えています（比が $0.53$〜$0.56$）。ただしこれは 30 本の標本で見た経験的な傾向で、定理ではありません。$U$ だけから測っているので、中間値をすべて含む成長因子はこれ以上です。 ### 4.2. 列対角優位な行列 ピボット選択が要らない重要な行列の族があります。狭義列対角優位とは、各列で対角成分の絶対値が、同じ列のほかの成分の絶対値の和より大きいことです。次の定理は、この性質が消去の途中で失われず、成長因子も小さく抑えられることを示します。 $A\in M_n(\mathbb{K})$ が**狭義列対角優位**、すなわちすべての $j$ で $|a_{jj}|\gt\sum_{i\ne j}|a_{ij}|$ とする。このとき（厳密な算術で） 1. ピボット選択なしの消去は最後まで進み、$A=LU$ の LU 分解がある。乗数は $|l_{ik}|\lt1$ で、部分ピボット選択は行を入れ替えない。 2. 消去の各段の残りのブロック（Schur 補行列）も狭義列対角優位である。 3. 各列の成分の絶対値の和は消去で増えず、成長因子は $\rho_n(A)\le2$ である。 1 段目を考えます。$A=\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{v}&B\end{pmatrix}$ とすると、第 1 列の狭義列対角優位から $|\alpha|\gt\sum_i|v_i|\ge0$ なので $\alpha\ne0$、$|l_i|=|v_i|/|\alpha|\lt1$ です。Schur 補行列の分解（ の 2 を $1\times1$ の左上ブロック $\alpha$ で使う）から、$A$ は $\begin{pmatrix}1&\boldsymbol{0}^{\mathsf{T}}\\\boldsymbol{v}/\alpha&I\end{pmatrix}\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{0}&S\end{pmatrix}$, $S=B-\boldsymbol{v}\boldsymbol{w}^{\mathsf{T}}/\alpha$ と書けます。 **2.** $S$ の第 $j$ 列（$A$ の第 $j+1$ 列に対応。以下 $b_{ij}$, $w_j$, $v_i$ は対応する成分）を見ます。$A$ のこの列の狭義列対角優位は $|b_{jj}|\gt|w_j|+\sum_{i\ne j}|b_{ij}|$ で、第 1 列のそれから $\sum_{i\ne j}|v_i|\lt|\alpha|-|v_j|$ です。よって $$ \sum_{i\ne j}|s_{ij}|\le\sum_{i\ne j}|b_{ij}|+\frac{|w_j|}{|\alpha|}\sum_{i\ne j}|v_i|\lt\bigl(|b_{jj}|-|w_j|\bigr)+|w_j|\Bigl(1-\frac{|v_j|}{|\alpha|}\Bigr)=|b_{jj}|-\frac{|w_j||v_j|}{|\alpha|}\le|s_{jj}| $$ です（最後は $s_{jj}=b_{jj}-v_jw_j/\alpha$ と三角不等式）。$S$ は狭義列対角優位なので、同じ議論を $S$ に繰り返せます。各段のピボットが $0$ でないので消去は最後まで進み、 の証明と同じ組み立て（$P=I$）で $A=LU$ を得ます。 **3.** $S$ の第 $j$ 列の絶対値の和は $\sum_i|s_{ij}|\le\sum_i|b_{ij}|+|w_j|\sum_i|v_i|/|\alpha|\le\sum_i|b_{ij}|+|w_j|$ で、これは $A$ の対応する列の絶対値の和です。各段で列の絶対値の和は増えないので、$k$ 段目の中間値について $|a^{(k)}_{ij}|\le\sum_p|a_{pj}|\lt2|a_{jj}|\le2\max|a_{pq}|$ です（狭義列対角優位から非対角の和は $|a_{jj}|$ 未満）。よって $\rho_n(A)\le2$ です。 列対角優位な行列は、偏微分方程式の上流差分、Markov 連鎖の遷移行列から作る行列などに現れます。行対角優位な行列は転置が列対角優位なので、$A^{\mathsf{T}}=LU$ から $A=U^{\mathsf{T}}L^{\mathsf{T}}$ の形の分解が同様に安定に作れます。 ## 5. 反復改良 計算した解 $\boldsymbol{x}_0$ の残差 $\boldsymbol{r}_0=\boldsymbol{b}-A\boldsymbol{x}_0$ を求め、補正方程式 $A\boldsymbol{d}_0=\boldsymbol{r}_0$ を解いて $\boldsymbol{x}_1=\boldsymbol{x}_0+\boldsymbol{d}_0$ とする操作を繰り返すのが**反復改良**です。補正方程式は、手元の（誤差のある）LU 分解で解きます。その誤差を $A$ の摂動 $E_k$ として表すと、次の定理が成り立ちます。 定理に出てくる 3 つの誤差は、それぞれ別の場所で生じるものです。$\boldsymbol{g}_k$ は残差 $\boldsymbol{b}-A\boldsymbol{x}_k$ を計算するときの誤差、$E_k$ は補正方程式を（精度の低い）LU 分解で解くことによる誤差、$\boldsymbol{f}_k$ は $\boldsymbol{x}_k+\boldsymbol{d}_k$ の足し算の誤差です。主張は、補正方程式の誤差が十分小さければ、1 回ごとに解の誤差が一定の割合で縮む、というものです。 $A\in M_n(\mathbb{K})$ を正則、$A\boldsymbol{x}=\boldsymbol{b}$ とし、ベクトルノルムとその誘導ノルムを 1 つ固定する。$\boldsymbol{x}_0$ から次のように $\boldsymbol{x}_k$ を作る。 $$ \boldsymbol{r}_k=\boldsymbol{b}-A\boldsymbol{x}_k+\boldsymbol{g}_k,\qquad(A+E_k)\boldsymbol{d}_k=\boldsymbol{r}_k,\qquad\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\boldsymbol{d}_k+\boldsymbol{f}_k. $$ ここで $\boldsymbol{g}_k$ は残差の計算誤差、$E_k$ は補正方程式を解く誤差、$\boldsymbol{f}_k$ は足し算の誤差である。$\|A^{-1}E_k\|\le\theta\lt1$ がすべての $k$ で成り立つとき、$\boldsymbol{e}_k=\boldsymbol{x}_k-\boldsymbol{x}$ は $$ \|\boldsymbol{e}_{k+1}\|\le\frac{\theta}{1-\theta}\|\boldsymbol{e}_k\|+\frac{\|A^{-1}\|\,\|\boldsymbol{g}_k\|}{1-\theta}+\|\boldsymbol{f}_k\| $$ を満たす。とくに $\boldsymbol{g}_k=\boldsymbol{f}_k=\boldsymbol{0}$ で $\theta\lt1/2$ なら、誤差は 1 回ごとに $\theta/(1-\theta)\lt1$ 倍以下に縮み、$\boldsymbol{x}_k\to\boldsymbol{x}$ である。 逆行列の摂動の補題（）を、正則行列 $I$（逆行列は $I$）と摂動 $A^{-1}E_k$ に使うと、$\|A^{-1}E_k\|\le\theta\lt1$ から $I+A^{-1}E_k$ は正則で $\|(I+A^{-1}E_k)^{-1}\|\le1/(1-\theta)$ です。$A+E_k=A(I+A^{-1}E_k)$ も正則で、$M_k=(A+E_k)^{-1}=(I+A^{-1}E_k)^{-1}A^{-1}$ です。$\boldsymbol{b}-A\boldsymbol{x}_k=-A\boldsymbol{e}_k$ なので $$ \boldsymbol{e}_{k+1}=\boldsymbol{e}_k+M_k(-A\boldsymbol{e}_k+\boldsymbol{g}_k)+\boldsymbol{f}_k=(I-M_kA)\boldsymbol{e}_k+M_k\boldsymbol{g}_k+\boldsymbol{f}_k $$ です。$I-M_kA=M_k\bigl((A+E_k)-A\bigr)=M_kE_k=(I+A^{-1}E_k)^{-1}A^{-1}E_k$ なので $\|I-M_kA\|\le\theta/(1-\theta)$、また $\|M_k\|\le\|A^{-1}\|/(1-\theta)$ です。三角不等式で主張の不等式を得ます。$\theta\lt1/2$ なら $\theta/(1-\theta)\lt1$ で、$\|\boldsymbol{e}_k\|\le(\theta/(1-\theta))^k\|\boldsymbol{e}_0\|\to0$ です。 この定理の読み方は次のとおりです。 - **縮小率** $\theta/(1-\theta)\approx\|A^{-1}E_k\|$ は、補正方程式を解く精度で決まります。分解を単位丸め $u_{\mathrm{f}}$ で行えば、 と成長因子の評価から $\|E_k\|\lesssim c_n\hat\rho_nu_{\mathrm{f}}\|A\|$ で、$\theta\approx c_n\hat\rho_n\kappa(A)u_{\mathrm{f}}$ です。 - **到達できる精度**は第 2 項・第 3 項、つまり残差と足し算の精度で決まります。残差を分解より高い精度 $u_{\mathrm{r}}$ で計算すれば、分解が低い精度でも、最終的な誤差は $u_{\mathrm{r}}$ の側で決まります。 これが**混合精度の反復改良**の考え方です。$O(n^3)$ の分解を速い低精度（binary32 や GPU の binary16）で行い、$O(n^2)$ の残差と更新だけを binary64 で行います。条件は $\kappa(A)u_{\mathrm{f}}$ が $1$ より十分小さいことです。 乱数種を固定した直交行列 $Q_1,Q_2$ と、$1$ から $1/\kappa$ まで等比に並べた特異値で $A=Q_1\Sigma Q_2^{\mathsf{T}}$（$\kappa_2(A)=\kappa$）を作り、$n=50$ で binary32 の LU 分解と binary64 の残差による反復改良を 4 回行います。1 行目は $\boldsymbol{x}_0$（binary32 の分解で解いたもの）から $\boldsymbol{x}_4$ までの相対誤差（$\infty$ ノルム、参照解は mpmath の 50 桁）、2 行目は最初の 2 回の縮小率、$\kappa u_{32}$、binary64 の `solve` で直接解いたときの相対誤差です。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 / mpmath 1.3.0 # （x86-64, OpenBLAS）／乱数種 20261006 import numpy as np import scipy.linalg as sla import mpmath as mp mp.mp.dps = 50 rng = np.random.default_rng(20261006) n, u32 = 50, 2.0 ** -24 Q1, _ = np.linalg.qr(rng.standard_normal((n, n))) Q2, _ = np.linalg.qr(rng.standard_normal((n, n))) b = rng.standard_normal(n) for kap in (1e2, 1e4): A = Q1 @ np.diag(np.geomspace(1, 1 / kap, n)) @ Q2.T # κ2 = kap xs = mp.lu_solve(mp.matrix(A.tolist()), mp.matrix(b.tolist())) xs = np.array([float(t) for t in xs]) # 参照解 f32 = sla.lu_factor(A.astype(np.float32)) # binary32 の LU x = sla.lu_solve(f32, b.astype(np.float32)).astype(np.float64) errs = [] for k in range(5): errs.append(np.abs(x - xs).max() / np.abs(xs).max()) r = b - A @ x # binary64 の残差 d = sla.lu_solve(f32, r.astype(np.float32)) x = x + d.astype(np.float64) x64 = np.linalg.solve(A, b) e64 = np.abs(x64 - xs).max() / np.abs(xs).max() rate = [errs[k + 1] / errs[k] for k in range(2)] print(f"κ={kap:.0e} " + " ".join(f"{e:.2e}" for e in errs)) print(f" 比 {rate[0]:.1e} {rate[1]:.1e} κu32 {kap * u32:.1e}" f" binary64 直接 {e64:.2e}") ``` ```text κ=1e+02 2.89e-06 3.84e-12 1.02e-15 1.59e-15 1.28e-15 比 1.3e-06 2.7e-04 κu32 6.0e-06 binary64 直接 2.39e-15 κ=1e+04 5.66e-05 1.99e-09 1.12e-13 5.25e-14 4.59e-14 比 3.5e-05 5.6e-05 κu32 6.0e-04 binary64 直接 1.55e-13 ``` **検算.** (a) 既知値: 特異値を指定して作ったので $\kappa_2(A)$ は構成どおり $10^2$, $10^4$ で、参照解は mpmath の 50 桁で解いたものです。(b) 別精度・別手法: 2〜3 回の改良で、binary64 の分解で直接解いた解と同じかそれ以上の精度に達します。(c) 率: 縮小率は $\kappa=10^4$ で $3.5\times10^{-5}$, $5.6\times10^{-5}$ と、 の予言する $\theta\approx\kappa u_{32}=6\times10^{-4}$ 以下です。$\kappa=10^2$ の 2 回目の比 $2.7\times10^{-4}$ は、誤差がすでに binary64 の到達精度（第 2・第 3 項）に達して止まったためで、縮小率ではありません。 ## 6. 演習 $A=\begin{pmatrix}2&1&1\\4&3&3\\8&7&9\end{pmatrix}$ の LU 分解（ピボット選択なし）を求め、 の $u_{11}u_{22}u_{33}=\det A$ を確かめよ。また部分ピボット選択を行うと、1 段目でどの行がピボットになるか。 1 段目の乗数は $l_{21}=2$, $l_{31}=4$ で、残りは $\begin{pmatrix}1&1\\3&5\end{pmatrix}$。2 段目の乗数は $l_{32}=3$ で、$u_{33}=5-3=2$ です。よって $$ L=\begin{pmatrix}1&0&0\\2&1&0\\4&3&1\end{pmatrix},\qquad U=\begin{pmatrix}2&1&1\\0&1&1\\0&0&2\end{pmatrix}, $$ $\det A=2\cdot1\cdot2=4$ です（$\det A_1=2$, $\det A_2=2\cdot3-1\cdot4=2=u_{11}u_{22}$ も確かめられます）。部分ピボット選択では第 1 列の絶対値最大の $8$ の第 3 行がピボットになり、乗数は $1/4$, $1/2$ です。 部分ピボット選択で $3\times3$ の行列の成長因子が上界 $4$ に達する例を作れ。 の $W_3$ の消去を手で追ってもよい。 $W_3=\begin{pmatrix}1&0&1\\-1&1&1\\-1&-1&1\end{pmatrix}$ です。1 段目は同点で第 1 行がピボット、乗数は $-1,-1$ で、残りのブロックは $\begin{pmatrix}1&2\\-1&2\end{pmatrix}$。2 段目も同点で第 2 行がピボット、乗数 $-1$ で $u_{33}=2+2=4$ です。$\max|a_{ij}|=1$ なので成長因子は $4=2^{3-1}$ です。各段で乗数の絶対値が $1$、かつ足し合わせる 2 つの成分の符号がそろうことが、上界に達する条件です。 $\kappa_\infty(A)=10^6$ の行列を binary32（$u_{32}\approx6\times10^{-8}$）の LU と binary64 の残差で反復改良する。 の縮小率を $\theta\approx\kappa u_{32}$ と見積もり、初期誤差 $10^{-1}$ から $10^{-8}$ に下げるのに何回の改良が要るかを見積もれ。$\kappa=10^8$ ではどうか。 $\theta\approx10^6\cdot6\times10^{-8}=0.06$ で、$\theta/(1-\theta)\approx0.064$、1 回でおよそ $1.2$ 桁縮みます。7 桁下げるには $7/1.2\approx5.9$ から約 $6$ 回です。それ以上改良しても、誤差は定理の第 2 項（binary64 の残差の誤差 $\|A^{-1}\|\|\boldsymbol{g}_k\|$）で止まり、最悪の場合の目安は $\kappa u_{64}\approx10^6\cdot1.1\times10^{-16}\approx10^{-10}$ です。$\kappa=10^8$ では $\theta\approx6$ で $\theta\lt1$ の仮定が成り立たず、定理は収束を保証しません。実際、binary32 の分解は $A$ を $\kappa u_{32}\gt1$ の相対誤差で表しており、補正の方向がほとんど情報をもたないので、反復改良は一般に収束しません。この場合は binary64 で分解します。 ## 参考文献 - N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — 第 8 章（三角方程式）、第 9 章（LU 分解と Gauss 消去、成長因子）、第 12 章（反復改良）。 - G. H. Golub, C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins University Press, 2013 — 第 3 章（一般の線形方程式）。 - L. N. Trefethen, D. Bau III, *Numerical Linear Algebra*, SIAM, 1997 — Lecture 20–23（Gauss 消去、ピボット選択、安定性）。 - J. H. Wilkinson, "Error analysis of direct methods of matrix inversion", *Journal of the ACM* 8 (1961), 281–330. [doi:10.1145/321075.321076](https://doi.org/10.1145/321075.321076) - L. N. Trefethen, R. S. Schreiber, "Average-case stability of Gaussian elimination", *SIAM Journal on Matrix Analysis and Applications* 11 (1990), 335–360. [doi:10.1137/0611023](https://doi.org/10.1137/0611023) - L. V. Foster, "The growth factor and efficiency of Gaussian elimination with rook pivoting", *Journal of Computational and Applied Mathematics* 86 (1997), 177–194.