This paper introduces biquaternion eigen-decomposition theory (via Peirce decomposition) with respect to a selected quaternion with a non-zero vector part. The eigen-decomposition allows evaluation of polynomials and power series with real coefficients as functions of quaternions. This extension of analytic functions to functions of quaternions requires only standard complex function evaluation. The theory also applies to quaternion rotations. The theory uses biquaternion calculations indicated by matrix methods via the algebraic isomorphism between Hamilton’s biquaternions and appropriate $$4\times 4$$ complex matrices. The isomorphism preserves algebraic structure. In particular, the left and right biquaternion multiplication by the selected quaternion maps to left and right matrix multiplication, respectively. This unifies the representation of the left and right quaternion multiplication as a linear map into a single matrix form. This matrix, as a linear operator, acts on matrices, so that the eigenvectors have matrix form that maps into the biquaternions. Use of an alternate quaternion basis results in a similarity transform of the representation matrix, preserving eigenvalues across change of basis. The similarity transform allows simple eigenvector calculation. The matrix for the selected quaternion has two identical, complex conjugate pairs of eigenvalues. Each pair corresponds to two complex conjugate pairs of eigenvector biquaternions, an idempotent pair and a nilpotent pair. Idempotent and nilpotent eigenvectors correspond to the commuting and non-commuting parts, respectively, of quaternion multiplication.