行列ノルム(作用素ノルム)を理解する — 1ノルム・2ノルム・∞ノルム

同じ連立一次方程式 $\bm{A}\bm{x} = \bm{b}$ を解いているのに、ある問題では入力データの下5桁が違うだけで答えが半分ずれ、別の問題ではびくともしない。この差はどこから来るのでしょうか。あるいは、反復法 $\bm{x}_{k+1} = \bm{M}\bm{x}_k + \bm{c}$ を回したとき、ある $\bm{M}$ では10回で収束し、ある $\bm{M}$ では発散して NaN になる。この境目はどこにあるのでしょうか。

どちらの問いにも、「その行列はベクトルを最大どれだけ引き伸ばすか」というたった一つの数が答えを与えます。それが行列ノルム、とくに作用素ノルム(誘導ノルム)です。行列を「数の表」ではなく「ベクトルを別のベクトルに送る変換」と見たとき、その変換の”最大ゲイン”を測る物差しが作用素ノルムです。

この一つの数が使われる場面は、思っている以上に広く分布しています。

  • 数値線形代数の誤差解析: 条件数 $\kappa(\bm{A}) = \|\bm{A}\|\,\|\bm{A}^{-1}\|$ は、入力の相対誤差が解の相対誤差に何倍に増幅されるかを与えます。有効桁が何桁失われるかが、この数だけで予測できます。
  • 反復法・力学系の安定性: 離散時間システム $\bm{x}_{k+1} = \bm{A}\bm{x}_k$ が原点に収束するかどうかは、$\bm{A}$ のスペクトル半径 $\rho(\bm{A})$ が 1 未満かどうかで決まり、その判定はゲルファントの公式を通じて行列ノルムと直結します。
  • 深層学習のリプシッツ定数: ニューラルネットの各線形層の作用素2ノルム(=最大特異値)を掛け合わせたものが、ネットワーク全体のリプシッツ定数の上界になります。敵対的頑健性やスペクトル正規化GANの理論は、この評価の上に立っています。
  • 制御工学の $\mathcal{H}_\infty$ 設計: 伝達関数行列の各周波数での最大特異値をとり、その上限を小さくするのがロバスト制御の基本方針です。

本記事の内容

  • ベクトルノルムから作用素ノルムが「誘導される」仕組みと、その直感的な意味
  • $\|\bm{A}\|_1$ が最大列和、$\|\bm{A}\|_\infty$ が最大行和になることの証明(最大値を達成する具体的な $\bm{x}$ を構成する)
  • $\|\bm{A}\|_2$ が最大特異値になることの、レイリー商からの導出
  • 劣乗法性、$\rho(\bm{A}) \le \|\bm{A}\|$、ゲルファントの公式、$\bm{A}^k \to \bm{O}$ の必要十分条件
  • 条件数 $\kappa(\bm{A})$ が誤差増幅率を与えることの摂動解析
  • Python による数値実験(単位球の直接最大化、単位円が写る楕円、$\bm{A}^k$ の推移、条件数と誤差の相関)

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

固有値・固有ベクトルと対称行列の直交対角化を知っていると、2ノルムの導出が一段とスムーズになります。

行列ノルムとは — 「最大の引き伸ばし率」を一つの数にする

行列は長さ1のベクトルを向きごとに違う倍率で伸ばし、その最大倍率が作用素ノルムになることを示す概念図

この図が本記事の全体像です。左の単位円上にある「長さ1のベクトル」たちに行列 $\bm{A}$ を掛けると、右のように長さも向きもバラバラに変わり、全体としては楕円になります。青や緑の矢印は伸びたり縮んだりまちまちですが、赤い矢印の向きだけは他のどの向きよりも大きく、6.71 倍に伸びています。この最大の伸び率という1つの数が作用素ノルムであり、以下ではこれを定義し、計算し、応用していきます。

ゴムのシートを両手で引っ張るところを想像してください。縦方向には 3 倍に伸び、横方向には 0.5 倍に縮む、といった具合に、方向によって伸び方が違います。「このシートはどれくらい伸びるシートですか」と聞かれたら、たいていの人は一番よく伸びる方向の伸び率、つまり 3 倍と答えるはずです。最悪ケースで語るのが安全側の評価だからです。

行列 $\bm{A}$ がベクトルに対してやっていることも、これとまったく同じです。$\bm{A}$ は入力ベクトル $\bm{x}$ を受け取り、$\bm{A}\bm{x}$ という別のベクトルを返します。方向によって伸び率 $\|\bm{A}\bm{x}\| / \|\bm{x}\|$ は違います。その最大値を $\bm{A}$ の大きさと定義しよう、というのが作用素ノルムの発想です。

もうひとつの見方は「増幅器のゲイン」です。オーディオアンプに信号を入れると出力が出てきます。入力の大きさを1に正規化したとき出力がどこまで大きくなるか、その最悪値がアンプのゲイン仕様です。行列も多入力多出力の線形増幅器だと思えば、作用素ノルムはそのゲインそのものです。この見方は、深層学習の層ごとの増幅率を評価するときにそのまま使えます。

ベクトルノルムの復習

行列ノルムを定義する前に、土台となるベクトルノルムを確認します。$\bm{x} = (x_1, \dots, x_n)^\top \in \mathbb{R}^n$ に対する $p$ ノルムは

$$ \begin{equation} \|\bm{x}\|_p = \left( \sum_{i=1}^{n} |x_i|^p \right)^{1/p}, \quad p \ge 1 \end{equation} $$

で定義されます。よく使うのは次の3つです。

$$ \|\bm{x}\|_1 = \sum_{i=1}^n |x_i|, \qquad \|\bm{x}\|_2 = \sqrt{\sum_{i=1}^n x_i^2}, \qquad \|\bm{x}\|_\infty = \max_{1 \le i \le n} |x_i| $$

$\|\bm{x}\|_\infty$ は $p \to \infty$ の極限として得られます。実際、最大成分を $|x_m|$ とすると $|x_m| \le \|\bm{x}\|_p \le n^{1/p} |x_m|$ が成り立ち、$n^{1/p} \to 1$ なので挟み撃ちで $\|\bm{x}\|_p \to |x_m|$ となります。

幾何学的には、$\|\bm{x}\|_p \le 1$ という「単位球」の形が $p$ ごとに違います。$p=2$ なら丸い球、$p=1$ なら頂点が座標軸上にある菱形(正軸体)、$p=\infty$ なら軸に平行な立方体です。この単位球の形の違いが、後で見る「最大列和」「最大行和」という公式の差を生みます。先に結論の一部を言ってしまうと、$p=1$ の単位球の頂点は基底ベクトル $\pm\bm{e}_j$、$p=\infty$ の単位球の頂点は成分が $\pm 1$ の符号ベクトルであり、最大の伸び率はまさにその頂点で達成されるのです。

p=1の菱形・p=2の円・p=∞の正方形という単位球の形の違いと、それぞれの頂点の位置

3つの単位球を並べると、同じ「長さ1の集合」でもノルムを変えると形がまったく違うことがわかります。$p=1$ では赤い頂点が座標軸上(基底ベクトル $\pm\bm{e}_j$)にあり、$p=\infty$ では頂点が $(\pm1,\pm1)$ という符号ベクトルにあります。$p=2$ だけは頂点がなく、どの向きも対等です。この「頂点があるかないか」の違いが、後で1ノルム・∞ノルムが単純な足し算の公式になり、2ノルムだけ SVD を要求する理由になります。

行列ノルムが満たすべき公理

行列全体の集合 $\mathbb{R}^{m \times n}$ もベクトル空間なので、そこに入るノルム $\|\cdot\|$ は次の3条件を満たす必要があります。

  1. 非負性・定値性: $\|\bm{A}\| \ge 0$ で、$\|\bm{A}\| = 0 \iff \bm{A} = \bm{O}$
  2. 斉次性: 任意のスカラー $c$ に対し $\|c\bm{A}\| = |c|\,\|\bm{A}\|$
  3. 三角不等式: $\|\bm{A} + \bm{B}\| \le \|\bm{A}\| + \|\bm{B}\|$

さらに、正方行列の積を扱うときには次の性質が非常に重要になります。

  1. 劣乗法性(submultiplicativity): $\|\bm{A}\bm{B}\| \le \|\bm{A}\|\,\|\bm{B}\|$

条件1〜3だけを満たすものは単に「行列ノルム」、4も満たすものを「乗法的行列ノルム」と呼び分けることもあります。作用素ノルムは 4 を自動的に満たします。しかも「$\bm{A}$ を掛けてから $\bm{B}$ を掛ける」を「ゲインの掛け算」で評価できるという、極めて実用的な性質です。これがあるからこそ $\|\bm{A}^k\| \le \|\bm{A}\|^k$ が言えて、べき乗の収束が議論できます。

ここまでで「最大の伸び率を測る」という方針と、ノルムが満たすべき条件が揃いました。次は、この方針を数式として正確に書き下します。

作用素ノルム(誘導ノルム)の定義

「最大の伸び率」をそのまま数式にします。$\bm{A} \in \mathbb{R}^{m \times n}$ に対し、定義域側に $\mathbb{R}^n$ のノルム $\|\cdot\|_p$、値域側に $\mathbb{R}^m$ のノルム $\|\cdot\|_p$ を置いて

$$ \begin{equation} \|\bm{A}\|_p \;=\; \max_{\bm{x} \neq \bm{0}} \frac{\|\bm{A}\bm{x}\|_p}{\|\bm{x}\|_p} \end{equation} $$

と定めます。これをベクトル $p$ ノルムから誘導された作用素ノルム、あるいは単に誘導ノルムと呼びます。「誘導された」というのは、行列ノルムを独立に定義したのではなく、既にあるベクトルノルムから自動的に決まる、という意味です。ベクトルノルムを変えれば行列ノルムも変わります。

3つの同値な表現

この定義式は次の2つと同値です。

$$ \|\bm{A}\|_p = \max_{\|\bm{x}\|_p = 1} \|\bm{A}\bm{x}\|_p = \max_{\|\bm{x}\|_p \le 1} \|\bm{A}\bm{x}\|_p $$

