In this paper we provide a detailed analysis of the preconditioned steepest descent (PSD) iteration solver for a convex splitting numerical scheme to the Cahn-Hilliard equation with variable mobility function. In more details, the convex-concave decomposition is applied to the energy functional, which in turn leads to an implicit treatment for the nonlinear term and the surface diffusion term, combined with an explicit update for the expansive concave term. In addition, the mobility function, which is solution-dependent, is explicitly computed, which ensures the elliptic property of the operator associated with the temporal derivative. The unique solvability of the numerical scheme is derived following the standard convexity analysis, and the energy stability analysis could also be carefully established. On the other hand, an efficient implementation of the numerical scheme turns out to be challenging, due to the coupling of the nonlinear term, the surface diffusion part, and a variable-dependent mobility elliptic operator. Since the implicit parts of the numerical scheme are associated with a strictly convex energy, we propose a preconditioned steepest descent iteration solver for the numerical implementation. Such an iteration solver consists of a computation of the search direction (involved with a Poisson-like equation), and a one-parameter optimization over the search direction, in which the Newton’s iteration becomes very powerful. In addition, a theoretical analysis is applied to the PSD iteration solver, and a geometric convergence rate is proved for the iteration. A few numerical examples are presented to demonstrate the robustness and efficiency of the PSD solver.