行列計算の認証と最適化の証明書：慣性・Rump の正定値性検証・Sturm 列・双対性

0. この記事の要点

コンピューターが出す数値には丸めの誤差があり、「たぶん正しい」だけでは数学の証明に使えません。この記事は、行列の計算結果から「この行列は確かに正定値（固有値がすべて正）である」「この最適化問題の最適値は確かに $\gamma$ 以上である」といった主張を、誤差を見積もったうえで確実に言える形にする方法を扱います。読み終えると、近似的な計算結果を、証明の一部として使える不等式に変える考え方がわかります。

• 対称行列を $A=LDL^{\mathsf{T}}$ と分解できれば、Sylvester の慣性法則により $A$ の正・負・零の固有値の個数は $D$ の対角成分の符号で決まります。有理数で厳密に分解すれば、それがそのまま証明書です。
• 浮動小数の Cholesky 分解が最後まで進めば、$\widetilde R^{\mathsf{T}}\widetilde R=A+\Delta A$, $|\Delta A|\le\gamma_{n+1}|\widetilde R^{\mathsf{T}}||\widetilde R|$ が成り立ちます。Rump の方法は、$A-cI$ の Cholesky 分解が成功すれば $A$ は正定値、となるシフト $c$ をこの評価から決めます。
• 三重対角行列では、LDLᵀ のピボットの符号の数え上げ（Sturm 列）が「$x$ 未満の固有値の個数」を与え、二分法で固有値を確実に包めます。
• 近似固有対の残差から、Weyl の不等式・Bauer–Fike の定理・Gershgorin の定理で固有値の位置を事後評価できます。
• LP と SDP の弱双対性により、近似的な双対解からでも、最適値の厳密な下界が計算できます。多項式の非負性の SOS 証明書（平方和の形に書けることを使った証明書）も、半正定値行列の検証に帰着します。

この記事で使う事前知識

• 正定値・半正定値: 対称行列の固有値がすべて正（半正定値は 0 以上）であること。→ 正定値
• Sylvester の慣性法則: 正則な行列 $S$ で $S^{\mathsf{T}}AS$ に変えても、正・負・零の固有値の個数は変わない、という法則。→ 行列の分解と標準形
• Weyl の不等式: 対称行列に別の対称行列 $B$ を足しても、固有値は $\|B\|$ 以下しか動かない、という評価。→ 行列の固有値評価と不等式
• Gershgorin の定理: 行列の固有値は、対角成分を中心とする円板のどれかに入る、という定理。→ 行列の固有値評価と不等式
• 丸め: 実数を浮動小数で表せる近くの数に置き換える操作。→ 丸め
• 浮動小数の誤差モデル: 1 回の演算で生じる相対誤差が単位丸め $u$ 以下、という見積もり。→ 浮動小数と区間演算

1. 動機

固有値の評価や最適化の結果を証明に使いたい場面では、「行列が正定値である」「この最適値は $\gamma$ 以上である」という主張が必要になります。浮動小数の固有値計算はこの主張を直接は保証しません。最小固有値が $10^{-15}$ と出ても、それが丸め誤差なのか本物なのかは区別できないからです。

本章では、行列について「計算機で確かめられる十分条件」を作ります。考え方は 2 つです。1 つは厳密な代数的証明書を作る方法で、有理数の LDLᵀ 分解や SOS 表示がそれです。もう 1 つは浮動小数の計算に後ろ向き誤差解析を組み合わせる方法で、Rump の正定値性検証がその代表です。どちらも第 1 章の誤差モデルと、線形代数の 2 つの定理、Sylvester の慣性法則と Weyl の不等式に帰着します。第 4 章の固有値の下界は、本章の慣性の判定を部品として使います。

2. 準備

$\mathbb{R}^{n\times n}$ の対称行列 $A$ の最小・最大の固有値を $\lambda_{\min}(A)$, $\lambda_{\max}(A)$、ユークリッドノルムから誘導される作用素ノルムを $\|\cdot\|$ と書きます。行列やベクトルの絶対値 $|A|$ は成分ごとにとり、不等号 $|E|\le F$ も成分ごとです。$A\succ0$ は正定値、$A\succeq0$ は半正定値を表します。

補題 2.1（成分ごとの評価からノルムの評価へ）✓

$|E|\le F$（$F$ は成分が非負）ならば $\|E\|\le\|F\|$ である。また $\boldsymbol{d}\in\mathbb{R}^n$ について $\|\boldsymbol{d}\boldsymbol{d}^{\mathsf{T}}\|=\|\boldsymbol{d}\|^2$ である。

証明（補題 2.1）

$\boldsymbol{x}\in\mathbb{R}^n$ について、成分ごとに $|E\boldsymbol{x}|\le|E||\boldsymbol{x}|\le F|\boldsymbol{x}|$ で、成分が非負のベクトルではノルムは成分について単調なので $\|E\boldsymbol{x}\|\le\|F|\boldsymbol{x}|\|\le\|F\|\,\|\boldsymbol{x}\|$ です。後半は $\boldsymbol{d}\boldsymbol{d}^{\mathsf{T}}\boldsymbol{x}=\langle\boldsymbol{d},\boldsymbol{x}\rangle\boldsymbol{d}$ と Cauchy–Schwarz の不等式から $\le$ が、$\boldsymbol{x}=\boldsymbol{d}$ で等号が出ます。

∎

