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

Question

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

假设不同的语言对指向浮点精度的处理略有不同，预计结果应该会有所不同，但到什么程度这是可以接受的？

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

虽然在某些计算中我可以看到其他计算中的值相似，但它们并不完全符合要求。

查看下面的步骤，其中将使用矩阵集的跟踪来确定是否存在解决方案不适定的点，我对 IDL 和 Python 都得到了以下结果:

这看起来相当不错(我可以这么说吗？)。

然后重新组织计算该迹线的(100，y 维度，x 维度)矩阵以计算最小二乘拟合，最终将产生创建均值的值。

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

预先感谢您的宝贵时间。

Answer 1

每个数字表示都会向计算过程注入(inject)两个量:
- 一些值(数字代表的主要数量)
- 一些错误(次要数量，作为副作用，由数字表示引起)

没有前者，任何人都无法进行计算(值(value)...)

没有人可以逃避后者，在实践中，这显然是计算过程流的最终结果错误的(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 工具中并且可以重复使用，不是吗？