<p>Quaternion optimization has attracted significant interest due to its broad applications, including color face recognition, video compression, and signal processing. Despite the growing literature on quadratic and matrix quaternion optimization, to the best of our knowledge, the study on quaternion polynomial optimization still remains blank. In this paper, we introduce the first investigation into this fundamental problem, and focus on the sphere-constrained homogeneous polynomial optimization over the commutative quaternion domain, which includes the best rank-one tensor approximation as a special case. Our study proposes a polynomial-time randomized approximation algorithm that employs tensor relaxation and random sampling techniques to tackle this problem. Theoretically, we prove an approximation ratio for the algorithm providing a worst-case performance guarantee.</p>