もう 1 つの道具は Weyl の不等式（定理 3.2［行列の固有値評価と不等式］）です。対称行列 $A$, $B$ について、大きい順に並べた $k$ 番目の固有値は $|\lambda_k(A+B)-\lambda_k(A)|\le\|B\|$ を満たします。とくに $\lambda_{\min}(A+B)\ge\lambda_{\min}(A)-\|B\|$ です。

3. LDLᵀ 分解と慣性による認証

対称行列 $A$ の正・負・零の固有値の個数（重複度込み）の組 $(n_+,n_-,n_0)$ を慣性といいます。

定理 3.1（LDLᵀ 分解による慣性の認証）✓

$A$ を実対称行列とし、対角成分が $1$ の下三角行列 $L$ と対角行列 $D=\operatorname{diag}(d_1,\ldots,d_n)$ で $A=LDL^{\mathsf{T}}$ と書けたとする。このとき $A$ の慣性は $(\#\{i:d_i\gt0\},\ \#\{i:d_i\lt0\},\ \#\{i:d_i=0\})$ である。とくに $A\succ0$ であることと、すべての $d_i\gt0$ であることは同値である。

証明（定理 3.1）

$L$ は行列式 $1$ で正則です。$S=(L^{-1})^{\mathsf{T}}$ とおくと $S^{\mathsf{T}}AS=L^{-1}A(L^{-1})^{\mathsf{T}}=D$ なので、Sylvester の慣性法則（定理 6.1［行列の分解と標準形］ の 2）により $A$ と $D$ の慣性は等しく、対角行列 $D$ の固有値は $d_i$ です。最後の主張は、正定値と固有値がすべて正であることの同値性（定理 6.3［スペクトル定理］ の 2）からです。

∎

$L$, $D$ はガウスの消去法で計算でき（行の入れ替えが要らない場合）、成分が有理数の行列なら、有理数の演算で厳密に求まります。こうして得た $(L,D)$ は、誰でも $LDL^{\mathsf{T}}=A$ を掛け算で確かめられる証明書です。たとえば

$$A=\begin{pmatrix}2&1&0\\1&2&1\\0&1&2\end{pmatrix}$$

の消去法のピボットは $d_1=2$, $d_2=2-\tfrac12=\tfrac32$, $d_3=2-\tfrac{1}{3/2}=\tfrac43$ ですべて正なので、$A\succ0$ です。

浮動小数で分解した場合は、得られるのは近くの行列の慣性です。それを $A$ に移すのが次の命題です。

命題 3.2（慣性の摂動による移転）✓

$A$, $\widetilde A$ を対称行列とし、$\widetilde A$ の固有値の絶対値の最小値が $\|A-\widetilde A\|$ より大きいとする。このとき $A$ と $\widetilde A$ の慣性は等しい。

証明（命題 3.2）

$\varepsilon=\|A-\widetilde A\|$ とおきます。Weyl の不等式から、大きい順の $k$ 番目の固有値について $|\lambda_k(A)-\lambda_k(\widetilde A)|\le\varepsilon\lt|\lambda_k(\widetilde A)|$ なので、$\lambda_k(A)$ は $0$ でなく $\lambda_k(\widetilde A)$ と同符号です。

∎

4. Cholesky 分解の誤差解析と Rump の方法

正定値性の検証では、もっと直接的な方法が使えます。Cholesky 分解 $A=R^{\mathsf{T}}R$（$R$ は対角成分が正の上三角行列）は、$A\succ0$ のときに限って最後まで計算できます。浮動小数で計算すると $\widetilde R^{\mathsf{T}}\widetilde R\ne A$ ですが、そのずれは評価できます。

アルゴリズムを固定します。$j=1,\ldots,n$ の順に、$i=1,\ldots,j$ について $s=a_{ij}$ から $s\leftarrow\mathrm{fl}(s-\mathrm{fl}(\tilde r_{ki}\tilde r_{kj}))$ を $k=1,\ldots,i-1$ の順に行い、$i\lt j$ なら $\tilde r_{ij}=\mathrm{fl}(s/\tilde r_{ii})$、$i=j$ なら $s\gt0$ を確かめて $\tilde r_{jj}=\mathrm{fl}(\sqrt s)$ とします。$s\le0$ になったら失敗として止めます。以下、アンダーフローとオーバーフローは起きないとし、$(n+1)u\lt\tfrac12$ とします。

次の補題のねらいは、引き算を重ねて 1 つの成分を求める 1 ステップで、各演算の小さなずれ（相対誤差が $u$ 以下）を「元の数にかかる係数 $1+\theta$」の形にまとめ直すことです。$\gamma_m$ は第 1 章で定めた、誤差が $m$ 回重なったときの係数のずれの上界です。証明では、各ステップの式を順に割って足し合わせ、$\theta$ を取り出します。

補題 4.1（1 ステップの丸め誤差）✓

$c,a_k,b_k,d\in\mathbb{F}$（$k=1,\ldots,m$）とし、$s_0=c$, $s_k=\mathrm{fl}(s_{k-1}-\mathrm{fl}(a_kb_k))$ とする。

1. $y=\mathrm{fl}(s_m/d)$ ならば $c=dy(1+\theta^{(0)})+\sum_{k=1}^ma_kb_k(1+\theta^{(k)})$ となる $|\theta^{(0)}|\le\gamma_{m+1}$, $|\theta^{(k)}|\le\gamma_{m+1}$ がある。
2. $s_m\gt0$, $y=\mathrm{fl}(\sqrt{s_m})$ ならば $c=y^2(1+\theta^{(0)})+\sum_{k=1}^ma_kb_k(1+\theta^{(k)})$ となる $|\theta^{(0)}|\le\gamma_{m+2}$, $|\theta^{(k)}|\le\gamma_{m+2}$ がある。

証明（補題 4.1）

定理 2.4［浮動小数と区間演算］ から $s_k=(s_{k-1}-a_kb_k(1+\varepsilon_k))(1+\delta_k)$（$|\varepsilon_k|,|\delta_k|\le u$）です。両辺を $\prod_{j\le k}(1+\delta_j)$ で割って足し合わせると

$$\frac{s_m}{\prod_{j=1}^m(1+\delta_j)}=c-\sum_{k=1}^ma_kb_k\frac{1+\varepsilon_k}{\prod_{j=1}^{k-1}(1+\delta_j)}$$

です。1 では $y=(s_m/d)(1+\delta_0)$ なので $s_m=dy/(1+\delta_0)$ を代入し、$1+\theta^{(0)}=\prod_{j=0}^m(1+\delta_j)^{-1}$（因子 $m+1$ 個）、$1+\theta^{(k)}=(1+\varepsilon_k)\prod_{j\lt k}(1+\delta_j)^{-1}$（因子 $k$ 個）とおけば、補題 3.2［浮動小数と区間演算］ から主張が出ます。2 では $y=\sqrt{s_m}(1+\delta_0)$ から $s_m=y^2/(1+\delta_0)^2$ で、$\theta^{(0)}$ の因子が $m+2$ 個になります。

∎

この定理は、浮動小数で得た $\widetilde R$ が、元の $A$ を少しだけ変えた行列 $A+\Delta A$ の Cholesky 因子になっている、と言っています。ずれ $\Delta A$ の大きさが、計算結果の $\widetilde R$ 自身で抑えられるのがポイントです。

定理 4.2（Cholesky 分解の後ろ向き誤差）✓

対称行列 $A\in\mathbb{F}^{n\times n}$ の浮動小数の Cholesky 分解が失敗せずに $\widetilde R$ を返したとする。このとき

$$\widetilde R^{\mathsf{T}}\widetilde R=A+\Delta A,\qquad|\Delta A|\le\gamma_{n+1}|\widetilde R|^{\mathsf{T}}|\widetilde R|,$$

さらに各 $j$ について $\sum_k\tilde r_{kj}^2\le a_{jj}/(1-\gamma_{n+1})$ が成り立つ。

証明（定理 4.2）

$i\lt j$ の成分には 補題 4.1 の 1 を $c=a_{ij}$, $m=i-1$, $d=\tilde r_{ii}$, $y=\tilde r_{ij}$, $a_kb_k=\tilde r_{ki}\tilde r_{kj}$ で使うと、$a_{ij}=\sum_{k\le i}\tilde r_{ki}\tilde r_{kj}(1+\theta^{(k)})$, $|\theta^{(k)}|\le\gamma_i\le\gamma_{n+1}$ です（$k=i$ の項が $dy$ の項）。対角成分には 2 を $m=j-1$ で使い、$a_{jj}=\sum_{k\le j}\tilde r_{kj}^2(1+\theta^{(k)})$, $|\theta^{(k)}|\le\gamma_{j+1}\le\gamma_{n+1}$ です。$\widetilde R$ は上三角なので $(\widetilde R^{\mathsf{T}}\widetilde R)_{ij}=\sum_{k\le\min(i,j)}\tilde r_{ki}\tilde r_{kj}$ で、差をとれば $|\Delta a_{ij}|\le\gamma_{n+1}\sum_k|\tilde r_{ki}||\tilde r_{kj}|$ です。対称性から $i\gt j$ も同じです。最後の不等式は、対角成分の式から $a_{jj}\ge(1-\gamma_{n+1})\sum_k\tilde r_{kj}^2$ となることによります。

∎

次の定理のねらいは、Cholesky 分解が成功したという計算結果から $A\succ0$ を結論することです。分解するのは $A$ そのものではなく、対角成分から $c$ を引いて固有値を全体に $c$ だけ下げた行列です。それでも分解が成功し、誤差の大きさ $\beta$ が $c$ より小さければ、元の最小固有値は $c-\beta$ 以上、つまり正だと言えます。証明は、分解の成功から対角が正であること、誤差のノルムを $\beta$ 以下で抑えること、Weyl の不等式で最小固有値に移すこと、の順で進みます。

定理 4.3（Rump の正定値性検証）✓

対称行列 $A\in\mathbb{F}^{n\times n}$ と浮動小数 $c\ge0$ をとり、

$$c\gt\beta:=\frac{\gamma_{n+1}}{1-\gamma_{n+1}}\operatorname{tr}A+u\max_ia_{ii}$$

とする。$\widetilde A$ を、$A$ の対角成分を $\mathrm{fl}(a_{ii}-c)$ に替えた行列とする。$\widetilde A$ の浮動小数の Cholesky 分解が失敗せずに終われば、$A\succ0$ であり、$\lambda_{\min}(A)\ge c-\beta$ である。

証明（定理 4.3）

分解の成功から、定理 4.2 により $\tilde a_{jj}\ge(1-\gamma_{n+1})\sum_k\tilde r_{kj}^2\gt0$ です。$\mathrm{fl}$ の単調性から $0\lt\tilde a_{jj}=\mathrm{fl}(a_{jj}-c)\le\mathrm{fl}(a_{jj})=a_{jj}$ です。

$\widetilde R^{\mathsf{T}}\widetilde R=\widetilde A+\Delta$ と書きます。$\widetilde R$ の第 $j$ 列を $\boldsymbol{r}_j$ とすると、Cauchy–Schwarz の不等式と 定理 4.2 から

$$(|\widetilde R|^{\mathsf{T}}|\widetilde R|)_{ij}=\sum_k|\tilde r_{ki}||\tilde r_{kj}|\le\|\boldsymbol{r}_i\|\|\boldsymbol{r}_j\|\le d_id_j,\qquad d_j=\sqrt{\tilde a_{jj}/(1-\gamma_{n+1})}$$

です。補題 2.1 を 2 回使うと

$$\|\Delta\|\le\gamma_{n+1}\|\boldsymbol{d}\boldsymbol{d}^{\mathsf{T}}\|=\gamma_{n+1}\|\boldsymbol{d}\|^2=\frac{\gamma_{n+1}}{1-\gamma_{n+1}}\sum_j\tilde a_{jj}\le\frac{\gamma_{n+1}}{1-\gamma_{n+1}}\operatorname{tr}A.$$

次に $F=A-cI-\widetilde A$ は対角行列で、定理 2.4［浮動小数と区間演算］ の 1 から $|f_{ii}|=|(a_{ii}-c)-\mathrm{fl}(a_{ii}-c)|\le u\,\tilde a_{ii}\le u\max_ia_{ii}$ です。以上から

$$A=\widetilde R^{\mathsf{T}}\widetilde R+cI-\Delta+F$$

で、$\widetilde R^{\mathsf{T}}\widetilde R\succeq0$ と Weyl の不等式により $\lambda_{\min}(A)\ge c-\|\Delta\|-\|F\|\ge c-\beta\gt0$ です。

∎

計算機で使うときは、$\beta$ を上向き丸めで計算し、$c$ をそれより大きくとります。シフト $c$ が要るのは、丸め誤差で $\lambda_{\min}$ が $\beta$ 程度動き得るからです。$c$ を引いても分解が通るほど $\lambda_{\min}(A)$ が大きければ、正定値性が証明されます。計算量は Cholesky 分解 1 回分です。Rump の原論文は、アンダーフローを含む場合にも対応した形でこの定理を述べています。

数直線上で見た Rump の方法。シフト c を β より大きくとり、A から c を引いた行列の分解が成功すれば、λ_min(A) は c−β 以上、つまり 0 より右にある。

例 4.4（Hilbert 行列の正定値性）

$H_n=(1/(i+j-1))_{i,j}$ は正定値ですが、最小固有値は $n=10$ で約 $10^{-13}$、$n=12$ で約 $10^{-16}$ と急速に小さくなります。成分は浮動小数ではないので、ここでは浮動小数に丸めた行列 $\mathrm{fl}(H_n)$ の正定値性を検証します（$H_n$ 自身については、丸めの差 $\|H_n-\mathrm{fl}(H_n)\|$ を Weyl の不等式で加えれば同じ議論ができます）。$\beta$ と $c=2\beta$ は nextafter による上向き丸めで計算します。

import math

def up(x):   return math.nextafter(x, math.inf)
def down(x): return math.nextafter(x, -math.inf)
u = 2.0 ** -53

def cholesky(A):
    """浮動小数の Cholesky 分解（本文の順序）。失敗したら None。"""
    n = len(A)
    R = [[0.0] * n for _ in range(n)]
    for j in range(n):
        for i in range(j + 1):
            s = A[i][j]
            for k in range(i):
                s = s - R[k][i] * R[k][j]
            if i < j:
                R[i][j] = s / R[i][i]
            elif s > 0:
                R[j][j] = math.sqrt(s)               # sqrt は IEEE 754 で正しく丸められる
            else:
                return None
    return R

def verify_pd(A):
    n = len(A)
    if min(A[i][i] for i in range(n)) <= 0:
        return False, None
    tr = 0.0
    for i in range(n):
        tr = up(tr + A[i][i])                        # tr A の上界
    t = up((n + 1) * u / down(1.0 - 2 * (n + 1) * u))   # γ/(1-γ) の上界
    beta = up(up(t * tr) + up(u * max(A[i][i] for i in range(n))))
    c = up(2.0 * beta)
    At = [row[:] for row in A]
    for i in range(n):
        At[i][i] = A[i][i] - c                       # fl(a_ii - c)
    if cholesky(At) is None:
        return False, None
    return True, down(c - beta)                      # λ_min(A) >= c - β

for n in range(8, 13):
    H = [[1.0 / (i + j + 1) for j in range(n)] for i in range(n)]
    print(n, *verify_pd(H))

8 True 2.131206743994237e-15
9 True 2.4209789186848627e-15
10 True 2.7162506501538186e-15
11 False None
12 False None

$n\le10$ では $\mathrm{fl}(H_n)\succ0$ が証明されました。$n=11$ では最小固有値（約 $3.4\times10^{-15}$）がシフト $c\approx6.0\times10^{-15}$ より小さく、検証は失敗します。失敗は「正定値でない」ことを意味しません。より高い精度で計算し直せば検証できます。

5. Sturm 列と二分法

対称な三重対角行列 $T$ の対角成分を $\alpha_1,\ldots,\alpha_n$、副対角成分を $\beta_1,\ldots,\beta_{n-1}$ とします。$x\in\mathbb{R}$ について

$$d_1(x)=\alpha_1-x,\qquad d_k(x)=\alpha_k-x-\frac{\beta_{k-1}^2}{d_{k-1}(x)}\quad(k=2,\ldots,n)$$

と定めます（途中で $d_{k-1}(x)=0$ になったら定義しない）。

定理 5.1（Sturm 列による固有値の数え上げ）✓

$d_1(x),\ldots,d_n(x)$ がすべて定義されて $0$ でないとき、$T$ の $x$ より小さい固有値の個数（重複度込み）は $\#\{k:d_k(x)\lt0\}$ に等しい。また $T-xI$ の左上 $k$ 次の主小行列式を $p_k(x)$（$p_0=1$）とすると $p_k(x)=d_1(x)\cdots d_k(x)$ であり、上の個数は列 $p_0,p_1(x),\ldots,p_n(x)$ の符号の変化の回数に等しい。

証明（定理 5.1）

$L$ を対角成分 $1$、副対角成分 $l_k=\beta_k/d_k$ の下二重対角行列、$D=\operatorname{diag}(d_k)$ とします。$LDL^{\mathsf{T}}$ の $(k,k)$ 成分は $d_k+l_{k-1}^2d_{k-1}=d_k+\beta_{k-1}^2/d_{k-1}=\alpha_k-x$、$(k+1,k)$ 成分は $l_kd_k=\beta_k$、それ以外の非対角成分は $0$ なので、$LDL^{\mathsf{T}}=T-xI$ です。定理 3.1 により $T-xI$ の負の固有値の個数は $\#\{k:d_k\lt0\}$ で、$T-xI$ の固有値は $T$ の固有値から $x$ を引いたものです。主小行列式については、左上 $k$ 次のブロックが $L_kD_kL_k^{\mathsf{T}}$（$L_k$, $D_k$ は左上 $k$ 次のブロック）なので $p_k=\det D_k=d_1\cdots d_k$ です。$d_k=p_k/p_{k-1}$ が負であることは、$p_{k-1}$ と $p_k$ の符号が変わることと同値です。

∎

この定理を使うと、$x$ を動かして個数の変わる点を二分法で追い、固有値を任意の幅の区間に閉じ込められます。$d_k(x)$ を区間演算で計算し、どの区間も $0$ を含まなければ、包含原理（定理 4.4［浮動小数と区間演算］）により符号は確定し、個数は厳密です。浮動小数で素朴に計算した場合でも、得られる個数はわずかに摂動した行列の正確な個数になることが知られています（Golub–Van Loan の第 8 章）。

x 未満の固有値の個数 N(x) は固有値で 1 段ずつ上がる階段関数。2 点 a, b で N(a)=1, N(b)=2 と数えられれば、[a, b) に 2 番目の固有値がちょうど 1 つある。中点で数え直すと区間は半分になる。

次の例は、直交多項式の記事の結果を使います。Legendre 多項式から作る三重対角の対称行列（ヤコビ行列）の固有値が、Gauss–Legendre 求積の節点になります。節点の位置を、Sturm 列の個数の数え上げで挟みます。

例 5.2（Legendre 多項式の零点の包含）

$[-1,1]$ の Legendre 多項式に対するヤコビ行列は $\alpha_k=0$, $\beta_k^2=k^2/(4k^2-1)$ で、定理 6.3［直交多項式］ によりその固有値は Gauss–Legendre 求積の節点です。$n=3$ では $\beta_1^2=\tfrac13$, $\beta_2^2=\tfrac4{15}$ です。有理数で厳密に計算すると、

• $x=\tfrac{77}{100}$: $d_1=-\tfrac{77}{100}$, $d_2=-\tfrac{7787}{23100}$, $d_3=-\tfrac{77}{100}+\tfrac{92400}{116805}\approx+0.0211$。負は 2 個。
• $x=\tfrac{78}{100}$: $d_1=-0.78$, $d_2\approx-0.3526$, $d_3\approx-0.0238$。負は 3 個。

よって最大の節点は $[0.77,0.78)$ にあります（実際 $\sqrt{3/5}=0.7745\ldots$）。計算は有理数で厳密なので、これは証明です。第 2 章の注意で述べた「節点を区間で包む」作業は、この方法でも行えます。

6. 事後誤差評価：Weyl・Gershgorin・Bauer–Fike

近似固有値・固有ベクトルが手元にあるとき、その誤差を事後に評価する定理を集めます。

定理 6.1（残差による固有値の評価）✓

$A$ を実対称行列、$\boldsymbol{x}\ne\boldsymbol{0}$、$\mu\in\mathbb{R}$ とし、$\boldsymbol{r}=A\boldsymbol{x}-\mu\boldsymbol{x}$ とおく。このとき $|\lambda-\mu|\le\|\boldsymbol{r}\|/\|\boldsymbol{x}\|$ を満たす $A$ の固有値 $\lambda$ が存在する。

証明（定理 6.1）

$A$ の正規直交固有基底 $\boldsymbol{v}_i$（固有値 $\lambda_i$）で $\boldsymbol{x}=\sum_ic_i\boldsymbol{v}_i$ と展開すると、$\boldsymbol{r}=\sum_i(\lambda_i-\mu)c_i\boldsymbol{v}_i$ なので

$$\|\boldsymbol{r}\|^2=\sum_i(\lambda_i-\mu)^2c_i^2\ge\min_i(\lambda_i-\mu)^2\sum_ic_i^2=\min_i(\lambda_i-\mu)^2\|\boldsymbol{x}\|^2$$

です。

∎

定理 6.2（Bauer–Fike の定理）✓

$A\in\mathbb{C}^{n\times n}$ が $A=V\Lambda V^{-1}$, $\Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n)$ と対角化できるとし、$\kappa(V)=\|V\|\|V^{-1}\|$ とおく。$\mu$ が $A+E$ の固有値ならば $\min_i|\mu-\lambda_i|\le\kappa(V)\|E\|$ である。

