qmatmul module

This module contains computational approaches for multiplication of quaternion-valued matrices. The default approach qmatmul_algo_numba_cuda is based on the fast algorithm proposed by authors (P. Klęsk, A. Cariow). 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.

Altogether, eight approaches (implementation variants) are provided. They cover several approaches based on NumPy library, thus supported by the underlying BLAS, but also several approaches employing Numba — a just-in-time compiler for Python targeting both CPU and GPU (CUDA). The design of CUDA computations for the proposed algorithm involves: six kernel functions with eleven invocations, suitable usage of tiling and shared memory, and few host-device memory transfers.

Usage example 1

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

\[\begin{split}{\scriptsize \begin{split} \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 \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} \\ &= \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{split} }\end{split}\]

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.

Dependencies

numpy: required for algebraic numerical computations.
numba and numba-cuda: required for just-in-time compilation of LLVM-targeted functions (decorated by @jit) and of CUDA kernels (decorated by @cuda.jit).
threadpoolctl: required for managing single- vs multi-threaded modes for NumPy-based CPU approaches.
To use the qmatmul module, NVIDIA CUDA drivers must be present in the operating system.

Link to project repository

https://github.com/pklesk/quaternions

qmatmul.dot(A, B, approach='algo_numba_cuda', extra_args={'verbose': False})[source]

The main function of qmatmul module — multiplies two quaternion-valued matrices.

Args:

A (numpy.ndarray[:, :, :]):: first factor — a matrix of quaternions.
B (numpy.ndarray[:, :, :]):: second factor — a matrix of quaternions.
approach (str):: choice of computational approach from {"naive_numba_st", "naive_numba_parallel", "direct_numpy_st", "direct_numpy_parallel", "direct_numba_cuda", "algo_numpy_st", "algo_numpy_parallel", "algo_numba_cuda"}, defaults to "algo_numba_cuda".
extra_args (dict):: dictionary of extra arguments to be passed to the invoked function that performs the chosen computational approach.

Returns:

numpy.ndarray[:, :, :]:: matrix-matrix product (quaternion-valued) of A and B.

Raises:

ValueError:: If input matrices are not 3-dimensional, their last dimension is not 4, or if their inner matrix dimensions do not match for multiplication.

qmatmul.qmatrand(M, N, range_min, range_max, dtype=<class 'numpy.float32'>, rounding=False)[source]: Helper function that generates a random three-dimensional numpy.ndarray of shape (M, N, 4), meant to represent a matrix of quaternions (the last dimension stores one real and three imaginary parts).

qmatmul.stack(E)[source]: (for internal use) For a three-dimensional numpy.ndarray of shape (R, S, 4) (matrix of quaternions), returns its stacked representation — the two-dimensional numpy.ndarray of shape (4 * R, S) with slices of imaginary parts stored as blocks of additional rows.

qmatmul.a44(A, a_blocks_signs=array([[1, -1, -1, -1], [1, 1, -1, 1], [1, 1, 1, -1], [1, -1, 1, 1]], dtype=int8), a_blocks_parts=array([[0, 1, 2, 3], [1, 0, 3, 2], [2, 3, 0, 1], [3, 2, 1, 0]], dtype=int8))[source]: (for internal use) For a three-dimensional numpy.ndarray of shape (M, N, 4) (matrix of quaternions), returns the two-dimensional transformation matrix numpy.ndarray of shape (4 * M, 4 * N) made of suitably permuted and signed blocks of real/imaginary parts.

qmatmul.a44_ubar(A, a_blocks_parts=array([[0, 1, 2, 3], [1, 0, 3, 2], [2, 3, 0, 1], [3, 2, 1, 0]], dtype=int8))[source]: (for internal use) For a three-dimensional numpy.ndarray of shape (M, N, 4) (matrix of quaternions), returns the two-dimensional transformation matrix with double diagonal block-symmetry — a numpy.ndarray of shape (4 * M, 4 * N) made of suitably permuted blocks of real/imaginary parts.