1つ目の同値性は、比 $\|\bm{A}\bm{x}\|_p / \|\bm{x}\|_p$ が $\bm{x}$ のスケールに依存しないことから従います。実際、$\bm{x} \neq \bm{0}$ を $\hat{\bm{x}} = \bm{x}/\|\bm{x}\|_p$ と正規化すると $\|\hat{\bm{x}}\|_p = 1$ であり、線形性より

$$ \frac{\|\bm{A}\bm{x}\|_p}{\|\bm{x}\|_p} = \left\| \bm{A}\frac{\bm{x}}{\|\bm{x}\|_p} \right\|_p = \|\bm{A}\hat{\bm{x}}\|_p $$

となるので、比の全体の集合と単位球面上の $\|\bm{A}\hat{\bm{x}}\|_p$ の集合は一致します。2つ目の同値性は、$\|\bm{x}\|_p \le 1$ の範囲では $\|\bm{A}\bm{x}\|_p \le \|\bm{A}(\bm{x}/\|\bm{x}\|_p)\|_p$($\|\bm{x}\|_p \le 1$ なので割ると大きくなる)となり、最大値は境界 $\|\bm{x}\|_p = 1$ で達成されるからです。

なお、$\sup$ ではなく $\max$ と書けるのは、単位球面 $\{\bm{x} : \|\bm{x}\|_p = 1\}$ が有界閉集合(コンパクト)であり、$\bm{x} \mapsto \|\bm{A}\bm{x}\|_p$ が連続関数だからです。ワイエルシュトラスの最大値定理より、最大値が実際に達成される $\bm{x}$ が存在します。この「達成する $\bm{x}$ が必ずある」という事実は、後で $\|\bm{A}\|_1$ や $\|\bm{A}\|_\infty$ の公式を証明するときに主役になります。公式を導く戦略は毎回同じで、

  1. すべての $\bm{x}$ に対して $\|\bm{A}\bm{x}\|_p \le C \|\bm{x}\|_p$ となる定数 $C$ を見つける(上界)
  2. 等号を達成する具体的な $\bm{x}^\star$ を作る(下界)

この2段構えで $\|\bm{A}\|_p = C$ を結論します。

p=1,2,∞それぞれの単位球面を一周させたときの伸び率の変化と、公式が与える最大値の水平線

$\bm{A} = \begin{pmatrix}3&0\\4&5\end{pmatrix}$ について、単位球面を一周しながら $\|\bm{A}\bm{x}\|_p$ を追いかけた図です。青い曲線は常に赤い破線(公式の値)以下にとどまり、赤い点でちょうど接しています。これが「上界を示し、等号を達成する $\bm{x}^\star$ を作る」という証明の2段構えの図示です。$p=1$ と $p=\infty$ の曲線が折れ線になり最大が尖った点で達成されるのに対し、$p=2$ の曲線はなめらかな山になっており、単位球の形の違いがそのまま伸び率の形に現れています。

定義から直ちに従う性質

(a) 基本不等式 任意の $\bm{x}$ に対して

$$ \begin{equation} \|\bm{A}\bm{x}\|_p \le \|\bm{A}\|_p \, \|\bm{x}\|_p \end{equation} $$

が成り立ちます。$\bm{x} = \bm{0}$ なら両辺 0 で成立、$\bm{x} \neq \bm{0}$ なら定義の最大値の性質からそのまま従います。実務ではこの形で使うことが圧倒的に多く、「行列を掛けると、ベクトルの長さは高々ノルム倍にしかならない」という安心を与えてくれます。

(b) 単位行列のノルムは 1 $\|\bm{I}\|_p = \max_{\|\bm{x}\|_p=1} \|\bm{x}\|_p = 1$。当たり前に見えますが、これは誘導ノルムを特徴づける重要な性質です。後で見るフロベニウスノルムは $\|\bm{I}\|_F = \sqrt{n}$ なので、フロベニウスノルムはどのベクトルノルムからも誘導されないことがここからわかります。

(c) 劣乗法性 $\bm{A} \in \mathbb{R}^{m\times n}$, $\bm{B} \in \mathbb{R}^{n \times \ell}$ に対し

$$ \|\bm{A}\bm{B}\|_p \le \|\bm{A}\|_p \|\bm{B}\|_p $$

証明は基本不等式を2回使うだけです。$\|\bm{x}\|_p = 1$ なる任意の $\bm{x}$ に対して

$$ \|\bm{A}\bm{B}\bm{x}\|_p = \|\bm{A}(\bm{B}\bm{x})\|_p \le \|\bm{A}\|_p \|\bm{B}\bm{x}\|_p \le \|\bm{A}\|_p \|\bm{B}\|_p \|\bm{x}\|_p = \|\bm{A}\|_p \|\bm{B}\|_p $$

最初の不等号で $\bm{B}\bm{x}$ をひとかたまりのベクトルと見て基本不等式を適用し、次の不等号で $\bm{B}\bm{x}$ の中身に再び基本不等式を適用しました。左辺の最大をとれば結論を得ます。

とくに $\bm{B} = \bm{A}$ と繰り返せば

$$ \|\bm{A}^k\|_p \le \|\bm{A}\|_p^k $$

が得られます。これが「$\|\bm{A}\| < 1$ なら $\bm{A}^k \to \bm{O}$」という収束判定の出発点です。

定義と一般論はここまでです。しかし定義式は最大化問題であり、このままでは計算できません。次からは、$p = 1, \infty, 2$ の各場合に対して、この最大化問題を閉じた式で解いていきます。

$\|\bm{A}\|_1$ は最大列和

まず $p = 1$ から始めます。$\ell_1$ の単位球は菱形で、その頂点は $\pm\bm{e}_1, \dots, \pm\bm{e}_n$ という基底ベクトルです。一方 $\bm{A}\bm{e}_j$ は $\bm{A}$ の第 $j$ 列そのものです。ということは「一番大きい列を選ぶだけ」で答えになりそうです。この直感が正しいことを、上界と達成の2段構えで示します。

主張: $\bm{A} = (a_{ij}) \in \mathbb{R}^{m\times n}$ に対して

$$ \begin{equation} \|\bm{A}\|_1 = \max_{1 \le j \le n} \sum_{i=1}^{m} |a_{ij}| \end{equation} $$

すなわち「各列の絶対値の和のうち最大のもの(最大列和)」です。

ステップ1(上界) 任意の $\bm{x} = (x_1,\dots,x_n)^\top$ に対して $\bm{A}\bm{x}$ の第 $i$ 成分は $\sum_j a_{ij}x_j$ なので

$$ \|\bm{A}\bm{x}\|_1 = \sum_{i=1}^m \left| \sum_{j=1}^n a_{ij} x_j \right| $$

三角不等式で絶対値を内側に入れます。

$$ \|\bm{A}\bm{x}\|_1 \le \sum_{i=1}^m \sum_{j=1}^n |a_{ij}| |x_j| $$

ここで有限和なので和の順序を入れ替えられます。$j$ を外に出すと

$$ \|\bm{A}\bm{x}\|_1 \le \sum_{j=1}^n |x_j| \left( \sum_{i=1}^m |a_{ij}| \right) = \sum_{j=1}^n |x_j| \, c_j $$

と書けます。ここで $c_j = \sum_i |a_{ij}|$ は第 $j$ 列の絶対値和です。$c_j \le \max_k c_k$ で上から押さえると

$$ \|\bm{A}\bm{x}\|_1 \le \left( \max_{k} c_k \right) \sum_{j=1}^n |x_j| = \left( \max_{k} c_k \right) \|\bm{x}\|_1 $$

したがって $\|\bm{A}\|_1 \le \max_k c_k$ が示せました。

ステップ2(達成) 最大列和を与える列番号を $j^\star = \arg\max_j c_j$ とし、$\bm{x}^\star = \bm{e}_{j^\star}$(第 $j^\star$ 成分だけ 1、他は 0)をとります。$\|\bm{e}_{j^\star}\|_1 = 1$ であり、$\bm{A}\bm{e}_{j^\star}$ は $\bm{A}$ の第 $j^\star$ 列なので

$$ \|\bm{A}\bm{e}_{j^\star}\|_1 = \sum_{i=1}^m |a_{i j^\star}| = c_{j^\star} = \max_k c_k $$

上界と一致するので、$\|\bm{A}\|_1 = \max_k c_k$ が確定しました。$\blacksquare$

証明のポイントは「$\ell_1$ ノルムでは、重みを1点に集中させるのが最も効率がよい」ことです。$\|\bm{x}\|_1 = 1$ という予算のもとで、$\bm{x}$ の重みを複数の列に分散させると、出力は各列の凸結合になってしまい、最良の列の値を超えられません。予算をすべて最良の列に注ぎ込むのが最適解です。この「予算を1点に集中」という構造は、$\ell_1$ 正則化がスパース解を生む理由とも通じています。

行列の各列の絶対値和を並べたヒートマップと棒グラフ、および基底ベクトルが最大値を達成することを示すヒストグラム

左の表で列ごとに絶対値を足すと $(6, 7, 3)$ となり、その最大値 7 がそのまま $\|\bm{A}\|_1$ です(中央の棒グラフの赤い柱)。右のヒストグラムは、$\|\bm{x}\|_1 = 1$ を満たすように重みを複数成分に分散させた3000通りの $\bm{x}$ に対する $\|\bm{A}\bm{x}\|_1$ の分布で、山全体が 7 の縦線より左に収まっています。重みを分散させると出力は各列の凸結合になり、最良の列の値 7 を超えられないことが実測でも確認でき、$\bm{x}^\star = \bm{e}_2$ に全振りするのが最適だとわかります。

上界の評価で一度も無駄がないことを確認できたので、次は同じ論法を $p=\infty$ に適用します。今度は単位球が立方体なので、最適な $\bm{x}$ の形も変わります。

$\|\bm{A}\|_\infty$ は最大行和

$\ell_\infty$ の単位球は各成分が $[-1,1]$ の立方体です。頂点は成分がすべて $\pm 1$ の符号ベクトルで、$2^n$ 個あります。出力の $\infty$ ノルムは「$\bm{A}\bm{x}$ の成分の最大絶対値」なので、どれか1行を狙い撃ちして最大化すればよいという構造になります。狙った行の内積を最大にするには、その行の各成分と同じ符号を $\bm{x}$ に持たせればよいはずです。

