将一段八度（mrdivide）代码转换为numpy答案

【问题标题】：Converting a piece of octave (mrdivide) code to numpy将一段八度（mrdivide）代码转换为numpy
【发布时间】：2019-08-30 13:05:44
【问题描述】：

我正在尝试将一段 octave 代码转换为 numpy，但我从 octave 和 numpy 得到了不同的结果。这是我的数据（实际上它们比下面给出的数据要大得多）：

A =

 Columns 1 through 6:

   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.908807141567176e-01   9.954090291002812e-01   9.908807141567179e-01
   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01
   9.908807141567176e-01   9.776578008378289e-01   9.567491788333946e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01
   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01
   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00
   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.608235911808946e-01   9.820064469806473e-01   9.954090291002812e-01
   9.776578008378289e-01   9.649505047327671e-01   9.448300707857176e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01
   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01
   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01
   9.567491788333946e-01   9.448300707857176e-01   9.259179407344914e-01   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01
   9.608235911808946e-01   9.567491788333946e-01   9.448300707857176e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01
   9.567491788333946e-01   9.608235911808946e-01   9.567491788333946e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01
   9.448300707857176e-01   9.567491788333946e-01   9.608235911808946e-01   9.448300707857176e-01   9.649505047327673e-01   9.776578008378289e-01

 Columns 7 through 12:

   9.776578008378289e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01   9.567491788333946e-01
   9.908807141567179e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.448300707857176e-01
   9.954090291002812e-01   9.448300707857176e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.259179407344914e-01
   9.608235911808946e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01   9.820064469806473e-01
   9.820064469806473e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.776578008378289e-01
   9.954090291002812e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.649505047327671e-01
   1.000000000000000e+00   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.448300707857176e-01
   9.567491788333946e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.608235911808946e-01   9.954090291002812e-01
   9.776578008378289e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.908807141567179e-01
   9.908807141567179e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.776578008378289e-01
   9.954090291002812e-01   9.608235911808946e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.567491788333946e-01
   9.448300707857176e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01   1.000000000000000e+00
   9.649505047327673e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.954090291002812e-01
   9.776578008378289e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.820064469806473e-01
   9.820064469806473e-01   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.608235911808946e-01

 Columns 13 through 15:

   9.608235911808946e-01   9.567491788333946e-01   9.448300707857176e-01
   9.567491788333946e-01   9.608235911808946e-01   9.567491788333946e-01
   9.448300707857176e-01   9.567491788333946e-01   9.608235911808946e-01
   9.776578008378289e-01   9.649505047327671e-01   9.448300707857176e-01
   9.820064469806473e-01   9.776578008378289e-01   9.649505047327673e-01
   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01
   9.649505047327673e-01   9.776578008378289e-01   9.820064469806473e-01
   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01
   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01
   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01
   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.954090291002812e-01   9.820064469806473e-01   9.608235911808946e-01
   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01
   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01
   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00

和

b =

  -1.024208397018539
  -1.055718555015945
  -1.066560607689640
  -1.187368188387253
  -1.258866007703282
  -1.305258462589997
  -1.321354530870290
  -1.333661132027602
  -1.421384660329320
  -1.478743779481671
  -1.498725636719488
  -1.385960967135295
  -1.479663779776475
  -1.541078471216082
  -1.562500000000000

在八度音阶中，我有x = b'/A。我找不到 / 对应的 numpy 函数。到目前为止，我从 numpy 尝试了x = np.dot(b.T,np.linalg.inv(A))，但结果与 octave 不同。

x = b'/A 的倍频结果为：

x =

 Columns 1 through 6:

  -5.642309525591432e+00   7.813412870218545e+00  -1.559855155426489e+02  -3.597241224212262e+01   2.201914551287831e+02  -3.100354445411479e+02

 Columns 7 through 12:

   7.253956488595386e+02   7.369595892720794e+01  -4.313273469816049e+02   6.064725968037579e+02  -9.855235323530542e+02  -4.111380598448122e+01

 Columns 13 through 15:

   2.334109900297194e+02  -3.269547704109582e+02   4.254317069619117e+02

numpy的结果是

x=np.array([[-5.642310487492068,    7.813414778371225,  -155.9855165364133, -35.9724138597885,  220.1914623928024,  -310.0354544342263, 725.3956532399461,  73.69596218669903,  -431.3273588509765, 606.472611254314,   -985.5235383154941, -41.11380770278629, 233.4109958125337,  -326.9547770833597, 425.4317096135928]])

如果有人可以帮助我找到与 numpy 八度音程相同的结果，我将不胜感激。或者比当前结果更接近于 Octave 结果。

【问题讨论】：

为了测试，请添加数据的python代码
我刚刚通过numpy.allclose() 运行了两个阵列。它们是一样的。
显式计算逆会导致糟糕的结果。请改用scypy.linalg.solve（NumPy 也有一个）。您将不得不将您的问题从xA=b 重写为Ax=b。
@KingStone，实际上代码很长。但是你可以通过save -6 matrix.mat A b将数据保存为八度，然后导入到numpy中。

标签： python numpy octave

【解决方案1】：

看起来 Matlab 和 numpy 使用不同的精度进行计算。我会尝试强制两者都使用一些浮点精度。

默认情况下，matlab 使用 16 位精度进行数值计算（请参阅https://www.mathworks.com/help/symbolic/increase-precision-of-numeric-calculations.html）。在此站点上还描述了如何更改精度。

numpy，就我而言，它使用 32 位精度（您可能会发现这篇文章很有用：how set numpy floating point accuracy?）。

【讨论】：

【解决方案2】：

获取两个数组（从您的帖子中复制和粘贴）

x = [-5.642309525591432e+00, 7.813412870218545e+00, -1.559855155426489e+02, -3.597241224212262e+01, 2.201914551287831e+02, -3.100354445411479e+02, 7.253956488595386e+02, 7.369595892720794e+01, -4.313273469816049e+02, 6.064725968037579e+02, -9.855235323530542e+02, -4.111380598448122e+01, 2.334109900297194e+02, -3.269547704109582e+02, 4.254317069619117e+02]
y = [-5.642310487492068,    7.813414778371225,  -155.9855165364133, -35.9724138597885,  220.1914623928024,  -310.0354544342263, 725.3956532399461,  73.69596218669903,  -431.3273588509765, 606.472611254314,   -985.5235383154941, -41.11380770278629, 233.4109958125337,  -326.9547770833597, 425.4317096135928]

并通过allclose 运行它们

import numpy as np
np.allclose(x, y) # True

它们（在some tolerance 内）相同。

【讨论】：

我需要更高的准确性
超过 1e-5？为什么？你知道浮点数本质上是有限的精度吗？
我知道，但如果可能的话，我需要更准确（更接近八度音阶）的结果。谢谢。
另外请记住，显示给您的值可能只是所保存值的近似值