New algorithm halves the cost of quaternion matrix multiplication

Quaternions, the four-dimensional cousins of the complex numbers, carry one real part and three imaginary parts, written q0 + q1i + q2j + q3k, with imaginary units obeying the rules i squared equals j squared equals k squared equals minus one, and ij = k while ji = minus k. That non-commutative structure makes them a natural tool for three-dimensional geometry and the computation of rotations, with documented applications in computer graphics, virtual reality, machine vision and robotics, in navigation systems and control theory, and in signal processing. Yet when engineers need to multiply large matrices whose entries are quaternions, mainstream programming languages offer no built-in quaternion type, and established libraries such as Eigen or Boost provide optimized tools for basic rotation and graphics tasks rather than for dense, large-scale quaternion linear algebra.

A new software package called qmatmul, released under the MIT license for Python 3.12 and above, tackles exactly that gap. Its authors describe a fast algorithm for quaternion matrix-matrix multiplication that cuts the number of elementary floating-point multiplications by a factor of two compared with the definition-based formula, regardless of the contents of the input matrices. The package ships with eight implementation variants spanning NumPy with its BLAS backends, Numba-compiled CPU kernels, and hand-written CUDA kernels for NVIDIA GPUs, all callable through a single function, qmm.dot. The code is available on GitHub and installable from PyPI, with documentation hosted for public use.

The mathematical trick behind the speedup begins with a stacked representation. A matrix of quaternions can be viewed as a three-dimensional array of real numbers whose last dimension of size four holds the real and three imaginary parts; stacking those slices produces a purely real block matrix of four times the height. Written this way, the multiplication decomposes into two subproducts. One involves a matrix with a doubly nested diagonal block symmetry, in which blocks of the form U, V, V, U repeat at two levels; the other involves a sparse matrix. Building on an earlier algebraic framework due to Cariow for exploiting block symmetries, the authors show that the symmetric subproduct can be evaluated using only one quarter of the multiplications the direct formula would demand, with the help of Hadamard transformations, which are sign-flipping matrices whose rows are mutually orthogonal, and block-wise permutations. The sparse term likewise needs just a quarter. Combining the two yields a total of 8MNP multiplications for matrices of sizes M by N and N by P, exactly half of the 16MNP required by the direct approach.

The operation counts extend to additions as well. The direct formula performs 16MNP multiplications and 16MNP minus 4MP additions, for a total of roughly 32MNP floating-point operations. The fast algorithm performs 8MNP multiplications plus a slightly larger number of additions, giving about 16MNP operations in total, again asymptotically half. The authors show that the ratio of total operations between the two approaches never exceeds two for any matrix dimensions, approaching that limit only as the matrices grow large. For small problems the lower-order terms matter: in a worked example multiplying five-by-three and three-by-two quaternion matrices, the fast algorithm actually needs slightly more additions than the direct formula, 448 versus 440.

To turn the algorithm into practical code, the team designed a CUDA implementation built from six kernel functions organized into eleven launches. Four kernels implement the stacking, Hadamard transformation, diagonal product and permutation operators; a fifth computes the final difference and negates the real part; a sixth unstacks the result back into quaternion format. The kernels use a tiling strategy that keeps data in fast on-chip shared memory, with a tile size of eight by eight threads found to saturate the streaming multiprocessors across three different GPU architectures. Thread indexing follows a column-major convention inside warps so that memory reads are coalesced, minimizing transactions to global memory. The implementation allocates ten intermediate device-side buffers, whose combined footprint for the largest tested case, three thousand by three thousand by three thousand in double precision, amounts to roughly 2.68 GiB of additional memory; the total footprint including inputs and output is 3.49 GiB, comfortably within typical VRAM capacities.

Benchmarking took place on three systems: a GeForce RTX 3090, a GeForce RTX 5090 and an RTX PRO 5000 Blackwell laptop GPU, each paired with a modern multi-core CPU. Problem sizes ranged from about one million to 2.7 times ten to the tenth power in the product of the matrix dimensions, with timings averaged over one hundred repetitions for smaller cases and ten for the largest, after a warm-up call to trigger just-in-time compilation. The single-threaded, definition-based Numba implementation served as the baseline and was, as expected, the slowest everywhere, despite low-level compilation through the LLVM toolchain; its CPU-parallelized variant ran faster by factors of roughly ten to seventy depending on problem size and processor.