証明（定理 6.2）

$\mu$ が $A$ の固有値なら左辺は $0$ です。そうでなければ $\mu I-\Lambda$ は正則で、

$$\mu I-A-E=V(\mu I-\Lambda)\bigl(I-(\mu I-\Lambda)^{-1}V^{-1}EV\bigr)V^{-1}$$

が特異なので、$M=(\mu I-\Lambda)^{-1}V^{-1}EV$ について $I-M$ は特異です。$(I-M)\boldsymbol{y}=\boldsymbol{0}$, $\boldsymbol{y}\ne\boldsymbol{0}$ なら $\|\boldsymbol{y}\|=\|M\boldsymbol{y}\|\le\|M\|\|\boldsymbol{y}\|$ なので $\|M\|\ge1$ です。$\|(\mu I-\Lambda)^{-1}\|=1/\min_i|\mu-\lambda_i|$ から

$$1\le\|M\|\le\frac{\kappa(V)\|E\|}{\min_i|\mu-\lambda_i|}$$

です。

∎

対称行列では $V$ を直交行列にとれるので $\kappa(V)=1$ で、固有値は $\|E\|$ しか動きません（Weyl の不等式と同じ結論）。非正規行列では $\kappa(V)$ が大きくなり、固有値は摂動に敏感になります。成分ごとの情報だけで固有値の位置を囲むのが Gershgorin の定理（定理 2.2［行列の固有値評価と不等式］）で、近似固有ベクトル行列 $X$ で $X^{-1}AX$ をほぼ対角にしてから使うと、小さな円板が得られます。

