特異な正方行列を解くためにlu分解を使用することについての知識はありますか？

Question

私はランクが不足している（ランク= 9）特異行列A（10 * 10）を持っており、Aの範囲空間にあるベクトルbを持っています。具体性のためにここに私の

array([[ 0.  , 0.  , 0.  , 0.86826141, 0.  , 
      0.  , 0.88788426, 0.  , 0.4089203 , 0.88134901], 
      [ 0.  , 0.  , 0.46416372, 0.  , 0.  , 
      0.  , 0.  , 0.  , 0.  , 0.  ], 
      [ 0.  , 0.  , 0.  , 0.  , 0.31303966, 
      0.  , 0.  , 0.  , 0.  , 0.  ], 
      [ 0.  , 0.  , 0.  , 0.  , 0.  , 
      0.  , 0.3155742 , 0.  , 0.64059294, 0.  ], 
      [ 0.  , 0.  , 0.  , 0.  , 0.51349938, 
      0.  , 0.  , 0.  , 0.53593509, 0.  ], 
      [ 0.  , 0.01252787, 0.  , 0.6870415 , 0.  , 
      0.  , 0.  , 0.  , 0.  , 0.  ], 
      [ 0.  , 0.  , 0.  , 0.  , 0.  , 
      0.  , 0.  , 0.  , 0.  , 0.  ], 
      [ 0.  , 0.  , 0.  , 0.16643105, 0.  , 
      0.  , 0.  , 0.  , 0.  , 0.  ], 
      [ 0.08626592, 0.  , 0.  , 0.  , 0.  , 
      0.  , 0.  , 0.  , 0.  , 0.66939531], 
      [ 0.43694586, 0.  , 0.  , 0.  , 0.  , 
      0.95941661, 0.  , 0.52936733, 0.79687149, 0.81463887]]) 

bはA.dot(np.ones(10))を使用して生成されています。今私は、この使用してLU分解を解決したかったし、そのために私はこれも

array([ nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]) 

を与えるlu_factorこの場合には正常に動作しているようだ

lu_fac=scipy.linalg.lu_factor(X) 
scipy.linalg.lu_solve(lu_fac,b) 

（それが実行時間を与えないいくつかの時間を以下でした"対角番号％dは正確にゼロです。特異行列"）。完全を期すためにここlu_factorからPLUを検証するためのコードは同じです：

L=np.tril(lu_fac[0]) 
np.fill_diagonal(L,1) 
U=np.triu(lu_fac[0]) 
perm=np.arange(10) 
ipiv=lu_factor[1] 
for i in range(10): 
    temp=perm[i] 
    perm[i]=perm[ipiv[i]] 
    perm[ipiv[i]]=temp 
np.allclose(X[perm,:],L.dot(U)) 

は今、私は私の行列が特異である知っていると私の問題に多くのソリューションは無限にあります。しかし、私は任意の解に興味があり、私はちょうどluの分解が失敗する理由を混乱させています。自由変数を0に設定して教えているような解決策を見つけることはできませんか？実行時間警告の対処は "対角番号％dは正確にゼロです。特異行列"。私はsvd/qrのアプローチでこれを解決するのには興味がありません。なぜluが特異行列に対して失敗するのか不思議です。どんな提案も大歓迎です。ありがとう。

Answer 1

0/lu_fac[0][9, 9] 

戻りnan - Uの最後の対角エントリは、ゼロです。したがって、このnanは9番目の変数の値になります。上記の方程式に代入され、当然残りはnanとして出てきます。 SciPyのLUコード、またはそれがラップするFortranコードは、ランクが不足している行列用に設計されていないため、決定できない変数の値を構成しません。

ランタイム警告の対処法「対角線番号％dはちょうどゼロです。特異行列」。

警告は明らかです：アルゴリズムが特異行列を検出しましたが、これは予想されません。また、実装が特異行列での使用を意図していないことも示しています。

は理論的だ

の範囲の空間であるベクトルbを持っています。実際には、浮動小数点演算に固有の誤差のために、ランク不足の行列の範囲空間に何かがあるかどうかは確かではありません。 b = A.dot(...)を計算し、Ax = bを解くことができます。浮動小数点数を操作するときにエラーが発生するため、解決策はありません。

ところで、PLU分解はすべての正方行列に存在すると述べましたが、SciPyは必ずしもそれを計算するように設計されているわけではありません。たとえば、

scipy.linalg.lu_factor(np.array([[0, 1], [0, 0]])) 

はNaNで行列を返します。解決策を見つけるためにしようとすると、あなたの応答のための要因U.

Answer 2

hereのように、行列は、rank(A11) + k >= rank([A11 A12]) + rank([A11 A21])の場合にのみLU分解を行います。あなたのケースでは、rank(A11) = 3,k = 5, およびrank([A11 A12]) + rank([A11 A21]) = 9です。ですから行列は条件を満足せず、LU分解をもたない。そのエントリため