qmatmul.a44_lbar(A)[source]: (for internal use) For a three-dimensional numpy.ndarray of shape (M, N, 4) (matrix of quaternions), returns the sparse transformation matrix numpy.ndarray of shape (4 * M, 4 * N) made of suitably placed blocks of real/imaginary parts.

qmatmul.i4_tilde(M, dtype)[source]: (for internal use) For a given size M, returns the block-like identity matrix numpy.ndarray of shape (4 * M, 4 * M) where the top-left block has negative signs.

qmatmul.qmatmul_algolike_numpy(A, B)[source]: (for internal use) Function provided only to check the compliance of computational outcomes of the algorithm (main expression with decomposition into two products).

qmatmul.c4_to_c(C4)[source]: (for internal use) For a two-dimensional numpy.ndarray of shape (4 * M, P) being a result of a matrix-matrix product, returns its unstacked version — the three-dimensional numpy.ndarray of shape (M, P, 4), i.e., a matrix of quaternions.

qmatmul.qmul_numba_float64(q1, q2)[source]: Returns the product of two single quaternions, each stored a 4-element long one-dimensional numpy.array of type float64.

qmatmul.qmul_numba_float32(q1, q2)[source]: Returns the product of two single quaternions, each stored a 4-element long one-dimensional numpy.array of type float32.

qmatmul.qmul_numba_int32(q1, q2)[source]: Returns the product of two single quaternions, each stored a 4-element long one-dimensional numpy.array of type int32.

qmatmul.qmatmul_naive_numba_st_float64(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="naive_numba_st") for input matrices of type float64.

qmatmul.qmatmul_naive_numba_st_float64_job(A, B)[source]: (for internal use) Actual computational job function for qmatmul_naive_numba_st_float64 - returns the result of dot(A, B, approach="naive_numba_st").

qmatmul.qmatmul_naive_numba_parallel_float64(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="naive_numba_parallel") for input matrices of type float64.

qmatmul.qmatmul_naive_numba_parallel_float64_job(A, B)[source]: (for internal use) Actual computational job function for qmatmul_naive_numba_parallel_float64 - returns the result of dot(A, B, approach="naive_numba_parallel").

qmatmul.qmatmul_naive_numba_st_float32(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="naive_numba_st") for input matrices of type float32.

qmatmul.qmatmul_naive_numba_st_float32_job(A, B)[source]: (for internal use) Actual computational job function for qmatmul_naive_numba_st_float32 - returns the result of dot(A, B, approach="naive_numba_st").

qmatmul.qmatmul_naive_numba_parallel_float32(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="naive_numba_parallel") for input matrices of type float32.

qmatmul.qmatmul_naive_numba_parallel_float32_job(A, B)[source]: (for internal use) Actual computational job function for qmatmul_naive_numba_parallel_float32 - returns the result of dot(A, B, approach="naive_numba_parallel").

qmatmul.qmatmul_direct_numpy_parallel(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="direct_numpy_parallel") for input matrices of type either float64 or float32.

qmatmul.qmatmul_direct_numpy_st(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="direct_numpy_st") for input matrices of type either float64 or float32.

qmatmul.had4(E4)[source]: (for internal use) For a stacked matrix of shape (4 * R, S), returns the result (of same shape) of Hadamard transform applied to it, with successive blocks representing: sum of sums, sum of differences, difference of sums, difference of differences.

qmatmul.matmuldiag(E4, F4, factor)[source]: (for internal use) For two stacked matrices of shape (4 * R, S) each, returns their diagonal product (same shape), with successive blocks being standard products of corresponding real and imaginary slices.

qmatmul.permute(E4, permutation)[source]: (for internal use) For a stacked matrix of shape (4 * R, S), returns its block-wise permuted version with order of blocks specified by the argument permutation.

qmatmul.qmatmul_algo_numpy_parallel(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="algo_numpy_parallel") for input matrices of type either float64 or float32.