主張:

$$ \begin{equation} \|\bm{A}\|_\infty = \max_{1 \le i \le m} \sum_{j=1}^{n} |a_{ij}| \end{equation} $$

すなわち「最大行和」です。

ステップ1(上界) $\|\bm{x}\|_\infty \le 1$ とすると、すべての $j$ で $|x_j| \le 1$ です。$\bm{A}\bm{x}$ の第 $i$ 成分は

$$ \left| \sum_{j=1}^n a_{ij}x_j \right| \le \sum_{j=1}^n |a_{ij}||x_j| \le \sum_{j=1}^n |a_{ij}| = r_i $$

と評価できます。1つ目の不等号は三角不等式、2つ目は $|x_j| \le 1$ を使いました。$r_i$ は第 $i$ 行の絶対値和です。すべての成分がこの形で押さえられるので、最大成分も

$$ \|\bm{A}\bm{x}\|_\infty = \max_i \left| \sum_j a_{ij}x_j \right| \le \max_i r_i $$

となり、$\|\bm{A}\|_\infty \le \max_i r_i$ を得ます。

ステップ2(達成) 最大行和を与える行番号を $i^\star = \arg\max_i r_i$ とし、その行の符号を並べたベクトル

$$ x^\star_j = \operatorname{sign}(a_{i^\star j}) = \begin{cases} 1 & (a_{i^\star j} \ge 0) \\ -1 & (a_{i^\star j} < 0)\end{cases} $$

をとります。$\|\bm{x}^\star\|_\infty = 1$ です($a_{i^\star j}$ がすべて 0 なら $r_{i^\star}=0$ で $\bm{A}=\bm{O}$ の自明な場合なので除外します)。このとき第 $i^\star$ 成分は

$$ \sum_{j=1}^n a_{i^\star j} x^\star_j = \sum_{j=1}^n a_{i^\star j} \operatorname{sign}(a_{i^\star j}) = \sum_{j=1}^n |a_{i^\star j}| = r_{i^\star} $$

となります。符号を合わせたので、すべての項が正の寄与をして打ち消し合いが起きない、というのがこの構成の狙いです。他の成分がこれより大きくなることはない(ステップ1より $\le \max_i r_i = r_{i^\star}$)ので、$\|\bm{A}\bm{x}^\star\|_\infty = r_{i^\star}$。上界と一致し、$\|\bm{A}\|_\infty = \max_i r_i$ が示されました。$\blacksquare$

行列の各行の絶対値和と、符号ベクトルを入力したときだけ最大行和8に到達することを示す棒グラフ

行ごとに絶対値を足すと $(3, 5, 8)$ で、最大行和 8 が $\|\bm{A}\|_\infty$ になります。右のグラフは、第3行の符号を並べた $\bm{x}^\star = (-1, 1, 1)$ と、ランダムに選んだ3つの $\bm{x} \in [-1,1]^3$ について $|(\bm{A}\bm{x})_i|$ を比較したものです。ランダムな入力では項どうしが打ち消し合って第3成分でも 2〜4 程度にしかならないのに対し、符号を揃えた $\bm{x}^\star$ だけが打ち消し合いをゼロにして 8 に到達していることが読み取れます。

最大列和の公式と最大行和の公式を見比べると、行と列がちょうど入れ替わっています。実際、$|a_{ij}|$ の行和と列和は転置で入れ替わるので

$$ \begin{equation} \|\bm{A}^\top\|_1 = \|\bm{A}\|_\infty, \qquad \|\bm{A}^\top\|_\infty = \|\bm{A}\|_1 \end{equation} $$

という双対関係が成り立ちます。これは $\ell_1$ と $\ell_\infty$ が互いに双対ノルム($1/p + 1/q = 1$ で $p=1, q=\infty$)であることの表れです。

1ノルムと∞ノルムは、いずれも「絶対値を足すだけ」で $O(mn)$ で計算できました。ところが最も自然に感じる2ノルム(ユークリッド距離での伸び率)は、そう簡単には行きません。次はそこに踏み込みます。

$\|\bm{A}\|_2$ は最大特異値 — レイリー商からの導出

$p=2$ の単位球は丸い球で、頂点というものがありません。したがって「頂点を調べれば終わり」という手が使えず、最大化問題を真面目に解く必要があります。しかし幸いなことに、$\|\bm{A}\bm{x}\|_2^2 = \bm{x}^\top \bm{A}^\top \bm{A} \bm{x}$ という二次形式に化けるので、対称行列の二次形式の最大化という古典的な問題に帰着します。その答えを与えるのがレイリー商の理論です。

二次形式への変形

$\|\bm{A}\bm{x}\|_2^2 = (\bm{A}\bm{x})^\top(\bm{A}\bm{x})$ を展開します。転置の性質 $(\bm{A}\bm{x})^\top = \bm{x}^\top\bm{A}^\top$ を使うと

$$ \|\bm{A}\bm{x}\|_2^2 = \bm{x}^\top \bm{A}^\top \bm{A} \bm{x} $$

ここで $\bm{G} = \bm{A}^\top\bm{A} \in \mathbb{R}^{n \times n}$ とおきます(複素行列なら $\bm{G} = \bm{A}^{H}\bm{A}$)。$\bm{G}$ には2つの良い性質があります。

(i) 対称性: $\bm{G}^\top = (\bm{A}^\top\bm{A})^\top = \bm{A}^\top(\bm{A}^\top)^\top = \bm{A}^\top\bm{A} = \bm{G}$。

(ii) 半正定値性: 任意の $\bm{x}$ に対して $\bm{x}^\top\bm{G}\bm{x} = \|\bm{A}\bm{x}\|_2^2 \ge 0$。

対称行列はスペクトル定理により直交行列で対角化できます。すなわち、正規直交な固有ベクトル $\bm{v}_1,\dots,\bm{v}_n$ と実固有値 $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_n$ が存在して

$$ \bm{G}\bm{v}_i = \lambda_i \bm{v}_i, \qquad \bm{v}_i^\top\bm{v}_j = \delta_{ij} $$

が成り立ちます。半正定値性から $\lambda_i = \bm{v}_i^\top\bm{G}\bm{v}_i \ge 0$、つまり固有値はすべて 0 以上です。

レイリー商の最大値

任意の $\bm{x}$ を固有ベクトルの線形結合 $\bm{x} = \sum_{i=1}^n c_i \bm{v}_i$ と展開します(正規直交基底なので必ずできます)。まず分母を計算すると、直交性から交差項が消えて

$$ \|\bm{x}\|_2^2 = \bm{x}^\top\bm{x} = \sum_{i}\sum_{j} c_i c_j \bm{v}_i^\top\bm{v}_j = \sum_{i=1}^n c_i^2 $$

次に分子です。$\bm{G}\bm{x} = \sum_i c_i \bm{G}\bm{v}_i = \sum_i c_i \lambda_i \bm{v}_i$ を使い、再び直交性で交差項を落とすと

$$ \bm{x}^\top\bm{G}\bm{x} = \left(\sum_j c_j \bm{v}_j\right)^\top \left(\sum_i c_i\lambda_i\bm{v}_i\right) = \sum_{i=1}^n \lambda_i c_i^2 $$

したがってレイリー商

$$ R(\bm{x}) = \frac{\bm{x}^\top\bm{G}\bm{x}}{\bm{x}^\top\bm{x}} = \frac{\sum_i \lambda_i c_i^2}{\sum_i c_i^2} $$

と書けます。この形は「重み $c_i^2 \ge 0$ による $\lambda_i$ の加重平均」です。加重平均は最大値 $\lambda_1$ を超えられません。

$$ R(\bm{x}) = \sum_i \lambda_i \frac{c_i^2}{\sum_k c_k^2} \le \lambda_1 \sum_i \frac{c_i^2}{\sum_k c_k^2} = \lambda_1 $$

しかも $\bm{x} = \bm{v}_1$(すなわち $c_1 = 1$、他は 0)とすれば $R(\bm{v}_1) = \lambda_1$ で等号が成立します。よって

$$ \max_{\bm{x}\neq\bm{0}} R(\bm{x}) = \lambda_1 = \lambda_{\max}(\bm{A}^\top\bm{A}) $$

レイリー商が単位ベクトルの向きに応じて最大固有値と最小固有値の間を動き、固有値の加重平均になっていることを示す図

左の図は $\bm{G} = \bm{A}^\top\bm{A}$($\bm{A} = \begin{pmatrix}3&0\\4&5\end{pmatrix}$ の場合)のレイリー商を、単位ベクトルの向きを 0° から 180° まで回しながらプロットしたものです。曲線は $\lambda_2 = 5$ と $\lambda_1 = 45$ の2本の破線の間だけを動き、$\bm{v}_1$ の向き(45°)で上の破線に接して最大、$\bm{v}_2$ の向き(135°)で下の破線に接して最小になります。右の図は同じことを重みの言葉で示したもので、$\bm{x}$ を $\bm{v}_1$ から離すほど $\lambda_1$ への重み $c_1^2$(赤)が減って $\lambda_2$ への重み(橙)に移り、加重平均である $R(\bm{x})$ が単調に下がっていきます。重みが $\lambda_1$ に全振りされたときだけ最大値 $\lambda_1$ が達成されるというのが、この節の証明の中身です。

特異値との接続

以上より

$$ \|\bm{A}\|_2^2 = \max_{\bm{x}\neq\bm{0}} \frac{\bm{x}^\top\bm{A}^\top\bm{A}\bm{x}}{\bm{x}^\top\bm{x}} = \lambda_{\max}(\bm{A}^\top\bm{A}) $$

両辺の平方根をとって

$$ \begin{equation} \|\bm{A}\|_2 = \sqrt{\lambda_{\max}(\bm{A}^\top\bm{A})} = \sigma_{\max}(\bm{A}) = \sigma_1 \end{equation} $$

