We derive the imaginary part of the potential NRQCD Hamiltonian up to order 1/m^4, when the typical momentum transfer between the heavy quarks is of the order of Lambda_{QCD} or greater, and the binding energy E much smaller than Lambda_{QCD}. We use this result to calculate the inclusive decay widths into light hadrons, photons and lepton pairs, up to O(mv^3 x (Lambda_{QCD}^2/m^2,E/m)) and O(mv^5) times a short-distance coefficient, for S- and P-wave heavy quarkonium states, respectively. We achieve a large reduction in the number of unknown non-perturbative parameters and, therefore, we obtain new model-independent QCD predictions. All the NRQCD matrix elements relevant to that order are expressed in terms of the wave functions at the origin and six universal non-perturbative parameters. The wave-function dependence factorizes and drops out in the ratio of hadronic and electromagnetic decay widths. The universal non-perturbative parameters are expressed in terms of gluonic field-strength correlators, which may be fixed by experimental data or, alternatively, by lattice simulations. Our expressions are expected to hold for most of the charmonium and bottomonium states below threshold. The calculations and methodology are explained in detail so that the evaluation of higher order NRQCD matrix elements in this framework should be straightforward. An example is provided.