locked
MatrixMultiply vs InnerProduct in loop vs double loop RRS feed

  • Question

  • I'm interested to know what the practical differences between the Variable.MatrixMultiply factor and other ways of achieving the same thing. In the code below, I've implemented a small example that uses Variable.MatrixMultiply, and also a version that uses Variable.InnerProduct in a for each block, and then another in a double for e

            public static void Test()
            {
                const int m = 4;
                const int n = 3;
                const int p = 2;
    
                Rand.Restart(0);
                var aInit = RandomGaussianArray(m, n);
                var cObs = RandomArray(m, p);
    
                var funcs = new Func<VariableArray2D<double>, VariableArray2D<double>, VariableArray2D<double>>[]
                {
                    Variable.MatrixMultiply, MatrixMultiplyDoubleLoop
                };
    
                foreach (var func in funcs)
                {
                    var mRange = new Range(m);
                    var nRange = new Range(n);
                    var pRange = new Range(p);
    
                    var a = Variable.Array<double>(mRange, nRange);
                    var b = Variable.Array<double>(nRange, pRange);
                    var tau = Variable.GammaFromShapeAndRate(1, 1);
    
                    a[mRange, nRange] = Variable.GaussianFromMeanAndPrecision(0, 1).ForEach(mRange, nRange);
                    b[nRange, pRange] = Variable.GaussianFromMeanAndPrecision(0, tau).ForEach(nRange, pRange);
    
                    a.InitialiseTo(Distribution<double>.Array(To2D(aInit)));
    
                    var noisy = Variable.Array(Variable.Array<double>(pRange), mRange);
    
                    var clean = func(a, b);
                    noisy[mRange][pRange] = Variable.GaussianFromMeanAndPrecision(clean[mRange, pRange], 10);
    
                    var engine = new InferenceEngine(new VariationalMessagePassing());
    
                    noisy.ObservedValue = cObs;
    
                    Console.WriteLine("Using " + func.Method.Name);
                    Console.WriteLine("Posterior for a:");
                    Console.WriteLine(engine.Infer(a));
                    Console.WriteLine("Posterior for b:");
                    Console.WriteLine(engine.Infer(b));
                    Console.WriteLine();
                }
    
                {
                    // Jagged version. Note that we transpose b here
                    var mRange = new Range(m);
                    var nRange = new Range(n);
                    var pRange = new Range(p);
    
                    var a = Variable.Array(Variable.Array<double>(nRange), mRange);
                    var b = Variable.Array(Variable.Array<double>(nRange), pRange);
                    var tau = Variable.GammaFromShapeAndRate(0.5, 1e-6);
    
                    a[mRange][nRange] = Variable.GaussianFromMeanAndPrecision(0, 1).ForEach(mRange, nRange);
                    b[pRange][nRange] = Variable.GaussianFromMeanAndPrecision(0, tau).ForEach(pRange, nRange);
    
                    a.InitialiseTo(Distribution<double>.Array(aInit));
    
                    var noisy = Variable.Array(Variable.Array<double>(pRange), mRange);
                    var clean = MatrixMultiplyInnerProduct(a, b);
                    noisy[mRange][pRange] = Variable.GaussianFromMeanAndPrecision(clean[mRange][pRange], 10);
    
                    var engine = new InferenceEngine(new VariationalMessagePassing());
    
                    noisy.ObservedValue = cObs;
    
                    Console.WriteLine("Using MatrixMultiplyInnerProduct");
                    Console.WriteLine("Posterior for a:");
                    Console.WriteLine(engine.Infer(a));
                    Console.WriteLine("Posterior for b:");
                    Console.WriteLine(engine.Infer(b));
                    noisy.ClearObservedValue();
                }
            }
    
            public static VariableArray<VariableArray<double>, double[][]> MatrixMultiplyInnerProduct(
                VariableArray<VariableArray<double>, double[][]> a,
                VariableArray<VariableArray<double>, double[][]> b)
            {
                // Note that b is transposed in this version
                var mRange = a.Range;
                var pRange = b.Range;
                var c = Variable.Array(Variable.Array<double>(pRange), mRange);
    
                using (Variable.ForEach(mRange))
                {
                    var v = Variable.Vector(a[mRange]);
    
                    using (Variable.ForEach(pRange))
                    {
                        c[mRange][pRange] = Variable.InnerProduct(b[pRange], v);
                    }
                }
    
                return c;
            }
    
    
            public static VariableArray2D<double> MatrixMultiplyDoubleLoop(
                VariableArray2D<double> a,
                VariableArray2D<double> b)
            {
                var mRange = a.Ranges[0];
                var nRange = a.Ranges[1];
                var pRange = b.Ranges[1];
                var c = Variable.Array<double>(mRange, pRange);
    
                using (Variable.ForEach(mRange))
                {
                    using (Variable.ForEach(pRange))
                    {
                        var products = Variable.Array<double>(nRange);
                        using (Variable.ForEach(nRange))
                        {
                            products[nRange] = a[mRange, nRange] * b[nRange, pRange];
                        }
    
                        c[mRange, pRange] = Variable.Sum(products);
                    }
                }
    
                return c;
            }
    
            public static double[][] RandomArray(int first, int second)
            {
                return Enumerable.Range(0, first).Select(
                    f => Enumerable.Range(0, second).Select(
                        t => (double)Rand.Int(10)).ToArray()).ToArray();
            }
    
            public static Gaussian[][] RandomGaussianArray(int first, int second)
            {
                return Enumerable.Range(0, first).Select(
                    f => Enumerable.Range(0, second).Select(
                        t => new Gaussian(Rand.Normal(), 1)).ToArray()).ToArray();
            }
    
            public static T[,] To2D<T>(T[][] jaggedArray)
            {
                int cols = jaggedArray.Max(row => row.Length);
                var darray = new T[jaggedArray.Length, cols];
                for (var i = 0; i < jaggedArray.Length; i++)
                {
                    for (var j = 0; j < cols; j++)
                    {
                        darray[i, j] = j < jaggedArray[i].Length ? jaggedArray[i][j] : default(T);
                    }
                }
    
                return darray;
            }
        }


    The output is as follows:

    Using MatrixMultiply
    Posterior for a:
    Compiling model...done.
    Iterating: 
    .........|.........|.........|.........|.........| 50
    [0,0] Gaussian(2.41, 0.007603)
    [0,1] Gaussian(6.44e-27, 0.5092)
    [0,2] Gaussian(1.793, 0.01583)
    [1,0] Gaussian(2.904, 0.007603)
    [1,1] Gaussian(-3.592e-27, 0.5092)
    [1,2] Gaussian(-1.557, 0.01583)
    [2,0] Gaussian(1.861, 0.007603)
    [2,1] Gaussian(1.188e-27, 0.5092)
    [2,2] Gaussian(0.1452, 0.01583)
    [3,0] Gaussian(2.819, 0.007603)
    [3,1] Gaussian(-2.556e-27, 0.5092)
    [3,2] Gaussian(-1.207, 0.01583)
    Posterior for b:
    [0,0] Gaussian(3.136, 0.003845)
    [0,1] Gaussian(1.792, 0.003845)
    [1,0] Gaussian(7.041e-27, 0.04819)
    [1,1] Gaussian(-1.427e-26, 0.04819)
    [2,0] Gaussian(0.7418, 0.01369)
    [2,1] Gaussian(-2.375, 0.01369)
    
    Using MatrixMultiplyDoubleLoop
    Posterior for a:
    Compiling model...done.
    Iterating: 
    .........|.........|.........|.........|.........| 50
    [0,0] Gaussian(2.746, 0.008259)
    [0,1] Gaussian(-1.401, 0.0127)
    [0,2] Gaussian(2.612e-43, 0.509)
    [1,0] Gaussian(2.619, 0.008259)
    [1,1] Gaussian(1.851, 0.0127)
    [1,2] Gaussian(1.291e-43, 0.509)
    [2,0] Gaussian(1.89, 0.008259)
    [2,1] Gaussian(0.09776, 0.0127)
    [2,2] Gaussian(1.404e-43, 0.509)
    [3,0] Gaussian(2.599, 0.008259)
    [3,1] Gaussian(1.507, 0.0127)
    [3,2] Gaussian(1.426e-43, 0.509)
    Posterior for b:
    [0,0] Gaussian(3.151, 0.003984)
    [0,1] Gaussian(1.439, 0.003984)
    [1,0] Gaussian(-0.1576, 0.01277)
    [1,1] Gaussian(2.779, 0.01277)
    [2,0] Gaussian(4.438e-43, 0.04823)
    [2,1] Gaussian(6.023e-44, 0.04823)
    
    Using MatrixMultiplyInnerProduct
    Posterior for a:
    Compiling model...done.
    Iterating: 
    .........|.........|.........|.........|.........| 50
    [0] [0] Gaussian(1.662e-58, 1)
        [1] Gaussian(2.09e-58, 1)
        [2] Gaussian(-5.647e-59, 1)
    [1] [0] Gaussian(2.72e-58, 1)
        [1] Gaussian(3.42e-58, 1)
        [2] Gaussian(-9.24e-59, 1)
    [2] [0] Gaussian(1.522e-58, 1)
        [1] Gaussian(1.914e-58, 1)
        [2] Gaussian(-5.171e-59, 1)
    [3] [0] Gaussian(2.582e-58, 1)
        [1] Gaussian(3.246e-58, 1)
        [2] Gaussian(-8.771e-59, 1)
    Posterior for b:
    [0] [0] Gaussian(1.847e-60, 1.999e-06)
        [1] Gaussian(2.322e-60, 1.999e-06)
        [2] Gaussian(-6.275e-61, 1.999e-06)
    [1] [0] Gaussian(1.38e-60, 1.999e-06)
        [1] Gaussian(1.735e-60, 1.999e-06)
        [2] Gaussian(-4.689e-61, 1.999e-06)

    The first two methods (Variable.MatrixMultiply and the double loop) get similar (but not identical) results. The version using Variable.InnerProduct gets completely different (poorly scaled) results.

    Are these differences down to differing schedules in the generated code, or are there actual inferential differences here? Is there any guidance for which method to choose?

    Thursday, June 2, 2016 2:42 PM

Answers

  • The first two methods should have the same set of fixed points.  Note that this model has parameter symmetries so there will generally be many equivalent fixed points.  Also, this model requires many iterations to converge.  The results you posted are for a non-converged inference.  When I run for 2000 iterations, I get similar results for the first two methods but with a parameter symmetry (permutation and sign flip).  The third method is quite different since it is using a different approximating family (VectorGaussian) and a non-conjugate factor (Variable.Vector).  In general, you are better off writing the model with the fewest number of factors, i.e. with MatrixMultiply.
    • Marked as answer by AndyMiles Friday, June 3, 2016 9:26 PM
    Friday, June 3, 2016 12:38 PM
    Owner