全固有値の認証。 実対称行列 $A$ について、$n$ 個の近似固有対に 定理 6.1 を使って得た区間が互いに交わらなければ、各区間は異なる固有値を含み、固有値は重複度込みで $n$ 個なので、各区間にちょうど 1 個ずつ入ります。残差 $\boldsymbol{r}$ は区間演算で、$\|\boldsymbol{r}\|$ は上向き丸め、$\|\boldsymbol{x}\|$ は下向き丸めで評価します。区間が交わる場合（近接・重複固有値）は、定理 5.1 や 命題 3.2 で個数を数えます。

7. LP・SDP の双対性と証明書

最適化問題には、値を最小にする「主問題」と、その値を下から押さえるための「双対問題」が組になっています。主問題の最適値を直接求めにくくても、双対問題の条件を満たす解が 1 つあれば、最適値はその値以上だと言えます。この節で使うのは、この弱双対性という性質だけです。

定義 7.1（線形計画の主問題と双対問題）↗

$A\in\mathbb{R}^{m\times n}$, $\boldsymbol{b}\in\mathbb{R}^m$, $\boldsymbol{c}\in\mathbb{R}^n$ に対し、主問題 (P) $\min\ \boldsymbol{c}^{\mathsf{T}}\boldsymbol{x}$ s.t. $A\boldsymbol{x}=\boldsymbol{b}$, $\boldsymbol{x}\ge\boldsymbol{0}$ と、双対問題 (D) $\max\ \boldsymbol{b}^{\mathsf{T}}\boldsymbol{y}$ s.t. $A^{\mathsf{T}}\boldsymbol{y}\le\boldsymbol{c}$ を考える。(P) の最適値を $p^\ast$ と書く。

