Ball matrix multiplication

The product $C = AB$ of an $M \times N$ matrix $A$ and an $N \times P$ matrix $B$ with ball entries $A = [A_m \pm A_r]$, $B = [B_m \pm B_r]$ must satisfy $\tilde A \tilde B \in C$ for all matrices $\tilde A \in A$, $\tilde B \in B$. Since

the entrywise bound

holds (with absolute values and inequalities taken entrywise), so it is enough to compute an enclosure of the midpoint product $A_m B_m$ and an upper bound for the radius product on the right. The classical algorithm computes each entry as a dot product in ball arithmetic (arb_dot, see Dot products of balls), which does exactly this entry by entry with $O(N)$ multiplications of $p$-bit numbers. For large matrices, both parts can instead be done with fast matrix multiplication, as observed by Rump for midpoint-radius interval matrices in floating-point arithmetic [Rum1999]; FLINT does this in arb_mat_mul_block, with exact integer products for the midpoints.

Midpoints: exact integer matrix products#

Every midpoint is a binary floating-point number $a = u \cdot 2^e$ with an odd integer $u$ of at most $p$ bits. If all entries in row $i$ of $A_m$ are multiples of $2^{s_i}$ and have absolute value less than $2^{s_i + h}$, then $2^{-s_i}$ times the row is a row of integers of at most $h$ bits, and similarly for the columns of $B_m$ with exponents $t_j$. The scaled product of integer matrices is then exact, and $(A_m B_m)_{ij} = 2^{s_i + t_j} (\hat A \hat B)_{ij}$.

For matrices whose entries in each row and column have similar magnitudes, which is the common case, there is a single block, and the computed midpoint of each entry has a single rounding error of at most half an ulp at $p$ bits: the result is as accurate as the classical algorithm, which also rounds each dot product once. With $k$ blocks, each entry accumulates at most $k$ roundings. Entries of very different magnitude in the same row (for example, a matrix with a huge first column) create several blocks, and a matrix with huge dynamic range degenerates to classical multiplication in blocks of 30. The bound $H$ of about $1.25p$ bits allows the integer entries to be somewhat longer than the precision, so that ordinary dynamic range within a row does not split blocks; the constant 192 accommodates low-precision inputs.

The integer matrix multiplication fmpz_mat_mul chooses its algorithm from the dimension and the entry size (fmpz_mat.h). At precisions of a few hundred bits, its multimodular method, which replaces the $p$-bit products by matrix products modulo word-size primes, gives most of the speedup over the classical algorithm.

Radii: products in double precision#

The radius term $A_r(|B_m| + B_r) + |A_m| B_r$ needs only an upper bound of moderate accuracy, so it is computed with matrices of nonnegative double numbers. The magnitudes $|A_m| + A_r$, $|B_m|$, $A_r$, $B_r$ are converted to mag_t upper bounds (30-bit mantissas), and _arb_mat_addmul_rad_mag_fast again splits the inner dimension into blocks in which the exponents of each row and column span at most 900 bits. In a block, each row and each column is scaled by a power of two so that all entries are representable as normal double numbers without rounding, and the product is computed with ordinary double arithmetic in cache-friendly blocks.

The scaled double products are multiplied by this factor and added, with rounding upwards, to the radii of $C$; blocks of at most 5 indices use mag_t arithmetic directly. The radius computation costs $O(MNP)$ floating-point operations in the worst case, but these are much cheaper than $p$-bit multiplications, and for exact inputs (all radii zero) it is skipped. If only one factor is exact, only one radius product is computed.

Algorithm selection#

arb_mat_mul uses the classical algorithm (threaded if several threads are available and the matrices are not tiny) when one of the dimensions is at most 60 (for $p \le 128$), 50 (for $p \le 512$) or 40, and arb_mat_mul_block otherwise. The block algorithm itself falls back to classical multiplication if some midpoint has more than 8000 bits, if both matrices have less than 10% nonzero entries and $H > 1024$, or if some exponent is too large for the fast integer exponent arithmetic.

For complex matrices, acb_mat_mul_reorder writes $A = A_1 + i A_2$ and $B = B_1 + i B_2$ and computes the real and imaginary parts $A_1 B_1 - A_2 B_2$ and $A_1 B_2 + A_2 B_1$ with four real block multiplications (two if one of the matrices is real, one if both are). acb_mat_mul uses it when the inner dimension $N$ and the number of columns $P$ of $B$ satisfy $\min(N, P) \ge \max(20, 5 + b/64)$, where $b < 8000$ is the minimum of $p$ and the largest number of bits of a midpoint; otherwise it multiplies classically (threaded if possible).

Use in FLINT#

• arb_mat_mul, arb_mat_mul_classical, arb_mat_mul_threaded, arb_mat_mul_block, arb_mat_sqr, _arb_mat_addmul_rad_mag_fast (arb_mat.h); acb_mat_mul, acb_mat_mul_reorder (acb_mat.h).
• gr_mat_mul over the contexts gr_ctx_init_real_arb and gr_ctx_init_complex_acb calls arb_mat_mul and acb_mat_mul (gr_mat.h).
• The same idea of exact integer products of blocks of similar magnitude is used for polynomials in Multiplying polynomials with ball coefficients.
• Preconditioned linear solving (Verified linear algebra) and the eigenvalue enclosures (Eigenvalue enclosures) spend most of their time in these products.