An iterative algorithm is developed for calculating nonlinear optical polarizabilities using a series of generalized sum rules that resemble the Lanczos algorithm and connect spectral moments of the driven single-electron density matrix to ground state charge distributions and bonding network, The size scaling and saturation of off-resonant polarizabilities (up to seventh order) of polyacetylene oligomers with up to 300 carbon atoms is analyzed in terms of collective electronic oscillators. Simple analytical expressions for size and bond-length alternation dependence of off-resonant polarizabilities are derived using a single-oscillator approximation.