次の定理で、$\bar{\boldsymbol{x}}$ は実行可能解の各成分の上界、$\boldsymbol{r}$ は双対問題の制約がどれだけ余っているか（負なら破れていること）を表します。

定理 7.2（LP の弱双対性と安全な下界）✓

(P) の任意の実行可能解 $\boldsymbol{x}$ が $\boldsymbol{x}\le\bar{\boldsymbol{x}}$ を満たすことが分かっているとする。任意の $\boldsymbol{y}\in\mathbb{R}^m$ について $\boldsymbol{r}=\boldsymbol{c}-A^{\mathsf{T}}\boldsymbol{y}$ とおくと

$$p^\ast\ \ge\ \boldsymbol{b}^{\mathsf{T}}\boldsymbol{y}+\sum_{j=1}^n\min(r_j,0)\,\bar x_j.$$

とくに $\boldsymbol{y}$ が (D) の実行可能解なら $\boldsymbol{r}\ge\boldsymbol{0}$ で、$p^\ast\ge\boldsymbol{b}^{\mathsf{T}}\boldsymbol{y}$（弱双対性。この場合 $\bar{\boldsymbol{x}}$ は不要）。

証明（定理 7.2）

実行可能解 $\boldsymbol{x}$ について $\boldsymbol{c}^{\mathsf{T}}\boldsymbol{x}=(A^{\mathsf{T}}\boldsymbol{y}+\boldsymbol{r})^{\mathsf{T}}\boldsymbol{x}=\boldsymbol{y}^{\mathsf{T}}\boldsymbol{b}+\sum_jr_jx_j$ です。$0\le x_j\le\bar x_j$ から $r_jx_j\ge\min(r_j,0)\bar x_j$ で、下限をとれば主張です。

