fredrikj.net / blog /

# High-precision linear algebra in Julia: BigFloat vs Arb

*July 31, 2018*

A few persons have asked me about the relative performance of Julia's native matrices with BigFloat elements and Arb matrices available through Nemo.jl. Unlike machine-precision matrices which build on BLAS technology, BigFloat matrices in Julia use generic code that has not been optimized specifically for BigFloat elements, so we can expect Arb's specialized C implementations to be faster. But by how much?

In this blog post, I provide some quick benchmarks done with Nemo 0.8.5 in Julia 0.6.4. We take the test matrix $A$ to be the DCT-II matrix of size $N$. This matrix can be constructed as follows as a Julia BigFloat matrix and as a Nemo Arb matrix:

setprecision(256) N=10; A = [j == 0 ? 1/sqrt(BigFloat(N)): cos(j*(k+0.5)*BigFloat(pi)/N) * sqrt(BigFloat(2)/N) for j in 0:N-1, k in 0:N-1];

R = ArbField(256) N=10; A = matrix(R, [j == 0 ? 1/sqrt(R(N)): cos(j*(k+0.5)*R(pi)/N) * sqrt(R(2)/N) for j in 0:N-1, k in 0:N-1]);

For these benchmarks, the precision is set to 256 bits, which is the default BigFloat precision in Julia. The relative performance will be slightly different with higher or lower precision. I only run the benchmarks with real matrices to keep things brief; complex matrices are about 3-4 times slower than real matrices with either BigFloat or Arb, so the relative performance figures should be similar.

## Matrix multiplication

We time computing `A * B` where `B = transpose(A)`.

$N$ | Julia BigFloat | Arb | Speedup |
---|---|---|---|

10 | 0.000310 s | 0.000112 s | 2.77 |

30 | 0.00940 s | 0.00243 s | 3.87 |

100 | 0.649 s | 0.0934 s | 6.95 |

300 | 23.645 s | 1.151 s | 20.54 |

1000 | 924.843 s | 24.606 s | 37.59 |

Not surprisingly, the clever matrix multiplication code in Arb (which uses FLINT behind the scenes) is leaps and bounds faster than Julia's generic code doing classical matrix multiplication over BigFloats.

## Inverse

We time computing `inv(A)` which is equivalent to solving the system $AX = I$ in both Arb and Julia.

The overhead inflicted by tracking error bounds in Arb is negligible for matrix multiplication.
However, linear solving switches from straightforward LU factorization to a significantly more
expensive algorithm
using preconditioning
to prevent the error bounds from blowing up exponentially with $N$.
Arb 2.14 introduced separate "approx" functions for
performing a few operations (including LU factorization and solving) *without error bounds*.
These functions may be useful as drop-in replacements for the usual numerical Julia functions
when heuristic error estimates are sufficient.
Wrappers for these functions are not yet implemented in Nemo, but it is easy
to write such functions manually in Julia. Here is a wrapper for in-place approximate solving,
and a corresponding approximate matrix inverse function:

function approx_solve!(z::arb_mat, x::arb_mat, y::arb_mat) r = ccall((:arb_mat_approx_solve, :libarb), Cint, (Ref{arb_mat}, Ref{arb_mat}, Ref{arb_mat}, Int), z, x, y, prec(base_ring(x))) r == 0 && error("Matrix cannot be inverted numerically") nothing end function approx_inv(x::arb_mat) y = one(parent(x)) approx_solve!(y, x, y) return y end

The table compares performance separately for the rigorous and approx versions:

$N$ | Julia BigFloat | Arb | Speedup | Arb (approx) | Speedup (approx) |
---|---|---|---|---|---|

10 | 0.000308 s | 0.000250 s | 1.23 | 0.000205 s | 1.50 |

30 | 0.00983 s | 0.00868 s | 1.13 | 0.00400 s | 2.46 |

100 | 0.786 s | 0.287 s | 2.74 | 0.111 s | 7.08 |

300 | 29.018 s | 4.372 s | 6.64 | 2.378 s | 12.20 |

1000 | 1641.659 s | 84.715 s | 19.38 | 51.168 s | 32.08 |

The rigorous algorithm in Arb is barely faster for small $N$ but eventually picks up steam. The approximate version performs even better.

## Matrix exponential

We time computing `exp(A)` with Arb (note that this gives the matrix exponential and not the elementwise exponential like `exp` applied to a normal Julia array).

Unfortunately, Julia's matrix exponential function `expm` does not work with `BigFloat` matrices,
and I'm not aware of any other Julia package that offers such a function except for mpmath via PyCall.
So, here is a totally lopsided comparison to `mpmath.expm` (also with 256-bit precision) to have
a point of reference.

$N$ | mpmath | Arb | Speedup |
---|---|---|---|

10 | 0.175 s | 0.00236 s | 74.1 |

30 | 3.762 s | 0.0456 s | 82.5 |

100 | 145.82 s | 1.334 s | 109.3 |

300 | 16.600 s | ||

1000 | 337.495 s |

## Conclusions

Arb matrices are likely to be slightly faster than Julia's native BigFloat matrices in general, just due to having much less garbage collection and memory allocation overhead.
In addition, Arb matrices are likely to be *much* faster for large dimensions due to using algorithms based on fast matrix multiplication.

Should you be using Nemo's Arb matrices instead of BigFloat matrices when 53-bit doubles don't cut it? If rigorous error bounds are useful, then the answer should be obvious! Otherwise, the short answer is that it depends. The benchmark results speak for themselves, but here are two more points to consider:

- In practice, Arb may perform worse if the error bounds in ball arithmetic are much more pessimistic than the actual numerical error so that higher precision has to be used. In such scenarios (assuming that rigorous error bounds are not necessary), a perfectly viable alternative to using BigFloats is to work with Arb matrices where one just sets the radii to zero and uses approx methods where available.
- Arb matrices currently provide relatively limited functionality (arithmetic, LU factorization, nonsingular solving, inverse, determinant, characteristic polynomial, exponential, and a few more functions) while the Julia standard library has a much wider range of functions (however, not all functions support BigFloats).