An ab initio density matrix algorithm for electronic structure computations of many-electron systems is proposed. The reduced single-electron density matrices are derived by mapping the density functional theory nonlinear optical response functions onto an effective multilevel system. These density matrices are then used as a zeroth order iteration into self-consistent equations whose solution should yield the exact energies and the complete set of (transition and diagonal) single-electron density matrices. Higher order (n electron) density matrices are not computed explicitly. The linear and nonlinear optical response functions may be obtained at a low computational cost. Application is made to constructing an exciton Hamiltonian for molecular aggregates using density matrices of isolated molecules, avoiding electronic structure calculations of the entire aggregate. (C) 2003 American Institute of Physics.