∎

この定理の価値は、$\boldsymbol{y}$ がソルバーの近似解で、わずかに実行不能でもよい点にあります。$\boldsymbol{r}$ と和を区間演算で計算すれば、右辺の下端は $p^\ast$ の厳密な下界です（Neumaier–Shcherbina）。強双対性 $p^\ast=\max(\mathrm{D})$ は、どちらかが実行可能で最適値が有限なら成り立ちますが（Farkas の補題から従います。本章では証明しません）、証明書の検証には弱双対性だけで足ります。

半正定値計画（SDP）では、変数がベクトルから対称行列に、$\boldsymbol{x}\ge\boldsymbol{0}$ が $X\succeq0$ に替わります。内積は $\langle C,X\rangle=\operatorname{tr}(CX)$ です。

補題 7.3（半正定値行列の内積）✓

$S$, $X$ が対称で $X\succeq0$ なら $\langle S,X\rangle\ge\lambda_{\min}(S)\operatorname{tr}X$ である。とくに $S,X\succeq0$ なら $\langle S,X\rangle\ge0$。

証明（補題 7.3）

$P=S-\lambda_{\min}(S)I\succeq0$ です。半正定値行列の平方根 $X^{1/2}$（補題 5.1［行列の分解と標準形］）をとると $\langle P,X\rangle=\operatorname{tr}(PX^{1/2}X^{1/2})=\operatorname{tr}(X^{1/2}PX^{1/2})\ge0$ です（$X^{1/2}PX^{1/2}\succeq0$ の対角成分は非負）。よって $\langle S,X\rangle\ge\lambda_{\min}(S)\langle I,X\rangle=\lambda_{\min}(S)\operatorname{tr}X$ です。

