How to repair a Parallel for()-loop not to return a different value than a serial for()-loop?

Question

I'm having an issue with the following code. The code works with no errors but I'm receiving different output values when using a parallel for loop vs a regular for loop. I need to get the parallel for loop working properly because I run this code thousands of times. Does anyone know why my parallel for loop is returning different outputs?

private object _lock = new object();

public double CalculatePredictedRSquared()
    {
        double press = 0, tss = 0, press2 = 0, press1 = 0;
        Vector<double> output = CreateVector.Dense(Enumerable.Range(0, 400).Select(i => Convert.ToDouble(i)).ToArray());
        List<double> input1 = new List<double>(Enumerable.Range(0, 400).Select(i => Convert.ToDouble(i)));
        List<double> input2 = new List<double>(Enumerable.Range(200, 400).Select(i => Convert.ToDouble(i)));

            Parallel.For(0, output.Count, i =>
            {
                ConcurrentBag<MultipleRegressionInfo> listMRInfoBag = new ConcurrentBag<MultipleRegressionInfo>(listMRInfo);
                ConcurrentBag<double> vectorArrayBag = new ConcurrentBag<double>(output);
                ConcurrentBag<double[]> matrixList = new ConcurrentBag<double[]>();

                lock (_lock)
                {
                    matrixList.Add(input1.Where((v, k) => k != i).ToArray());
                    matrixList.Add(input2.Where((v, k) => k != i).ToArray());
                }

                var matrixArray2 = CreateMatrix.DenseOfColumnArrays(matrixList);
                var actualResult = vectorArrayBag.ElementAt(i);
                var newVectorArray = CreateVector.Dense(vectorArrayBag.Where((v, j) => j != i).ToArray());
                var items = FindBestMRSolution(matrixArray2, newVectorArray);
                double estimate1 = 0;

                if (items != null)
                {
                    lock (_lock)
                    {
                        var y = 0d;
                        var independentCount = matrixArray2.RowCount;
                        var dependentCount = newVectorArray.Count;

                        if (independentCount == dependentCount)
                        {
                            var populationCount = independentCount;
                            y = newVectorArray.Average();

                            for (int l = 0; l < matrixArray2.ColumnCount; l++)
                            {
                                var avg = matrixArray2.Column(l).Average();
                                y -= avg * items[l];
                            }
                        }

                        for (int m = 0; m < 2; m++)
                        {
                            var coefficient = items[m];

                            if (m == 0)
                            {
                                estimate1 += input1.ElementAt(i) * coefficient;
                            }
                            else
                            {
                                estimate1 += input2.ElementAt(i) * coefficient;
                            }
                        }

                        estimate1 += y;
                    }
                }
                else
                {
                    lock (_lock)
                    {
                        estimate1 = 0;
                    }
                }

                lock (_lock)
                {
                    press1 += Math.Pow(actualResult - estimate1, 2);
                }
            });

            for (int i = 0; i < output.Count; i++)
            {
                List<double[]> matrixList = new List<double[]>();
                matrixList.Add(input1.Where((v, k) => k != i).ToArray());
                matrixList.Add(input2.Where((v, k) => k != i).ToArray());
                var matrixArray = CreateMatrix.DenseOfColumnArrays(matrixList);
                var actualResult = output.ElementAt(i);
                var newVectorArray = CreateVector.Dense(output.Where((v, j) => j != i).ToArray());
                var items = FindBestMRSolution(matrixArray, newVectorArray);
                double estimate = 0;

                if (items != null)
                {
                    var y = CalculateYIntercept(matrixArray, newVectorArray, items);
                    for (int m = 0; m < 2; m++)
                    {
                        var coefficient = items[m];

                        if (m == 0)
                        {
                            estimate += input1.ElementAt(i) * coefficient;
                        }
                        else
                        {
                            estimate += input2.ElementAt(i) * coefficient;
                        }
                    }
                }
                else
                {
                    estimate = 0;
                }

                press2 += Math.Pow(actualResult - estimate, 2);
            }

            tss = CalculateTotalSumOfSquares(vectorArray.ToList());
            var test1 = 1 - (press1 / tss);
            var test2 = 1 - (press2 / tss);
}

