A = {{5, 4, 1}, {4, 7, 3}, {1, 3, 5}};
 B = {{3, 2, -1}, {2, 5, 1}, {-1, 1, 3}};
 Xsam = RandomVariate[MultinormalDistribution[{0, 0, 0}, B], 1000];
 Ysam = RandomVariate[MultinormalDistribution[{0, 0, 0}, A], 1000];

M0 = Array[Subscript[a, ##] &, {3, 3}];

r[X_,Y_,M_]:=Exp[X.M.Transpose[Y]];

pX=1/1000;
pY=1/1000;
r1=r[Xsam,Ysam,M0]

g[M0_]:=Take[RecurrenceTable[{b[n + 1] == pY/(Transpose[r1].a[n]),
a[n + 1] == pX/(r1.b[n]), a[1] == Exp[-ConstantArray[1, Length [Xsam]]],
b[1] == Exp[-ConstantArray[1, Length [Xsam]]]}, {a, b}, {n, 1,
10}], {10}]