∎

定理 7.4（SDP の安全な下界）✓

主問題 $\min\langle C,X\rangle$ s.t. $\langle A_i,X\rangle=b_i$ ($i=1,\ldots,m$), $X\succeq0$ の任意の実行可能解が $\operatorname{tr}X\le\bar t$ を満たすとする。$\boldsymbol{y}\in\mathbb{R}^m$ について $S=C-\sum_iy_iA_i$ とし、$\lambda_{\min}(S)\ge\underline\lambda$ が分かっているとする。このとき最適値 $p^\ast$ について $p^\ast\ge\boldsymbol{b}^{\mathsf{T}}\boldsymbol{y}+\min(\underline\lambda,0)\,\bar t$ である。

証明（定理 7.4）

実行可能解 $X$ について $\langle C,X\rangle-\boldsymbol{b}^{\mathsf{T}}\boldsymbol{y}=\langle C,X\rangle-\sum_iy_i\langle A_i,X\rangle=\langle S,X\rangle$ です。補題 7.3 から $\langle S,X\rangle\ge\lambda_{\min}(S)\operatorname{tr}X$ で、$\lambda_{\min}(S)\ge0$ なら右辺は $0$ 以上、負なら $\underline\lambda\,\bar t$ 以上です。

∎

$\underline\lambda$ は、$S-\underline\lambda I$ に 定理 4.3 を使って検証します（Jansson–Chaykin–Keil）。SDP の典型的な応用が、多項式の非負性の証明です。

命題 7.5（SOS 証明書）✓

多項式 $p(x)=\sum_{k=0}^{2d}p_kx^k$ と単項式のベクトル $\boldsymbol{z}=(1,x,\ldots,x^d)^{\mathsf{T}}$ について、$p(x)=\boldsymbol{z}^{\mathsf{T}}Q\boldsymbol{z}$ を恒等的に満たす対称行列 $Q\succeq0$ があれば、すべての実数 $x$ で $p(x)\ge0$ である。$Q\succ0$ なら $p(x)\ge\lambda_{\min}(Q)(1+x^2+\cdots+x^{2d})\gt0$ である。

証明（命題 7.5）

$\boldsymbol{z}^{\mathsf{T}}Q\boldsymbol{z}\ge\lambda_{\min}(Q)\|\boldsymbol{z}\|^2$ と $\|\boldsymbol{z}\|^2=1+x^2+\cdots+x^{2d}\ge1$ からです。$Q=LDL^{\mathsf{T}}$ なら $p=\sum_id_i\bigl((L^{\mathsf{T}}\boldsymbol{z})_i\bigr)^2$ と平方和で書けます。

∎

例 7.6（4 次多項式の正値性の証明書）

$p(x)=x^4-4x^3+6x^2-3x+1=(x-1)^4+x$ の正値性は一目ではわかりません（$x\lt0$ で第 2 項が負）。$\boldsymbol{z}=(1,x,x^2)$ について $\boldsymbol{z}^{\mathsf{T}}Q\boldsymbol{z}=p$ となる $Q$ は、係数比較から $q_{00}=1$, $2q_{01}=-3$, $2q_{02}+q_{11}=6$, $2q_{12}=-4$, $q_{22}=1$ で、$t=q_{02}$ が自由に選べます。$t=\tfrac12$ として有理数で LDLᵀ 分解します。

from fractions import Fraction as Fr

def ldlt(A):
    """有理数行列の LDL^T 分解（ピボット選択なし）。ピボットが 0 なら None。"""
    n = len(A)
    L = [[Fr(int(i == j)) for j in range(n)] for i in range(n)]
    d = [Fr(0)] * n
    for j in range(n):
        d[j] = A[j][j] - sum(L[j][k] ** 2 * d[k] for k in range(j))
        if d[j] == 0:
            return None
        for i in range(j + 1, n):
            L[i][j] = (A[i][j] - sum(L[i][k] * L[j][k] * d[k] for k in range(j))) / d[j]
    return L, d

