矩阵乘法的精度误差

Question

在我的程序中编写矩阵乘法时，出现精度错误（大型矩阵的结果不准确）。

这是我的代码。当前对象的数据逐行存储在扁平数组中。其他矩阵 B 将数据存储在展平数组中，一列又一列（因此我可以使用指针算术）。

protected double[,] multiply (IMatrix B)
{
    int columns = B.columns;
    int rows = Rows;
    int size = Columns;

    double[,] result = new double[rows,columns];
    for (int row = 0; row < rows; row++)
    {
       for (int col = 0; col < columns; col++)
       {
           unsafe
           {
               fixed (float* ptrThis = data)
               fixed (float* ptrB = B.Data)
               {
                   float* mePtr = ptrThis + row*rows;
                   float* bPtr = ptrB + col*columns;
                   double value = 0.0;
                   for (int i = 0; i < size; i++)
                   {
                       value += *(mePtr++) * *(bPtr++);
                   }
                   result[row, col] = value;
               }
           }
       }
    }
}

实际上，代码有点复杂：我对几个块进行乘法操作（因此我不是从 0 到 size，而是从 localStart 到 localStop），然后对结果矩阵求和。

我的问题：对于一个大矩阵，我得到精度误差：

NUnit.Framework.AssertionException: Error at (0,1)
    expected: <6.4209571409444209E+18>
     but was: <6.4207619776304906E+18>

有什么想法吗？

Answer 1

也许您所要做的就是使用 Kahan 求和。但您永远不会指望得到确切特定的浮点数学结果。

Answer 2

事实证明这只是......一个错误。最终结果是：

float* mePtr = ptrThis + row*rows;
float* bPtr = ptrB + col*columns;

我的行的正确索引器是：

float* mePtr = ptrThis + row * size;
float* bPtr = ptrB + col * size;

抱歉，这里的答案并不是很奇特。但感谢您的帮助！

Answer 3

我最初说过，您应该将浮点型转换为双精度型。然而，正如您所指出的，这会破坏您的算法。

您可以尝试：

value += (double)*(mePtr++) * (double)*(bPtr++);

您的代码现在存在的问题是乘法是以float精度完成的，然后添加到double。首先转换为 double 会在一定程度上有所帮助。

使用中间 double 变量可能会更清楚 - 但这取决于您。

如果这不能给您带来所需的准确性，那么您需要考虑使用decimal而不是double。但是，这可能会导致性能下降，因此首先进行一些基准测试。

Answer 4

哼，它并没有真正解决你的问题，但在 NUnit 中，你可以允许有精度误差并选择这个 epsilon 的值

Answer 5

首先，在任何地方使用 double 而不是 float。

Answer 6

至少，你应该自始至终都使用双打。浮动非常不精确。

Answer 7

这是一种称为“矩阵蠕变”的现象，如果您不始终对矩阵进行标准化，这种现象会在矩阵操作过程中逐渐发生。