为什么在 python 中求解 Xc=y 的不同方法会给出不同的解,而它们不应该给出不同的解?

2023-12-14

我试图解决线性系统Xc=y那是方形的。我知道解决这个问题的方法有:

  1. 使用逆c=<X^-1,y>
  2. 使用高斯消去法
  3. 使用伪逆

据我所知,这些似乎与我认为的基本事实不符。

  1. 首先通过将 30 次多项式拟合到频率为 5 的余弦来生成真值参数。所以我有y_truth = X*c_truth.
  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)时误差会累积(很多)?老实说,我不知道,但所有这些假设对我来说似乎都很愚蠢。如果有人发现我的错误,请随时指出。否则,我可能不懂线性代数之类的东西……这更令人担忧。

另外,如果我能得到有关此工作的建议,我将不胜感激。我是否:

  1. 将间隔的大小增加到不小于 1(数量级)?
  2. 我可以使用的最大多项式大小是多少?
  3. 不同的语言...?或者提高精度?

任何建议表示赞赏!


问题是浮点精度。您将 0 到 1 之间的数字求 30 次方,然后将它们相加。如果您使用无限精度算术来执行此操作,这些方法将恢复输入。对于浮点运算,精度损失意味着你不能。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

为什么在 python 中求解 Xc=y 的不同方法会给出不同的解,而它们不应该给出不同的解? 的相关文章

  • 为什么从 Pandas 1.0 中删除了日期时间?

    我在 pandas 中处理大量数据分析并每天使用 pandas datetime 最近我收到警告 FutureWarning pandas datetime 类已弃用 并将在未来版本中从 pandas 中删除 改为从 datetime 模块
  • Django 的内联管理:一个“预填充”字段

    我正在开发我的第一个 Django 项目 我希望用户能够在管理中创建自定义表单 并向其中添加字段当他或她需要它们时 为此 我在我的项目中添加了一个可重用的应用程序 可在 github 上找到 https github com stephen
  • 如何用python脚本控制TP LINK路由器

    我想知道是否有一个工具可以让我连接到路由器并关闭它 然后从 python 脚本重新启动它 我知道如果我写 import os os system ssh l root 192 168 2 1 我可以通过 python 连接到我的路由器 但是
  • 安装了 32 位的 Python,显示为 64 位

    我需要运行 32 位版本的 Python 我认为这就是我在我的机器上运行的 因为这是我下载的安装程序 当我重新运行安装程序时 它会将当前安装的 Python 版本称为 Python 3 5 32 位 然而当我跑步时platform arch
  • 删除flask中的一对一关系

    我目前正在使用 Flask 开发一个应用程序 并且在删除一对一关系中的项目时遇到了一个大问题 我的模型中有以下结构 class User db Model tablename user user id db Column db String
  • Pandas 日期时间格式

    是否可以用零后缀表示 pd to datetime 似乎零被删除了 print pd to datetime 2000 07 26 14 21 00 00000 format Y m d H M S f 结果是 2000 07 26 14
  • 使用 kivy textinput 的 'input_type' 属性的问题

    您好 我在使用 kivy 的文本输入小部件的 input type 属性时遇到问题 问题是我制作了两个自定义文本输入 其中一个称为 StrText 其中设置了 input type text 然后是第二个文本输入 名为 NumText 其
  • Python zmq SUB 套接字未接收 MQL5 Zmq PUB 套接字

    我正在尝试在 MQL5 中设置一个 PUB 套接字 并在 Python 中设置一个 SUB 套接字来接收消息 我在 MQL5 中有这个 include
  • 将 python2.7 与 Emacs 24.3 和 python-mode.el 一起使用

    我是 Emacs 新手 我正在尝试设置我的 python 环境 到目前为止 我已经了解到在 python 缓冲区中使用 python mode el C c C c将当前缓冲区的内容加载到交互式 python shell 中 显然使用了什么
  • 从Python中的字典列表中查找特定值

    我的字典列表中有以下数据 data I versicolor 0 Sepal Length 7 9 I setosa 0 I virginica 1 I versicolor 0 I setosa 1 I virginica 0 Sepal
  • Numpy - 根据表示一维的坐标向量的条件替换数组中的值

    我有一个data多维数组 最后一个是距离 另一方面 我有距离向量r 例如 Data np ones 20 30 100 r np linspace 10 50 100 最后 我还有一个临界距离值列表 称为r0 使得 r0 shape Dat
  • pip 列出活动 virtualenv 中的全局包

    将 pip 从 1 4 x 升级到 1 5 后pip freeze输出我的全局安装 系统 软件包的列表 而不是我的 virtualenv 中安装的软件包的列表 我尝试再次降级到 1 4 但这并不能解决我的问题 这有点类似于这个问题 http
  • 使用特定颜色和抖动在箱形图上绘制数据点

    我有一个plotly graph objects Box图 我显示了箱形 图中的所有点 我需要根据数据的属性为标记着色 如下所示 我还想抖动这些点 下面未显示 Using Box我可以绘制点并抖动它们 但我不认为我可以给它们着色 fig a
  • 如何断言 Unittest 上的可迭代对象不为空?

    向服务提交查询后 我会收到一本字典或一个列表 我想确保它不为空 我使用Python 2 7 我很惊讶没有任何assertEmpty方法为unittest TestCase类实例 现有的替代方案看起来并不正确 self assertTrue
  • 如何在 Windows 命令行中使用参数运行 Python 脚本

    这是我的蟒蛇hello py script def hello a b print hello and that s your sum sum a b print sum import sys if name main hello sys
  • Python ImportError:无法导入名称 __init__.py

    我收到此错误 ImportError cannot import name life table from cdc life tables C Users tony OneDrive Documents Retirement retirem
  • 使用for循环时如何获取前一个元素? [复制]

    这个问题在这里已经有答案了 可能的重复 Python 循环内的上一个和下一个值 https stackoverflow com questions 1011938 python previous and next values inside
  • 模拟pytest中的异常终止

    我的多线程应用程序遇到了一个错误 主线程的任何异常终止 例如 未捕获的异常或某些信号 都会导致其他线程之一死锁 并阻止进程干净退出 我解决了这个问题 但我想添加一个测试来防止回归 但是 我不知道如何在 pytest 中模拟异常终止 如果我只
  • Scipy Sparse:SciPy/NumPy 更新后出现奇异矩阵警告

    我的问题是由大型电阻器系统的节点分析产生的 我基本上是在设置一个大的稀疏矩阵A 我的解向量b 我正在尝试求解线性方程A x b 为了做到这一点 我正在使用scipy sparse linalg spsolve method 直到最近 一切都
  • Django-tables2 列总计

    我正在尝试使用此总结列中的所有值文档 https github com bradleyayers django tables2 blob master docs pages column headers and footers rst 但页