qmatmul.qmatmul_algo_numpy_st(A, B, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="algo_numpy_st") for input matrices of type either float64 or float32.

qmatmul.qmatmul_direct_numba_cuda_float64(A, B, tile_size=8, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="direct_numba_cuda_float64") for matrices of type float64.

qmatmul.stack_numba_cuda_job_float64(E, E4)[source]: (CUDA kernel for internal use) For a three-dimensional device array float64[:, :, :] of shape (R, S, 4) named E (matrix of quaternions), prepares its stacked representation — the two-dimensional device array float64[:, :] of shape (4 * R, S) named E4 with slices of imaginary parts stored as blocks of additional rows, by making each thread write a single entry of E4.

qmatmul.a44_numba_cuda_job_float64(A, A44)[source]: (CUDA kernel for internal use) For a three-dimensional device array float64[:, :, :] of shape (M, N, 4) named A (matrix of quaternions), prepares the two-dimensional transformation matrix float64[:, :] named A44 of shape (4 * M, 4 * N) made of suitably permuted and signed blocks of real/imaginary parts, by making each thread write a single entry of A44.

qmatmul.matmul_numba_cuda_job_float64(A44, B4, C4)[source]: (CUDA kernel for internal use) For a two-dimensional device array float64[:, :] of shape (4 * M, 4 * N) named A44 (transformation matrix) and a two-dimensional device array float64[:, :] of shape (4 * N, P) named B4, computes their matrix-matrix product named C4 being a two-dimensional device array float64[:, :] of shape (4 * M, P), using shared memory and tiling, with each thread responsible for computing a single row-column inner product being an entry in C4.

qmatmul.c4_to_c_numba_cuda_job_float64(C4, C)[source]: (CUDA kernel for internal use) For a two-dimensional device array float64[:, :] of shape (4 * M, P) named C4 being a result of a matrix-matrix product, prepares its unstacked version — the three-dimensional device array float64[:, :, :] of shape (M, P, 4) named C, by making each thread write a single entry of C4.

qmatmul.qmatmul_direct_numba_cuda_float32(A, B, tile_size=8, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="direct_numba_cuda_float32") for matrices of type float32.

qmatmul.stack_numba_cuda_job_float32(E, E4)[source]: (CUDA kernel for internal use) For a three-dimensional device array float32[:, :, :] of shape (R, S, 4) named E (matrix of quaternions), prepares its stacked representation — the two-dimensional device array float32[:, :] of shape (4 * R, S) named E4 with slices of imaginary parts stored as blocks of additional rows, by making each thread write a single entry of E4.

qmatmul.a44_numba_cuda_job_float32(A, A4)[source]: (CUDA kernel for internal use) For a three-dimensional device array float32[:, :, :] of shape (M, N, 4) named A (matrix of quaternions), prepares the two-dimensional transformation matrix float32[:, :] named A44 of shape (4 * M, 4 * N) made of suitably permuted and signed blocks of real/imaginary parts, by making each thread write a single entry of A44.

qmatmul.matmul_numba_cuda_job_float32(A44, B4, C4)[source]: (CUDA kernel for internal use) For a two-dimensional device array float32[:, :] of shape (4 * M, 4 * N) named A44 (transformation matrix) and a two-dimensional device array float32[:, :] of shape (4 * N, P) named B4, computes their matrix-matrix product named C4 being a two-dimensional device array float32[:, :] of shape (4 * M, P), using shared memory and tiling, with each thread responsible for computing a single row-column inner product being an entry in C4.

qmatmul.c4_to_c_numba_cuda_job_float32(C4, C)[source]: (CUDA kernel for internal use) For a two-dimensional device array float32[:, :] of shape (4 * M, P) named C4 being a result of a matrix-matrix product, prepares its unstacked version — the three-dimensional device array float32[:, :, :] of shape (M, P, 4) named C, by making each thread write a single entry of C4.

