PythonでXc = yを解決するためのさまざまな方法が、なぜそうしてはいけないときに異なる解決策を与えるのですか？

Question

私は正方形の線形システムXc=yを解決しようとしていました。これを解決するために、私が知っている方法がある：それは私の知る限り、これらのドン」を言うことができるように思わ擬似逆

を使用してガウスの消去

を使用して逆c=<X^-1,y>

を使用して

1. 私は地面の真実と思ったものにマッチする。
   1. まず、周波数30の余弦に次数30の多項式を当てはめることによって、真理値パラメータを生成します。だから私はy_truth = X*c_truthを持っています。上記の3つの方法が、私はそれを試してみましたが、方法が一致していないようだし、それはケースであるべき理由を私は見ていない真実に

   と一致した場合

2. は、それから私は確認してください。

   私は完全に実行可能な再現性のあるコードを生成：

   import numpy as np 
   from sklearn.preprocessing import PolynomialFeatures 
   
   ## some parameters 
   degree_target = 25 
   N_train = degree_target+1 
   lb,ub = -2000,2000 
   x = np.linspace(lb,ub,N_train) 
   ## generate target polynomial model 
   freq_cos = 5 
   y_cos = np.cos(2*np.pi*freq_cos*x) 
   c_polyfit = np.polyfit(x,y_cos,degree_target)[::-1] ## needs to me reverse to get highest power last 
   ## generate kernel matrix 
   poly_feat = PolynomialFeatures(degree=degree_target) 
   K = poly_feat.fit_transform(x.reshape(N_train,1)) # generates degree 0 first 
   ## get target samples of the function 
   y = np.dot(K,c_polyfit) 
   ## get pinv approximation of c_polyfit 
   c_pinv = np.dot(np.linalg.pinv(K), y) 
   ## get Gaussian-Elminiation approximation of c_polyfit 
   c_GE = np.linalg.solve(K,y) 
   ## get inverse matrix approximation of c_polyfit 
   i = np.linalg.inv(K) 
   c_mdl_i = np.dot(i,y) 
   ## check rank to see if its truly invertible 
   print('rank(K) = {}'.format(np.linalg.matrix_rank(K))) 
   ## comapre parameters 
   print('--c_polyfit') 
   print('||c_polyfit-c_GE||^2 = {}'.format(np.linalg.norm(c_polyfit-c_GE))) 
   print('||c_polyfit-c_pinv||^2 = {}'.format(np.linalg.norm(c_polyfit-c_pinv))) 
   print('||c_polyfit-c_mdl_i||^2 = {}'.format(np.linalg.norm(c_polyfit-c_mdl_i))) 
   print('||c_polyfit-c_polyfit||^2 = {}'.format(np.linalg.norm(c_polyfit-c_polyfit))) 
   ## 
   print('--c_GE') 
   print('||c_GE-c_GE||^2 = {}'.format(np.linalg.norm(c_GE-c_GE))) 
   print('||c_GE-c_pinv||^2 = {}'.format(np.linalg.norm(c_GE-c_pinv))) 
   print('||c_GE-c_mdl_i||^2 = {}'.format(np.linalg.norm(c_GE-c_mdl_i))) 
   print('||c_GE-c_polyfit||^2 = {}'.format(np.linalg.norm(c_GE-c_polyfit))) 
   ## 
   print('--c_pinv') 
   print('||c_pinv-c_GE||^2 = {}'.format(np.linalg.norm(c_pinv-c_GE))) 
   print('||c_pinv-c_pinv||^2 = {}'.format(np.linalg.norm(c_pinv-c_pinv))) 
   print('||c_pinv-c_mdl_i||^2 = {}'.format(np.linalg.norm(c_pinv-c_mdl_i))) 
   print('||c_pinv-c_polyfit||^2 = {}'.format(np.linalg.norm(c_pinv-c_polyfit))) 
   ## 
   print('--c_mdl_i') 
   print('||c_mdl_i-c_GE||^2 = {}'.format(np.linalg.norm(c_mdl_i-c_GE))) 
   print('||c_mdl_i-c_pinv||^2 = {}'.format(np.linalg.norm(c_mdl_i-c_pinv))) 
   print('||c_mdl_i-c_mdl_i||^2 = {}'.format(np.linalg.norm(c_mdl_i-c_mdl_i))) 
   print('||c_mdl_i-c_polyfit||^2 = {}'.format(np.linalg.norm(c_mdl_i-c_polyfit))) 
   

   を、私は結果を得る：

   rank(K) = 4 
   --c_polyfit 
   ||c_polyfit-c_GE||^2 = 4.44089220304006e-16 
   ||c_polyfit-c_pinv||^2 = 1.000000000000001 
   ||c_polyfit-c_mdl_i||^2 = 1.1316233165135605e-06 
   ||c_polyfit-c_polyfit||^2 = 0.0 
   --c_GE 
   ||c_GE-c_GE||^2 = 0.0 
   ||c_GE-c_pinv||^2 = 1.0000000000000007 
   ||c_GE-c_mdl_i||^2 = 1.1316233160694804e-06 
   ||c_GE-c_polyfit||^2 = 4.44089220304006e-16 
   --c_pinv 
   ||c_pinv-c_GE||^2 = 1.0000000000000007 
   ||c_pinv-c_pinv||^2 = 0.0 
   ||c_pinv-c_mdl_i||^2 = 0.9999988683985006 
   ||c_pinv-c_polyfit||^2 = 1.000000000000001 
   --c_mdl_i 
   ||c_mdl_i-c_GE||^2 = 1.1316233160694804e-06 
   ||c_mdl_i-c_pinv||^2 = 0.9999988683985006 
   ||c_mdl_i-c_mdl_i||^2 = 0.0 
   ||c_mdl_i-c_polyfit||^2 = 1.1316233165135605e-06 
   

   なぜそれがあるの？それは機械の精密なものですか？それとも、度が大きい（1より大きい）ときに誤差が累積する（たくさん）ためですか？正直なところ私は知らないけど、これらの仮説はすべて私にはばかげているようだ。誰かが私の間違いを見つけられる場合は、それを指摘してください。そうでなければ、私は線形代数や何かを理解していないかもしれません...より多くの心配です。

   また、私はこれがうまくいくための提案を得ることができれば、非常に感謝します。 Do I：

   1. 間隔のサイズを1より大きくしないでください（大きさで）。
   2. 私が使用できる多項式の最大次数はどれくらいですか？
   3. 異なる言語...？または精度を上げる？

   どのような提案もありがとうございます！

Answer 1

問題は浮動小数点精度です。あなたは、0と1との間の数字を30番目の力に上げてからそれらを加算します。無限精度演算でこれを行っていた場合、メソッドは入力を回復します。浮動小数点演算では、精度の低下はできないことを意味します。