を得ます。ここで特異値 $\sigma_i = \sqrt{\lambda_i}$ は $\bm{A}^\top\bm{A}$ の固有値の平方根として定義されるものであり、特異値分解 $\bm{A} = \bm{U}\bm{\Sigma}\bm{V}^\top$ の対角成分に並ぶ量です。この理由から $\|\bm{A}\|_2$ はスペクトルノルムとも呼ばれます。最大化を達成する $\bm{x}$ は $\bm{v}_1$、つまり最大特異値に対応する右特異ベクトルであり、そのとき出力は $\bm{A}\bm{v}_1 = \sigma_1 \bm{u}_1$ という左特異ベクトル方向を向きます。

この最大特異値の公式の代償は計算量です。1ノルム・∞ノルムが $O(mn)$ で済むのに対し、2ノルムは特異値分解(またはべき乗法)が必要で $O(mn\min(m,n))$ 程度かかります。大規模行列で「とりあえずノルムを押さえたい」だけなら、後述の $\|\bm{A}\|_2 \le \sqrt{\|\bm{A}\|_1\|\bm{A}\|_\infty}$ のような安価な上界が実用的です。

なお、$\bm{A}^\top\bm{A}$ と $\bm{A}\bm{A}^\top$ は非零固有値を共有するので、$\|\bm{A}\|_2 = \|\bm{A}^\top\|_2$ が成り立ちます。転置の関係式で1ノルムと∞ノルムが入れ替わったのに対し、2ノルムは転置で不変です。

3つの主要な誘導ノルムが出揃いました。ここで一度立ち止まって、「誘導ノルムではないが実務で頻繁に使われるノルム」との違いを確認しておきます。

誘導ノルムではないもの — フロベニウスノルム

行列を単に $mn$ 次元のベクトルだと思って、成分の二乗和の平方根をとったものがフロベニウスノルムです。

$$ \|\bm{A}\|_F = \sqrt{\sum_{i=1}^m\sum_{j=1}^n a_{ij}^2} = \sqrt{\operatorname{tr}(\bm{A}^\top\bm{A})} = \sqrt{\sum_{i=1}^{r}\sigma_i^2} $$

最後の等号は、トレースが固有値の和であること($\operatorname{tr}(\bm{A}^\top\bm{A}) = \sum_i\lambda_i = \sum_i \sigma_i^2$)から従います。$r = \operatorname{rank}(\bm{A})$ です。

フロベニウスノルムは計算が極めて安く、微分も素直なので、機械学習の重み減衰(weight decay)や行列補完の目的関数に頻出します。しかし $\|\bm{I}_n\|_F = \sqrt{n}$ なので、性質 (b) よりどのベクトルノルムからも誘導されません。それでも劣乗法性 $\|\bm{A}\bm{B}\|_F \le \|\bm{A}\|_F\|\bm{B}\|_F$ は満たします(各成分にコーシー・シュワルツを適用すれば示せます)。

特異値表現から、2ノルムとの関係が直ちに読み取れます。

$$ \|\bm{A}\|_2 = \sigma_1 \le \sqrt{\sigma_1^2 + \dots + \sigma_r^2} = \|\bm{A}\|_F \le \sqrt{r\,\sigma_1^2} = \sqrt{r}\,\|\bm{A}\|_2 $$

左の不等号は「最大の1項 ≤ 全部の和」、右は「各 $\sigma_i \le \sigma_1$」から従います。つまり両者は $\sqrt{r}$ 倍以内で一致します。ランク1行列では $r=1$ なので $\|\bm{A}\|_2 = \|\bm{A}\|_F$ となり、完全に一致します。

「$\|\bm{A}\|_F$ を小さくすれば $\|\bm{A}\|_2$ も小さくなる」ので、リプシッツ定数を抑えたいときにフロベニウスノルムで代用するのは安全側です。ただし $\sqrt{r}$ 倍の緩みがあり、高ランクの層では過度に強い制約になってしまいます。スペクトル正規化がわざわざ最大特異値を推定しに行くのは、この緩みを避けるためです。

特異値の棒グラフに2ノルム・フロベニウスノルム・√r倍の線を重ねた図と、比がランクとともに変化する図

左は $6\times6$ のランダム行列の特異値です。$\|\bm{A}\|_2 = 4.029$ は最大の棒1本の高さ、$\|\bm{A}\|_F = 5.602$ は全部の棒を二乗和で束ねた高さで、確かに $\|\bm{A}\|_2 \le \|\bm{A}\|_F \le \sqrt{6}\|\bm{A}\|_2 = 9.868$ に収まっています。右はランクを 1 から 20 まで変えて比 $\|\bm{A}\|_F/\|\bm{A}\|_2$ を測ったもので、ランク1では実測が厳密に 1(両ノルムが完全一致)、ランクが上がるほど比が増えますが上限 $\sqrt{r}$ には届きません。理論の上下界がどちらも破られず、しかも実測は上限よりかなり手前にいることが確認できます。

複数のノルムが出てきたところで、それらの間にどんな不等式が成り立つのかを整理しておきましょう。ノルム同士の橋渡しができると、計算しやすいノルムで計算しにくいノルムを評価できます。

ノルム同士の関係と有用な不等式

有限次元では、すべてのノルムは互いに同値(定数倍で挟める)です。行列ノルムでも具体的な定数が知られています。$\bm{A} \in \mathbb{R}^{m\times n}$ に対して代表的なものを挙げます。

$$ \frac{1}{\sqrt{n}}\|\bm{A}\|_\infty \le \|\bm{A}\|_2 \le \sqrt{m}\,\|\bm{A}\|_\infty, \qquad \frac{1}{\sqrt{m}}\|\bm{A}\|_1 \le \|\bm{A}\|_2 \le \sqrt{n}\,\|\bm{A}\|_1 $$

とくに実用上よく使うのが次の幾何平均による評価です。

$$ \begin{equation} \|\bm{A}\|_2 \le \sqrt{\|\bm{A}\|_1 \, \|\bm{A}\|_\infty} \end{equation} $$

証明はこうです。まず後で示す $\rho(\bm{M}) \le \|\bm{M}\|$(任意の誘導ノルム)を $\bm{M} = \bm{A}^\top\bm{A}$ に、$\infty$ ノルムで適用します。

$$ \|\bm{A}\|_2^2 = \lambda_{\max}(\bm{A}^\top\bm{A}) = \rho(\bm{A}^\top\bm{A}) \le \|\bm{A}^\top\bm{A}\|_\infty $$

ここで劣乗法性を使って積を分解し、さらに転置の関係式を使うと

$$ \|\bm{A}^\top\bm{A}\|_\infty \le \|\bm{A}^\top\|_\infty \|\bm{A}\|_\infty = \|\bm{A}\|_1 \|\bm{A}\|_\infty $$

両辺の平方根をとれば幾何平均による評価が得られます。$\blacksquare$

この評価の実用的な価値は、$O(mn)$ で計算できる2つの量から、計算コストの高い $\|\bm{A}\|_2$ の上界が即座に得られることです。巨大な疎行列に対して安全側の評価がほしいとき、SVD を回さずに済みます。

600個のランダム行列で2ノルムと幾何平均の上界を比較した散布図と、その比のヒストグラム

サイズも大きさもばらばらなランダム行列 600 個で、横軸に $\sqrt{\|\bm{A}\|_1\|\bm{A}\|_\infty}$、縦軸に真の $\|\bm{A}\|_2$ をとった散布図です。点は1つ残らず $y=x$ の破線より下側にあり、不等式が破られていないことが見て取れます。右のヒストグラムを見ると比の中央値は 0.666、最大でも 0.977 で、上界は「常に安全側だが、行列によっては 1.5 倍ほど過大評価する」程度の実用的な緩さだとわかります。

さて、ここまでは「行列がベクトルをどれだけ伸ばすか」という話でした。一方、行列の”大きさ”を測るもう一つの自然な量に、固有値の最大絶対値があります。この2つはどう関係するのでしょうか。

スペクトル半径 $\rho(\bm{A})$ との関係

正方行列 $\bm{A} \in \mathbb{C}^{n\times n}$ の固有値を $\lambda_1,\dots,\lambda_n$(重複込み)とするとき、

$$ \rho(\bm{A}) = \max_{1\le i \le n} |\lambda_i| $$

スペクトル半径と呼びます。複素平面上で、すべての固有値を含む原点中心の最小円の半径です。

固有ベクトル方向では $\bm{A}\bm{v} = \lambda\bm{v}$ なので伸び率はちょうど $|\lambda|$ です。作用素ノルムは全方向の最大伸び率でしたから、固有ベクトル方向だけを見たものより大きいはずです。これを正確に述べます。

$\rho(\bm{A}) \le \|\bm{A}\|$ の証明

主張: 任意の誘導ノルム $\|\cdot\|$ に対して $\rho(\bm{A}) \le \|\bm{A}\|$。

証明: $\lambda$ を $|\lambda| = \rho(\bm{A})$ なる固有値、$\bm{v} \neq \bm{0}$ を対応する固有ベクトルとします。実行列でも固有値は複素になり得るので、複素ベクトル空間で考えます($\mathbb{C}^n$ 上のベクトルノルムと、そこから誘導される行列ノルムを使います)。$\bm{A}\bm{v} = \lambda\bm{v}$ の両辺のノルムをとると、斉次性より

$$ \|\bm{A}\bm{v}\| = \|\lambda\bm{v}\| = |\lambda|\,\|\bm{v}\| $$

一方、基本不等式から $\|\bm{A}\bm{v}\| \le \|\bm{A}\|\,\|\bm{v}\|$。この2つを合わせると

$$ |\lambda|\,\|\bm{v}\| \le \|\bm{A}\|\,\|\bm{v}\| $$

$\bm{v}\neq\bm{0}$ より $\|\bm{v}\| > 0$ なので、両辺を $\|\bm{v}\|$ で割って $|\lambda| \le \|\bm{A}\|$。$\lambda$ は最大絶対値の固有値だったので $\rho(\bm{A}) \le \|\bm{A}\|$。$\blacksquare$

この不等式は一般に等号ではありません。対称行列(より一般に正規行列 $\bm{A}^H\bm{A} = \bm{A}\bm{A}^H$)では $\|\bm{A}\|_2 = \rho(\bm{A})$ となりますが、非正規行列では大きな差が出ます。たとえば

$$ \bm{A} = \begin{pmatrix} 3 & 0 \\ 4 & 5 \end{pmatrix} $$

