为什么python中解决Xc = y的不同方法在不应该的时候给出不同的解决方案?
问题描述:
我试图解决方形的线性系统Xc=y
。我知道解决这个的方法是:为什么python中解决Xc = y的不同方法在不应该的时候给出不同的解决方案?
- 使用逆
c=<X^-1,y>
- 用高斯消元法
- 使用伪逆
似乎只要我可以告诉大家,这些别不符合我认为的基本事实。
- 首先通过拟合30次多项式到频率为5的余弦来生成真值参数。所以我有
y_truth = X*c_truth
。 - 然后我检查,如果上述三种方法相匹配的真相
我尝试过,但似乎该方法不匹配,我不明白为什么这应该是这样的。
我公司生产的完全可运行的重复性代码:
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)时,错误累积(很多)?老实说,我不知道,但所有这些假设对我来说都很愚蠢。如果有人能够发现我的错误,请随时指出。否则,我可能不会理解线性代数或其他......这更令人担忧。
此外,如果我能得到这个工作的建议,它将不胜感激。我是否:
- 增加间隔的大小不小于1(大小)?
- 我可以使用的最大多项式大小度数是多少?
- 不同的语言......?或者提高精度?
任何建议表示赞赏!
作为一个有趣的侧面问题,如果我要用GD或SGD训练我的多项式模型,如果我将模型限制在'[-1,1]'中,那么我是否会有相同的数值错误,或者这只是特别的(伪)反转矩阵? – Pinocchio
另一个跟进问题。我已经多次提出建议将间隔改为'[-5,5]'。但是,现在我的问题是,为什么不出现相反的问题,然后取代'1.0 + eps_num = NaN'的'1.0 + eps = 1.0'?如果'big_number = number^30'?就像为什么数值问题似乎在提高到30的幂数时小于1的数字更敏感? – Pinocchio