Make c# matrix code faster

Question

Hi all,

Working on some matrix code, I'm concerned of performance issues.

here's how it works : I've a IMatrix abstract class (with all matrices operations etc), implemented by a ColumnMatrix class.

abstract class IMatrix
{
    public int Rows {get;set;}
    public int Columns {get;set;}
    public abstract float At(int row, int column);
}

class ColumnMatrix : IMatrix
{
    private data[];

    public override float At(int row, int column)
    {
        return data[row + columns * this.Rows];
    }
}

This class is used a lot across my application, but I'm concerned with performance issues. Testing only read for a 2000000x15 matrix against a jagged array of the same size, I get 1359ms for array access agains 9234ms for matrix access :

public void TestAccess()
{
    int iterations = 10;

    int rows = 2000000;
    int columns = 15;
    ColumnMatrix matrix = new ColumnMatrix(rows, columns);
    for (int i = 0; i < rows; i++)
        for (int j = 0; j < columns; j++)
            matrix[i, j] = i + j;

    float[][] equivalentArray = matrix.ToRowsArray();

    TimeSpan totalMatrix = new TimeSpan(0);
    TimeSpan totalArray = new TimeSpan(0);

    float total = 0f;
    for (int iteration = 0; iteration < iterations; iteration++)
    {
        total = 0f;
        DateTime start = DateTime.Now;
        for (int i = 0; i < rows; i++)
            for (int j = 0; j < columns; j++)
                total = matrix.At(i, j);
        totalMatrix += (DateTime.Now - start);

        total += 1f; //Ensure total is read at least once.
        total = total > 0 ? 0f : 0f;

        start = DateTime.Now;
        for (int i = 0; i < rows; i++)
            for (int j = 0; j < columns; j++)
                total = equivalentArray[i][j];
        totalArray += (DateTime.Now - start);
    }

    if (total < 0f)
        logger.Info("Nothing here, just make sure we read total at least once.");

    logger.InfoFormat("Average time for a {0}x{1} access, matrix : {2}ms", rows, columns, totalMatrix.TotalMilliseconds);
    logger.InfoFormat("Average time for a {0}x{1} access, array : {2}ms", rows, columns, totalArray.TotalMilliseconds);
    Assert.IsTrue(true);
}

So my question : how can I make this thing faster ? Is there any way I can make my ColumnMatrix.At faster ? Cheers !

Answer 1

If a two-dimensional array performs so much better, which don't you use a two-dimensional array for your class's internal storage, rather than the one-dimensional one with the overhead of calculating the index?

Answer 2

For every access of an element you do a multiplication: row + columns * this.Rows. You might see if internally you could also use a 2 dimensional array

You also gain extra overhead that the thing is abstracted away in a class. You are doing an extra method call everytime you access an element in the matrix

Answer 3

1. Remove abstract class IMatrix. This is wrong because it's not interface and calling overridden methods is slower than calling final (aka non-modifier methods).
2. You could use unsafe code (pointers) to get elements of the array without array-bounds-checks (faster, but more work and unsafe)

Answer 4

The array code you've written can be optimized easily enough as it's clear that you're accessing memory sequentially. This means the JIT compiler will probably do a better job at converting it to native code and that will result in better performance.
Another thing you're not considering is that inlining is still hit and miss so if your At method (why not using an indexer property, by the way?) is not inlined you'll suffer a huge performance hit due to the use of call and stack manipulation. Finally you should consider sealing the ColumnMatrix class because that will make the optimization much easier for the JIT compiler (call is definitely better than callvirt).

Answer 5

As you are using DateTime.Now to measure the performance, the result is quite random. The resolution of the clock is something like 1/20 second, so instead of measuring the actual time, you are measuring where in the code the clock happens to tick.

You should use the Stopwatch class instead, which has much higher resolution.

Answer 6

Change to this:

interface IMatrix
{
    int Rows {get;set;}
    int Columns {get;set;}
    float At(int row, int column);
}

class ColumnMatrix : IMatrix
{
    private data[,];

    public int Rows {get;set;}
    public int Columns {get;set;}

    public float At(int row, int column)
    {
        return data[row,column];
    }
}

You're better off with the interface than the abstract class - if you need common functions of it add extension methods for the interface.

Also a 2D matrix is quicker than either the jagged one or your flattened one.

Answer 7