の固有値は 3 と 5 なので $\rho(\bm{A}) = 5$ ですが、$\bm{A}^\top\bm{A} = \begin{pmatrix}25 & 20\\ 20 & 25\end{pmatrix}$ の固有値は 45 と 5 なので $\|\bm{A}\|_2 = \sqrt{45} = 3\sqrt{5} \approx 6.708$ です。固有ベクトル方向 $(0,1)^\top$ では 5 倍にしかなりませんが、$(1,1)^\top/\sqrt{2}$ の方向では 6.708 倍に伸びます。固有ベクトルが直交していないと、非固有ベクトル方向で”より伸びる”ことがあるというのが差の正体です。

複素平面上の固有値とスペクトル半径の円・作用素ノルムの円を、非正規行列と対称行列で比較した図

複素平面に固有値(赤点)を打ち、半径 $\rho(\bm{A})$ の赤い破線円と半径 $\|\bm{A}\|_2$ の青い円を重ねた図です。左の非正規行列では 5.000 と 6.708 の2つの円のあいだにはっきり隙間があり、$\rho < \|\bm{A}\|_2$ が図として見えています。右の対称行列では2つの円が完全に重なり、$\rho = \|\bm{A}\|_2 = 6.236$ となって隙間が消えます。固有値は必ず青い円の内側にある($\rho \le \|\bm{A}\|$)が、非正規行列ではその内側にかなり余裕を残す、というのがこの不等式の実像です。

スペクトル半径はノルムではない

$\rho$ が行列ノルムの公理を満たさないことも確認しておきます。

$$ \bm{A} = \begin{pmatrix} 0 & 1 \\ 0 & 0\end{pmatrix}, \qquad \bm{B} = \begin{pmatrix} 0 & 0 \\ 1 & 0\end{pmatrix} $$

はどちらもべき零行列で $\rho(\bm{A}) = \rho(\bm{B}) = 0$ ですが、$\bm{A} \neq \bm{O}$ なので定値性が破れています。さらに $\bm{A}+\bm{B} = \begin{pmatrix}0&1\\1&0\end{pmatrix}$ の固有値は $\pm 1$ で $\rho(\bm{A}+\bm{B}) = 1 > 0 = \rho(\bm{A}) + \rho(\bm{B})$ となり、三角不等式も破れます。$\rho$ は「大きさ」らしく見えて、実はノルムではないのです。

ゲルファントの公式

$\rho(\bm{A}) \le \|\bm{A}\|$ には隙間がありますが、べき乗をとって $k$ 乗根を計算すると、その隙間が消えていきます

$$ \begin{equation} \rho(\bm{A}) = \lim_{k\to\infty} \|\bm{A}^k\|^{1/k} \end{equation} $$

これがゲルファントの公式(スペクトル半径公式)で、任意の劣乗法的行列ノルムで成立します。ノルムの選び方によらず同じ極限になる、というのが驚くべき点です。

証明の骨格を2つの向きに分けて説明します。

下からの評価 $\bm{A}\bm{v} = \lambda\bm{v}$ を繰り返すと $\bm{A}^k\bm{v} = \lambda^k\bm{v}$ なので、$\lambda^k$ は $\bm{A}^k$ の固有値です。よって $\rho(\bm{A})^k = \rho(\bm{A}^k) \le \|\bm{A}^k\|$、両辺の $k$ 乗根をとって

$$ \rho(\bm{A}) \le \|\bm{A}^k\|^{1/k} \quad (\forall k) $$

つまり数列 $\|\bm{A}^k\|^{1/k}$ は常に $\rho(\bm{A})$ 以上です。

上からの評価 任意の $\varepsilon > 0$ に対し、$\tilde{\bm{A}} = \bm{A}/(\rho(\bm{A})+\varepsilon)$ とおくと $\rho(\tilde{\bm{A}}) < 1$ です。後述の定理より $\tilde{\bm{A}}^k \to \bm{O}$、よって十分大きな $k$ で $\|\tilde{\bm{A}}^k\| \le 1$、すなわち $\|\bm{A}^k\| \le (\rho(\bm{A})+\varepsilon)^k$。$k$ 乗根をとれば $\|\bm{A}^k\|^{1/k} \le \rho(\bm{A})+\varepsilon$ です。$\varepsilon$ は任意なので、上下から挟んでゲルファントの公式が従います。$\blacksquare$

さらに、次の”逆向き”の事実も重要です。任意の $\varepsilon>0$ に対して、$\|\bm{A}\|_\star \le \rho(\bm{A}) + \varepsilon$ となる誘導ノルム $\|\cdot\|_\star$ が存在する。 構成は簡単で、$\bm{A}$ をジョルダン標準形 $\bm{J}$ に相似変換し、$\bm{D}_t = \operatorname{diag}(1, t, t^2, \dots, t^{n-1})$ で $\bm{D}_t^{-1}\bm{J}\bm{D}_t$ を作ると、ジョルダンブロックの超対角成分が $t$ になります。$t$ を小さくすれば行和は $\rho(\bm{A})+\varepsilon$ 以下にできます。そして $\|\bm{X}\|_\star := \|\bm{S}^{-1}\bm{X}\bm{S}\|_\infty$ は、相似変換された座標系でのベクトル $\infty$ ノルムから誘導される立派な作用素ノルムです。

具体例で見てみます。$\bm{J} = \begin{pmatrix}0.9 & 10 \\ 0 & 0.9\end{pmatrix}$ は $\rho = 0.9$ ですが $\|\bm{J}\|_\infty = 10.9$ です。$\bm{D} = \operatorname{diag}(1, 0.001)$ で相似変換すると

$$ \bm{D}^{-1}\bm{J}\bm{D} = \begin{pmatrix}0.9 & 0.01 \\ 0 & 0.9\end{pmatrix} $$

となり、$\infty$ ノルムは 0.91 まで下がります。「非正規性による見かけ上のゲインは、座標を取り替えれば小さくできる」というのがこの構成の意味です。

相似変換のパラメータtを小さくするとジョルダン型行列の∞ノルムがスペクトル半径0.9に近づく様子

$\bm{D}_t = \operatorname{diag}(1, t)$ の $t$ を 1 から 0.001 まで小さくしていくと、$\|\bm{D}_t^{-1}\bm{J}\bm{D}_t\|_\infty$ は $0.9 + 10t$ に沿って 10.900 → 1.900 → 1.000 → 0.910 と減り、赤い破線 $\rho = 0.9$ に漸近します。右の棒グラフでは $t \le 0.001$ で初めてノルムが 1 を下回り(緑)、この座標系を選べば「ノルムが1未満だから収束する」という単純な論法が使えるようになることがわかります。元の座標のままでは $\|\bm{J}\|_\infty = 10.9 > 1$ で収束を何も言えなかったことと対比すると、ノルムの選び方が判定力を左右することが実感できます。

$\bm{A}^k \to \bm{O}$ の必要十分条件

これらを使うと、力学系の安定性を完全に特徴づけられます。

定理: $\bm{A}^k \to \bm{O}$ ($k\to\infty$) $\iff$ $\rho(\bm{A}) < 1$。

($\Leftarrow$) $\rho(\bm{A}) < 1$ なら、$\varepsilon = (1-\rho(\bm{A}))/2 > 0$ をとると、上で構成した誘導ノルムで $\|\bm{A}\|_\star \le \rho(\bm{A}) + \varepsilon = (1+\rho(\bm{A}))/2 =: q < 1$ となります。劣乗法性から $\|\bm{A}^k\|_\star \le q^k \to 0$。有限次元ではすべてのノルムが同値なので、成分ごとにも 0 に収束します。

($\Rightarrow$) 対偶を示します。$\rho(\bm{A}) \ge 1$ なら、$|\lambda| \ge 1$ なる固有値と固有ベクトル $\bm{v}\neq\bm{0}$ があり、$\bm{A}^k\bm{v} = \lambda^k\bm{v}$ の大きさは $|\lambda|^k\|\bm{v}\| \ge \|\bm{v}\| > 0$ なので 0 に収束しません。$\blacksquare$

ここで注意したいのは、$\|\bm{A}\| < 1$ は収束の十分条件にすぎないことです。$\rho(\bm{A}) < 1$ でも $\|\bm{A}\|_2 \gg 1$ の場合があり、そのとき $\|\bm{A}^k\|$ は最初のうち増大してから減衰します。これを過渡的増幅(transient growth)と呼び、流体の遷移や非正規な制御系で実務的な問題になります。「漸近的には安定だが、短期的には数十倍に増幅される」システムは、実際には使い物にならないことがあるのです。この現象は後の Python 実験で目に見える形で確認します。

同じ議論はノイマン級数にも直結します。$\rho(\bm{M}) < 1$ のとき $\bm{I} - \bm{M}$ は正則で

$$ (\bm{I}-\bm{M})^{-1} = \sum_{k=0}^{\infty} \bm{M}^k $$

が収束し、$\|\bm{M}\| < 1$ なる誘導ノルムがあれば $\|(\bm{I}-\bm{M})^{-1}\| \le 1/(1-\|\bm{M}\|)$ という定量評価まで得られます。ヤコビ法・ガウス-ザイデル法の収束判定(対角優位なら収束する、という定理)は、まさに $\|\bm{M}\|_\infty < 1$ を示すことで得られます。

安定性の話が片付いたので、最後の大きな応用に進みます。連立一次方程式を解いたとき、答えはどれくらい信用してよいのか、という問いです。

条件数 $\kappa(\bm{A})$ — 誤差増幅率の正体

センサから得た右辺 $\bm{b}$ には必ず測定誤差が乗っています。$\bm{A}\bm{x} = \bm{b}$ を解いて得た $\bm{x}$ は、その誤差をどれだけ受け継ぐのでしょうか。「入力の相対誤差 → 出力の相対誤差」の増幅率を評価するのが、ここでの目的です。

右辺の摂動に対する誤差評価

$\bm{A}$ を正則な正方行列とし、真の系 $\bm{A}\bm{x} = \bm{b}$ と摂動系 $\bm{A}(\bm{x}+\delta\bm{x}) = \bm{b} + \delta\bm{b}$ を考えます。両者を引き算すると

