Is it this what you are looking for?
bmix = {b1, a1, b3, a3};
amix = {a2, b2, a4, b4};
tT = Table[t[i, j], {i, 4}, {j, 4}];
tT // MatrixForm
lsg1 = Flatten[Solve[bmix == tT. amix // Thread, {b1, b2, b3, b4}]] //
FullSimplify;
r1 = {b1, b2, b3, b4} /. lsg1;
aa = {a1, a2, a3, a4};
sM = Table[Coefficient[r1[[i]], aa[[j]]], {i, 4}, {j, 4}];
sM // MatrixForm
r1 - sM.aa // FullSimplify