qmatmul 1.0.0

Fast multiplications of quaternion-valued matrices.

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

We present an algorithm for fast multiplication of matrices whose elements are quaternions - hypercomplex numbers consisting of one real and three imaginary parts. The number of elementary floating-point multiplications involved in the algorithm is reduced twice with respect to the definition-based formula, regardless of the input matrices. This is owed to a suitable representation and decomposition into two products, one of which takes advantage of certain diagonal symmetry properties, the other of sparsity.

The qmatmul package is suitable for Python's ecosystem. Altogether, we provide 8 implementation variants of matrix-matrix multiplication for quaternion-valued inputs. The variants cover several approaches based on NumPy, thus supported by BLAS, but also several approaches employing Numba - a just-in-time compiler targeting both CPU and GPU (CUDA). Our design of CUDA computations for the proposed algorithm involves: 6 kernel functions with 11 invocations, tiling and shared memory, and few host-device memory transfers.

Selected kernels - flows of CUDA computations

Speed-ups

Installation

pip install TODO

Note: for further usage, NVIDIA CUDA drivers must be present in the operating system.

Usage example 1

Suppose one would like to multiply the following matrices of quaternions:

$$ \begin{aligned} {\tiny \begin{pmatrix} 2i +j & -1 +i -j + 4k & 5 -i +3j \ 1 +4i -4j -4k & -5 +3i +3j +4k & 5 +3i +3k \ -4 +i -4j + 4k & -i -2j +3k & i -5j +k \ 1 +i +4j + 2k& -1 -i +2j -4k & 2 +2i -3j -4k \ -2 -i +j -k & 5 -4i -3j -3k & 2 -2i -3k \end{pmatrix} } {\cdot} {\tiny \begin{pmatrix} -3 -4i + 2j -4k & -3 -i +3j -4k \ 3 - 4i +5j & 5 +i +2j -5k\ -2 -4i -2j -4k & -2 -i -4j +2k \end{pmatrix} } {\tiny =} {\tiny \begin{pmatrix} 4 -53i -39j +15k & 16 -6i -17j +52k \ 1 -3i +4j +7k & -30 +3i +30j +56k \ 44 +23i -16j -38k & 40 -4i -5j +7k \ -34 -14i +29j -17k & -40 -62i -10j -21k \ -14 -18i +21j -26k & 18 -41j -18k \end{pmatrix}. } \end{aligned} $$

With qmatmul module installed, one can write:

import qmatmul as qmm
import numpy as np
import time

print("QMATMUL EXAMPLE...")
A = np.array([
    [[ 0,  2,  1,  0], [-1,  1, -1,  4], [ 5, -1,  3,  0]],
    [[ 1,  4, -4, -4], [-5,  3,  3,  4], [ 5,  3,  0,  3]],
    [[-4,  1, -4,  4], [ 0, -1, -2, -3], [ 0,  1, -5,  1]],
    [[ 1,  1,  4,  2], [-1, -1,  2, -4], [ 2,  2, -3, -4]],
    [[-2, -1,  1, -1], [ 5, -4, -3, -3], [ 2, -2,  0, -3]]
    ])
B = np.array([
    [[-3, -4,  2, -4], [-3, -1,  3, -4]], 
    [[ 3, -4,  5,  0], [ 5,  1,  2, -5]],
    [[-2, -4, -2, -4], [-2, -1, -4,  2]]
    ])
t1 = time.time()
C = qmm.dot(A, B)
t2 = time.time()
print(f"RESULT -> C:")
print(C)
print(f"QMATMUL EXAMPLE DONE. TIME OF qmm.dot: {t2 - t1:.6f} s.")

Running the code above produces the following output:

QMATMUL EXAMPLE...
RESULT -> C:
[[[  4. -53. -39.  15.]
  [ 16.  -6. -17.  52.]]

 [[  1.  -3.   4.   7.]
  [-30.   3.  30.  56.]]

 [[ 44.  53.   8. -56.]
  [ 10.   8. -11. -23.]]

 [[-34. -14.  29. -17.]
  [-40. -62. -10. -21.]]

 [[-14. -18.  21. -26.]
  [ 18.   0. -41. -18.]]]
QMATMUL EXAMPLE DONE. TIME OF qmm.dot: 0.001703 s.

Usage example 2 (large arguments)

In the example below, two large random matrices with quaternions are multiplied.

import qmatmul as qmm
import numpy as np
import time

print("QMATMUL EXAMPLE (LARGE ARGUMENTS)...")
M, N, P = 1000, 3000, 2000
np.random.seed(0)
A = np.random.rand(M, N, 4) # M x N matrix of quaternions
B = np.random.rand(N, P, 4) # N x P matrix of quaternions
t1 = time.time()
C = qmm.dot(A, B)
t2 = time.time()
print(f"RESULT FRAGMENT -> C[:3, :3]:")
print(C[:3, :3])
print(f"QMATMUL EXAMPLE (LARGE ARGUMENTS) DONE. TIME OF qmm.dot: {t2 - t1:.6f} s.")

The result is computed fast:

QMATMUL EXAMPLE (LARGE ARGUMENTS)...
RESULT FRAGMENT -> C[:3, :3]:
[[[-1528.6768062   1482.01579334  1482.64352966  1469.29588132]
  [-1474.39555984  1485.26884228  1485.81638515  1486.03433938]
  [-1459.21558118  1487.12054919  1468.20437822  1461.65795584]]

 [[-1502.50516591  1465.24326325  1494.5907814   1503.74503685]
  [-1472.32677282  1493.41728185  1487.01751106  1506.41597882]
  [-1460.42240679  1488.1077977   1474.34164258  1504.36179871]]

 [[-1558.07311493  1475.34714296  1475.731701    1495.41197899]
  [-1528.83559572  1476.28138065  1474.62854662  1504.35168932]
  [-1508.66904595  1501.06439708  1459.07714887  1471.97545486]]]
QMATMUL EXAMPLE (LARGE ARGUMENTS) DONE. TIME OF qmm.dot: 0.151431 s.

Choice of approach

An additional optional argument approach, e.g., qmm.dot(A, B, approach="..."), allows the user to select one of the following eight computational approaches:

approach	target, mode	description
"naive_numba_st"	CPU, single-threaded	naive single-threaded implementation of definition-based formula consisting of three nested loops, low-level compiled via Numba and LLVM
"naive_numba_parallel"	CPU, multi-threaded	as above, but parallelized over CPU cores
"direct_numpy_st"	CPU, single-threaded	direct implementation of formula based on the transformation matrix and stacked representation (followed by unstack), using NumPy/BLAS
"direct_numpy_parallel"	CPU, multi-threaded	as above, but allowing for CPU parallelization supported by NumPy/BLAS
"direct_numba_cuda"	GPU, multi-threaded	direct implementation of formula based on the transformation matrix and stacked representation (followed by unstack), using CUDA, compiled via Numba and LLVM to PTX/SASS
"algo_numpy_st"	CPU, single-threaded	implementation of the proposed fast algorithm, using NumPy/BLAS
"algo_numpy_parallel"	CPU, multi-threaded	as above, but allowing for CPU parallelization supported by NumPy/BLAS
"algo_numba_cuda"	GPU, multi-threaded	implementation of the proposed fast algorithm, using CUDA, compiled via Numba and LLVM to PTX/SASS

The default setting is "algo_numba_cuda".

Documentation

Developer documentation of the project is accessible at: https://pklesk.github.io/quaternions.

License

This project is licensed under the MIT License.

Acknowledgments and credits

• NumPy: the fundamental package for scientific computing with Python.
• Numba: a high-performance just-in-time Python compiler.