Python：矩阵乘法的C连续数组和F连续数组之间的舍入误差

Question

---

问题 ：

我偶然发现了一个对行业产生深远（且非常有影响力）影响的问题，并且似乎没有在任何地方记录：

似乎在 python 中，矩阵乘法（使用@或np.matmul）在 C 连续数组和 F 连续数组之间给出了不同的答案：

import platform
print(platform.platform()) # Linux-5.10.0-23-cloud-amd64-x86_64-with-glibc2.31
import numpy as np
print(np.__version__) # 1.23.5

np.random.seed(0) # Will work with most of seed, I believe
M, N = 10, 5

X_c = np.random.normal(size=M*N).reshape(M,N)
X_f = np.asfortranarray(X_c)
assert (X_c == X_f).all()

p = np.random.normal(size=N)
assert (X_c @ p == X_f @ p).all() # FAIL

据您所知，此问题是否记录在任何地方？

---

后果示例：

在某些情况下，这些错误可能会变得很大：

示例：在sklearn.linear_model.LinearRegression类中，fit在 F 连续的近奇异矩阵 X（通常来自 pandas DataFrame）的情况下，该方法将导致非常错误的参数。这可能会导致预测到处都是并且 R2 为负值。

Answer 1

概括

\n

在我看来，您注意到的问题是浮点运算的基本限制造成的；尤其

\n

\n
• 浮点计算的顺序（在给定情况下：求和）很重要，并且
\n
• 因此，当给定完全相同的输入 \xc2\xad\xe2\x80\x93 时，选择不同顺序的不同算法可能会导致略有不同的结果，其中任何一个结果都不应被视为比另一个“更正确”。
\n

\n

为什么求和顺序很重要

\n

考虑下面代码中的示例，其中我们想要对数组的所有元素求和a。我们将以两种不同的顺序执行此操作：sum_1从左到右对所有值求和，sum_2逐步对所有相邻的值对求和，直到只剩下一个值。从数学上讲，这应该没有什么区别，因为求和是实数的关联运算。然而，正如我们从最后一行的输出中看到的那样，它与计算机上的浮点数不关联：

\n

import numpy as np\n\neps_half = np.finfo(float).eps / 2\na = np.asarray([0.5, 0.5, eps_half, eps_half], dtype=float)\n\nsum_1 = ((a[0] + a[1]) + a[2]) + a[3]\nsum_2 =  (a[0] + a[1]) + (a[2] + a[3])\nprint(sum_1, "versus", sum_2)\n# >>> 1.0 versus 1.0000000000000002\n

\n

这里发生了什么？请注意 的使用eps，它被定义为1.0 和下一个大于 1.0 的最小可表示浮点数之间的差。对于该float类型，至少在我的系统上，该值是 ca。2e-16，即 0.0000000000000002。该定义意味着，如果我们添加一个小于eps1.0 的值，结果仍然是 1.0！

\n

在上面构建的示例中，后者发生在 for 的计算中，sum_1但不发生在 的计算中sum_2，从它们的计算步骤可以看出：

\n

\n
• 步骤sum_1：\n
  \n
  1. 0.5 + 0.5 \xe2\x86\x92 1.0 添加 a[0] 和 a[1]：按预期工作
  \n
  2. 1.0 + eps/2 \xe2\x86\x92 1.0 将 a[2] 添加到先前的总和中：eps/2 太小，无法将结果表示为大于 1.0 的自身浮点值，因此总和仍为 1.0
  \n
  3. 1.0 + eps/2 \xe2\x86\x92 1.0 将 a[3] 添加到先前的总和中：与步骤 2 中的问题相同
  \n
  \n
\n
• 步骤sum_2：\n
  \n
  1. 0.5 + 0.5 \xe2\x86\x92 1.0 添加 a[0] 和 a[1]：按预期工作
  \n
  2. eps/2 + eps/2 \xe2\x86\x92 eps 添加 a[2] 和 a[3]：按预期工作
  \n
  3. 1.0 + eps \xe2\x86\x92 1.0000000000000002 添加步骤 1 和 2 的结果：按预期工作，因为 eps 足够大，结果可以表示为其自身大于 1.0 的浮点值
  \n
  \n
\n

\n

这与矩阵向量积有什么关系

\n

你可能会问自己，你的观察结果与我上面的例子有什么关系。最有可能发生的情况如下：

\n

对于矩阵向量乘积的结果向量中的每个元素，如代码中计算的那样，本质上我们计算两个向量的点积：(1) 给定矩阵的相应行和 (2) 给定向量。根据定义，点积将给定向量的相应元素相乘，然后将这些中间结果相加。现在，在根据相同的输入值计算矩阵向量乘积时，唯一可能给我们带来不同结果的选择可能在于这个求和；或者更确切地说，它的顺序，作为每个点积的最后一步。事情是这样的：

\n

\n
1. 正如上面我的例子中那样，我们总是看到浮点求和的抵消效应；这不是对 1.0 和 求和的特殊问题eps，而是对大浮点值和小浮点值求和的一般问题。
\n
2. 在数组求和中，没有哪个求和顺序比任何其他顺序“更正确”或“更不正确”；事实上，求和结果与“真实结果”（即以无限精度计算的结果）的接近程度取决于数组中值的顺序。（某些求和顺序通常会比其他求和顺序产生更稳健的结果，但在这一点上讨论这一点就太过分了。）
\n
3. 选择哪种求和顺序取决于我们的矩阵向量乘法算法的实现，并且它可能会选择不同的方式，例如出于速度原因以及给定数据的不同内存布局。在您的情况下，您的矩阵数据肯定有不同的内存布局（C 与 Fortran 顺序），因此最有可能发生的情况是：对于 C 顺序数组，选择与 Fortran 顺序不同的求和顺序一，因此会发生不同的抵消效应，导致（略有）不同的结果。
\n

\n

这一切有什么影响

\n

我尝试总结一些应该作为关键要点的要点：

\n

\n
• 浮点数有一定的局限性，导致它们在实践中的表现与实数在数学理论中的表现不同。浮点运算的局限性是众所周知的，不应感到太意外。如果它们确实令人惊讶，那么以下几乎“经典”且非常广泛的阅读可能会有所帮助：D. Goldberg：每个计算机科学家应该了解浮点算术（我承认，我自己从未完全阅读过它）。
\n
• 在你的例子中，效果相当小。在我的机器上，如果我计算代码的两个矩阵向量乘积版本之间的最大绝对差（使用相同的随机种子），我得到：\n

  print(np.abs(X_c @ p - X_f @ p).max())\n# >>> 2.220446049250313e-16\n

  \n\xe2\x80\xa6 非常小eps（顺便说一句，它是上面的精确值）。这并不意味着应该忽略这种不精确性，例如，它们可能会因渐进计算而放大。
\n
• 为了处理这种不精确性，比较浮点计算的结果通常倾向于使用诸如np.allclose()精确相等比较之类的函数。
\n

\n

我承认，我有点怀疑你的最终结论；也就是说，你看到的巨大错误sklearn.linear_model.LinearRegression有相同的原因，但也许是上述放大效应的结果。也许您应该提供一个最小的可重复示例以供进一步调查。

\n