Trinity Operator
作用素を回す
ITERATE · DIAGNOSE · SEE
Trinity-Infinity の三篇が扱うのは、次の一行を繰り返す作用素です。
x ← D Q x + (I − D) p
三篇は、Q を巡回置換に、混合係数 a を 0 と 1 の間に限っています。その仮定の下では作用素ノルム ‖DQ‖₂ が 1 未満になり、Banach の不動点定理がそのまま使えます。ですが、収束を決めているのはそこではありません。
任意の初期値から一意の不動点へ収束する ⟺ ρ(A) < 1
誤差が初手から単調に減衰する ⟸ ‖A‖₂ < 1
ρ ≤ ‖A‖₂ なので、論文の条件は必要以上に強いことになります。等号が成り立つのは A が正規行列のときで、三篇の設定（Q が置換、a が一様）はちょうどその場合です。だから区別が表に出てきませんでした。
これは新しい数学ではありません。非正規行列の過渡的増幅も、収束条件が ρ(A) < 1 であることも、数値線形代数の標準的な事実です。ここで言うのは一点だけ —— 三篇が置いた仮定は、必要以上に強い。
数値を書き換えると、その場で計算し直します。このページは外部へ一切リクエストを出しません。固有値・特異値・連立一次方程式は、すべてこのページの中で解いています。
壊れた議論を、その場で組み直す
ここまでは、論文の条件が必要以上に強いことを測っていました。ここからはその続きで、足りなくなった分を作り直します。
三篇の議論は「‖A‖₂ < 1 だから縮小写像、よって Banach の不動点定理」と進みます。ρ < 1 ≤ ‖A‖₂ の領域では、結論は正しいまま、議論だけが壊れます。ユークリッド距離で測るかぎり縮小写像ではないので、定理をそのまま当てられません。
やることは一つです。縮小写像になる距離を、その場で作る。
P − BᵀPB = I B = A/γ、 ρ < γ < 1
‖Ax‖_P ≤ κ‖x‖_P κ = γ√(1 − 1/λmax(P)) < γ < 1
‖x‖_P = √(xᵀPx) と置いたときの話です。P はリアプノフ（Stein）方程式の解で、下のスライダを動かすとその場で解き直します。κ はこのノルムでの ‖A‖_P の上界ではなく、達成される値そのものです（定義から取り直して突き合わせています）。
‖·‖_P はユークリッド距離と同値なので (ℝⁿ, d_P) は完備で、Banach の不動点定理が字義どおり当たります。ユークリッド距離に戻すと ‖xₖ − x*‖ ≤ √(λmax/λmin)·κᵏ·‖x₀ − x*‖ になり、この係数が上で測った過渡的増幅にあたります。誤差はいったん増えてよい。増えたうえで κ のレートに乗る、という形です。
これも新しい数学ではありません。ρ に任意に近いノルムが存在することは Ostrowski / Householder の古典的な結果で、リアプノフ方程式でそれを構成するのも標準的な手順です。ここでやっているのは、存在証明で終わっている構成を実際に走らせ、出したものが主張どおりかを当たることだけです。
γ は目盛りです。ρ に寄せるほど κ は小さくなりますが、係数は大きくなります。「どんな ρ < 1 でも κ < 1 の距離が取れる」が内容であって、「良い定数が取れる」ではありません。ここでも外部へのリクエストは出していません。n²×n² の連立一次方程式も、コレスキー分解も、このページの中で解いています。
何を計算しているか
A = D Q（D は a を並べた対角行列）、b = (I − D) p と置くと、反復は x ← Ax + b です。これが任意の初期値から収束するのは、ρ(A) < 1 のとき、かつそのときに限ります。ρ はスペクトル半径、つまり固有値の絶対値の最大です。
不動点は (I − A)x = b を解いて求めています。反復を回して近づけたのではなく、直接解いた値です。
固有値は、Householder 変換で上 Hessenberg 形にしてから、複素数のままシフト付き QR 法を回して取っています。実数のシフトだけだと複素固有値の対で止まるためです。特異値は AᵀA の固有値の平方根で、対称なので Jacobi 回転で取っています。連立一次方程式は部分ピボット付きの Gauss 消去です。外部のライブラリは使っていません。
縮小になる距離のほうは、B = A/γ と置いて P − BᵀPB = I を解いています。行優先の vec では vec(BᵀPB) = (Bᵀ ⊗ Bᵀ) vec(P) になるので、クロネッカー積で n²×n² の連立一次方程式に落として、上と同じ Gauss 消去で解いています。P が正定値であることはコレスキー分解が通ることで確かめ、λmax・λmin は Jacobi 回転で取っています。‖A‖_P は、y = Lᵀx と置き直したときの作用素 Lᵀ A L⁻ᵀ の 2-ノルムです。
この実装は NumPy と突き合わせて検査してあります。大きさ 1 から 8 までの行列 728 個 —— 三角行列、強い非正規行列、複素固有値を持つもの、固有値が重なるもの、特異行列を含みます —— について、ρ と ‖A‖₂ の相対誤差は 6.7×10⁻¹⁵ と 1.7×10⁻¹⁵ を超えませんでした。記録した値は verification/trinity_fixtures.json にあり、node verification/check_trinity.js で走ります。
縮小になる距離のほうも、同じやり方で突き合わせています。証書 23 件 —— 論文の設定、反例、三角行列、対称行列、直交行列、大きさ 1 から 8 まで —— について、κ の相対差は 3.9×10⁻¹⁶、係数と P そのものは 4.2×10⁻¹⁰ を超えませんでした。P のほうが桁を落とすのは、γ を ρ に寄せると解く連立一次方程式が悪条件になるからです。それでも κ の精度は保たれます。記録した値は verification/certificate_fixtures.json にあります。
数値計算である以上、限界はあります。固有値が高度に縮退した行列では、ρ の桁数は表示ほど信用できません。ρ が 1 のごく近くにあるとき、「収束する」の判定は浮動小数点の誤差の中に入ります。「境界 ρ = 1」の見本は、まさにその場所を見るためのものです。
作用素そのものの実装と、論文の数値の再現は trinity-operator にあります（Python、MIT）。三篇の PDF と、その検証は trinity-infinity にあります。