We derive and implement a new way of solving coupled cluster equations with lower computational scaling. Our method is based on the decomposition of both amplitudes and two electron integrals, using a combination of tensor hypercontraction and canonical polyadic decomposition. While the original theory scales as O(N6) with respect to the number of basis functions, we demonstrate numerically that we achieve sub-millihartree difference from the original theory with O(N4) scaling. This is accomplished by solving directly for the factors that decompose the cluster operator. The proposed scheme is quite general and can be easily extended to other many-body methods.