Consider the following code

import sympy as sym

A = sym.MatrixSymbol('A', 5, 3)
E = sym.MatrixSymbol('E', 5, 3)
sym.expand((A*A.T + E*E.T)*(A*A.T + E*E.T))

Sympy just returns (A*A.T + E*E.T)**2. Is it possible to force expansion of the square?