In
y_hat = gamma + np.einsum("jik,kl->ij", X, beta) + u
gamma and u are pymc distribution objects but np.einsum will return an ndarray. So to make y_hat an object in the pm.Model() context you should use
y_hat = pm.Deterministic("y_hat", gamma + np.einsum("jik,kl->ij", X, beta) + u)
However, if this doesn’t work (I haven’t checked), you may need to use a tensor algebra for the np.einsum calculation, such as this example.