# MatrixMultiply vs InnerProduct in loop vs double loop

• ### 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