public Vector<double> CalculateWithQR(Matrix<double> x, Vector<double> y)
    {
        Vector<double> result = null;

            result = MultipleRegression.QR(x, y);

            for (int i = 0; i < result.Count; i++)
            {
                var value = result.ElementAt(i);

                if (Double.IsNaN(value) || Double.IsInfinity(value))
                {
                    return null;
                }
            }

        return result;
    }

    public Vector<double> CalculateWithNormal(Matrix<double> x, Vector<double> y)
    {
        Vector<double> result = null;

            result = MultipleRegression.NormalEquations(x, y);

            for (int i = 0; i < result.Count; i++)
            {
                var value = result.ElementAt(i);

                if (Double.IsNaN(value) || Double.IsInfinity(value))
                {
                    return null;
                }
            }

        return result;
    }

    public Vector<double> CalculateWithSVD(Matrix<double> x, Vector<double> y)
    {
        Vector<double> result = null;

            result = MultipleRegression.Svd(x, y);

            for (int i = 0; i < result.Count; i++)
            {
                var value = result.ElementAt(i);

                if (Double.IsNaN(value) || Double.IsInfinity(value))
                {
                    return null;
                }
            }

        return result;
    }

    public Vector<double> FindBestMRSolution(Matrix<double> x, Vector<double> y)
    {
        Vector<double> result = null;

            result = CalculateWithNormal(x, y);

            if (result != null)
            {
                return result;
            }
            else
            {
                result = CalculateWithSVD(x, y);

                if (result != null)
                {
                    return result;
                }
                else
                {
                    result = CalculateWithQR(x, y);

                    if (result != null)
                    {
                        return result;
                    }
                }
            }

        return result;
    }

public double CalculateTotalSumOfSquares(List<double> dependentVariables)
    {
        double tts = 0;

            for (int i = 0; i < dependentVariables.Count; i++)
            {
                tts += Math.Pow(dependentVariables.ElementAt(i) - dependentVariables.Average(), 2);
            }

        return tts;
    }

Actual Output (Updated results):

test1 = 137431.12889999992 (parallel for loop)

test2 = 7.3770258447689254E- (regular for loop)

---

Epilogue: How to setup an MCVE-compliant testing

This may be a fair way to prepare an indeed fully reproducible setup of an MCVE-code + A/B/C/... DataSET-s, put inside a ready-to-run [IDE & Testing Sandbox, hyperlinked here][1], so that Community members can click a re-run button and focus on root-cause analysis, not on decoding and re-engineering the heaps of incomplete SLOCs.

If this runs for the O/P, it will run for other Community Members, whom the O/P has asked for an answer or help.

Try it online!

My new version of the code:

public double CalculatePredictedRSquared()
    {
        Vector<double> output = CreateVector.Dense(Enumerable.Range(0, 400).Select(i => Convert.ToDouble(i)).ToArray());
        List<double> input1 = new List<double>(Enumerable.Range(0, 400).Select(i => Convert.ToDouble(i)));
        List<double> input2 = new List<double>(Enumerable.Range(200, 400).Select(i => Convert.ToDouble(i)));
        double tss = CalculateTotalSumOfSquares(output.ToList());

        IEnumerable<int> range = Enumerable.Range(0, output.Count);
        var query = range.Select(i => DoIt(i, output, input1, input2));
        var result = 1 - (query.Sum() / tss);
        return result;
    }

    public double DoIt(int i, Vector<double> output, List<double> input1, List<double> input2)
    {
        List<double[]> matrixList = new List<double[]>
                {
                    input1.Where((v, k) => k != i).ToArray(),
                    input2.Where((v, k) => k != i).ToArray()
                };
        var matrixArray = CreateMatrix.DenseOfColumnArrays(matrixList);
        var actualResult = output.ElementAt(i);
        var newVectorArray = CreateVector.Dense(output.Where((v, j) => j != i).ToArray());
        var items = FindBestMRSolution(matrixArray, newVectorArray);
        double estimate = 0;

        if (items != null)
        {
            var y = CalculateYIntercept(matrixArray, newVectorArray, items);
            for (int m = 0; m < 2; m++)
            {
                var coefficient = items[m];

                if (m == 0)
                {
                    estimate += input1.ElementAt(i) * coefficient;
                }
                else
                {
                    estimate += input2.ElementAt(i) * coefficient;
                }
            }
        }
        else
        {
            estimate = 0;
        }

       return Math.Pow(actualResult - estimate, 2);
    }

