You can simply perform a cross product of matrics.
can be expressed using Einstein notation as
Code implementation of this problem using NumPy:
numpy.einsum('az,bz,cz,dz -> abcd', A, B, C, D)
so in your case, you would use
numpy.einsum('az,bz->ab', P.U, P.U)
or, in your 3-matrix case
numpy.einsum('az,bz,cz->abc', P.U, P.U, P.U)
sktensor.ktensor.ktensor also have a method totensor() that does exactly this:
np.allclose(np.einsum('az,bz->ab', P.U, P.U), P.totensor())
Hope this answer helps.