We study several approaches to orbital optimization in selected configuration interaction plus perturbation theory (SCI+PT) methods, and test them on the ground and excited states of three molecules using the semistochastic heatbath configuration interaction (SHCI) method. We discuss the ways in which the orbital optimization problem in SCI resembles and differs from that in complete active space self-consistent field (CASSCF). Starting from natural orbitals, these approaches divide into three classes of optimization methods according to how they treat coupling between configuration interaction (CI) coefficients and orbital parameters, namely uncoupled, fully coupled, and quasi-fully coupled methods. We demonstrate that taking the coupling into account is crucial for fast convergence and recommend two quasi-fully coupled methods for such applications: accelerated diagonal Newton and Broyden-Fletcher-Goldfarb-Shanno (BFGS).