行列の安定性と条件数・特異値・逆行列

線形計算の分野で大きな比重を占める連立方程式 Ax=b の解について注意点をまとめてみました（Aは正則な正方行列、x と b はベクトル、サイズは全て同じ）。

Ax=b 自体は

  LuSolver lu = LuSolver.Create(A);  // A を因子分解する（A の性質により因子分解法が自動的に選択されます）.
  VectorDenseDouble x = lu.Solve(b); // 因子分解の結果(luが保持します)を用いて連立方程式を解きます.

と簡単に解くことができます。

逆行列

浮動小数の数値計算では常に誤差の影響を考慮する必要があります。上記 Ax=b を解いたつもりで得られた結果の x^* を元に残差 E = Ax^*-b を計算してみたら |E| >> 0 といったことが、実は、珍しくありません。

特に、前記のように因子分解法を用いず、A の逆行列 A^-1 を一旦計算してから x^* = A^-1b とする方法は誤差等の点でまずいやり方です。理由は、もとの行列が帯行列やスパース構造であっても、逆行列は一般に密な行列（帯行列の LU 分解は帯行列の構造を保持します）になってしまうため、メモリー増加、計算回数増加による誤差の蓄積、等が顕著になってしまうためです。

現時点で Mol は、逆行列そのものを保持する機会は少ないと判断して、逆行列の計算をサポートしていません。
もちろん、行列 I を単位行列（全対角要素が１、その他の要素はゼロ）として、下の例で X = lu.Solve(I) とすれば X は Y の逆行列（Y^-1）になります。

例えば、統計処理などでよく P = X・Y^-1・Z のような計算が必要な場合があります。この様な場合、Y^-1 のみに注目せず、 Y^-1・Z が一気に計算できることに注目します（X・Y^-1 は転置すれば同じ事になります）。


  LuSolver lu = LuSolver.Create(Y);         // Y を因子分解する.
  MatrixDenseGeneralDouble W = lu.Solve(Z); // W = Y^-1Z （Y・W = Z だから…）
  P = X*W;                                  // P = X・(Y^-1・Z)

特異値

以下では、x は十分 Ax=b を満たしているとして、連立方程式解の精度については問題にしません。
それでも、計算された x が「望んでいる x か？」ということを考えます。

行列 A やベクトル b （の片方、または両方）が実験値や測定値で、そもそもの A や b に誤差が混入している場合、「その誤差が真の解 x にどの程度影響するのか？」を考察します。

行列 A の安定性には A の条件数（C：Condition Number）が大きく影響します。
※条件数（C≧1）は連立方程式（LuSolver）や特異値分解（SvdSolver）の過程で簡単に計算できます。

  b に誤差が入っている場合： ∥Δx∥・∥x∥^-1    ≦ C・∥Δb∥・∥b∥^-1
  A に誤差が入っている場合： ∥Δx∥・∥x+Δx∥^-1 ≦ C・∥ΔA∥・∥A∥^-1

従って、条件数 C が大きいと誤差（∥Δb∥や∥ΔA∥）が小さい場合でも、解の相対誤差（∥Δx∥・∥x∥^-1）が非常に大きくなる可能性があります。例えば b が測定値の場合、わずかな測定誤差（Δb）によって∥Δx∥・∥x∥^-1が大きく変動した結果、x は正確に計算できたとしても、望む x とは異なってしまう場合があります。

※上式の導出は比較的簡単ですが、詳細は文献などを参照してください（例えば岩波講座情報科学-18 森正武他著「数値計算」等）。
※∥b∥や∥A∥はそれぞれベクトルと行列のノルム（正の実数）です。
※ノルムにはいくつかの定義があります。条件数の計算もノルムに依存するので値は唯一とはなりません。

ベクトルと行列のノルム

ノルムはベクトルや行列の「大きさ」（数の絶対値に相当）を表現するもので、以下の定義があります。

   VectorDenseDouble        v = new VectorDenseDouble(2);
   v[1] = 1; v[2] = 2;
   // ベクトルノルムの定義：
   //    1-Norm         ： ∥v∥₁ =  Σ_i｜v_i｜
   //    2-Norm         ： ∥v∥₂ = (Σ_i｜v_i｜²)^1/2 
   //    ユークリッド-Norm： ∥v∥_E = (Σ_i｜v_i｜²)^1/2  ... ベクトルでは 2-Norm と同じです。
   W("ベクトル v の 1-Norm:         " + v.Norm1()); // ==> 3
   W("ベクトル v の 2-Norm:         " + v.Norm2()); // ==> 2.23606797749979
   W("ベクトル v の ユークリッド-Norm:" + v.NormE()); // ==> 2.23606797749979

   MatrixDenseGeneralDouble A = new MatrixDenseGeneralDouble(2, 2);
   A[1, 1] = 11; A[1, 2] = 12;
   A[2, 1] = 21; A[2, 2] = 22;
   // 行列ノルムの定義：
   //    1-Norm：          ∥A∥₁ = Max_j(Σ_i(|A_ij|)) 
   //    無限-Norm：        ∥A∥_∞ = Max_i(Σ_j(|A_ij|)) 
   //    ユークリッド-Norm： ∥A∥_E = (Σ_iΣ_j|A_ij|)^1/2
   W("行列 A の 1-Norm:         " + A.Norm1()); // ==> 34
   W("行列 A の 無限-Norm:       " + A.NormI()); // ==> 43
   W("行列 A の ユークリッド-Norm:" + A.NormE()); // ==> 34.4963766213207

条件数の定義

