We discuss a new theoretical framework for modeling molecular electron densities. Our approach decomposes the total density into contributions from basis function products and then approximates each product using constrained least-squares approximation in a tailored local basis of functions with adjustable non-linear parameters. We show how to solve directly for the expansion coefficients and Lagrange multipliers and present an iterative method to optimize the non-linear parameters. Example products from the Dunning cc-pVTZ basis set are discussed.