随机推荐

  • Objective-c 上的指针

    据我了解 如果我错了 请纠正我 int x count 10 int hello hello count x hello 这里变量 x 和 count 被声明为整数类型 此外 变量 count 的值被指定为 10 hello 是一个指向整数
  • R将布局对象的网格单位转换为原生

    我的问题有点与使用 R 中的网格将单位从 npc 转换为本地单位 我试图找出 ggplot2 对象中某些绘图元素的位置 轴 主图等 我找到了以下代码 rm list ls library ggplot2 library grid libra
  • 多次调用setcontentview?

    如果我的布局相同但资源发生变化 我可以多次调用 setcontentview 吗 例如 如果图像在 2 个 imageview 小部件中交换 这实际上是我的应用程序中发生的所有情况 您可以多次切换 setContentView 然而 我在实
  • 无法读取 servlet 中的表单字段[重复]

    这个问题在这里已经有答案了 嘿 我对 servlet 环境很陌生 在这里 我尝试将一个表单发布到我的 servlet 如下所示
  • 如果图像被裁剪/调整大小,Camera Intrinsics 将如何变化?

    我有一个来自 Realsense 相机的录制相机 ROS 包文件 所记录设置的相机内部结构已经知道 图像的初始分辨率为848 480 由于相机视场中有一些视觉障碍 我想裁剪掉图像的顶部 这样我正在使用的视觉 SLAM 算法就不会检测到它 由
  • VBA一次性删除所有幻灯片

    我找到了一段代码 可以一张一张地删除除活动幻灯片 索引1 之外的所有ppt幻灯片 但是 任何人都可以帮助我重写这段代码 以便一键操作该代码 我不想循环播放每张幻灯片 因为大约有 300 张幻灯片需要删除 这是我的代码 Sub Deletes
  • 如何使用smack 4.1发送信息查询包到xmpp服务器?

    如何向xmpp服务器发送信息查询包 换句话说 如何向服务器发送 来查询一些信息
  • 同时读取多个文件是个好主意吗?

    我们公司的一台服务器有32个CPU 我们有1000 个非常大的文件需要处理 我不确定同时读取 32 个文件是否是一个好主意 这样所有内核也可以同时执行独立计算 谁能简单解释一下硬盘的工作原理 如果同时读取32个文件 读取速度会不会变慢 谢谢
  • 如何获取构造函数的代码引用?

    我有以下代码 my coderef MyModule MyTool gt new 但当我尝试时 coderef gt 我收到错误 Not a CODE reference 如何引用构造函数 不调用它 并稍后运行引用的代码 The 是标量解引
  • 使用 SortDescriptor AND Predicate 对 NSMutableArray 进行排序可能吗?

    我有一个 Restaurant 类型的数组 其中有一个 Rating 的 NSSet 评级有一个 ID 和一个值 我想按 ID 为 01 的评级从高到低对餐厅数组进行排序 类似于下面的内容 但是如何在originalArray上一起使用谓词
  • 当输入无效值时,用户必须在 C 中重试新值

    对于我的硬件分配 我必须创建一个程序 根据用户输入输出基于星号的三角形 我已经让我的程序正常工作 只要用户输入一个整数 就会输出正确的三角形 但我的问题是 当输入无效值时 如何使用户必须重新尝试提交值 我查看了论坛 但没有找到类似的问题 i
  • 在摘要页面上显示用户输入

    我在产品旁边有几个数量框 以便用户可以输入他们想要的某种产品的数量 该设置使用 Jquery 逐步进行 第一步由复选框组成 第二步由数量输入组成 我需要帮助 最后一步显示已选择的内容 全部提交收到电子邮件给我 带复选框的第 1 步已完成 在
  • Meteor:观察回调中的 Meteor.call() 不执行

    是否有可能从内部调用服务器方法observeMeteor 中的回调 我整理了一个重现该问题的示例 即Meteor call 从回调中调用myCursor observe 不执行 当从观察回调中调用时 Meteor method它本身也不会回
  • 右值引用的对象内部的对象也是右值引用的吗?

    右值引用的对象内部的对象也是右值引用的吗 struct A struct B A a2 template
  • NullPointerException 接近警报

    我正在创建一个应用程序 它使用 GPS 跟踪用户位置 将经度和纬度存储在 SQLite 数据库中并添加接近警报 参考 http www javacodegeeks com 2011 01 android proximity alerts t
  • MongoDB 正则表达式查询:为什么这不起作用?

    请参阅倒数第二个输入 注意 我正在使用http try mongodb org gt person fullname Derek Litz fullname Derek Litz gt db people save person ok gt
  • 在一个查询中选择多对多表

    我有桌子 RUBRIC RubricID RubricName AUTOR AutorID FirstName LastName BOOK BookID BookName book photo BOOKAUTOR BookID AutorI
  • 从 R 数据帧的行中提取 JSON 数据

    我有一个数据框 其中列的值参数是Json数据 Parameters 1 a 0 b 10 2 11 5 22 1 2 a 3 b 4 0 6 2 3 3 我想提取每行的参数并将它们作为列附加到数据框A B1 B2 and B3 我该怎么做
  • 模板文字语法在 IE11 中不起作用

    当使用 use strict 指令时 反引号字符在 IE11 中不会被识别为有效字符 而它在其他浏览器 例如 Chrome 中工作 考虑到 IE11 即使在 Windows 10 用户中仍然广泛使用 对此行为的解释是什么 use stric
  • 为什么在 python 中求解 Xc=y 的不同方法会给出不同的解,而它们不应该给出不同的解?

    我试图解决线性系统Xc y那是方形的 我知道解决这个问题的方法有 使用逆c