t = Fr(1, 2)
Q = [[Fr(1), Fr(-3, 2), t], [Fr(-3, 2), 6 - 2 * t, Fr(-2)], [t, Fr(-2), Fr(1)]]
L, d = ldlt(Q)
print([str(x) for x in d])                      # ['1', '11/4', '2/11']
print([[str(x) for x in row] for row in L])     # [['1','0','0'], ['-3/2','1','0'], ['1/2','-5/11','1']]

ピボット $1,\tfrac{11}4,\tfrac2{11}$ はすべて正なので、定理 3.1 により $Q\succ0$ です。命題 7.5 の証明の式に代入すると、

$$p(x)=\Bigl(1-\tfrac32x+\tfrac12x^2\Bigr)^2+\tfrac{11}4\Bigl(x-\tfrac5{11}x^2\Bigr)^2+\tfrac2{11}x^4$$

という平方和の表示が得られ、展開すれば誰でも検算できます。実際の応用では $t$ のような自由度が多く、SDP ソルバーで $Q$ を探し、有理数に丸めてから厳密に検証するか、定理 4.3 で検証します。

8. 演習

演習 8.1易

$A=\begin{pmatrix}4&0.1&0\\0.1&2&0.1\\0&0.1&1\end{pmatrix}$ の固有値を、$A=\operatorname{diag}(4,2,1)+E$ と Weyl の不等式で包め。Gershgorin の円板と比べよ。

解答

$E$ は副対角が $0.1$ の三重対角行列で、固有値は $0,\pm0.1\sqrt2$ なので $\|E\|=0.1\sqrt2\approx0.1415$ です。Weyl の不等式から、大きい順に $\lambda_1\in[3.858,4.142]$, $\lambda_2\in[1.858,2.142]$, $\lambda_3\in[0.858,1.142]$ です。Gershgorin の円板（実軸上の区間）は $[3.9,4.1]$, $[1.8,2.2]$, $[0.9,1.1]$ で、第 1・第 3 行ではこちらが狭く、第 2 行では Weyl が狭くなります。この例では円板が互いに交わらないので、Gershgorin の定理からも各円板にちょうど 1 個の固有値があることが言えます。Weyl の評価の利点は、円板が重なる場合でも、各区間に入る固有値の番号（大きい順で何番目か）まで決まることです。

演習 8.2標準

LP $\min\ x_1+x_2$ s.t. $x_1+2x_2=2$, $\boldsymbol{x}\ge\boldsymbol{0}$ について、(1) $\bar{\boldsymbol{x}}$ を求めよ。(2) わずかに実行不能な近似双対解 $y=0.5000001$ から 定理 7.2 の下界を計算し、最適値と比べよ。

解答

(1) $x_1=2-2x_2\le2$, $x_2=(2-x_1)/2\le1$ なので $\bar{\boldsymbol{x}}=(2,1)$ です。(2) $\boldsymbol{r}=\boldsymbol{c}-A^{\mathsf{T}}y=(1-y,\,1-2y)=(0.4999999,\,-0.0000002)$ で、$2y+\min(r_2,0)\cdot\bar x_2=1.0000002-0.0000002=1$ です。実行可能解 $\boldsymbol{x}=(0,1)$ の目的関数値は $1$ なので最適値は $1$ で、下界は最適値に一致します。$y$ をそのまま使った $\boldsymbol{b}^{\mathsf{T}}y=1.0000002$ は下界ではない（最適値を超える）ことに注意してください。

演習 8.3標準

例 7.6 で $t=0$ とした $Q$ は半正定値でないことを、LDLᵀ 分解で示せ。

解答

$Q=\begin{pmatrix}1&-\frac32&0\\-\frac32&6&-2\\0&-2&1\end{pmatrix}$ です。$d_1=1$, $l_{21}=-\tfrac32$, $l_{31}=0$, $d_2=6-\tfrac94=\tfrac{15}4$, $l_{32}=-2/\tfrac{15}4=-\tfrac8{15}$, $d_3=1-\tfrac{64}{225}\cdot\tfrac{15}4=1-\tfrac{16}{15}=-\tfrac1{15}$ です。負のピボットがあるので 定理 3.1 により $Q$ は負の固有値を 1 つもち、半正定値ではありません。同じ多項式でも、Gram 行列の選び方によって証明書になったりならなかったりするので、SDP で $t$ を探す必要があります。

参考文献

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002 — 第 8 章（三角系の誤差解析）、第 10 章（Cholesky 分解）。
• S. M. Rump, “Verification of positive definiteness”, BIT Numerical Mathematics 46 (2006). doi:10.1007/s10543-006-0056-1
• G. H. Golub, C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013 — 第 4 章（Cholesky 分解）、第 8 章（対称固有値問題、Sturm 列と二分法）。
• F. L. Bauer, C. T. Fike, “Norms and exclusion theorems”, Numerische Mathematik 2 (1960), 137–141.
• A. Neumaier, O. Shcherbina, “Safe bounds in linear and mixed-integer linear programming”, Mathematical Programming 99 (2004), 283–296.
• C. Jansson, D. Chaykin, C. Keil, “Rigorous error bounds for the optimal value in semidefinite programming”, SIAM Journal on Numerical Analysis 46 (2007), 180–200.