qmatmul.qmatmul_algo_numba_cuda_float64(A, B, tile_size=8, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="algo_numba_cuda") for input matrices of type float64.

qmatmul.had4_numba_cuda_job_float64(E4, H4E4)[source]: (CUDA kernel for internal use) For a two-dimensional device array float64[:, :] of shape (4 * R, S), computes the Hadamard transform applied to it, resulting in a two-dimensional device array float64[:, :] named H4E4, with each thread responsible for computing four entries of H4E4 representing: a sum of sums, a sum of differences, a difference of sums, a difference of differences (for suitable block pairs).

qmatmul.matmuldiag_numba_cuda_job_float64(E4, F4, factor, G4)[source]: (CUDA kernel for internal use) For a two-dimensional device array float64[:, :] of shape (4 * R, S) named E4 and a two-dimensional device array float64[:, :] of shape (4 * S, T) named F4, computes their diagonal matrix-matrix product named G4 (successive blocks being standard products of corresponding real and imaginary slices) being a two-dimensional device array float64[:, :] of shape (4 * R, T), using shared memory and tiling, with each thread responsible for computing a single row-column inner product being an entry in G4.

qmatmul.permute_numba_cuda_job_float64(E4, permutation, E4_permuted)[source]: (CUDA kernel for internal use) For a two-dimensional device array float64[:, :] of shape (4 * R, S) named E4, prepares its block-wise permuted version named E4_permuted (of same shape and type) with order of blocks specified by the argument permutation given as int8[:], using shared memory and tiling, with each thread responsible for writing four entries of E4_permuted.

qmatmul.matsub_numba_cuda_job_float64(C4_left, C4_right, C4)[source]: (CUDA kernel for internal use) For two two-dimensional device arrays float64[:, :] named C4_left and C4_right, of shape (4 * M, P) each, computes their difference named C4 (of same shape and type), using shared memory and tiling, with each thread responsible for computing and writing a single entry in C4.

qmatmul.qmatmul_algo_numba_cuda_float32(A, B, tile_size=8, verbose=False)[source]: (to be invoked by the main function qmatmul.dot) Returns the result of dot(A, B, approach="algo_numba_cuda") for input matrices of type float32.

qmatmul.had4_numba_cuda_job_float32(E4, H4E4)[source]: (CUDA kernel for internal use) For a two-dimensional device array float32[:, :] of shape (4 * R, S), computes the Hadamard transform applied to it, resulting in a two-dimensional device array float32[:, :] named H4E4, with each thread responsible for computing four entries of H4E4 representing: a sum of sums, a sum of differences, a difference of sums, a difference of differences (for suitable block pairs).

qmatmul.matmuldiag_numba_cuda_job_float32(E4, F4, factor, G4)[source]: (CUDA kernel for internal use) For a two-dimensional device array float32[:, :] of shape (4 * R, S) named E4 and a two-dimensional device array float32[:, :] of shape (4 * S, T) named F4, computes their diagonal matrix-matrix product named G4 (successive blocks being standard products of corresponding real and imaginary slices) being a two-dimensional device array float32[:, :] of shape (4 * R, T), using shared memory and tiling, with each thread responsible for computing a single row-column inner product being an entry in G4.

qmatmul.permute_numba_cuda_job_float32(E4, permutation, E4_permuted)[source]: (CUDA kernel for internal use) For a two-dimensional device array float32[:, :] of shape (4 * R, S) named E4, prepares its block-wise permuted version named E4_permuted (of same shape and type) with order of blocks specified by the argument permutation given as int8[:], using shared memory and tiling, with each thread responsible for writing four entries of E4_permuted.

qmatmul.matsub_numba_cuda_job_float32(C4_left, C4_right, C4)[source]: (CUDA kernel for internal use) For two two-dimensional device arrays float32[:, :] named C4_left and C4_right, of shape (4 * M, P) each, computes their difference named C4 (of same shape and type), using shared memory and tiling, with each thread responsible for computing and writing a single entry in C4.