You can use Parallel programming for speed up your algorithm. You can compile this code, and compare the performance for normal matrix equations (MultiplyMatricesSequential function) and parallel matrix equations (MultiplyMatricesParallel function). You have implemented compare functions of performance of this methods (in Main function).

You can compile this code under Visual Studio 2010 (.NET 4.0)

namespace MultiplyMatrices
{
    using System;
    using System.Collections.Generic;
    using System.Collections.Concurrent;
    using System.Diagnostics;
    using System.Linq;
    using System.Threading;
    using System.Threading.Tasks;

    class Program
    {
        #region Sequential_Loop
        static void MultiplyMatricesSequential(double[,] matA, double[,] matB,
                                                double[,] result)
        {
            int matACols = matA.GetLength(1);
            int matBCols = matB.GetLength(1);
            int matARows = matA.GetLength(0);

            for (int i = 0; i < matARows; i++)
            {
                for (int j = 0; j < matBCols; j++)
                {
                    for (int k = 0; k < matACols; k++)
                    {
                        result[i, j] += matA[i, k] * matB[k, j];
                    }
                }
            }
        }
        #endregion

        #region Parallel_Loop

        static void MultiplyMatricesParallel(double[,] matA, double[,] matB, double[,] result)
        {
            int matACols = matA.GetLength(1);
            int matBCols = matB.GetLength(1);
            int matARows = matA.GetLength(0);

            // A basic matrix multiplication.
            // Parallelize the outer loop to partition the source array by rows.
            Parallel.For(0, matARows, i =>
            {
                for (int j = 0; j < matBCols; j++)
                {
                    // Use a temporary to improve parallel performance.
                    double temp = 0;
                    for (int k = 0; k < matACols; k++)
                    {
                        temp += matA[i, k] * matB[k, j];
                    }
                    result[i, j] = temp;
                }
            }); // Parallel.For
        }

        #endregion


        #region Main
        static void Main(string[] args)
        {
            // Set up matrices. Use small values to better view 
            // result matrix. Increase the counts to see greater 
            // speedup in the parallel loop vs. the sequential loop.
            int colCount = 180;
            int rowCount = 2000;
            int colCount2 = 270;
            double[,] m1 = InitializeMatrix(rowCount, colCount);
            double[,] m2 = InitializeMatrix(colCount, colCount2);
            double[,] result = new double[rowCount, colCount2];

            // First do the sequential version.
            Console.WriteLine("Executing sequential loop...");
            Stopwatch stopwatch = new Stopwatch();
            stopwatch.Start();

            MultiplyMatricesSequential(m1, m2, result);
            stopwatch.Stop();
            Console.WriteLine("Sequential loop time in milliseconds: {0}", stopwatch.ElapsedMilliseconds);

            // For the skeptics.
            OfferToPrint(rowCount, colCount2, result);

            // Reset timer and results matrix. 
            stopwatch.Reset();
            result = new double[rowCount, colCount2];

            // Do the parallel loop.
            Console.WriteLine("Executing parallel loop...");
            stopwatch.Start();
            MultiplyMatricesParallel(m1, m2, result);
            stopwatch.Stop();
            Console.WriteLine("Parallel loop time in milliseconds: {0}", stopwatch.ElapsedMilliseconds);
            OfferToPrint(rowCount, colCount2, result);

            // Keep the console window open in debug mode.
            Console.WriteLine("Press any key to exit.");
            Console.ReadKey();
        }


        #endregion

        #region Helper_Methods

        static double[,] InitializeMatrix(int rows, int cols)
        {
            double[,] matrix = new double[rows, cols];

            Random r = new Random();
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                {
                    matrix[i, j] = r.Next(100);
                }
            }
            return matrix;
        }

        private static void OfferToPrint(int rowCount, int colCount, double[,] matrix)
        {
            Console.WriteLine("Computation complete. Print results? y/n");
            char c = Console.ReadKey().KeyChar;
            if (c == 'y' || c == 'Y')
            {
                Console.WindowWidth = 180;
                Console.WriteLine();
                for (int x = 0; x < rowCount; x++)
                {
                    Console.WriteLine("ROW {0}: ", x);
                    for (int y = 0; y < colCount; y++)
                    {
                        Console.Write("{0:#.##} ", matrix[x, y]);
                    }
                    Console.WriteLine();
                }

            }
        }

        #endregion
    }

}