Погрешность при умножении матриц

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;
               }
           }
       }
    }
}

На самом деле, код немного сложнее: я делаю умножение для нескольких кусков (поэтому вместо того, чтобы иметь значение i от 0 до размера, я иду от localStart к localStop), затем суммирую полученные матрицы.

Моя проблема: для большой матрицы я получаю ошибку точности:

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

Есть идеи?

Answer 1

Возможно, все, что вам нужно сделать, это использовать Суммирование Кахана .Но вы можете никогда не ожидать, что получит точно определенный результат с математикой с плавающей точкой.

Answer 2

Оказывается, это была просто ... ошибка.Закончилось это тем, что вместо:

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

правильными индексаторами для моих строк были:

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

Извините за это, на самом деле не очень интересный ответ.Но спасибо за помощь!

Answer 3

Я изначально заявил, что вы должны конвертировать floats в doubles. Однако, как вы указали, это нарушит ваш алгоритм.

Вы можете попробовать:

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

Проблема с вашим кодом в его нынешнем виде заключается в том, что умножение выполняется с точностью float, а затем добавляется к double. В некоторой степени поможет приведение к double.

Возможно, было бы понятнее использовать промежуточные double переменные - но это ваше дело.

Если это не дает вам желаемой точности, вам нужно рассмотреть возможность использования decimal вместо double. Однако это может привести к снижению производительности, поэтому сначала сделайте несколько тестов.

Answer 4

Это явление называется «Матричный ползучесть», которое происходит постепенно во время манипуляций с матрицей, если вы не постоянно нормализуете свои матрицы.

Answer 5

По крайней мере, вы должны использовать двойные повсюду. Поплавки очень неточные.

Answer 6

В качестве отправной точки, используйте double везде вместо float.

Answer 7

Подол, на самом деле это не решит вашу проблему, но в NUnit вы можете разрешить ошибку точности и выбрать значение этого эпсилона