比较使用两种不同语言执行的数学运算之间的数值结果
Comparing numerical results between mathematical operations performed with two different languages
我目前正在尝试将算法从 IDL 移植到 Python3，在比较结果时我遇到了以下问题：我如何查看数字并确定我是否有效地重现了结果？
假设不同的语言对指向浮点精度的处理略有不同，预计结果应该会有所不同，但到什么时候这是可以接受的？
在下图中，我将使用 IDL 和 python 生成的数据集的平均值来说明我的观点：
虽然在某些计算中我可以看到其他计算中的值相似，但它们并没有完全达到目标。
查看下面的步骤，其中一组矩阵的踪迹将用于确定是否存在解不适定的点，我对 IDL 和 Python:
看起来不错（我可以这么说吗？）。
然后重新组织从中计算此轨迹的（100，y 维，x 维）矩阵以计算最小二乘拟合，该拟合最终将产生产生均值的值。
我正在使用这些比较来寻找关于 python 版本需要更改和改进的线索，因此我感谢任何可以引导我朝这个方向发展的想法。
提前感谢您的宝贵时间。
每个数字表示在计算过程中注入两个量：
- some value（数字代表的初级数量）
- 一些错误（次要数量，如side-effect，由数字表示引起）
没有前者就无法计算（值...）
没有人能逃避后者，在实践中可见，computation-process 流的最终结果错误（更好的累积不确定性水平）的 (i-) 责任。
正如 IEEE-754(-1985) 所提倡的那样，没有多少策略可以应对嵌入在简化 "common" 表示中的主要错误（不确定性）。
然而，
在许多科学领域，
其中结果精度很快就会下降的数值方法
是不够的。 . .那么？
无论是天文学，还是行星际飞行动力学计算，在某些情况下，IEEE-754 编号很快就无法提供可接受的服务。
此处计算工具提供了几种方案供选择。
>>> import decimal
>>>
>>> with decimal.localcontext() as locCTX:
... for aPREC in range( 20, 31 ):
... locCTX.prec = aPREC
... ( pure_dec_LSQ_5DoF( locCTX, dec_fmin_x0_SEARCH_TRIM_TO_BE_PRECISE, decX, decY ), pure_dec_RESi( locCTX, dec_fmin_x0_SEARCH_TRIM_TO_BE_PRECISE, decX, decY ) )
...
(Decimal('0.038471115298826195147'), (Decimal('0.023589050081780503'), Decimal('-0.082605913918299990'), Decimal('0.150647690402532134'), Decimal('-0.091630826566012630')))
(Decimal('0.0384711152988261953165'), (Decimal('0.0235890500817804889'), Decimal('-0.0826059139182999933'), Decimal('0.1506476904025321349'), Decimal('-0.0916308265660126301')))
(Decimal('0.03847111529882619531420'), (Decimal('0.02358905008178048823'), Decimal('-0.08260591391829999331'), Decimal('0.15064769040253213501'), Decimal('-0.09163082656601263007')))
(Decimal('0.038471115298826195324048'), (Decimal('0.023589050081780488368'), Decimal('-0.082605913918299993309'), Decimal('0.150647690402532135021'), Decimal('-0.091630826566012630071')))
(Decimal('0.0384711152988261953231489'), (Decimal('0.0235890500817804883582'), Decimal('-0.0826059139182999933087'), Decimal('0.1506476904025321350199'), Decimal('-0.0916308265660126300707')))
(Decimal('0.03847111529882619532322276'), (Decimal('0.02358905008178048835950'), Decimal('-0.08260591391829999330863'), Decimal('0.15064769040253213501998'), Decimal('-0.09163082656601263007070')))
(Decimal('0.038471115298826195323213788'), (Decimal('0.023589050081780488359358'), Decimal('-0.082605913918299993308625'), Decimal('0.150647690402532135019974'), Decimal('-0.091630826566012630070702')))
(Decimal('0.0384711152988261953232136753'), (Decimal('0.0235890500817804883593541'), Decimal('-0.0826059139182999933086251'), Decimal('0.1506476904025321350199740'), Decimal('-0.0916308265660126300707023')))
(Decimal('0.03847111529882619532321367314'), (Decimal('0.02358905008178048835935336'), Decimal('-0.08260591391829999330862505'), Decimal('0.15064769040253213501997413'), Decimal('-0.09163082656601263007070231')))
(Decimal('0.038471115298826195323213665675'), (Decimal('0.023589050081780488359353229'), Decimal('-0.082605913918299993308625043'), Decimal('0.150647690402532135019974132'), Decimal('-0.091630826566012630070702306')))
(Decimal('0.0384711152988261953232136649869'), (Decimal('0.0235890500817804883593532187'), Decimal('-0.0826059139182999933086250437'), Decimal('0.1506476904025321350199741307'), Decimal('-0.0916308265660126300707023064')))
Python乐享其成几乎-无限精度数学 ，因此最简单的向前迈出的一步是 re-design python 侧的算法，因此它纯粹包含这种 almost-non-degrading 精度数学和你突然站在 safer-side 上，不管 IDL 原件在哪里。
以这种精度 non-degrading 方式给出一个 re-formulated 所有计算步骤，结果值得花时间：
def pure_dec_LSQ_5DoF( decCTX, Xopt, decX_measured, decY_measured ): # [PERF] ~ 2400 [us] @ .prec = 14
return decCTX.add( decCTX.add( decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[0], Xopt[2] ) ), Xopt[3] ), decY_measured[0] ), decimal.Decimal( 2 ) ), # ~ 2800 [us] @ .prec = 28
decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[1], Xopt[2] ) ), Xopt[3] ), decY_measured[1] ), decimal.Decimal( 2 ) ) # ~ 7700 [us] @ .prec = 100
),
""" [0] [4] [1] [2] [3] _measured[i] ~ [X1,Y1], ...
| | | | |
| | | | | Xopt[0,1,2,3,4] ~ [a,b,c,d,f]
| | | | | | | | | |
+----------------------|--------------------|--------------------------|------------|----------------------------+ | | | |
| +--------------------------|------------|------------------------------+ | | |
| +------------|--------------------------------+ | |
| +----------------------------------+ |
+-------------------------------------------------------------------------------------------------+
"""
decCTX.add( decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[2], Xopt[2] ) ), Xopt[3] ), decY_measured[2] ), decimal.Decimal( 2 ) ), # ~ 1340 [ms] @ .prec = 1000
decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[3], Xopt[2] ) ), Xopt[3] ), decY_measured[3] ), decimal.Decimal( 2 ) ) #
)
)
如果需要 ~ 14 位精度，只需花费 ~ 2.4 [ms] 每一步，
如果需要 ~ 28 位精度，只需花费 ~ 2.8 [ms] 每个这样的步骤，
如果需要 ~100 位精度，只需花费 ~ 7.7 [ms] 每个这样的步骤，
如果需要 1000 位精度，只需花费 ~ 1.3 [ s] 每个这样的步骤，
还不错，
是吗？
# [PERF] ~ 2400 [us] @ .prec = 14
# ~ 2800 [us] @ .prec = 28
# ~ 7700 [us] @ .prec = 100
# ~ 1340 [ms] @ .prec = 1000
这一切都已经包含在 python 工具中并且很酷 re-use，不是吗？
真正的问题不是不同的实现是否给你相似的值，这可能会让你觉得它们是正确的。
真正的问题是这些值是否有意义！
我目前正在尝试将算法从 IDL 移植到 Python3，在比较结果时我遇到了以下问题：我如何查看数字并确定我是否有效地重现了结果？
假设不同的语言对指向浮点精度的处理略有不同，预计结果应该会有所不同，但到什么时候这是可以接受的？
在下图中，我将使用 IDL 和 python 生成的数据集的平均值来说明我的观点：
虽然在某些计算中我可以看到其他计算中的值相似，但它们并没有完全达到目标。
查看下面的步骤，其中一组矩阵的踪迹将用于确定是否存在解不适定的点，我对 IDL 和 Python:
看起来不错（我可以这么说吗？）。
然后重新组织从中计算此轨迹的（100，y 维，x 维）矩阵以计算最小二乘拟合，该拟合最终将产生产生均值的值。
我正在使用这些比较来寻找关于 python 版本需要更改和改进的线索，因此我感谢任何可以引导我朝这个方向发展的想法。
提前感谢您的宝贵时间。
每个数字表示在计算过程中注入两个量：
- some value（数字代表的初级数量）
- 一些错误（次要数量，如side-effect，由数字表示引起）
没有前者就无法计算（值...）
没有人能逃避后者，在实践中可见，computation-process 流的最终结果错误（更好的累积不确定性水平）的 (i-) 责任。
正如 IEEE-754(-1985) 所提倡的那样，没有多少策略可以应对嵌入在简化 "common" 表示中的主要错误（不确定性）。
然而，
在许多科学领域，
其中结果精度很快就会下降的数值方法
是不够的。 . .那么？
无论是天文学，还是行星际飞行动力学计算，在某些情况下，IEEE-754 编号很快就无法提供可接受的服务。
此处计算工具提供了几种方案供选择。
>>> import decimal
>>>
>>> with decimal.localcontext() as locCTX:
... for aPREC in range( 20, 31 ):
... locCTX.prec = aPREC
... ( pure_dec_LSQ_5DoF( locCTX, dec_fmin_x0_SEARCH_TRIM_TO_BE_PRECISE, decX, decY ), pure_dec_RESi( locCTX, dec_fmin_x0_SEARCH_TRIM_TO_BE_PRECISE, decX, decY ) )
...
(Decimal('0.038471115298826195147'), (Decimal('0.023589050081780503'), Decimal('-0.082605913918299990'), Decimal('0.150647690402532134'), Decimal('-0.091630826566012630')))
(Decimal('0.0384711152988261953165'), (Decimal('0.0235890500817804889'), Decimal('-0.0826059139182999933'), Decimal('0.1506476904025321349'), Decimal('-0.0916308265660126301')))
(Decimal('0.03847111529882619531420'), (Decimal('0.02358905008178048823'), Decimal('-0.08260591391829999331'), Decimal('0.15064769040253213501'), Decimal('-0.09163082656601263007')))
(Decimal('0.038471115298826195324048'), (Decimal('0.023589050081780488368'), Decimal('-0.082605913918299993309'), Decimal('0.150647690402532135021'), Decimal('-0.091630826566012630071')))
(Decimal('0.0384711152988261953231489'), (Decimal('0.0235890500817804883582'), Decimal('-0.0826059139182999933087'), Decimal('0.1506476904025321350199'), Decimal('-0.0916308265660126300707')))
(Decimal('0.03847111529882619532322276'), (Decimal('0.02358905008178048835950'), Decimal('-0.08260591391829999330863'), Decimal('0.15064769040253213501998'), Decimal('-0.09163082656601263007070')))
(Decimal('0.038471115298826195323213788'), (Decimal('0.023589050081780488359358'), Decimal('-0.082605913918299993308625'), Decimal('0.150647690402532135019974'), Decimal('-0.091630826566012630070702')))
(Decimal('0.0384711152988261953232136753'), (Decimal('0.0235890500817804883593541'), Decimal('-0.0826059139182999933086251'), Decimal('0.1506476904025321350199740'), Decimal('-0.0916308265660126300707023')))
(Decimal('0.03847111529882619532321367314'), (Decimal('0.02358905008178048835935336'), Decimal('-0.08260591391829999330862505'), Decimal('0.15064769040253213501997413'), Decimal('-0.09163082656601263007070231')))
(Decimal('0.038471115298826195323213665675'), (Decimal('0.023589050081780488359353229'), Decimal('-0.082605913918299993308625043'), Decimal('0.150647690402532135019974132'), Decimal('-0.091630826566012630070702306')))
(Decimal('0.0384711152988261953232136649869'), (Decimal('0.0235890500817804883593532187'), Decimal('-0.0826059139182999933086250437'), Decimal('0.1506476904025321350199741307'), Decimal('-0.0916308265660126300707023064')))
Python乐享其成几乎-无限精度数学 ，因此最简单的向前迈出的一步是 re-design python 侧的算法，因此它纯粹包含这种 almost-non-degrading 精度数学和你突然站在 safer-side 上，不管 IDL 原件在哪里。
以这种精度 non-degrading 方式给出一个 re-formulated 所有计算步骤，结果值得花时间：
def pure_dec_LSQ_5DoF( decCTX, Xopt, decX_measured, decY_measured ): # [PERF] ~ 2400 [us] @ .prec = 14
return decCTX.add( decCTX.add( decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[0], Xopt[2] ) ), Xopt[3] ), decY_measured[0] ), decimal.Decimal( 2 ) ), # ~ 2800 [us] @ .prec = 28
decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[1], Xopt[2] ) ), Xopt[3] ), decY_measured[1] ), decimal.Decimal( 2 ) ) # ~ 7700 [us] @ .prec = 100
),
""" [0] [4] [1] [2] [3] _measured[i] ~ [X1,Y1], ...
| | | | |
| | | | | Xopt[0,1,2,3,4] ~ [a,b,c,d,f]
| | | | | | | | | |
+----------------------|--------------------|--------------------------|------------|----------------------------+ | | | |
| +--------------------------|------------|------------------------------+ | | |
| +------------|--------------------------------+ | |
| +----------------------------------+ |
+-------------------------------------------------------------------------------------------------+
"""
decCTX.add( decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[2], Xopt[2] ) ), Xopt[3] ), decY_measured[2] ), decimal.Decimal( 2 ) ), # ~ 1340 [ms] @ .prec = 1000
decCTX.power( decCTX.subtract( decCTX.fma( Xopt[0], decCTX.power( Xopt[4], decCTX.fma( Xopt[1], decX_measured[3], Xopt[2] ) ), Xopt[3] ), decY_measured[3] ), decimal.Decimal( 2 ) ) #
)
)
如果需要 ~ 14 位精度，只需花费 ~ 2.4 [ms] 每一步，
如果需要 ~ 28 位精度，只需花费 ~ 2.8 [ms] 每个这样的步骤，
如果需要 ~100 位精度，只需花费 ~ 7.7 [ms] 每个这样的步骤，
如果需要 1000 位精度，只需花费 ~ 1.3 [ s] 每个这样的步骤，
还不错，
是吗？
# [PERF] ~ 2400 [us] @ .prec = 14
# ~ 2800 [us] @ .prec = 28
# ~ 7700 [us] @ .prec = 100
# ~ 1340 [ms] @ .prec = 1000
这一切都已经包含在 python 工具中并且很酷 re-use，不是吗？
真正的问题不是不同的实现是否给你相似的值，这可能会让你觉得它们是正确的。
真正的问题是这些值是否有意义！