The inverse robot kinematics problem has been extensively studied by many workers, but still some problems related to the complexity and strong nonlinear of the inverse kinematics process need suitable heuristic and adhoc techniques and simplifications. A novel method based on uniqueness domains notion was proposed. With using the boundary confirmed by robot’s Jacobian matrix determinant equal to zero, the joint space of the robot was divided into uniqueness domains with the same number of solutions as the inverse kinematics, and the boundary of each uniqueness domain was used as a constraint condition. Then the inverse kinematics solution in the uniqueness domain was transformed into the constrained optimization of the CMA-ES algorithm, the initial mean points of the CMA-ES algorithm in the uniqueness domain were optimized by using the characteristics of the uniform distribution of the good point set. The application of the presented method was described in detail by solving the inverse kinematics multiple solution of the 6R industrial robot, and comparing with the numerical method on Qianjiang No.1 industrial robot and the KUKA humanoid manipulator. The simulation results showed that under the precondition of accuracy requirement, the proposed method had a faster solution speed. For the industrial 6R robot, the average solution time of CMA-ES algorithm was about 5.1ms/time, and that of numerical method was about 7.5ms/time, and KUKA humanoid manipulator, the average solution time of inverse solutions was about 18.9ms/time, and the average solution time of numerical method was about 54.8ms/time. The presented CMA-ES algorithm stabilized the position tracking accuracy of both robots at 10-6mm level.