矩阵乘法的精度误差答案

【问题标题】：Precision error on matrix multiplication矩阵乘法的精度误差
【发布时间】：2010-04-27 12:22:01
【问题描述】：

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

这是我的代码。当前对象的数据逐行存储在展平数组中。其他矩阵 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;
               }
           }
       }
    }
}

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

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

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

有什么想法吗？

【问题讨论】：

也许使用precision 标签应该会自动打开那些经常出现在这些问题中的关于浮点数学的网页？

标签： c# matrix floating-point double precision

【解决方案1】：

也许您所要做的就是使用Kahan summation。但是您可以通过never expect 获得准确浮点数学的特定结果。

【讨论】：

【解决方案2】：

原来这只是……一个错误。结束了，而不是：

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

我的行的正确索引器是：

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

很抱歉，这里不是很花哨的答案。但是感谢您的帮助！

【讨论】：

【解决方案3】：

我最初声明您应该将floats 转换为doubles。但是，正如您指出的那样，这会破坏您的算法。

你可以试试：

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

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

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

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

【讨论】：

是的，但是将 mePtr 和 bPtr 转换为双精度，我还可以使用 mePtr++ 还是只使用数组的一半元素？ b.data 是 float[]，如果 mePtr 是指向 b.data[0] 的双指针，mePtr++ 不会是右边 8 个字节，所以 b.data[2] 而不是 b.data[1] ?
@Wam - 在这种情况下，您必须在计算中将值加倍 - 我会更新答案

【解决方案4】：

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

【讨论】：

【解决方案5】：

作为起点，在任何地方都使用double，而不是float。

【讨论】：

我的数组已经是 float[]，无法对此进行任何更改。在某些时候出现内存问题，整个应用程序使用 float[] 等。
那我不确定你能做多少。你抛弃了精确性；再多的聪明也无法挽回。

【解决方案6】：

至少，您应该始终使用双打。浮点数非常不精确。

【讨论】：

嗯，浮点数不如双精度，这是肯定的。但是所有（嗯，几乎所有，我的统计数据都不准确）浮点数都是不精确的。
@High，双精度数比浮点数多 29 位。这使它们的精确度提高了大约 1 亿倍。
@Marcelo -- 确实如此。我们中的许多人仍然在进行认真的数字运算，不得不担心长时间运行计算中的精度损失。将浮点数更改为双精度数只会推迟它不会消除它的问题。当然，有时它会将问题推迟到计算完成之后。

【解决方案7】：

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

【讨论】：