scales as the 6th power of the basis-set size, significantly worse than the asymptotic 3rd power scaling of Hartree-Fock (HF) theory. Here we introduce a new natural-orbital-functional theory, OP-NOFT, in which the basic variables are the NSOs, their occupation numbers, and their joint occupation probabilities (OP). The latter allow us to represent the 2-DM accurately and transcend the limitations of the 1-DM theories. Its general form contains single-NSO through 4-NSO joint-occupation probabilities and scales as the 5th power of the basis-set size. Its simplest formulation, for seniority 0, OP-NOFT-0, corresponds to doubly-occupied configuration interaction (DOCI) [@Weinhold1967]. It contains only single- and 2-natural-orbital (NO) OPs and retains the 3rd power scaling of HF energy-functional minimization, albeit with a higher prefactor. It describes the dissociation of simple diatomic molecules and multi-atom chains with accuracy comparable to that of DOCI, which uses a compact basis of Slater determinants (SD) but retains exponential scaling. OP-NOFT-0 is powerful at high correlation, i.e. for static correlation at intermediate and large interatomic separations where HF fails due to the multi-reference character of the ground-state wavefunction. There, OP-NOFT-0 outperforms HF, DFT and quantum-chemistry methods such as (single-reference) coupled cluster with single, double and perturbative triple electron-hole excitations (CCSD(T)), 