The explicit split-operator algorithm is often used for solving the linear and nonlinear time-dependent Schrödinger equations. However, when applied to certain nonlinear time-dependent Schrödinger equations, this algorithm loses time reversibility and second-order accuracy, which makes it very inefficient. Here, we propose to overcome the limitations of the explicit split-operator algorithm by abandoning its explicit nature. We describe a family of high-order implicit split-operator algorithms that are norm-conserving, time-reversible, and very efficient. The geometric properties of the integrators are proven analytically and demonstrated numerically on the local control of a two-dimensional model of retinal. Although they are only applicable to separable Hamiltonians, the implicit split-operator algorithms are, in this setting, more efficient than the recently proposed integrators based on the implicit midpoint method.