The headline result is the ratio of execution times between the direct and algorithmic versions of the same backend, which isolates the effect of the algorithm from hardware quirks. For single-threaded NumPy computations on the largest problems, the fast algorithm ran about 1.8 times faster, matching the theoretical prediction almost perfectly. The CUDA-based algorithm proved consistently best for large problems, with speed-ups against the naive baseline reaching roughly 15,300 times in single precision and 6,550 times in double precision on the RTX 5090 system. At small problem sizes, however, the picture reverses: the direct CUDA approach, which issues only four kernels and about eleven API calls compared with eleven kernels and twenty-seven calls for the fast algorithm, wins because per-call dispatch overhead dominates when there is little actual computation to do. The crossover appears at problem sizes around ten to the ninth. Device memory tells a similar story, with the fast algorithm using about forty-four percent more total memory than the direct implementation when all dimensions are equal, a price paid for avoiding the direct approach’s structural redundancy.

Numerical accuracy was examined as a complement to the speed results. In double precision, both approaches agree with a reference computation to machine precision, with relative errors around ten to the minus sixteen. In single precision, errors grow with the inner dimension of the product, as expected from the accumulation of independent rounding errors, but stay below two times ten to the minus six across all tested sizes. Contrary to the plausible worry that the Hadamard transformations, scaling and final subtraction might amplify rounding error through cancellation, the fast algorithm was consistently more accurate than the direct formula across every tested configuration and every error metric, never worse in any case. A worst-case tolerance-normalized error remained below one throughout, though the safety margin narrowed with size, and the authors note that extrapolating the growth trend, any fixed absolute tolerance would eventually be exceeded at sufficiently large, untested sizes, an intrinsic property of floating-point summation rather than a flaw of either formula.

To demonstrate a realistic use case, the authors built a proof-of-concept Quaternion Convolutional Neural Network, a class of networks in which multi-channel inputs such as color images can be represented as matrices of quaternions, so that inter-channel structure is preserved rather than degraded by simple summations over channels. Their example network, with two convolutional layers and a dense layer processing batches of thirty-two sixty-four by sixty-four four-channel images, uses the qmatmul package for every layer and applies a component-wise quaternion ReLU activation. The per-layer timings reproduce the crossover pattern: the first convolution, with a small inner dimension, shows no advantage, the second convolution runs about thirty-one percent faster with the algorithm, and the dense layer, with an inner dimension of over a quarter million, runs nearly twice as fast, approaching the asymptotic two-fold gain. A correctness check against the direct approach on the dense-layer output yielded a relative Frobenius error of 1.58 times ten to the minus five, confirming the draft implementation is numerically sound.

Additional comparisons with CuPy, which dispatches to the highly optimized cuBLAS library, showed that for single precision the CuPy-backed fast algorithm outperformed both Numba CUDA versions across all tested shapes, running 4.2 times faster than the Numba algorithm on the largest case. For double precision and skewed, thin rectangular shapes resembling the dense layer of the QCNN example, the Numba-based implementations held the advantage instead. Kernel-level profiling shows the two matrix-multiplication stages dominating at large sizes, reaching over eighty percent of total execution time, while non-computational overheads shrink to a small remainder. The authors see high potential impact for the package, noting that its reliance on NumPy arrays and just-in-time compilation keeps it embedded entirely in the Python ecosystem without requiring integration with the native CUDA C++ API, and they point to quaternion neural networks as a probable future research direction, asking whether such networks can outperform traditional ones at similar parameter counts and whether quaternion matrix products can support the forward and backward computations of their training.

Subject of Research: Fast algorithm and CUDA implementations for quaternion-valued matrix multiplication

Article Title: qmatmul: Fast multiplication of quaternion-valued matrices—algorithm and its implementations for sequential and CUDA computations

Article References: Klȩsk, P., & Cariow, A. (2026). qmatmul: Fast multiplication of quaternion-valued matrices—algorithm and its implementations for sequential and CUDA computations. SoftwareX, 36, Article 103159. https://doi.org/10.1016/j.softx.2026.103159

Image Credits: AI Generated

DOI: 10.1016/j.softx.2026.103159

Keywords: quaternions, matrix multiplication, CUDA, GPU computing, Numba, NumPy, BLAS, linear algebra, Hadamard transform, quaternion neural networks, scientific software, Python