Speeding up Math.NET declaration for symmetric matrices

Question

I'm building a reasonably huge matrix with this constructor:

var M = Matrix<double>.Build.Dense(N, N, (i, j) => SomeRoutine(i, j));

N is big and SomeRoutine is slow, so I'm trying to optimize things here and there. I noticed that for any i, j holds SomeRoutine(i, j) == SomeRoutine(j, i), i.e. M is symmetric, therefore one can define only an upper (or lower) triangle, reducing the number of calls to SomeRoutine from N^2 to N(N+1)/2, which is nice.

Here's my approach to this optimization.

var arr = new double[N, N];
for (int i = 0; i < N; i++)
{
  for (int j = 0; j < N; j++)
  {
    arr[i, j] = (i <= j) ? SomeRoutine(i, j) : arr[j, i];
  }
}
var M = Matrix<double>.Build.DenseOfArray(arr);

Doesn't seem very elegant to me. Is there any way I could approach the same optimization idea while retaining the lambda-style declaration? Or maybe I should write some sort of wrapper that would mask for loops?

Answer 1

I think it would be better to split this into two: Firstly, calculate the lower triangle, then in a separate loop assign the values from the lower triangle to the upper triangle. It's likely to be more performant. You might also be able to use Parallel.For for the outer loop (during calculations) to speed it up.

Something like this:

public static void Main()
{
    var test = CreateLowerMatrix(5);
    CopyLowerToUpperMatrix(test, 5);
}

public static double[,] CreateLowerMatrix(int n)
{
    var result = new double[n,n];

    Parallel.For(0, n, r => {
        for (int c = r; c < n; ++c)
            result[r, c] = calc(r, c); });

    return result;
}

public static void CopyLowerToUpperMatrix(double[,] matrix, int n)
{
    for (int r = 1; r < n; ++r)
        for (int c = r + 1; c < n; ++c)
            matrix[c, r] = matrix[r, c];
}

public static double calc(int x, int y)
{
    return x * y;
}

I doubt that it would be useful to parallel up the CopyLowerToUpperMatrix().

Answer 2

If you want something short you can use this. Just note it's not really recommended since it's not really readable, and lambdas with side effects are usually frowned upon. Just imagine the race conditions that would occur if the matrix building function decided to do things in parallel.

I would recommend building the matrix from an array, like in the question or other answers, but if you do decide to use this, just make sure it's well commented.

var cache = new double?[N, N];
Func<int, int, double> WrapAndCacheSomeRoutine = (i, j) => 
    cache[i,j] ?? cache[j,i] ?? (cache[i,j] = SomeRoutine(i, j));
var M = Matrix<double>.Build.Dense(N, N, WrapAndCacheSomeRoutine);