We propose a structure-exploiting conjugate-gradient (CG)–type algorithm for solving the matrix equation \(AXB = C\) over the generalized quaternions. The proposed method is developed from an idea of operating the linear map \(\mathcal {K}(X) = AXB\) directly on the matrix space. This enables a matrix-free Krylov subspace implementation that avoids the explicit construction of the associated large-scale Kronecker matrix. By exploiting the intrinsic component-wise structure of quaternion matrices, the procedure performs all computations through matrix–matrix multiplications, leading to significant reductions in computational complexity and memory requirements. Finite-step convergence of the algorithm is established under suitable positive-definiteness assumptions. The proposed framework naturally includes real- and complex-valued matrix equations, the Hamilton quaternion case, and other quaternion algebras as special cases. Numerical experiments confirm the efficiency, robustness, and scalability of the proposed algorithm compared with existing methods.