Answer 1

This whole thing is a dog's breakfast; you should abandon that attempt at parallelism entirely.

Start over. Here's what I want you to do. I want you to write a method DoIt that returns double and takes a int i and whatever other state is required to do a single iteration of the loop.

You will then rewrite your method as follows:

public double CalculatePredictedRSquared()
{
    Vector<double> output = whatever;
    // Whatever other state you need here
    IEnumerable<int> range = Enumerable.Range(0, output.Count);
    var query = range.Select(i => DoIt(i, whatever_other_state));
    return query.Sum();
}

Got it? DoIt is the thing that is inside your loop right now. It must take in i, and output and whatever other vectors you need to pass into it. It must only compute a double -- in this case, the square of the estimate error -- and return that double.

It must be pure: It must not read or write any non-local variable, it must not call any non-pure method, and it must give exactly the same results when given the same inputs, every time. Pure methods are the easiest methods to write, to read, to understand, to test and to parallelize; always try to write pure methods when doing math computations.

Write test cases for DoIt, and test the heck out of it. It's a pure method; you should be able to write lots of test cases. Similarly test any of the pure methods called by DoIt.

Once you are satisfied that DoIt is both correct and pure, then the magic happens. Just change it to:

range.AsParallel().Select...

Then compare the parallel and non-parallel versions. They should produce the same result; if not, then something was impure. Figure out what it was.

Then, verify that the parallel version was faster. If not, then you have failed to do enough work in DoIt to justify parallelism; see https://en.wikipedia.org/wiki/Amdahl%27s_law for details.

Answer 2

A few things:

lock (_lock)
{
   matrixList.Add(input1.Where((v, k) => k != i).ToArray());
   matrixList.Add(input2.Where((v, k) => k != i).ToArray());
}

You're adding items to a collection that is already thread-safe by design, so no need to lock. While List is not thread-safe, it should be OK to read from it concurrently. From the documentation:

It is safe to perform multiple read operations on a List, but issues can occur if the collection is modified while it’s being read. To ensure thread safety, lock the collection during a read or write operation. To enable a collection to be accessed by multiple threads for reading and writing, you must implement your own synchronization.

Also note that matrixList is stored in a local variable; in this case, the collection can't be called from multiple threads because the entire body of the delegate is guaranteed to run on the same thread - it won't be the case that half of the body of the Parallel.For loop will be run on thread A and that the other half will be run on thread B, for example.

Similarly, there's no reason to lock while making changes to estimate1 because it can't possibly be modified from other threads.

Disclaimer: there's no guarantee as to the degree of parallelism of the Parallel.For loop overall. There's not even a guarantee that it will run in parallel at all.

press1 and press2, however, are not local variables, so you do need to synchronize these somehow. (It would be better if you found some way to avoid locking every time, though, because that'll at least partially kill the point of multithreading).

Perhaps most critically, ConcurrentBag is an unordered collection. You don't show all of the operations you're doing on your matrices, but if you're doing matrix multiplication anywhere this could easily cause wrong results. There is no guarantee that matrix multiplication will commute. While A * B = B * A for integers, this is not true in general for matrices. It's quite possible that your logic is subtly dependent on operations occurring in a particular order (and they won't because ConcurrentBag is unordered).