Yes, they are fully factorized. If you don't want this, you need to rewrite the model to use inner products of Vectors (Variable<Vector> types and Variable.Inner), like so:
Range N = data.Range0.Named("N");
Range D = data.Range1.Named("D");
Range d = new Range(10).Named("d");
Variable<Gaussian> priorMu = Variable.New<Gaussian>().Named("PriorMu");
Variable<Gamma> priorPrec = Variable.New<Gamma>().Named("PriorPrec");
priorMu.ObservedValue = Gaussian.FromMeanAndPrecision(0.0, 1.0);
priorPrec.ObservedValue = Gamma.FromShapeAndScale(2.0, 1 / 2.0);
// Mixing matrix
VariableArray<Vector> W = Variable.Array<Vector>(D).Named("W");
var Alpha = Variable.WishartFromShapeAndScale(2, PositiveDefiniteMatrix.IdentityScaledBy(d.SizeAsInt, .5)).Named("Alpha");
W[D] = Variable.VectorGaussianFromMeanAndPrecision(Vector.Zero(d.SizeAsInt), Alpha).ForEach(D);
// Latent variables
var Z = Variable.Array<Vector>(N).Named("Z");
Z[N] = Variable.VectorGaussianFromMeanAndPrecision(Vector.Zero(d.SizeAsInt), PositiveDefiniteMatrix.Identity(d.SizeAsInt)).ForEach(N);
// Multiply
var ZtimesW = Variable.Array<double>(N, D).Named("ZTimesW");
ZtimesW[N, D] = Variable.InnerProduct(Z[N], W[D]);
// Bias
VariableArray<double> mu = Variable.Array<double>(D).Named("mu");
mu[D] = Variable.Random<double, Gaussian>(priorMu).ForEach(D);
VariableArray2D<double> ZWplusMu = Variable.Array<double>(N, D).Named("ZWplusMu");
ZWplusMu[N, D] = ZtimesW[N, D] + mu[D];
// Observation noise
VariableArray<double> pi = Variable.Array<double>(D).Named("pi");
pi[D] = Variable.Random<double, Gamma>(priorPrec).ForEach(D);
// The observation
data[N, D] = Variable.GaussianFromMeanAndPrecision(ZWplusMu[N, D], pi[D]);
Also see
this post.
