# Cholesky 分解と帯行列・疎行列の直接法：後退安定性・慣性・fill-in の経路定理 > 正定値対称行列の Cholesky 分解の存在と一意性、ピボット選択なしでの後退安定性と解の後退誤差を証明し、LDLᵀ 分解と慣性、帯行列で帯幅が保たれること、疎行列の格納形式と消去グラフ、fill-in の経路定理を扱う。並べ替えによる fill の違いを SciPy で確かめる。 > https://rikai.mugen-giken.com/mathematics/numerical-linear-algebra/cholesky-and-sparse-direct ## 0. この記事の要点 連立一次方程式 $A\boldsymbol{x}=\boldsymbol{b}$ の係数行列が、応用で非常によく現れる「対称で、しかも正定値と呼ばれる性質をもつ」場合には、解くための手間を半分にでき、誤差の心配もほとんどなくなります。この記事は、その解き方（Cholesky 分解）がなぜ安全なのかを示します。後半では、ほとんどの成分が零である巨大な行列を解くとき、途中で零だった成分が非零に変わってしまう現象を、グラフの上の道として予測する方法を扱います。 - 実対称行列 $A$ が[正定値](/mathematics/linear-algebra/spectral-theorem#def-positive-definite)であることと、対角が正の上三角行列 $R$ で $A=R^{\mathsf{T}}R$ と書けることは同値で、$R$ は一意です（**Cholesky 分解**）。演算量は [LU 分解](/mathematics/numerical-linear-algebra/lu-and-pivoting#def-lu)の半分の $\tfrac13n^3$ です。 - Cholesky 分解は**ピボット選択なしで後退安定**（丸め誤差の影響が、$A$ をわずかに動かしたことと同じになる）です。計算した $\hat R$ の列のノルムは $\sqrt{a_{jj}}$ 程度に抑えられて成長が起きず、$\|\hat R^{\mathsf{T}}\hat R-A\|_2\le cn^2u\|A\|_2$ の形の評価が成り立ちます。解も $(A+\Delta A)\hat{\boldsymbol{x}}=\boldsymbol{b}$, $\|\Delta A\|_2\le c'n^2u\|A\|_2$ を満たします。 - 対称行列の $LDL^{\mathsf{T}}$ 分解の $D$ の符号は慣性（正・負・零の固有値の個数）を与えます。ただし不定値の行列では、対角からピボットを選ぶだけでは安定になりません。 - 帯行列の LU・Cholesky 分解は帯幅を保ち、演算量は $O(np^2)$ です。疎行列は CSR・CSC（非零成分だけを並べて持つ格納法）で格納し、行列ベクトル積は $O(\mathrm{nnz})$ です。 - **経路定理**: 数値的な打ち消しがなければ、Cholesky 因子の $(i,j)$ 成分が非零であることと、$A$ のグラフで $i$ と $j$ を、途中の頂点がすべて $\min(i,j)$ より小さい番号の道で結べることは同値です。番号の付け方（並べ替え）で fill-in（もとは零だった成分が、分解の途中で非零になること）は桁違いに変わります。 ## 1. 動機 最小二乗法の正規方程式 $A^{\mathsf{T}}A\boldsymbol{x}=A^{\mathsf{T}}\boldsymbol{b}$、有限要素法の剛性行列、Gauss 過程の共分散行列など、応用に現れる連立一次方程式の多くは係数行列が**正定値対称**です。この場合、[Gauss 消去と LU 分解](/mathematics/numerical-linear-algebra/lu-and-pivoting)の章で見た成長因子とピボット選択の心配が消え、計算量も半分になります。それを保証するのが Cholesky 分解の理論です。 もう 1 つの主題は**疎行列**です。偏微分方程式を格子で離散化すると、$n=10^6$ 元でも 1 行の非零は数個しかありません。密行列として $n^2$ 個の数を持つことも、$\tfrac13n^3$ の演算をすることもできないので、非零だけを格納し、分解の途中で新しく生じる非零（**fill-in**）を抑えます。fill-in がどこに生じるかは、行列の数値によらずグラフだけで決まる、というのが経路定理です。 以下、行列は実行列、$\|\cdot\|_2$ はユークリッドノルムとその誘導ノルム、$|A|$ と不等号は成分ごととします。実行列に限るのは、3 節で引用する後ろ向き誤差の定理が実対称行列についての主張だからです。複素 Hermite 行列でも、転置を共役転置に替えれば同じ証明が通ります。 ## 2. Cholesky 分解 $A\in M_n(\mathbb{R})$ を対称行列とする。$A$ が正定値（）であることと、対角成分がすべて正の上三角行列 $R$ で $$ A=R^{\mathsf{T}}R $$ と書けることは同値である。このとき $R$ は一意である。さらに $A=\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{w}&B\end{pmatrix}$ と分けると、$R$ の第 1 行は $(\sqrt\alpha,\ \boldsymbol{w}^{\mathsf{T}}/\sqrt\alpha)$ で、$R$ から第 1 行と第 1 列を除いた行列は Schur 補行列 $S=B-\boldsymbol{w}\boldsymbol{w}^{\mathsf{T}}/\alpha$ の Cholesky 因子である。 **⇐.** $R$ の対角成分が正なので $\det R\gt0$ で $R$ は正則です。$\boldsymbol{x}\ne\boldsymbol{0}$ なら $R\boldsymbol{x}\ne\boldsymbol{0}$ で、$\boldsymbol{x}^{\mathsf{T}}A\boldsymbol{x}=\|R\boldsymbol{x}\|_2^2\gt0$ です。 **⇒.** $n$ についての帰納法です。$n=1$ は $A=(\alpha)$, $\alpha=\boldsymbol{e}_1^{\mathsf{T}}A\boldsymbol{e}_1\gt0$ で $R=(\sqrt\alpha)$ です。$n\ge2$ とし、上のように分けます。$\alpha=\boldsymbol{e}_1^{\mathsf{T}}A\boldsymbol{e}_1\gt0$ です。第 1 行・第 1 列を最後に移す置換行列 $Q$ で $M=Q^{\mathsf{T}}AQ=\begin{pmatrix}B&\boldsymbol{w}\\\boldsymbol{w}^{\mathsf{T}}&\alpha\end{pmatrix}$ とすると、$\boldsymbol{x}^{\mathsf{T}}M\boldsymbol{x}=(Q\boldsymbol{x})^{\mathsf{T}}A(Q\boldsymbol{x})$ なので $M$ も正定値です。右下の $1\times1$ ブロック $\alpha$ は正定値なので、Schur 補行列の定理（ の 3）により、$M$ が正定値であることから $S=B-\boldsymbol{w}\alpha^{-1}\boldsymbol{w}^{\mathsf{T}}$ は正定値です。帰納法の仮定から $S=R'^{\mathsf{T}}R'$（$R'$ は対角が正の上三角）で、 $$ R=\begin{pmatrix}\sqrt\alpha&\boldsymbol{w}^{\mathsf{T}}/\sqrt\alpha\\\boldsymbol{0}&R'\end{pmatrix} $$ とおくと $R^{\mathsf{T}}R=\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{w}&\boldsymbol{w}\boldsymbol{w}^{\mathsf{T}}/\alpha+R'^{\mathsf{T}}R'\end{pmatrix}=A$ です。 **一意性と最後の主張.** 対角が正の上三角行列 $R$ で $A=R^{\mathsf{T}}R$ なら、$(1,1)$ 成分から $r_{11}^2=\alpha$、$r_{11}\gt0$ なので $r_{11}=\sqrt\alpha$、第 1 行の比較から $r_{11}r_{1j}=a_{1j}$ で第 1 行は上の形に決まります。残りのブロックの比較から、$R$ の右下のブロック $R''$ は $R''^{\mathsf{T}}R''=B-\boldsymbol{w}\boldsymbol{w}^{\mathsf{T}}/\alpha=S$ を満たします。$n$ についての帰納法で $R''$ は $S$ の Cholesky 因子として一意で、$R$ も一意です。 最後の主張を繰り返すと、$k$ 段目の Schur 補行列を $A^{(k)}=(a^{(k)}_{ij})_{i,j\ge k}$（$A^{(1)}=A$）として $$ r_{kk}=\sqrt{a^{(k)}_{kk}},\qquad r_{kj}=a^{(k)}_{kj}/r_{kk}\ (j\gt k),\qquad a^{(k+1)}_{ij}=a^{(k)}_{ij}-\frac{a^{(k)}_{ik}a^{(k)}_{kj}}{a^{(k)}_{kk}} $$ です。$A^{(k)}$ はすべて正定値なので、平方根の中身はつねに正です。逆に、対称行列にこの計算を行って途中で $a^{(k)}_{kk}\le0$ になれば、$A$ は正定値ではありません。**Cholesky 分解は、正定値性の最も安い判定法**でもあります。 演算量は、$k$ 段目（$m=n-k$）で対称な Schur 補行列の上三角 $m(m+1)/2$ 成分を乗算・減算 2 回ずつで更新するので、主要項は $\sum_mm^2\approx\tfrac13n^3$ で、LU 分解の $\tfrac23n^3$ の半分です。対称性を使って片側だけを計算するからです。 ## 3. 後退誤差解析 浮動小数点数の Cholesky 分解の誤差は、既存の定理で評価されています。算法は、列 $j=1,\ldots,n$ の順に、$i=1,\ldots,j$ について $s=a_{ij}$ から $s\leftarrow\mathrm{fl}(s-\mathrm{fl}(\hat r_{ki}\hat r_{kj}))$ を $k=1,\ldots,i-1$ の順に行い、$i\lt j$ なら $\hat r_{ij}=\mathrm{fl}(s/\hat r_{ii})$、$i=j$ なら $s\gt0$ を確かめて $\hat r_{jj}=\mathrm{fl}(\sqrt s)$ とするものです（$s\le0$ なら失敗として止める）。オーバーフロー・アンダーフローは起きず、$(3n+1)u\lt\tfrac12$ とします。 ここで $u$ は 1 回の演算で生じる相対誤差の上限（単位丸め。binary64 では $2^{-53}$）、$\gamma_k=ku/(1-ku)$ は $k$ 回の演算で積み上がる誤差をまとめて書くための量です。以下の評価は、丸め誤差の影響が「$A$ をほんの少し動かしたもの」に置き換えられる、という形で読んでください。 > **Cholesky 分解の後ろ向き誤差**（）: 対称行列 $A\in\mathbb{F}^{n\times n}$ の分解が失敗せずに $\hat R$ を返せば、$\hat R^{\mathsf{T}}\hat R=A+\Delta A$, $|\Delta A|\le\gamma_{n+1}|\hat R|^{\mathsf{T}}|\hat R|$ で、各 $j$ について $\sum_k\hat r_{kj}^2\le a_{jj}/(1-\gamma_{n+1})$。 後半の不等式が要点です。$\hat R$ の第 $j$ 列の 2 乗和は $a_{jj}$ 程度に抑えられ、LU 分解のような成長が起こりません。これをノルムの評価に直します。 上の仮定のもとで、分解が失敗せずに $\hat R$ を返せば $$ \hat R^{\mathsf{T}}\hat R=A+\Delta A,\qquad\|\Delta A\|_2\le\frac{\gamma_{n+1}}{1-\gamma_{n+1}}\operatorname{tr}A\le\frac{n\gamma_{n+1}}{1-\gamma_{n+1}}\|A\|_2. $$ 既存の定理（）の $|\Delta A|\le\gamma_{n+1}|\hat R|^{\mathsf{T}}|\hat R|$ と、成分ごとの評価からノルムの評価への補題（）から $\|\Delta A\|_2\le\gamma_{n+1}\||\hat R|^{\mathsf{T}}|\hat R|\|_2$ です。誘導ノルムの劣乗法性と $\|X^{\mathsf{T}}\|_2=\|X\|_2$（ の 3）、$\|X\|_2\le\|X\|_F$（ の 3）から $$ \||\hat R|^{\mathsf{T}}|\hat R|\|_2\le\||\hat R|\|_2^2\le\||\hat R|\|_F^2=\sum_j\sum_k\hat r_{kj}^2\le\frac{\operatorname{tr}A}{1-\gamma_{n+1}} $$ です。最後に $a_{jj}=\boldsymbol{e}_j^{\mathsf{T}}A\boldsymbol{e}_j\le\|A\|_2$（Cauchy–Schwarz の不等式と誘導ノルムの定義）から $\operatorname{tr}A\le n\|A\|_2$ です。 $n\gamma_{n+1}/(1-\gamma_{n+1})\approx n(n+1)u$ なので、$\|\Delta A\|_2\le cn^2u\|A\|_2$ の形です。ピボット選択も成長因子もなしに、行列だけで決まる評価が得られます。 上の $\hat R$ を使い、前進代入 $\hat R^{\mathsf{T}}\hat{\boldsymbol{y}}=\boldsymbol{b}$ と後退代入 $\hat R\hat{\boldsymbol{x}}=\hat{\boldsymbol{y}}$ を三角方程式の算法（）で解くと、 $$ (A+\Delta A)\hat{\boldsymbol{x}}=\boldsymbol{b},\qquad|\Delta A|\le\gamma_{3n+1}|\hat R|^{\mathsf{T}}|\hat R|,\qquad\|\Delta A\|_2\le\frac{n\gamma_{3n+1}}{1-\gamma_{n+1}}\|A\|_2 $$ が成り立つ。 三角方程式の後退誤差（）から $(\hat R^{\mathsf{T}}+\Delta_1)\hat{\boldsymbol{y}}=\boldsymbol{b}$, $|\Delta_1|\le\gamma_n|\hat R|^{\mathsf{T}}$、$(\hat R+\Delta_2)\hat{\boldsymbol{x}}=\hat{\boldsymbol{y}}$, $|\Delta_2|\le\gamma_n|\hat R|$ です。代入すると $$ \boldsymbol{b}=(\hat R^{\mathsf{T}}+\Delta_1)(\hat R+\Delta_2)\hat{\boldsymbol{x}}=(A+\Delta A_0+\hat R^{\mathsf{T}}\Delta_2+\Delta_1\hat R+\Delta_1\Delta_2)\hat{\boldsymbol{x}} $$ （$\hat R^{\mathsf{T}}\hat R=A+\Delta A_0$）で、成分ごとに $|\Delta A|\le(\gamma_{n+1}+2\gamma_n+\gamma_n^2)|\hat R|^{\mathsf{T}}|\hat R|$ です。$\gamma$ の計算規則（ の 2）から $2\gamma_n+\gamma_n^2\le\gamma_{2n}$、$\gamma_{n+1}+\gamma_{2n}\le\gamma_{3n+1}$ です。ノルムの評価は の証明と同じ計算です。 したがって[条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-condition-number)との積で、前進誤差はおよそ $cn^2u\,\kappa_2(A)$ 以下です（）。正定値の系では Cholesky 分解がもっとも速く、しかも安定な直接法です。 乱数種を固定した直交行列 $Q$ と $1$ から $10^{-6}$ まで等比に並べた固有値で正定値行列 $A=Q\Lambda Q^{\mathsf{T}}$ を作り、`numpy.linalg.cholesky`（下三角 $L=\hat R^{\mathsf{T}}$ を返す）の残差 $\|LL^{\mathsf{T}}-A\|_2/\|A\|_2$、$nu$、 の評価 $\gamma_{n+1}\operatorname{tr}A/((1-\gamma_{n+1})\|A\|_2)$ を並べます。後半は Hilbert 行列 $H_n=(1/(i+j-1))$ を binary64 に丸めた行列の分解が成功するかと、$H_n$ の最小固有値（mpmath の 50 桁）です。BLAS のスレッド数で和の順序が変わり、表の残差は先頭の桁から変わる（たとえば $n=400$ が $2.6\times10^{-16}$ から $1.8\times10^{-16}$ になる）ので、1 スレッドに固定しています。どちらでも評価の内側にあることは変わりません。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 / mpmath 1.3.0 # （x86-64, OpenBLAS, 1 スレッド）／乱数種 20261006 import os os.environ["OPENBLAS_NUM_THREADS"] = "1" # 和の順序を固定する import numpy as np import scipy.linalg as sla import mpmath as mp u = 2.0 ** -53 g = lambda k: k * u / (1 - k * u) rng = np.random.default_rng(20261006) for n in (50, 100, 200, 400): Q, _ = np.linalg.qr(rng.standard_normal((n, n))) A = Q @ np.diag(np.geomspace(1, 1e-6, n)) @ Q.T # 正定値, κ2=1e6 A = (A + A.T) / 2 L = np.linalg.cholesky(A) # A = L L^T res = np.linalg.norm(L @ L.T - A, 2) / np.linalg.norm(A, 2) bnd = g(n + 1) / (1 - g(n + 1)) * np.trace(A) / np.linalg.norm(A, 2) print(f"{n:>3} {res:.1e} {n * u:.1e} {bnd:.1e}") mp.mp.dps = 50 for n in range(10, 16): try: np.linalg.cholesky(sla.hilbert(n)) ok = "成功" except np.linalg.LinAlgError: ok = "失敗" lmin = min(mp.eigsy(mp.hilbert(n), eigvals_only=True)) print(n, ok, f"λmin={float(lmin):.1e}") ``` ```text 50 9.2e-17 5.6e-15 2.3e-14 100 1.2e-16 1.1e-14 8.6e-14 200 1.5e-16 2.2e-14 3.3e-13 400 2.6e-16 4.4e-14 1.3e-12 10 成功 λmin=1.1e-13 11 成功 λmin=3.4e-15 12 成功 λmin=1.0e-16 13 成功 λmin=3.2e-18 14 失敗 λmin=9.9e-20 15 失敗 λmin=3.0e-21 ``` **検算.** (a) 既知値: 構成に使った固有値から $\operatorname{tr}A$ と $\|A\|_2=1$ は既知で、評価の列はそれで計算した値です。Hilbert 行列の最小固有値は有理数の行列を 50 桁で対角化したものです。(b) 別の物差し: 残差はどの $n$ でも評価の $10^{-2}$〜$10^{-4}$ 倍、$nu$ の $10^{-2}$ 倍程度で、 の内側にあります。(c) 増え方: 評価は $n^2$ に比例して（$n$ を 2 倍するごとに約 4 倍）増えますが、実際の残差は $n$ を 8 倍しても 3 倍程度しか増えません。最悪評価と典型的な誤差の差です。Hilbert 行列では $\lambda_{\min}$ が $u\|H_n\|\approx2\times10^{-16}$ を下回った後もしばらく分解は成功し（丸めた行列 $\mathrm{fl}(H_n)$ の分解で、成功しても $H_n$ の正定値性の証明にはなりません）、$n=14$ で失敗します。[精度保証つきの線形代数](/mathematics/computer-assisted-proofs/certified-linear-algebra)の章の Rump の検証が $n=11$ で失敗するのは、丸め誤差の分の余裕 $c\gt\beta$ を固有値に要求するからで、こちらの「成功」とは意味が違います。 ## 4. LDLᵀ 分解と慣性 正定値でない対称行列にも、平方根を使わない対称な分解があります。慣性が分かれば、たとえば行列が正定値かどうか（固有値がすべて正か）も分かります。次の命題は、その個数が分解の対角成分 $D$ の符号を数えるだけで読み取れることを言います。 $A\in M_n(\mathbb{R})$ を対称行列とし、先頭主小行列式がすべて $0$ でないとする。このとき単位下三角行列 $L$ と対角行列 $D=\operatorname{diag}(d_1,\ldots,d_n)$ で $A=LDL^{\mathsf{T}}$ と一意に書け、$d_k=\det A_k/\det A_{k-1}$（$\det A_0=1$）である。$A$ の慣性（正・負・零の固有値の個数）は $D$ の正・負の対角成分の個数に等しく、零固有値はない。 LU 分解の存在と一意性（）から $A=LU$ で、$U$ の対角は $u_{kk}=\det A_k/\det A_{k-1}\ne0$ です。$D=\operatorname{diag}(u_{11},\ldots,u_{nn})$、$M=D^{-1}U$（対角が $1$ の上三角）とおくと $A=LDM$ です。$A=A^{\mathsf{T}}=M^{\mathsf{T}}(DL^{\mathsf{T}})$ は、単位下三角 $M^{\mathsf{T}}$ と対角が $0$ でない上三角 $DL^{\mathsf{T}}$ の積なので、LU 分解の一意性から $M^{\mathsf{T}}=L$ で、$A=LDL^{\mathsf{T}}$ です。分解が $A=L'D'L'^{\mathsf{T}}$ と別にあれば、$D'L'^{\mathsf{T}}$ は対角が $d'_1,\ldots,d'_n$ の上三角で、 の 2 ⇒ 1 の計算から $d'_1\cdots d'_k=\det A_k\ne0$ なので $d'_k\ne0$ です。よって $L'(D'L'^{\mathsf{T}})$ も LU 分解なので $L'=L$, $D'L'^{\mathsf{T}}=DL^{\mathsf{T}}$ から $D'=D$ です。慣性の主張は、$LDL^{\mathsf{T}}$ 分解による慣性の認証（。Sylvester の慣性法則 による）そのもので、$d_k\ne0$ なので零固有値はありません。 正定値なら $d_k\gt0$ で、$R=D^{1/2}L^{\mathsf{T}}$ が Cholesky 因子です。ところが不定値のときは、$d_k$ が小さくなって成長が起こりえます。 $A=\begin{pmatrix}\varepsilon&1\\1&\varepsilon\end{pmatrix}$（$0\lt\varepsilon\ll1$）は固有値が $\varepsilon\pm1$ の不定値行列で、$\kappa_2(A)=(1+\varepsilon)/(1-\varepsilon)\approx1$ と条件がよいです。対称性を保つには行と列を同時に入れ替える（対角成分どうしを入れ替える）しかありませんが、対角は両方 $\varepsilon$ なので、どちらを選んでも $d_1=\varepsilon$, $l_{21}=1/\varepsilon$, $d_2=\varepsilon-1/\varepsilon$ となり、$|L||D||L^{\mathsf{T}}|$ の成分は $1/\varepsilon$ に達します。$A=\begin{pmatrix}0&1\\1&0\end{pmatrix}$ では $LDL^{\mathsf{T}}$ 分解そのものが存在しません。 不定値の対称行列では、$D$ に $1\times1$ と $2\times2$ のブロックを許す $PAP^{\mathsf{T}}=LDL^{\mathsf{T}}$ が標準です。$2\times2$ のピボット $\begin{pmatrix}a_{kk}&a_{kr}\\a_{rk}&a_{rr}\end{pmatrix}$ を、成長が定数倍に抑えられる基準（閾値 $\alpha=(1+\sqrt{17})/8$）で選ぶ方法が Bunch–Kaufman（1977）で、LAPACK の `sytrf`（SciPy の `scipy.linalg.ldl`）が使っています。成長因子の評価と後退安定性は Higham 第 2 版 第 11 章にあります。この記事の証明はこれを使いません。 ## 5. 帯行列 $a_{ij}=0$（$i\gt j+p$ または $j\gt i+q$）のとき、$A$ は**下帯幅** $p$、**上帯幅** $q$ をもつといいます。三重対角行列は $p=q=1$ です。 $A\in M_n(\mathbb{R})$ が下帯幅 $p$・上帯幅 $q$ をもち、LU 分解 $A=LU$（）をもつとする。 1. $L$ は下帯幅 $p$、$U$ は上帯幅 $q$ をもつ。消去の演算量は $n(2pq+p)$ 以下である。 2. $A$ が正定値対称で $p=q$ なら、Cholesky 因子 $R$ は上帯幅 $p$ をもち、分解の演算量は $n(p+1)(2p+1)$ 以下である。 1. 先頭主小行列式が $0$ でないので、ピボット選択なしの消去は各段のピボット $u_{kk}=\det A_k/\det A_{k-1}\ne0$ で進み、 の証明の組み立て（$P=I$）で LU 分解を作ります。LU 分解は一意（）なので、これが $A=LU$ です。1 段目で $A=\begin{pmatrix}\alpha&\boldsymbol{w}^{\mathsf{T}}\\\boldsymbol{v}&B\end{pmatrix}$ と分けると、帯の仮定から $\boldsymbol{v}$ の非零は最初の $p$ 成分、$\boldsymbol{w}$ の非零は最初の $q$ 成分に限られます。$L$ の第 1 列 $(1,\boldsymbol{v}/\alpha)$ と $U$ の第 1 行 $(\alpha,\boldsymbol{w}^{\mathsf{T}})$ は帯の中にあります。Schur 補行列 $S=B-\boldsymbol{v}\boldsymbol{w}^{\mathsf{T}}/\alpha$ で $B$ から変わるのは左上の $p\times q$ ブロック（$A$ の添字で $2\le i\le p+1$, $2\le j\le q+1$）だけで、そこでは $i-j\le p-1$, $j-i\le q-1$ なので、$S$ も同じ帯幅をもちます。$S$ について同じ議論を繰り返せば、$L$ と $U$ の各列・各行が帯の中にあることがわかります。各段の演算は、割り算 $p$ 回以下と $p\times q$ 個以下の成分の更新（各 2 回）なので、$n$ 段で $n(2pq+p)$ 以下です。 2. の $R$ について、$D=\operatorname{diag}(r_{11},\ldots,r_{nn})$ とすると $A=(R^{\mathsf{T}}D^{-1})(DR)$ は LU 分解（単位下三角 × 対角が正の上三角）なので、一意性から $U=DR$ で、1 から $U$ の上帯幅は $p$、したがって $R$ もそうです。2 節の算法で $\hat r_{ij}$（$j-p\le i\le j$）を計算するとき、内積の項は $\max(1,j-p)\le k\lt i$ の $p$ 個以下なので、1 成分あたり $2p+1$ 回以下、1 列あたり $p+1$ 成分で、全体で $n(p+1)(2p+1)$ 以下です。 $p$ が $n$ に比べて小さければ、$O(n^3)$ が $O(np^2)$ に下がります。1 次元の境界値問題の差分行列（三重対角）は $O(n)$ で解けます。ただし部分ピボット選択を行うと、行の入れ替えで $U$ の上帯幅は $p+q$ まで広がりえます（$L$ の各列の非零は $p$ 個以下のまま）。正定値対称や列対角優位（）ならピボット選択は要らず、帯幅は保たれます。
→ A R
## 6. 疎行列の直接法 ### 6.1. 格納形式と行列ベクトル積 非零成分の個数を $\mathrm{nnz}(A)$ と書きます。 $A\in M_{m,n}(\mathbb{R})$ の非零成分を、次のいずれかの形で格納する（添字は $0$ 始まりとする）。 1. **COO 形式**: 長さ $\mathrm{nnz}$ の 3 つの配列 `row`, `col`, `val` で、$a_{\mathtt{row}[t],\mathtt{col}[t]}=\mathtt{val}[t]$。 2. **CSR 形式**（圧縮行格納）: 非零を行の順に並べた長さ $\mathrm{nnz}$ の配列 `val` と列番号の配列 `colind`、長さ $m+1$ の配列 `rowptr` で、第 $i$ 行の非零は `val[rowptr[i]:rowptr[i+1]]`、その列番号は `colind` の同じ範囲。 3. **CSC 形式**（圧縮列格納）: CSR の行と列の役割を入れ替えたもの（$A^{\mathsf{T}}$ の CSR と同じ配列）。 COO は組み立て（要素ごとに足し込む有限要素法など）に、CSR は行ごとの演算と行列ベクトル積に、CSC は列ごとの演算と直接法（SuperLU・CHOLMOD）に向いています。SciPy の `scipy.sparse` はこの 3 つを `coo_array`, `csr_array`, `csc_array` などで持ち、`tocsr()` などで相互に変換できます。 CSR 形式の $A\in M_{m,n}(\mathbb{R})$ と $\boldsymbol{x}\in\mathbb{R}^n$ について、$y_i=\sum_{t=\mathtt{rowptr}[i]}^{\mathtt{rowptr}[i+1]-1}\mathtt{val}[t]\,x_{\mathtt{colind}[t]}$ で $\boldsymbol{y}=A\boldsymbol{x}$ が求まり、乗算と加算はそれぞれ $\mathrm{nnz}(A)$ 回である。記憶量は $2\,\mathrm{nnz}(A)+m+1$ 語である。 第 $i$ 行の非零はちょうど `val` の `rowptr[i]` から `rowptr[i+1]-1` までなので、$y_i=\sum_ja_{ij}x_j$ の非零の項だけを足していることになり、正しい値です。各非零 $t$ について乗算 1 回と加算 1 回（$y_i=0$ から始めて足し込む）を行い、各 $t$ はちょうど 1 つの行に属するので、合計はそれぞれ $\mathrm{nnz}(A)$ 回です。記憶量は `val`, `colind` が $\mathrm{nnz}$ 語ずつ、`rowptr` が $m+1$ 語です。 反復法（Krylov 部分空間法）が大規模な疎行列に向くのは、1 反復の主な仕事がこの $O(\mathrm{nnz})$ の積だからです。直接法の問題は、分解の因子が $A$ より密になることです。 ### 6.2. 消去グラフと fill-in の経路定理 対称な非零パターンをグラフで表します。 対称行列 $A\in M_n(\mathbb{R})$ の**グラフ** $G(A)$ を、頂点集合 $\{1,\ldots,n\}$、辺集合 $\{\{i,j\}:i\ne j,\ a_{ij}\ne0\}$ の無向グラフとする。**消去グラフ**の辺集合 $E_k$（頂点集合 $\{k,\ldots,n\}$ の上）を、$E_1=G(A)$ の辺集合とし、 $$ E_{k+1}=\bigl\{\{i,j\}\in E_k:i,j\gt k\bigr\}\cup\bigl\{\{i,j\}:i\ne j,\ \{i,k\}\in E_k,\ \{j,k\}\in E_k\bigr\} $$ で定める。すなわち頂点 $k$ を消し、その隣接点どうしをすべて辺で結ぶ。$E_k$ に属するが $G(A)$ の辺でないものを **fill-in**（の位置）という。
k a b c d → a b c d
この操作は、Schur 補行列の公式 $a^{(k+1)}_{ij}=a^{(k)}_{ij}-a^{(k)}_{ik}a^{(k)}_{kj}/a^{(k)}_{kk}$ で、$a^{(k)}_{ik}$ と $a^{(k)}_{kj}$ がともに非零なら $(i,j)$ 成分が非零になりうる、ということをグラフで書いたものです。消去グラフを $k$ ごとに作らなくても、fill-in の位置は元のグラフ $G(A)$ だけから読み取れます。上の図の d は k と隣り合っていないので、k を消しても d につながる辺は増えません。 次の定理は、これをさらに進めて、どこが非零になるかを消去を実行せずに当てられることを言います。歩道とは、辺をたどって頂点から頂点へ進む道で、同じ頂点を何度通ってもよいものです。直感的には、頂点 $i$ と $j$ が結ばれるのは、自分より番号の小さい頂点だけを経由して行ける場合に限る、ということです。 $A\in M_n(\mathbb{R})$ を対称行列とし、$E_k$ を の消去グラフの辺集合とする。 1. $k\le i,j$, $i\ne j$ について、$\{i,j\}\in E_k$ であることと、$G(A)$ に $i$ から $j$ への歩道 $i=v_0,v_1,\ldots,v_m=j$（$m\ge1$、隣り合う頂点は $G(A)$ の辺で結ばれ、頂点の重複を許す）で、途中の頂点 $v_1,\ldots,v_{m-1}$ がすべて $k$ より小さいものがあることは同値である。 2. $A$ が正定値で、$A^{(k)}$ をその $k$ 段目の Schur 補行列、$R$ を Cholesky 因子とする。$k\le i,j$, $i\ne j$ で $a^{(k)}_{ij}\ne0$ なら $\{i,j\}\in E_k$ であり、$i\lt j$ で $r_{ij}\ne0$ なら $\{i,j\}\in E_i$ である。 3. 2 で、さらに消去の途中で打ち消しが起こらない（各段の $a^{(k+1)}_{ij}=a^{(k)}_{ij}-a^{(k)}_{ik}a^{(k)}_{kj}/a^{(k)}_{kk}$ の右辺の 2 項のうち少なくとも一方が非零なら左辺も非零）とすると、$i\lt j$ について $$ r_{ij}\ne0\iff G(A)\text{ に }i\text{ から }j\text{ への、途中の頂点がすべて }i\text{ より小さい歩道がある}. $$ 証明の方針です。1 は、頂点を 1 つ消すごとに、歩道が経由してよい頂点の範囲（$k$ より小さい番号）が 1 つ広がる、という見方で $k$ についての帰納法を行います。2 と 3 は、更新式の右辺の項が非零かどうかを見るだけです。 1. $k$ についての帰納法です。$k=1$ では、途中の頂点が $1$ より小さい歩道は途中の頂点のない長さ 1 の歩道、つまり $G(A)$ の辺だけなので、主張は $E_1$ の定義そのものです。 $k$ で正しいとし、$i,j\ge k+1$, $i\ne j$ をとります。 ($\Rightarrow$) $\{i,j\}\in E_{k+1}$ とします。$\{i,j\}\in E_k$ なら、帰納法の仮定から途中の頂点が $k$ より小さい（したがって $k+1$ より小さい）歩道があります。そうでなければ $\{i,k\},\{k,j\}\in E_k$ で、帰納法の仮定から $i$ から $k$ へ、$k$ から $j$ への、途中の頂点が $k$ より小さい歩道があります。つなげると、途中の頂点が $k$ 以下、つまり $k+1$ より小さい $i$ から $j$ への歩道です。 ($\Leftarrow$) 途中の頂点がすべて $k+1$ より小さい歩道 $i=v_0,\ldots,v_m=j$ があるとします。途中の頂点に $k$ が現れなければ、途中の頂点はすべて $k$ より小さく、帰納法の仮定から $\{i,j\}\in E_k$ で、$i,j\gt k$ なので $\{i,j\}\in E_{k+1}$ です。$k$ が現れるなら、最初に現れる位置を $s$、最後に現れる位置を $t$ とします。$v_0,\ldots,v_s$ は $i$ から $k$ への歩道で途中の頂点は $k$ でも $k+1$ 以上でもないので $k$ より小さく、帰納法の仮定から $\{i,k\}\in E_k$ です（$i\ne k$）。同様に $v_t,\ldots,v_m$ から $\{k,j\}\in E_k$ です。定義から $\{i,j\}\in E_{k+1}$ です。 2. $k$ についての帰納法で、$a^{(k)}_{ij}\ne0$（$i\ne j$）$\Rightarrow\{i,j\}\in E_k$ を示します。$k=1$ は $G(A)$ の定義です。$a^{(k+1)}_{ij}=a^{(k)}_{ij}-a^{(k)}_{ik}a^{(k)}_{kj}/a^{(k)}_{kk}\ne0$ なら、$a^{(k)}_{ij}\ne0$（このとき $\{i,j\}\in E_k$ で、$i,j\gt k$ なので $E_{k+1}$ に残る）か、$a^{(k)}_{ik}\ne0$ かつ $a^{(k)}_{kj}\ne0$（このとき $\{i,k\},\{k,j\}\in E_k$ で、定義から $\{i,j\}\in E_{k+1}$）です。$R$ については、 の最後の主張を繰り返し使うと $r_{ij}=a^{(i)}_{ij}/\sqrt{a^{(i)}_{ii}}$（$i\lt j$）なので、$r_{ij}\ne0\Rightarrow a^{(i)}_{ij}\ne0\Rightarrow\{i,j\}\in E_i$ です。 3. 打ち消しがないという仮定のもとでは、2 の帰納法の各段で逆向きも成り立ちます（$\{i,j\}\in E_{k+1}$ なら右辺の一方の項が非零で、仮定から左辺も非零）。よって $a^{(k)}_{ij}\ne0\iff\{i,j\}\in E_k$ で、$r_{ij}\ne0\iff\{i,j\}\in E_i$ です。これに 1 を $k=i=\min(i,j)$ で使えば主張です。 経路定理は、fill-in が**数値ではなくグラフと番号の付け方だけで決まる**ことを言っています。そこで、分解の前に頂点の番号を付け替えて（$PAP^{\mathsf{T}}$ を分解して）fill-in を減らすことができます。対称な置換は正定値性を保つので、Cholesky 分解ではピボット選択を気にせず並べ替えを選べます。 対角と第 1 行・第 1 列だけが非零の $n\times n$ 正定値行列（たとえば $a_{11}=n$, $a_{ii}=1$, $a_{1i}=a_{i1}=1/2$）を考えます。$G(A)$ は頂点 $1$ を中心とする星形です。$2\le i\lt j$ について歩道 $i,1,j$ は途中の頂点 $1\lt i$ をもつので、 の 1 から $\{i,j\}\in E_i$ で、打ち消しがなければ $R$ は上三角全体が非零、$\mathrm{nnz}(R)=n(n+1)/2$ です。中心を最後の番号 $n$ に付け替えると、$i\lt j\lt n$ を結ぶ歩道は必ず $n$ を通り、途中の頂点 $n$ は $i$ より大きいので fill-in はなく、$R$ の非零は対角と最後の列だけの $2n-1$ 個です。次数の大きい頂点を後回しにするのが、最小次数順序の考え方です。 $m\times m$ の格子の 5 点差分行列（$n=m^2$。格子点を行ごとに番号付けする自然な順序）は、上下帯幅 $m$ の帯行列です。 から因子は帯の外に非零をもたず、経路定理によれば帯の中はほとんど埋まります（格子の同じ行の点を経由する歩道で、途中の番号はどれも小さくできるから）。したがって $\mathrm{nnz}(L)+\mathrm{nnz}(U)\approx2nm=2n^{3/2}$ で、$A$ の $\mathrm{nnz}\approx5n$ よりずっと多くなります。 まず小さい格子（$m=5$）で、経路定理から計算した $L$ の非零の位置（下三角、対角を除く）と、SuperLU（`scipy.sparse.linalg.splu`）にピボット選択と並べ替えを止めさせた（`permc_spec="NATURAL"`, `diag_pivot_thresh=0`）分解の $L$ の非零の位置を照合します。正定値対称で並べ替えなしの LU 分解の $L$ は、Cholesky 因子の転置と同じ非零パターンをもちます（ の 2 の証明と同じく $L=R^{\mathsf{T}}D^{-1}$）。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64） import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla m = 5 # 5×5 格子, n = 25 T = sp.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(m, m), dtype=float) I = sp.identity(m, dtype=float) A = (sp.kron(I, T) + sp.kron(T, I)).tocsc() n = A.shape[0] adj = [set(A[:, j].nonzero()[0]) - {j} for j in range(n)] def fill_path(i, j): # i から j へ, 途中の頂点がすべて < min(i,j) lo, seen, stack = min(i, j), {i}, [i] while stack: for w in adj[stack.pop()]: if w == j: return True if w < lo and w not in seen: seen.add(w) stack.append(w) return False pred = {(i, j) for j in range(n) for i in range(j + 1, n) if fill_path(i, j)} # L の (i, j), i > j lu = spla.splu(A, permc_spec="NATURAL", diag_pivot_thresh=0.0) Lc = lu.L.tocoo() got = {(i, j) for i, j in zip(Lc.row, Lc.col) if i > j} idn = np.arange(n) print("置換なし", (lu.perm_r == idn).all() and (lu.perm_c == idn).all()) print("辺", sum(map(len, adj)) // 2, "経路定理", len(pred), "splu", len(got), "一致", pred == got) ``` ```text 置換なし True 辺 40 経路定理 104 splu 104 一致 True ``` **検算.** (a) 既知値: $5\times5$ 格子の辺は $2\cdot5\cdot4=40$ 本です。帯の中の下三角の位置は $\sum_{j}\min(5,24-j)=110$ 個で、そのうち最初の格子行の中で隣接しない 6 個（$(3,1)$ など）だけは、途中の頂点を自分より小さい番号にとれないので埋まらず、$110-6=104$ 個です。(b) 別手法: 経路定理（グラフの探索）で数えた位置と、SuperLU が数値の分解で作った $L$ の非零の位置は、集合として完全に一致します（`True`）。行・列の置換が恒等であることも確かめています。離散 Laplace 行列は消去の各段で同符号の項どうしを足すので、打ち消しが起こりません。(c) 次の表で、格子を大きくしたときの増え方を見ます。 $n^{3/2}$ になる理由は数え方だけの話です。帯幅が $m=\sqrt n$ で、各行の帯の中に約 $m$ 個の位置があり、それがほぼ埋まるので、全体で $n\times m=n^{3/2}$ 程度、という計算です。 次に、$m=16,32,64,128$ の格子で、SuperLU の列の並べ替え（`permc_spec`）を自然な順序（`NATURAL`）、$A^{\mathsf{T}}+A$ の最小次数順序（`MMD_AT_PLUS_A`）、列近似最小次数順序（`COLAMD`）に変えて $\mathrm{nnz}(L)+\mathrm{nnz}(U)$ を比べます。列は $m$, $n$, $\mathrm{nnz}(A)$、3 つの順序の $\mathrm{nnz}(L)+\mathrm{nnz}(U)$、自然な順序の値 $/n^{3/2}$、最小次数順序の値 $/(n\log_2n)$ です。時間は環境に依存するので載せません。 ```python # Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1（x86-64） import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla def laplace2d(m): # m×m 格子の 5 点差分（n = m^2） T = sp.diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(m, m), dtype=float) I = sp.identity(m, dtype=float) return (sp.kron(I, T) + sp.kron(T, I)).tocsc() specs = ("NATURAL", "MMD_AT_PLUS_A", "COLAMD") for m in (16, 32, 64, 128): A = laplace2d(m) nnz = [] for s in specs: lu = spla.splu(A, permc_spec=s, diag_pivot_thresh=0.0, options=dict(SymmetricMode=True)) nnz.append(lu.L.nnz + lu.U.nnz) n = m * m print(f"{m:>3} {n:>6} {A.nnz:>6} " + " ".join(f"{k:>8}" for k in nnz) + f" {nnz[0] / n**1.5:.2f} {nnz[1] / (n * np.log2(n)):.2f}") ``` ```text 16 256 1216 8222 4192 5450 2.01 2.05 32 1024 4992 65598 23016 37334 2.00 2.25 64 4096 20224 524414 126532 220624 2.00 2.57 128 16384 81408 4194558 664094 1221240 2.00 2.90 ``` **検算.** (a) 既知値: $\mathrm{nnz}(A)=n+4m(m-1)$ で、$m=16$ なら $256+960=1216$ です。自然な順序の $\mathrm{nnz}(L)+\mathrm{nnz}(U)$ は帯幅 $m$ の帯をほぼ埋めた $2nm$ に一致し、比 $\mathrm{nnz}/n^{3/2}$ は $2.00$ です。(b) 別手法: 最小次数順序は自然な順序の $1/2$〜$1/6$、COLAMD はその中間で、$m=128$ では自然な順序の約 $1/6$ と $1/3$ です。(c) 増え方: 格子幅 $m$ を 2 倍（$n$ を 4 倍）するごとに、自然な順序は $8$ 倍（$n^{3/2}$）、最小次数順序は $5.2$〜$5.5$ 倍で、$n\log_2n$ の増え方（約 $4.7$〜$5.0$ 倍）よりやや速く、比が $2.05$ から $2.90$ にゆっくり増えています。最小次数順序は発見的な方法で、増え方の定理はありません（次の注意）。 fill-in を最小にする並べ替えを求める問題は NP 困難です（Yannakakis 1981）。実用の方法は発見的です。**逆 Cuthill–McKee**（RCM）は幅優先探索で帯幅を小さくし、帯行列・スカイライン形式の分解に向きます（`scipy.sparse.csgraph.reverse_cuthill_mckee`）。**最小次数順序**は、消去グラフでその時点の次数が最小の頂点を次に消す方法で、近似最小次数（AMD、Amestoy–Davis–Duff 1996）が標準です。**入れ子分割**は、グラフを小さな分離集合で 2 つに分け、分離集合を最後に番号付けすることを再帰的に繰り返します（George 1973、METIS）。詳しくは Davis『Direct Methods for Sparse Linear Systems』第 7 章、George–Liu『Computer Solution of Large Sparse Positive Definite Systems』を参照してください。 2 次元の $m\times m$ 格子（$n=m^2$）を入れ子分割で番号付けすると、Cholesky 因子の非零は $O(n\log n)$、分解の演算量は $O(n^{3/2})$ で、後者はどの順序でも改善できない（下界がある）ことが知られています（George 1973、Hoffman–Martin–Rose 1973）。自然な順序の $O(n^{3/2})$ 個・$O(n^2)$ 回と比べて大きな差です。3 次元の格子では非零 $O(n^{4/3})$、演算量 $O(n^2)$ で、直接法が反復法に比べて不利になる理由の 1 つです。評価はグラフの分割の議論が長いので、この記事では扱いません。 非対称な疎行列の LU 分解（SuperLU、UMFPACK）では、数値的安定性のためのピボット選択と、fill-in を抑える並べ替えが衝突します。SuperLU は列の並べ替えを分解の前に決め、行のピボット選択には**しきい値ピボット選択**を使います。対角成分が列の最大の $\tau$ 倍以上（`diag_pivot_thresh` $=\tau\in[0,1]$）なら対角を選び、そうでなければ最大の成分を選ぶので、乗数は $|l_{ij}|\le1/\tau$ に抑えられ、成長因子の評価は部分ピボット選択より弱くなります。$\tau=1$ が部分ピボット選択、$\tau=0$ はピボット選択なしです。Li「An overview of SuperLU」（2005）、Davis 第 6 章を参照してください。 ## 7. 演習 $A=\begin{pmatrix}4&2&2\\2&5&3\\2&3&6\end{pmatrix}$ の Cholesky 因子 $R$ を求め、$A$ が正定値であることを確かめよ。 $r_{11}=2$, $r_{12}=2/2=1$, $r_{13}=2/2=1$。Schur 補行列は $\begin{pmatrix}5-1&3-1\\3-1&6-1\end{pmatrix}=\begin{pmatrix}4&2\\2&5\end{pmatrix}$ で、$r_{22}=2$, $r_{23}=1$、最後に $5-1=4$ から $r_{33}=2$ です。 $$ R=\begin{pmatrix}2&1&1\\0&2&1\\0&0&2\end{pmatrix} $$ で対角がすべて正なので、 により $A$ は正定値です。 $n\ge4$ の巡回グラフ（辺 $\{i,i+1\}$, $i=1,\ldots,n-1$ と $\{1,n\}$）をグラフにもつ正定値行列を自然な番号で Cholesky 分解するとき、経路定理から fill-in の位置を求めよ。打ち消しはないとする。 $i\lt j$ で $\{i,j\}$ が辺でない組について、途中の頂点が $i$ より小さい歩道を探します。$i\ge2$, $j=n$ なら、$i,i-1,\ldots,1,n$ は途中の頂点 $i-1,\ldots,1$ がすべて $i$ より小さいので、$r_{in}\ne0$ です。それ以外（$2\le i\lt j\lt n$, $j\ne i+1$）では、$i$ から $j$ に行くには $i+1$ 側を通るか（途中に $i$ 以上の頂点が現れる）、$1$ と $n$ を経由するか（途中に $n\gt i$ が現れる）しかないので、fill-in はありません。よって fill-in は最後の列の $(i,n)$, $2\le i\le n-2$ の $n-3$ 個です（$(n-1,n)$ は元の辺）。$R$ は上二重対角に最後の列を加えた形になります。 $n=10^6$、上下帯幅 $p=100$ の正定値対称な帯行列について、 の演算量の上界と、密行列としての Cholesky 分解の演算量 $\tfrac13n^3$ を比べよ。因子の記憶量（語数）も比べよ。 帯の評価は $n(p+1)(2p+1)=10^6\cdot101\cdot201\approx2.0\times10^{10}$ 回、密行列では $\tfrac13\cdot10^{18}\approx3.3\times10^{17}$ 回で、約 $1.6\times10^7$ 倍の差です。記憶量は帯の因子が $n(p+1)\approx10^8$ 語（倍精度で約 0.8 GB）、密な上三角は $n(n+1)/2\approx5\times10^{11}$ 語（約 4 TB）です。 ## 参考文献 - N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — 第 10 章（Cholesky 分解）、第 11 章（対称不定値行列と Bunch–Kaufman）。 - G. H. Golub, C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins University Press, 2013 — 第 4 章（特殊な構造の線形方程式：正定値・帯行列・対称不定値）。 - T. A. Davis, *Direct Methods for Sparse Linear Systems*, SIAM, 2006 — 第 2 章（格納形式）、第 4 章（Cholesky 分解と消去木）、第 7 章（並べ替え）。 - A. George, J. W. H. Liu, *Computer Solution of Large Sparse Positive Definite Systems*, Prentice-Hall, 1981 — 消去グラフ、帯・包絡形式、入れ子分割。 - D. J. Rose, R. E. Tarjan, G. S. Lueker, "Algorithmic aspects of vertex elimination on graphs", *SIAM Journal on Computing* 5 (1976), 266–283. — fill-in の経路による特徴づけ。 - A. George, "Nested dissection of a regular finite element mesh", *SIAM Journal on Numerical Analysis* 10 (1973), 345–363.