条件数 C は以下のように定義されます。

  C = ∥A∥∥A^-1∥   …　連立方程式解法における因子分解 LuSolver.Solve(A) の結果として得られます。

また、行列（サイズ n x n）A の最大特異値 σ_n と最小特異値 σ₁ の商として

  C = σ₁/σ_n   　 …　特異値分解 SvdSolver.Solve(A) の結果を用いて定義通りに計算できます。

とも定義されます。

行列の特異値分解

A を (m × n) の一般行列として、A の特異値分解とは

　A = UΣV^H

の形に分解することです（V^Hは行列 V の複素共役転置）。ここで、U (m × m) と V (n × n) は、ユニタリ行列 (A が複素行列の場合) または直交行列 (Aが実行列の場合)です。 Σ は対角行列で、対角成分に特異値 σ_i を持つ m × n の対角行列になります。

 σ₁ ≥ σ₂ ≥ ... ≥ σ_{min(m, n)} ≥ 0

※ σ_i は A の特異値。また、非ゼロの特異値の個数は A の階数(ランク)となります。
※ σ_i は実数です。また対称（エルミート）行列 A^HA の固有値の平方根でもあります。従って直交行列の特異値は全て 1.0 となります。

連立方程式の解と条件数

岩波講座情報科学-18 森正武他著「数値計算」にある例題を Mol を使用して解いてみます。

  int N = 2;
  MatrixDenseGeneralDouble A = new MatrixDenseGeneralDouble(2, 2);
  A[1,1] = 1; A[1,2] =  1;
  A[2,1] = 1; A[2,2] =  1.0001;

  LuSolver lu = LuSolver.Create(A);
  W("Lu :条件数=" + lu.ConditionNumberI(A.NormI()));         // ==> 40004.0001000044
  MatrixDenseBandDouble B  // SvdSolver.Solve()は行列を破壊するので Clone() します。
       = SvdSolver.Solve((MatrixDenseGeneralDouble)A.Clone());
  W("Svd:条件数=" + (B[1,1] / B[N,N]));                      // ==> 40002.0000749119
  
  VectorDenseDouble b = new VectorDenseDouble(N);
  b[1] = 2; b[2] = 2.0001;
  VectorDenseDouble x = lu.Solve(b);
  T.P("摂動前の解", x);                                     // ==> 1.000    1.000
  T.P("(Ax-b).NormE=" + (A * x - b).NormE());               // ==> 0
  // 摂動
  b[2] =  2.0002;
  x = lu.Solve(b);
  T.P("摂動後の解", x);                                     // ==> 0.000    2.000
  T.P("(Ax-b).NormE=" + (A * x - b).NormE());               // ==> 0

b[2]に非常に小さな摂動 0.0001 を加えたことで解 x の変動は殆ど上限に近い影響が出ていることは以下のように確かめることができます。

1-ノルムからプログラミングするまでもなく簡単に ∥Δx∥・∥x∥^-1 = 1 ∥Δb∥・∥b∥^-1 = 0.25x10^-4 となります。従って ∥Δx∥・∥x∥^-1≒C・∥Δb∥・∥b∥^-1。

上記の計算は２次元ですので逆行列のノルムも簡単に計算でき、∥A∥∥A^-1∥≒C を確認できます。

次は非常に不安定な行列として有名なヒルベルト行列です。

  int N = 5;
  MatrixDenseSymmetricDouble A = new MatrixDenseSymmetricDouble(N);
  VectorDenseDouble b = new  VectorDenseDouble(N);
  for (int i = 1; i <= N; ++i)
  {
     b[i] = i;
     for(int j=i;j<=N;++j) A[i,j] = 1.0/((double)(i+j-1));
  }
  LuSolver lu = LuSolver.Create(A);
  double c =  lu.ConditionNumber1(A.Norm1());
  W("条件数=" + c);                                      // ==> 943656.000001442
  VectorDenseDouble x = lu.Solve(b);
  T.P("摂動前の解", x);                                  // ==> 125.000  -2880.000  14490.000  -24640.000  13230.000
  T.P("(Ax-b).NormE=" + (A * x - b).NormE());            // ==> 1.13686837721616E-12

  double normb = b.Norm1();
  double normx = x.Norm1();

  // 摂動
  b[1] += 0.0001;
  b[3] -= 0.0001;
  b[5] += 0.0001;
  VectorDenseDouble x2 = lu.Solve(b);
  T.P("摂動後の解", x2);                                  // ==> 124.961  -2879.400  14487.837  -24637.200  13228.803
  T.P("(Ax-b).NormE=" + (A * x2 - b).NormE());            // ==> 7.92355866838669E-13

  double normdb = 0.0003;
  double normdx = (x2 - x).Norm1();
  T.P("相対誤差=" + (normdx / normx));                    // ==> 0.00012281224600415
  T.P("最大誤差=" + (c * (normdb / normb)));              // ==> 18.8731200000288

最初の例より条件数は一桁大きくなっていますが、x の相対誤差はそれほどでもありません。

条件数は誤差の取り得る（可能性の）最大値を示しているので、大きければ必ず x が大きく変動するとは限りません。
しかし、条件数が大きい行列は、連立方程式に限らず、多くの計算過程で誤差の影響が大きくなる可能性があるので注意が必要です。

ヒルベルト行列はサイズが大きくなるに従って条件数は加速度的に大きくなり、直ぐに手におえなくなります。

 条件数=29070279.0017475 (N=6)
 条件数=35352585302606.9 (N=10)

例で示したように条件数の計算は簡単ですので、時には b （または A）に小さな摂動を与えて上記のように解析してみることも大変重要です。

※現時点で疎な行列の条件数計算はサポートされていません。