>> X_{ijk} = \sum_{l} A_{il}*B_{jl}*C_{kl} > > (A[:,newaxis,newaxis]*B[newaxis,:,newaxis]*C[newaxis,newaxis,:]).sum(axis=-1) Thanks for the quick solution and practical exercise in broadcasting! :-) However, this creates a temporary 4-array, right? Is there a way of avoiding this memory requirement? - Hagen