$$ \bm{A}\,\delta\bm{x} = \delta\bm{b} \quad \Longrightarrow \quad \delta\bm{x} = \bm{A}^{-1}\delta\bm{b} $$

誤差だけの方程式が得られました。ここに基本不等式を適用します。

$$ \|\delta\bm{x}\| \le \|\bm{A}^{-1}\|\,\|\delta\bm{b}\| $$

一方、元の式 $\bm{b} = \bm{A}\bm{x}$ にも同じ不等式を使うと $\|\bm{b}\| \le \|\bm{A}\|\,\|\bm{x}\|$、変形して

$$ \frac{1}{\|\bm{x}\|} \le \frac{\|\bm{A}\|}{\|\bm{b}\|} $$

この2つを掛け合わせます。左辺同士・右辺同士を掛けると

$$ \begin{equation} \frac{\|\delta\bm{x}\|}{\|\bm{x}\|} \le \|\bm{A}\|\,\|\bm{A}^{-1}\| \, \frac{\|\delta\bm{b}\|}{\|\bm{b}\|} \end{equation} $$

この係数を条件数と呼びます。

$$ \kappa(\bm{A}) = \|\bm{A}\|\,\|\bm{A}^{-1}\| $$

この不等式は「右辺の相対誤差は、解の相対誤差として高々 $\kappa(\bm{A})$ 倍に増幅される」と読めます。$\kappa(\bm{A}) \approx 10^t$ なら、有効数字が約 $t$ 桁失われる、という実務的な目安になります。倍精度(有効約16桁)で $\kappa = 10^{12}$ の問題を解けば、答えの有効数字は4桁程度しか残りません。

条件数の基本性質

(a) 常に $\kappa(\bm{A}) \ge 1$。誘導ノルムの劣乗法性より $1 = \|\bm{I}\| = \|\bm{A}\bm{A}^{-1}\| \le \|\bm{A}\|\|\bm{A}^{-1}\| = \kappa(\bm{A})$。等号は直交行列などで達成されます。$\kappa = 1$ の行列は「誤差を一切増幅しない」最良の行列です。

(b) 2ノルムでは特異値の比。$\bm{A} = \bm{U}\bm{\Sigma}\bm{V}^\top$ なら $\bm{A}^{-1} = \bm{V}\bm{\Sigma}^{-1}\bm{U}^\top$ であり、$\bm{\Sigma}^{-1}$ の最大対角成分は $1/\sigma_n$ です。よって

$$ \kappa_2(\bm{A}) = \frac{\sigma_{\max}(\bm{A})}{\sigma_{\min}(\bm{A})} = \frac{\sigma_1}{\sigma_n} $$

幾何学的には、単位球が写った楕円体の最長軸と最短軸の比、つまり楕円体の”つぶれ具合”です。つぶれた楕円体ほど、逆変換で誤差が大きく拡大されます。

(c) スケール不変。$\kappa(c\bm{A}) = \|c\bm{A}\|\|(c\bm{A})^{-1}\| = |c|\|\bm{A}\| \cdot |c|^{-1}\|\bm{A}^{-1}\| = \kappa(\bm{A})$。行列全体を何倍しても条件数は変わりません。「行列の要素が小さいから危ない」のではなく、「方向によるゲインの差が大きいから危ない」のです。

(d) 特異行列への距離。2ノルムでは $\min\{\|\bm{E}\|_2 : \bm{A}+\bm{E} \text{ が特異}\} = \sigma_n$ が成り立ちます(エッカート・ヤングの定理)。これを $\|\bm{A}\|_2$ で割ると

$$ \frac{\text{特異行列までの相対距離}}{1} = \frac{\sigma_n}{\sigma_1} = \frac{1}{\kappa_2(\bm{A})} $$

つまり条件数が大きい ⟺ 特異行列のすぐ近くにいる。この解釈が最も本質的です。

行列側の摂動

右辺だけでなく $\bm{A}$ 自身が $\bm{A}+\delta\bm{A}$ と揺れる場合も、同様の評価が得られます。$(\bm{A}+\delta\bm{A})(\bm{x}+\delta\bm{x}) = \bm{b}$ を展開し、$\bm{A}\bm{x} = \bm{b}$ を引くと

$$ \bm{A}\,\delta\bm{x} + \delta\bm{A}(\bm{x}+\delta\bm{x}) = \bm{0} $$

$\delta\bm{x}$ について解くと $\delta\bm{x} = -\bm{A}^{-1}\delta\bm{A}(\bm{x}+\delta\bm{x})$。ノルムをとって基本不等式を2回使うと

$$ \|\delta\bm{x}\| \le \|\bm{A}^{-1}\|\,\|\delta\bm{A}\|\,\|\bm{x}+\delta\bm{x}\| $$

両辺を $\|\bm{x}+\delta\bm{x}\|$ で割り、右辺の分母分子に $\|\bm{A}\|$ を掛けて整理すると

$$ \frac{\|\delta\bm{x}\|}{\|\bm{x}+\delta\bm{x}\|} \le \kappa(\bm{A}) \frac{\|\delta\bm{A}\|}{\|\bm{A}\|} $$

やはり増幅率は $\kappa(\bm{A})$ です。行列の丸め誤差も、右辺の測定誤差も、同じ係数で効いてくるわけです。

手計算でわかる具体例

この誤差評価の不等式がどれだけタイトかを、小さな例で確かめます。

$$ \bm{A} = \begin{pmatrix} 1 & 1 \\ 1 & 1.0001 \end{pmatrix}, \qquad \bm{b} = \begin{pmatrix} 2 \\ 2.0001\end{pmatrix} $$

行列式は $1 \times 1.0001 – 1\times 1 = 0.0001$ なので

$$ \bm{A}^{-1} = \frac{1}{0.0001}\begin{pmatrix}1.0001 & -1 \\ -1 & 1\end{pmatrix} = \begin{pmatrix}10001 & -10000 \\ -10000 & 10000\end{pmatrix} $$

解は $\bm{x} = (1, 1)^\top$ です(実際 $1+1=2$、$1+1.0001 = 2.0001$)。ここで右辺の第2成分だけを $0.0001$ だけ動かし、$\bm{b}’ = (2, 2.0002)^\top$ にしてみます。$\delta\bm{b} = (0, 0.0001)^\top$ なので

$$ \delta\bm{x} = \bm{A}^{-1}\delta\bm{b} = \begin{pmatrix}10001\cdot 0 + (-10000)\cdot 0.0001 \\ -10000\cdot 0 + 10000 \cdot 0.0001\end{pmatrix} = \begin{pmatrix}-1 \\ 1\end{pmatrix} $$

新しい解は $\bm{x}’ = (0, 2)^\top$ です。右辺を 0.0001 動かしただけで、解が $(1,1)$ から $(0,2)$ へ完全に変わってしまいました。

相対量で見ましょう。1ノルムで測ると

$$ \frac{\|\delta\bm{b}\|_1}{\|\bm{b}\|_1} = \frac{0.0001}{4.0001} \approx 2.50\times 10^{-5}, \qquad \frac{\|\delta\bm{x}\|_1}{\|\bm{x}\|_1} = \frac{2}{2} = 1 $$

増幅率は $1 / (2.50\times10^{-5}) \approx 4.00\times10^{4}$。一方、条件数は $\|\bm{A}\|_1 = \max(2, 2.0001) = 2.0001$、$\|\bm{A}^{-1}\|_1 = \max(20001, 20000) = 20001$ より

$$ \kappa_1(\bm{A}) = 2.0001 \times 20001 \approx 4.0004\times10^4 $$

増幅率 $4.0001\times10^4$ が条件数 $4.0004\times10^4$ とほぼ一致しました。条件数による誤差評価の不等式は、悪い方向に摂動を入れるとほぼ等号で達成されるのです。逆に言えば、条件数は「最悪ケースでどこまで悪くなるか」の正しい見積もりを与えます。

幾何的にも納得できます。$\bm{A}$ の2行 $(1,1)$ と $(1, 1.0001)$ が定める2直線はほぼ平行で、交点の位置は直線をわずかに動かすだけで大きく滑ります。条件数が大きいとは、この「ほぼ平行」の度合いを数値化したものにほかなりません。

ほぼ平行な2直線と、右辺の微小な変化で交点が(1,1)から(0,2)へ滑る様子を拡大表示した図

左は3本の直線をそのまま描いたもので、傾きの差が $10^{-4}$ しかないため完全に重なって見分けがつきません。それでも交点(解)は $(1,1)$ と $(0,2)$ という遠く離れた2点になっています。右はその仕掛けを見せるために、第1式の直線からの縦方向のずれを1万倍に拡大した図で、緑(元の第2式)と赤(右辺を $0.0001$ 動かした第2式)はほとんど同じ傾きのまま平行移動しているだけなのに、ずれがゼロになる位置=交点が $x_1 = 1$ から $x_1 = 0$ まで動いてしまうことが読み取れます。これが条件数 $4\times10^4$ の正体です。

ここまでの理論をすべて数値実験で確かめます。定義通りの最大化が公式と一致するか、楕円の絵で2ノルムが見えるか、$\bm{A}^k$ が理論どおり振る舞うか、条件数が本当に誤差を予測するかを順に見ていきます。

Python実装 1: 定義通りの最大化と公式の一致

まず、作用素ノルムの定義式を素朴に「単位球からたくさんサンプリングして最大の伸び率を探す」方法で近似し、最大列和・最大行和・最大特異値の公式と一致するかを確かめます。使う行列は

$$ \bm{A} = \begin{pmatrix} 1 & -2 & 0 \\ 3 & 1 & -1 \\ -2 & 4 & 2 \end{pmatrix} $$

です。列和は $(6, 7, 3)$、行和は $(3, 5, 8)$ なので、公式によれば $\|\bm{A}\|_1 = 7$、$\|\bm{A}\|_\infty = 8$ になるはずです。

import numpy as np

rng = np.random.default_rng(0)
A = np.array([[1., -2., 0.],
              [3.,  1., -1.],
              [-2., 4.,  2.]])
N = 200000  # サンプル数

