pearson_bechmark.cs

using System.Runtime.InteropServices;
using System.Runtime.Intrinsics;
using BenchmarkDotNet.Attributes;
using BenchmarkDotNet.Jobs;
using BenchmarkDotNet.Running;

#pragma warning disable CS8618

var summary = BenchmarkRunner.Run<MainTest>();

// =256
// | Method         | Mean      | Error    | StdDev   | Ratio | RatioSD | Allocated | Alloc Ratio |
// |--------------- |----------:|---------:|---------:|------:|--------:|----------:|------------:|
// | pearson_simple | 330.82 ns | 6.616 ns | 8.125 ns |  1.00 |    0.03 |         - |          NA |
// | pearson_SIMD   |  56.33 ns | 0.672 ns | 0.629 ns |  0.17 |    0.00 |         - |          NA |

// =1024
// | Method         | Mean       | Error   | StdDev  | Ratio | Allocated | Alloc Ratio |
// |--------------- |-----------:|--------:|--------:|------:|----------:|------------:|
// | pearson_simple | 1,303.1 ns | 2.67 ns | 2.49 ns |  1.00 |         - |          NA |
// | pearson_SIMD   |   183.6 ns | 0.43 ns | 0.38 ns |  0.14 |         - |          NA |

// =512K
// | Method         | Mean     | Error   | StdDev  | Ratio | Allocated | Alloc Ratio |
// |--------------- |---------:|--------:|--------:|------:|----------:|------------:|
// | pearson_simple | 677.9 us | 3.29 us | 2.57 us |  1.00 |         - |          NA |
// | pearson_SIMD   | 101.4 us | 0.34 us | 0.30 us |  0.15 |         - |          NA |

[SimpleJob(RuntimeMoniker.Net90, baseline: true)]
// [WarmupCount(2)]
// [IterationCount(2)]
[MemoryDiagnoser]
public class MainTest
{
    const    int      N    = 512*1024 + 3;
    readonly double[] arr1 = Enumerable.Range(0, N).Select(_ => Random.Shared.NextDouble()).ToArray();
    readonly double[] arr2 = Enumerable.Range(0, N).Select(_ => Random.Shared.NextDouble()).ToArray();

    [Benchmark(Baseline = true)]
    public double pearson_simple()
    {
        return pearson_simple(arr1, arr2);
    }

    [Benchmark]
    public double pearson_SIMD()
    {
        return pearson_SIMD(arr1, arr2);
    }

    static double pearson_simple(double[] x, double[] y)
    {
        ArgumentNullException.ThrowIfNull(x);
        ArgumentNullException.ThrowIfNull(y);
        if (x.Length != y.Length)
            throw new ArgumentException("Arrays must have the same length");

        var sumX = 0.0;
        var sumY = 0.0;
        for (var i = 0; i < x.Length; i++)
        {
            sumX += x[i];
            sumY += y[i];
        }

        var sumSquaredXY = 0.0;
        var sumSquaredX  = 0.0;
        var sumSquaredY  = 0.0;
        for (var i = 0; i < x.Length; i++)
        {
            sumSquaredXY += x[i] * y[i];
            sumSquaredX  += x[i] * x[i];
            sumSquaredY  += y[i] * y[i];
        }

        var r1 = (x.Length * sumSquaredXY - sumX * sumY);
        var r2 = Math.Sqrt((x.Length * sumSquaredX - sumX * sumX) * (x.Length * sumSquaredY - sumY * sumY));
        return r1 / r2;
    }

    static double pearson_SIMD(double[] x, double[] y)
    {
        ArgumentNullException.ThrowIfNull(x);
        ArgumentNullException.ThrowIfNull(y);
        if (x.Length != y.Length)
            throw new ArgumentException("Arrays must have the same length");

        var vx = MemoryMarshal.Cast<double, Vector256<double>>(x);
        var vy = MemoryMarshal.Cast<double, Vector256<double>>(y);

        var vectorSumX         = Vector256<double>.Zero;
        var vectorSumY         = Vector256<double>.Zero;
        var vectorSumSquaredXY = Vector256<double>.Zero;
        var vectorSumSquaredX  = Vector256<double>.Zero;
        var vectorSumSquaredY  = Vector256<double>.Zero;
        for (var i = 0; i < vx.Length; i++)
        {
            var vectorX = vx[i];
            var vectorY = vy[i];
            vectorSumX         += vectorX;
            vectorSumY         += vectorY;
            vectorSumSquaredXY += vectorX * vectorY;
            vectorSumSquaredX  += vectorX * vectorX;
            vectorSumSquaredY  += vectorY * vectorY;
        }

        var sumX         = vectorSumX[0]         + vectorSumX[1]         + vectorSumX[2]         + vectorSumX[3];
        var sumY         = vectorSumY[0]         + vectorSumY[1]         + vectorSumY[2]         + vectorSumY[3];
        var sumSquaredXY = vectorSumSquaredXY[0] + vectorSumSquaredXY[1] + vectorSumSquaredXY[2] + vectorSumSquaredXY[3];
        var sumSquaredX  = vectorSumSquaredX[0]  + vectorSumSquaredX[1]  + vectorSumSquaredX[2]  + vectorSumSquaredX[3];
        var sumSquaredY  = vectorSumSquaredY[0]  + vectorSumSquaredY[1]  + vectorSumSquaredY[2]  + vectorSumSquaredY[3];
        for (var i = vx.Length * 4; i < x.Length; i++)
        {
            sumX         += x[i];
            sumY         += y[i];
            sumSquaredXY += x[i] * y[i];
            sumSquaredX  += x[i] * x[i];
            sumSquaredY  += y[i] * y[i];
        }

        var r1 = (x.Length * sumSquaredXY - sumX * sumY);
        var r2 = Math.Sqrt((x.Length * sumSquaredX - sumX * sumX) * (x.Length * sumSquaredY - sumY * sumY));
        return r1 / r2;
    }
}