  We consider the problem of constructing a canonical polyadic (CP)
decomposition for a tensor network, rather than a single tensor. We illustrate
how it is possible to reduce the complexity of constructing an approximate CP
representation of the network by leveraging its structure in the course of the
CP factor optimization. The utility of this technique is demonstrated for the
order-4 Coulomb interaction tensor approximated by 2 order-3 tensors via an
approximate generalized square-root (SQ) factorization, such as density fitting
or (pivoted) Cholesky. The complexity of constructing a 4-way CP decomposition
is reduced from $\mathcal{O}(n^4 R_\text{CP})$ (for the non-approximated
Coulomb tensor) to $\mathcal{O}(n^3 R_\text{CP})$ for the SQ-factorized tensor,
where $n$ and $R_\text{CP}$ are the basis and CP ranks, respectively. This
reduces the cost of constructing the CP approximation of 2-body interaction
tensors of relevance to accurate many-body electronic structure by up to 2
orders of magnitude for systems with up to 36 atoms studied here. The full
4-way CP approximation of the Coulomb interaction tensor is shown to be more
accurate than the known approaches utilizing CP-decomposed SQ factors (also
obtained at the $\mathcal{O}(n^3 R_\text{CP})$ cost), such as the algebraic
pseudospectral and tensor hypercontraction approaches. The CP decomposed SQ
factors can also serve as a robust initial guess for the 4-way CP factors.
