In this paper, a quaternion boundary element method (BEM) is proposed to solve the magnetostatic problem. The present quaternion-valued BEM is developed by discretizing the quaternion-valued boundary integral equation (BIE). The quaternion-valued BIE can be seen as an extension of the generalized complex variable BIE in 3-D space. In other words, quaternion algebra is an extension of the complex variable in 3-D space. To derive quaternion-valued BIEs, the quaternion-valued Stokes’ theorem is utilized. The quaternion-valued BIEs are noted for singularity, which exists in the sense of the Cauchy principal value (CPV). An analytical scheme is developed to evaluate the CPV by introducing a simple quaternion-valued harmonic function. For the domain points close to the boundary, some sorts of analogous, nearly singular, so-called “numerical boundary layer” phenomena appear and are remedied by using a similar analytic evaluation. The quaternion BEM features the oriented surface element, combining the unit outward normal vector with the ordinary surface element. It is noted that all derivations are done in quaternion algebra. In addition, quaternion algebra is more flexible than vector algebra in solving some 3-D problems from the point view of algebraic space. In deriving BIEs for exterior fields, the conditions at infinity for the quaternion-valued functions are carefully examined. Later, a magnetic sphere in a uniform magnetic field is considered. This problem is a magnetostatic problem of coupled exterior and interior magnetostatic fields. Finally, we apply the present approach to solve the magnetostatic problem. By comparing with exact solutions, the validity of the present approach is checked.