# --- p=2: 単位球面上を一様サンプリング(正規分布を正規化) ---
X2 = rng.normal(size=(N, 3))
X2 /= np.linalg.norm(X2, axis=1, keepdims=True)
ratio2 = np.linalg.norm(X2 @ A.T, axis=1)          # ||Ax||_2 / 1

# --- p=1: 単位球(菱形)の表面 = 単純体の facet を符号付きでサンプリング ---
G = rng.exponential(size=(N, 3)); G /= G.sum(1, keepdims=True)
X1 = G * rng.choice([-1., 1.], size=(N, 3))
ratio1 = np.abs(X1 @ A.T).sum(1) / np.abs(X1).sum(1)

# --- p=inf: 単位球(立方体)の内部をサンプリング ---
Xi = rng.uniform(-1, 1, size=(N, 3))
ratioi = np.abs(Xi @ A.T).max(1) / np.abs(Xi).max(1)

print("p=1   サンプル最大 %.6f  公式(最大列和) %.6f" % (ratio1.max(), np.linalg.norm(A, 1)))
print("p=2   サンプル最大 %.6f  公式(最大特異値) %.6f" % (ratio2.max(), np.linalg.norm(A, 2)))
print("p=inf サンプル最大 %.6f  公式(最大行和) %.6f" % (ratioi.max(), np.linalg.norm(A, np.inf)))

実行すると、$p=1$ でサンプル最大 6.987 に対し公式 7.000、$p=2$ で 5.403994 に対し 5.404062、$p=\infty$ で 7.991 に対し 8.000 となります。どのサンプル最大値も公式の値を決して超えず、下からじわじわ近づいている点が重要です。これは定義が「最大値」である以上当然の振る舞いですが、公式が上界として正しいことの数値的な裏づけになっています。$p=2$ の一致がとくに良いのは、球面上では最大値が滑らかな峰になっていて近傍でも値が落ちにくいのに対し、$p=1,\infty$ では最大値が菱形・立方体の尖った頂点でしか達成されず、ランダム点がその1点を引き当てにくいためです。

p=1,2,∞それぞれで20万点をサンプリングした伸び率のヒストグラムと、公式の値・サンプル最大値の縦線

3つのヒストグラムはいずれも右端で赤い破線(公式の値)にぶつかって切れており、公式の値を超えるサンプルが1点も存在しないことが視覚的に確認できます。$p=2$ では分布が 5.4041 のすぐ手前まで密に詰まっていて、サンプル最大 5.4040 が公式とほぼ一致します。一方 $p=1$ と $p=\infty$ では右端に向かって度数が急激に減っており(縦軸は対数)、頂点付近のごく狭い領域を引き当てる確率が低いために 6.9873 / 7.9908 と、わずかに届いていない様子が見て取れます。

そこで次に、証明で構成した「最大値を達成するベクトル」を直接作って、公式の値をぴったり実現できることを確認します。

import numpy as np

A = np.array([[1., -2., 0.],
              [3.,  1., -1.],
              [-2., 4.,  2.]])

# p=1: 最大列和を与える列の基底ベクトル e_j
j = int(np.argmax(np.abs(A).sum(axis=0)))
e = np.zeros(3); e[j] = 1.0
print("p=1   j*=%d, x*=%s -> ||Ax||_1 = %.4f" % (j, e, np.abs(A @ e).sum()))

# p=inf: 最大行和を与える行の符号ベクトル
i = int(np.argmax(np.abs(A).sum(axis=1)))
x_inf = np.sign(A[i])
print("p=inf i*=%d, x*=%s -> ||Ax||_inf = %.4f" % (i, x_inf, np.abs(A @ x_inf).max()))

# p=2: 最大特異値に対応する右特異ベクトル v1
U, s, Vt = np.linalg.svd(A)
v1 = Vt[0]
print("p=2   v1=%s -> ||Av1||_2 = %.6f (sigma_1 = %.6f)"
      % (np.round(v1, 4), np.linalg.norm(A @ v1), s[0]))
print("特異値:", np.round(s, 6))
print("A^T A の固有値:", np.round(np.linalg.eigvalsh(A.T @ A), 6))

出力は $j^\star = 1$(0始まりなので第2列)で $\|\bm{A}\bm{e}_1\|_1 = 7.0000$、$i^\star = 2$(第3行)で符号ベクトル $(-1, 1, 1)$ が $\|\bm{A}\bm{x}\|_\infty = 8.0000$、$\bm{v}_1 = (-0.5242, 0.7656, 0.3730)$ が $\|\bm{A}\bm{v}_1\|_2 = 5.404062 = \sigma_1$ を与えます。証明で「こう作れば等号になる」と主張したベクトルが、数値的にも寸分違わず公式の値を実現していることがわかります。$\bm{A}^\top\bm{A}$ の固有値は $(0.66228, 10.13384, 29.20388)$ で、その最大値の平方根 $\sqrt{29.20388} = 5.404062$ が確かに $\sigma_1$ と一致しており、最大特異値の公式の導出が裏づけられました。

数値が合ったところで、今度は2ノルムの幾何学的な意味を目で見てみます。

Python実装 2: 単位円が写る楕円と $\|\bm{A}\|_2$

$\|\bm{A}\|_2$ は「単位円が $\bm{A}$ で楕円に写されるとき、その長半径」です。$\rho(\bm{A})$ との違いも同じ絵の中で見えるように、非正規行列

$$ \bm{A} = \begin{pmatrix} 3 & 0 \\ 4 & 5\end{pmatrix} $$

(固有値 3, 5 で $\rho = 5$、特異値 $3\sqrt{5} \approx 6.708$ と $\sqrt{5}\approx 2.236$)を使います。

import numpy as np
import matplotlib, matplotlib.pyplot as plt

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

A = np.array([[3., 0.], [4., 5.]])
th = np.linspace(0, 2*np.pi, 400)
circle = np.vstack([np.cos(th), np.sin(th)])   # 単位円 (2, 400)
ellipse = A @ circle                            # 写った楕円

U, s, Vt = np.linalg.svd(A)
v1, v2 = Vt[0], Vt[1]                           # 右特異ベクトル
w, V = np.linalg.eig(A)
v_eig = V[:, np.argmax(np.abs(w))]              # 最大固有値の固有ベクトル

fig, ax = plt.subplots(1, 2, figsize=(12, 5.5))
ax[0].plot(circle[0], circle[1], color="steelblue", lw=2, label="単位円 $\\|x\\|_2=1$")
for v, c, lb in [(v1, "crimson", "$v_1$(最大に伸びる方向)"),
                 (v2, "darkorange", "$v_2$(最小に伸びる方向)"),
                 (v_eig, "seagreen", "固有ベクトル(伸び率5)")]:
    ax[0].arrow(0, 0, v[0], v[1], head_width=0.06, color=c, length_includes_head=True, label=lb)
ax[0].set_title("入力側:単位円と特徴的な方向"); ax[0].set_aspect("equal"); ax[0].legend(fontsize=9)
ax[0].grid(alpha=.3); ax[0].set_xlim(-1.6, 1.6); ax[0].set_ylim(-1.6, 1.6)

ax[1].plot(ellipse[0], ellipse[1], color="steelblue", lw=2, label="$Ax$ が描く楕円")
for v, c, lb in [(v1, "crimson", "$Av_1$:長さ$\\sigma_1$=6.708"),
                 (v2, "darkorange", "$Av_2$:長さ$\\sigma_2$=2.236"),
                 (v_eig, "seagreen", "$Av$:長さ5(固有値)")]:
    y = A @ v
    ax[1].arrow(0, 0, y[0], y[1], head_width=0.25, color=c, length_includes_head=True, label=lb)
ax[1].set_title("出力側:楕円の長半径が $\\|A\\|_2$"); ax[1].set_aspect("equal")
ax[1].legend(fontsize=9, loc="lower right"); ax[1].grid(alpha=.3)
plt.tight_layout(); plt.show()

単位円が行列Aで楕円に写り、長半径が最大特異値6.708、短半径が2.236になる様子

左が入力側の単位円、右が $\bm{A}$ で写された楕円です。読み取れることは3つあります。第一に、楕円の長半径がちょうど $\sigma_1 = 6.708$、短半径が $\sigma_2 = 2.236$ になっており、$\|\bm{A}\|_2$ が「最も伸びる方向の伸び率」であることが視覚的に確認できます。第二に、$\bm{v}_1$ と $\bm{v}_2$ は入力側で直交しており、出力側でも $\bm{A}\bm{v}_1 \perp \bm{A}\bm{v}_2$(楕円の主軸)になっています。これは特異値分解が「直交系を直交系に写す」変換であることの図示です。第三に、緑の固有ベクトル方向は長さ 5 にしかならず、赤の $\bm{v}_1$ 方向(6.708)に負けています。$\rho(\bm{A}) = 5 < 6.708 = \|\bm{A}\|_2$ という不等式が、そのまま矢印の長さの差として見えているわけです。

この「固有値では捉えきれない増幅」が時間発展でどんな悪さをするかを、次に見ます。

Python実装 3: $\bm{A}^k$ のノルム推移と収束条件

$\bm{A}^k \to \bm{O}$ の判定基準は $\rho(\bm{A}) < 1$ でした。しかし非正規行列では $\|\bm{A}\|_2 \gg 1$ になり得るので、収束するまでに大きな山を作ることがあります。3つの行列で比較します。

import numpy as np
import matplotlib, matplotlib.pyplot as plt

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

mats = {
    "非正規 ρ=0.9(過渡的増幅あり)": np.array([[0.9, 10.], [0., 0.9]]),
    "対角 ρ=0.9(増幅なし)":        np.array([[0.9, 0.], [0., 0.4]]),
    "対角 ρ=1.05(発散)":            np.array([[1.05, 0.], [0., 0.4]]),
}
K = 120
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
for name, M in mats.items():
    P = np.eye(2); vals = []
    for k in range(1, K+1):
        P = P @ M; vals.append(np.linalg.norm(P, 2))
    vals = np.array(vals); ks = np.arange(1, K+1)
    ax[0].semilogy(ks, vals, lw=2, label=name)
    ax[1].plot(ks, vals**(1/ks), lw=2, label=name)
    rho = max(abs(np.linalg.eigvals(M)))
    ax[1].axhline(rho, ls="--", lw=1, color="gray")

