Symbolic Matrix/Operator Multiplication in Mathematica

Question

Let A,Band C be three 4x4 matrices defined as A=x*P+y*Q-R, C=x^(-s)(L0+sL1+s^(2)L2+s^(3)L3). Assuming the actual form of P, Q, R,L0,L1,L2,L3 are not known. Can we compute symbolically, while preserving the order of multiplication (as matrices do not commute in general), the following quantity in Mathematica:$ABA C^TA^TB^TA^T$ (T=Transpose)?

Edit: $x,y,s$ are scalars and not matrices.

Answer 1

You could use TensorExpand:

r = TensorExpand[
    A . B . A . Transpose[C0] . Transpose[A] . Transpose[B] . Transpose[A],
    Assumptions -> (x|y|s) ∈ Reals
] /. Transpose[a_, {2, 1}] :> Transpose[a];
Short[r, 10]

x^(4-s) P.B.P.Transpose[L0].Transpose[P].Transpose[B].Transpose[P]+x^(3-s) y P.B.P.Transpose[L0].Transpose[P].Transpose[B].Transpose[Q]-x^(3-s) P.B.P.Transpose[L0].Transpose[P].Transpose[B].Transpose[R]+x^(3-s) y P.B.P.Transpose[L0].Transpose[Q].Transpose[B].Transpose[P]+<<478>>+s^3 x^-s R.B.R.Transpose[L3].Transpose[R].Transpose[B].Transpose[R]

Answer 2

How about

Q = Array[q, {4, 4}];
R = Array[r, {4, 4}];
B = Array[b, {4, 4}];
P = Array[p, {4, 4}];
A = x*P + y*Q - R
L0 = Array[l0, {4, 4}];
L1 = Array[l1, {4, 4}];
L2 = Array[l2, {4, 4}];
L3 = Array[l3, {4, 4}];
C0 = x^(-s)*(L0 + s*L1 + s^(2)*L2 + s^(3)*L3);
res=A . B . A . Transpose@C0 . Transpose@A . Transpose@B . Transpose@A;
  
Dimensions[res]
  (* {4,4} *)

Result is too large to display fully