¿Por qué los diferentes métodos para resolver Xc = y en python dan una solución diferente cuando no deberían?

Estaba tratando de resolver un sistema lineal Xc=y que era cuadrado. Los métodos que conozco para resolver esto son:

  1. usando inversa c=
  2. utilizando la eliminación de Gauss
  3. utilizando el pseudoinverso

Por lo que puedo decir, parece que no coinciden con lo que pensé que sería la verdad fundamental.

  1. Primero genere los parámetros de verdad ajustando un polinomio de grado 30 a un coseno con frecuencia 5. Entonces tengo y_truth = X*c_truth .
  2. Luego verifico si los tres métodos anteriores coinciden con la verdad

Lo intenté, pero los métodos no parecen coincidir y no veo por qué debería ser así.

Produje código reproducible completamente ejecutable:

 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) )) 

y me sale el resultado:

 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 

Por que es ¿Es una máquina de precisión? ¿O es porque el error se acumula (mucho) cuando el grado es grande (mayor que 1)? Sinceramente, no lo sé, pero todas esas hipótesis me parecen tontas. Si alguien puede detectar mi error, no dude en señalarlo. De lo contrario, podría no entender el álgebra lineal o algo … lo cual es mucho más preocupante.

Además, si puedo obtener sugerencias para que esto funcione, sería muy apreciado. Yo:

  1. ¿Aumentar el tamaño del intervalo para que no sea menor que 1 (en magnitud)?
  2. ¿Cuál es el mayor tamaño de polinomio que puedo usar?
  3. lenguaje diferente…? ¿O boost la precisión?

Cualquier sugerencia es apreciada!

El problema es la precisión de punto flotante. Estás aumentando los números entre cero y uno a la potencia 30, y luego sumándolos. Si estuviera haciendo esto con aritmética de precisión infinita, los métodos recuperarían las entradas. Con la aritmética de punto flotante, la pérdida de precisión significa que no se puede.