ax[0].axhline(1.0, color="k", ls=":", lw=1)
ax[0].set_xlabel("べき乗の回数 $k$"); ax[0].set_ylabel("$\\|A^k\\|_2$(対数軸)")
ax[0].set_title("$A^k$ のノルムの推移"); ax[0].legend(fontsize=9); ax[0].grid(alpha=.3)
ax[1].set_xlabel("べき乗の回数 $k$"); ax[1].set_ylabel("$\\|A^k\\|_2^{1/k}$")
ax[1].set_title("ゲルファントの公式:$k$乗根はスペクトル半径へ(破線)")
ax[1].legend(fontsize=9); ax[1].grid(alpha=.3)
plt.tight_layout(); plt.show()

3つの行列についてA^kのノルムの推移と、そのk乗根がスペクトル半径に収束する様子

左のグラフで、非正規行列 $\begin{pmatrix}0.9&10\\0&0.9\end{pmatrix}$ の挙動が際立ちます。$k=1$ ですでに $\|\bm{A}\|_2 = 10.08$、そこからさらに増えて $k=9$ で最大 38.75 に達し、$\|\bm{A}^k\|_2$ が 1 を下回るのはようやく $k=63$ です。$\rho = 0.9$ なので最終的には $k=120$ で $0.0043$ まで落ちますが、「固有値を見て安定と判断したのに、実運用の 60 ステップの間ずっと増幅されている」という事態がここで起きています。一方、同じ $\rho = 0.9$ でも対角行列は最初から単調減少します。$\rho$ は漸近的な減衰率しか教えてくれず、短期的な振る舞いは作用素ノルムが決めるというのが要点です。$\rho = 1.05$ の場合は指数的に増え、$k=120$ で 348.9 に達しています。

右のグラフはゲルファントの公式の可視化です。3本とも $\|\bm{A}^k\|_2^{1/k}$ がそれぞれのスペクトル半径(破線)に収束していきます。対角行列では $k=1$ から既にぴったり $\rho$ ですが、非正規行列では $k=1$ で 10.08、$k=10$ で 1.442、$k=100$ で 0.965 と、非常にゆっくりとしか 0.9 に近づきません。$k$ 乗根が非正規性を”洗い流す”のに時間がかかる様子が読み取れます。

安定性の話が済んだので、最後に条件数と数値誤差の関係を実験で確かめます。

Python実装 4: 条件数は本当に誤差を予言するか

条件数による誤差評価の不等式が実際の数値計算でどれだけ機能するかを、条件数を人工的に制御した行列で検証します。$20\times20$ の行列を、直交行列 $\bm{U}, \bm{V}$ と対数的に並べた特異値から $\bm{A} = \bm{U}\bm{\Sigma}\bm{V}^\top$ として作れば、条件数を $10^0$ から $10^8$ まで自由に設定できます。

import numpy as np
import matplotlib, matplotlib.pyplot as plt

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

rng = np.random.default_rng(42)
n, eps = 20, 1e-8          # 右辺に与える相対摂動の大きさ
conds, errs, bounds = [], [], []
for _ in range(300):
    k = 10 ** rng.uniform(0, 8)                       # 目標条件数
    U, _ = np.linalg.qr(rng.normal(size=(n, n)))
    V, _ = np.linalg.qr(rng.normal(size=(n, n)))
    s = np.logspace(0, -np.log10(k), n)               # 特異値を対数等間隔に
    A = U @ np.diag(s) @ V.T
    x = rng.normal(size=n); b = A @ x                  # 真の解と右辺
    db = rng.normal(size=n)
    db *= eps * np.linalg.norm(b) / np.linalg.norm(db) # 相対 eps の摂動
    x2 = np.linalg.solve(A, b + db)
    rel_b = np.linalg.norm(db) / np.linalg.norm(b)
    rel_x = np.linalg.norm(x2 - x) / np.linalg.norm(x)
    kappa = np.linalg.cond(A, 2)
    conds.append(kappa); errs.append(rel_x); bounds.append(kappa * rel_b)

conds, errs, bounds = map(np.array, (conds, errs, bounds))
print("誤差/上界 の中央値 %.4f, 最大 %.4f, 上界超過 %d 件"
      % (np.median(errs/bounds), (errs/bounds).max(), (errs > bounds).sum()))
print("log-log 相関係数 %.4f" % np.corrcoef(np.log10(conds), np.log10(errs))[0, 1])

出力は「誤差/上界の中央値 0.0836、最大 0.9765、上界超過 0 件」「log-log 相関係数 0.9888」です。300 回すべてで理論の不等式が守られており、しかも最大では上界の 97.7% まで到達していることがわかります。ランダムな摂動方向では平均的に上界の 1 割弱にとどまりますが、最悪ケースはきちんと上界に張り付きます。相関係数 0.989 は、条件数の対数と誤差の対数がほぼ直線関係にあること、つまり「$\kappa$ が 10 倍になれば誤差も概ね 10 倍になる」ことを示しています。

この関係を散布図にして、上界の線と重ねてみます(前のブロックで作った conds, errs, eps をそのまま使います)。

import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8.5, 5.5))
ax.loglog(conds, errs, "o", ms=4, alpha=.55, color="steelblue", label="ランダムな摂動での実測誤差")
kk = np.logspace(0, 8, 50)
ax.loglog(kk, kk * eps, "r--", lw=2, label="理論上界 $\\kappa_2(A)\\cdot 10^{-8}$")
ax.axhline(1.0, color="gray", ls=":", lw=1)
ax.text(1.5, 1.3, "相対誤差100% = 解が完全に無意味", fontsize=9, color="gray")
ax.set_xlabel("条件数 $\\kappa_2(A)$"); ax.set_ylabel("解の相対誤差 $\\|\\delta x\\|_2/\\|x\\|_2$")
ax.set_title("条件数が誤差増幅率を決める(右辺の相対摂動 $10^{-8}$)")
ax.legend(fontsize=10); ax.grid(alpha=.3, which="both")
plt.tight_layout(); plt.show()

条件数と解の相対誤差の両対数散布図、および理論上界を表す傾き1の破線

散布図は赤い破線(理論上界)の下側にきれいに収まりつつ、傾き 1 の直線に沿って並びます。$\kappa_2 < 10^2$ の領域では相対誤差の中央値が $1.8\times10^{-8}$ で入力の摂動とほぼ同じですが、$\kappa_2 > 10^6$ になると中央値が $2.9\times10^{-3}$、つまり約 5 桁分の精度が失われています。実験設定の摂動が $10^{-8}$ なので、$\kappa_2 = 10^8$ 付近では相対誤差が 1 のオーダー、すなわち解が完全に意味を失う領域に入ります。「有効桁数の損失 $\approx \log_{10}\kappa$」という経験則が、そのまま図に現れているわけです。

なお、この実験では倍精度演算そのものの丸め誤差($\approx 10^{-16}$)も入っていますが、与えた摂動 $10^{-8}$ の方が 8 桁大きいので、観測される誤差はほぼ摂動由来です。逆に摂動をゼロにして丸め誤差だけを見れば、同じ図が $\kappa_2 \cdot 10^{-16}$ の線に沿って現れます。条件数は「どんな解法を使っても超えられない、問題そのものの難しさ」を表しているのです。

まとめ

本記事では、ベクトルノルムから誘導される作用素ノルムを、定義・計算公式・応用の3方向から見てきました。

  • 定義: $\|\bm{A}\|_p = \max_{\bm{x}\neq\bm{0}} \|\bm{A}\bm{x}\|_p/\|\bm{x}\|_p$ は「最大の引き伸ばし率」。単位球面のコンパクト性から最大値は必ず達成され、その達成ベクトルを構成することで閉じた公式が得られます。
  • 3つの公式: $\|\bm{A}\|_1$ は最大列和(達成ベクトルは基底 $\bm{e}_{j^\star}$)、$\|\bm{A}\|_\infty$ は最大行和(達成ベクトルは行の符号ベクトル)、$\|\bm{A}\|_2$ は最大特異値(達成ベクトルは右特異ベクトル $\bm{v}_1$)。前2つは $O(mn)$、2ノルムだけは SVD が必要です。
  • 性質: 誘導ノルムは $\|\bm{I}\|=1$ と劣乗法性 $\|\bm{A}\bm{B}\|\le\|\bm{A}\|\|\bm{B}\|$ を満たします。フロベニウスノルムは $\|\bm{I}\|_F=\sqrt{n}$ なので誘導ノルムではありませんが、劣乗法性は持ち、$\|\bm{A}\|_2 \le \|\bm{A}\|_F \le \sqrt{r}\|\bm{A}\|_2$ で挟めます。
  • スペクトル半径: 任意の誘導ノルムで $\rho(\bm{A}) \le \|\bm{A}\|$。$\rho$ 自体はノルムではありません。ゲルファントの公式 $\rho(\bm{A}) = \lim_k \|\bm{A}^k\|^{1/k}$ が両者を結び、$\bm{A}^k\to\bm{O} \iff \rho(\bm{A})<1$ という完全な判定条件を与えます。ただし非正規行列では収束前に大きな過渡的増幅が起こり得ます(実験では $\rho=0.9$ でも $k=9$ で 38.7 倍)。
  • 条件数: $\kappa(\bm{A}) = \|\bm{A}\|\|\bm{A}^{-1}\| \ge 1$ が相対誤差の増幅率を与え、$\kappa_2 = \sigma_1/\sigma_n$ は「特異行列までの相対距離の逆数」と解釈できます。数値実験では 300 例すべてで理論上界が守られ、最悪ケースは上界の 97.7% に達しました。

行列ノルムは、線形代数の抽象論と数値計算の現場をつなぐ蝶番のような概念です。ここを押さえておくと、反復解法の収束証明、最適化アルゴリズムのステップ幅設計(勾配法の $1/L$ はヘッセ行列の作用素ノルム)、深層学習のリプシッツ定数評価、ロバスト制御の $\mathcal{H}_\infty$ ノルムまで、同じ道具で読み解けるようになります。次は特異値分解の中身をもう一段深く追いかけるか、条件数と数値安定性の具体的なアルゴリズム設計に進むのがおすすめです。

次のステップとして、以下の記事も参考にしてください。