Precision error on matrix multiplication

Question

Hello all,

Coding a matrix multiplication in my program, I get precision errors (inaccurate results for large matrices).

Here's my code. The current object has data stored in a flattened array, row after row. Other matrix B has data stored in a flattened array, column after column (so I can use pointer arithmetic).

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

Actually, the code is a bit more complicated : I do the multiply thing for several chunks (so instead of having i from 0 to size, I go from localStart to localStop), then sum up the resulting matrices.

My problem : for a big matrix I get precision error :

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

Any idea ?

Answer 1

Hem, it doesn't really solve your problem but in NUnit, you can allow to have a precision error and choose the value of this epsilon

Answer 2

As a starting point, use double everywhere instead of float.

Answer 3

I originally stated that you should convert the floats to doubles. However, as you point out that will break your algorithm.

You could try:

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

A problem with your code as it now stands is that the multiplication is being done in float precision then added to a double. Casting to double first will help to some extent.

It might be clearer to use intermediate double variables - but that's up to you.

If this doesn't give you the desire accuracy then you'll need to consider using decimal instead of double. However, this may result in a performance hit so do some benchmarks first.

Answer 4

At the very least, you should be using doubles throughout. Floats are very imprecise.

Answer 5

For the kind of long-running iterative matrix-based computations I work on that's very good agreement. To be more helpful:

• What size are your arrays ? Or, how many f-p operations are you executing ?

• How did you get your expected value ?

Answer 6

Perhaps all you have to do is use Kahan summation. But you can never expect to get exactly a specific result with floating-point math.

Answer 7

This is a phenomenon called "Matrix Creep" which happens gradually during matrix manipulations if you don't consistently normalize your matrices.

Answer 8

Turns out it was just ... a bug. Ended up that instead of having :

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

The correct indexers for my rows were :

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

Sorry for that, not really fancy answer here. But thanks for the help !