比较用两种不同语言执行的数学运算之间的数值结果

Question

我目前正在尝试将算法从IDL移植到Python 3，在比较结果时，我遇到了以下问题：如何查看数字并确定我是否有效地再现结果？

假设不同的语言在指向浮点精度方面略有不同，预计结果应该有点不同，但到目前为止这是可以接受的？

在下图中，我将使用IDL和python生成的数据集的平均值来说明我的观点：

虽然在某些计算中，我可以看到其他值相似，但它们并没有达到标准。

查看下面的步骤，其中将使用将要使用的一组矩阵的跟踪来确定是否存在解决方案生成的点，我对IDL和Python都有以下结果：

哪个看起来很不错（我可以这么说吗？）。

然后重新组织计算该迹线的（100，y维，x维）矩阵以计算最小二乘拟合，该最小二乘拟合最终将产生将产生该装置的值。

我正在使用这些比较来找到有关python版本需要更改和改进的线索，因此我感谢任何可以引导我朝这个方向发展的想法。

我要感谢你的时间。

Answer 1

每个数字表示向计算过程注入两个量： - 一些值（由数字表示的主要数量） - 一些错误（次要数量，作为副作用，由数字表示引起）

没有前者，没有人可以计算（价值......）

没有人可以逃避后者，实际上可见，因为（i-）负责计算过程流程的最终结果错误（更好的累积不确定性水平）。

正如IEEE-754（-1985）所倡导的那样，没有太多策略来处理嵌入在简化的“常见”表示中的主要错误（不确定性）。

然而， 有很多科学领域， 其中数值方法很快就会降低结果精度 不够。 。 。所以？

无论是天文学还是行星际飞行动力学计算，有些情况下，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很高兴享受它几乎无限的精确数学，所以最简单的一步就是在python方面重新设计算法，这样就完全包含了这种几乎没有降级的精度数学，你突然站在更安全，独立于IDL原件的位置。

鉴于以这种精确的非降级方式重新制定了所有计算步骤，结果值得花时间：

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工具中并且很酷，可以重复使用，不是吗？

Answer 2

真正的问题不是如果不同的实现给你相似的值，这可能会让你觉得他们是对的。

真正的问题是价值是否有意义！