SIMD vectorization strategies for group-by operations on multiple, very large data arrays

As said in the comments, vectorization is hard for such aggregation use case. However, it does not mean SIMD is completely useless for your problem. Try this version (untested).

The main idea, this version saves 50% of random loads/stores spent updating the accumulators. It interleaves the accumulators in memory, uses 128-bit load/add/store instructions, and splits the result back into 2 C# arrays after it consumed all input values.

static unsafe void aggregateSse2( double* accumulators, double* a, double* b, byte* dimension, int count )
{
    Debug.Assert( count >= 0 );
    double* aEnd = a + ( count & ( ~1 ) );
    while( a < aEnd )
    {
        // Load accumulator corresponding to the first bucket
        double* accPointer = accumulators + ( 2u * dimension[ 0 ] );
        Vector128<double> acc = Sse2.LoadAlignedVector128( accPointer );

        // Load 2 values from each input array.
        // BTW, possible to use AVX and unroll by 4 instead of 2, using GetLow/GetHigh to extract the 16-byte pieces.
        // Gonna save a bit of loads at the cost of more shuffles, might be slightly faster overall.
        Vector128<double> va = Sse2.LoadVector128( a );
        Vector128<double> vb = Sse2.LoadVector128( b );

        // Increment accumulator with the first value of each array, store back to RAM
        acc = Sse2.Add( acc, Sse2.UnpackLow( va, vb ) );
        Sse2.StoreAligned( accPointer, acc );

        // Load accumulator corresponding to the second bucket.
        // Potentially it might be the same pointer, can't load both in advance.
        accPointer = accumulators + ( 2u * dimension[ 1 ] );
        acc = Sse2.LoadAlignedVector128( accPointer );
        a += 2;
        b += 2;
        dimension += 2;

        // Increment accumulator with the second value of each array, store back to RAM
        acc = Sse2.Add( acc, Sse2.UnpackHigh( va, vb ) );
        Sse2.StoreAligned( accPointer, acc );
    }

    if( 0 != ( count & 1 ) )
    {
        // The input size was odd number, one item left at these pointers.

        // Load a scalar from first input array into lower lane of a vector
        Vector128<double> vec = Sse2.LoadScalarVector128( a );

        // Load the accumulator corresponding to the bucket
        double* accPointer = accumulators + ( 2u * dimension[ 0 ] );
        Vector128<double> acc = Sse2.LoadAlignedVector128( accPointer );

        // Load scalar from second input array into higher lane of that vector
        vec = Sse2.LoadHigh( vec, b );

        // Increment accumulator and store back to RAM
        acc = Sse2.Add( acc, vec );
        Sse2.StoreAligned( accPointer, acc );
    }
}

static unsafe void splitAccumulators( double* values, double* weights, double* accumulators, int numberOfValues )
{
    double* end = accumulators + numberOfValues * 2;
    while( accumulators < end )
    {
        Vector128<double> vec = Sse2.LoadAlignedVector128( accumulators );
        accumulators += 2;
        Sse2.StoreScalar( values, vec );
        values++;
        Sse2.StoreHigh( weights, vec );
        weights++;
    }
}

/// <summary>Align pointer by 16 bytes, rounding up.</summary>
[MethodImpl( MethodImplOptions.AggressiveInlining )]
static unsafe void* roundUpBy16( void* pointer )
{
    if( Environment.Is64BitProcess )  // This branch is missing from JIT output BTW, it's free.
    {
        long a = (long)pointer;
        a = ( a + 15L ) & ( -16L );
        return (void*)a;
    }
    else
    {
        int a = (int)pointer;
        a = ( a + 15 ) & ( -16 );
        return (void*)a;
    }
}

[SkipLocalsInit] // Otherwise the runtime gonna zero-initialize the stack allocated buffer, very slowly with `push 0` instructions in a loop.
public static (double[], double[]) AggregateDataSse2( double[] data, double[] weight, byte[] dimension )
{
    Debug.Assert( data.Length == weight.Length && data.Length == dimension.Length );

    const int numberOfValues = 0x100;
    unsafe
    {
        // The buffer is about 4kb RAM, fits in L1D cache.
        // Allocating 2 extra doubles (16 extra bytes) to align the pointer.
        double* accumulators = stackalloc double[ ( numberOfValues * 2 ) + 2 ];
        // Align by 16 bytes
        accumulators = (double*)roundUpBy16( accumulators );
        // Clear accumulators with zeros, let's hope the implementation of that standard library method is good.
        new Span<double>( accumulators, numberOfValues * 2 ).Fill( 0 );

        // Process the input data
        fixed( double* a = data )
        fixed( double* b = weight )
        fixed( byte* dim = dimension )
            aggregateSse2( accumulators, a, b, dim, data.Length );

        // Split the result into 2 arrays
        double[] totalValue = new double[ numberOfValues ];
        double[] totalWeight = new double[ numberOfValues ];
        fixed( double* values = totalValue )
        fixed( double* weights = totalWeight )
            splitAccumulators( values, weights, accumulators, numberOfValues );

        return (totalValue, totalWeight);
    }
}

