For accurate ab initio description of Rydberg excited states, this study suggests generating appropriate diffuse basis functions by cheap variational optimization of virtual orbitals of the corresponding ion core. By following this approach, dozens of converged correlated lithium Rydberg states, namely, all the states up to 24 2S, 25 2P, 14 2D, 16 2F and 16 2G, not yet achieved via other ab initio approaches, could be obtained at the EOM-CCSD level of theory with compact and mostly state-selective contracted Gaussian basis sets. Despite its small size and Gaussian character, the optimized basis leads to highly accurate excitation energies that differ merely in the order of meV from the reference state-of-the-art explicitly correlated Gaussian method and even surpass Full-CI results on the Slater basis by an order of magnitude.