It only uses SSE2 because it doesn’t quite need wider ones, but still, should save non-trivial amount of instructions and RAM transactions compared to your scalar version. I would expect some measurable improvement on all computers.

[Sorry if this should be a comment, not an anwser, but not sure how to format code within a comment]

A version of the answer by Soonts that uses interleaved inputs and outputs is below. However, it ends-up looking equivalent to the Vector2 struct, except using doubles instead of floats. Interestingly, this version is slower than the original answer.

Below are the results for the original answer (SSE2), the optimised answer (SSE2a), my "interleaved" version (SSE2i) and the Vector2 (float-based) results (for comparison).

|                           Type |             Method | Seed |      Size |     Mean |   Error |  StdDev |
|------------------------------- |------------------- |----- |---------- |---------:|--------:|--------:|
|  AggregatorBenchmarkDoubleSSE2 | AggregateBenchmark |   42 | 100000000 | 288.3 ms | 2.25 ms | 2.10 ms |
| AggregatorBenchmarkDoubleSSE2a | AggregateBenchmark |   42 | 100000000 | 275.6 ms | 5.31 ms | 5.45 ms |
| AggregatorBenchmarkDoubleSSE2i | AggregateBenchmark |   42 | 100000000 | 280.3 ms | 2.69 ms | 2.52 ms |
|     AggregatorBenchmarkVector2 | AggregateBenchmark |   42 | 100000000 | 259.9 ms | 1.48 ms | 1.31 ms |

        static unsafe void aggregateSse2i(double* accumulators, double* a, byte* dimension, int count)
        {
            Debug.Assert(count >= 0);

            double* aEnd = a + 2 * count;
            while (a < aEnd)
            {
                // Load accumulator corresponding to the first bucket
                double* accPointer = accumulators + (2u * *dimension);
                Vector128<double> acc = Sse2.LoadAlignedVector128(accPointer);
                dimension++;

                // Load pair of values from input array.
                // BTW, possible to use AVX and unroll by 4 instead of 2, using GetLow/GetHigh to extract the 16-byte pieces.
                // Gonna save a bit of loads at the cost of more shuffles, might be slightly faster overall.
                Vector128<double> va = Sse2.LoadVector128(a);
                a += 2;

                // Increment accumulator with the first value of each array, store back to RAM
                acc = Sse2.Add(acc, va);
                Sse2.StoreAligned(accPointer, acc);
            }
        }

        static unsafe void splitAccumulators2(double* values, double* accumulators, int numberOfValues)
        {
            double* end = accumulators + numberOfValues * 2;
            while (accumulators < end)
            {
                Vector128<double> vec = Sse2.LoadAlignedVector128(accumulators);
                accumulators += 2;
                Sse2.Store(values, vec);
                values += 2;
            }
        }

        [SkipLocalsInit] // Otherwise the runtime gonna zero-initialize the stack allocated buffer, very slowly with `push 0` instructions in a loop.
        public static double[] AggregateDataSse2i(double[] data, byte[] dimension)
        {
            Debug.Assert(data.Length == 2 * dimension.Length);

            const int numberOfValues = 0x100;
            unsafe
            {
                // The buffer is about 4kb RAM, fits in L1D cache.
                // Allocating 2 extra doubles (16 extra bytes) to align the pointer.
                double* accumulators = stackalloc double[(numberOfValues * 2) + 2];
                // Align by 16 bytes
                accumulators = (double*)roundUpBy16(accumulators);
                // Clear accumulators with zeros, let's hope the implementation of that standard library method is good.
                new Span<double>(accumulators, numberOfValues * 2).Fill(0);

                // Process the input data
                fixed (double* a = data)
                fixed (byte* dim = dimension)
                    aggregateSse2i(accumulators, a, dim, dimension.Length);

                return new Span<double>(accumulators, numberOfValues * 2).ToArray();
            }
        }