BallArithmetic

Documentation for BallArithmetic.

In this package we use the tecniques first introduced in Ref. [1], following the more recent work Ref. [2] to implement a rigorous matrix product in mid-radius arithmetic.

This allows to implement numerous algorithms developed by Rump, Miyajima, Ogita and collaborators to obtain a posteriori guaranteed bounds.

The main object are BallMatrices, i.e., midpoint matrices equipped with non-negative radii that provide rigorous entrywise enclosures.

Sylvester equations

sylvester_miyajima_enclosure provides a componentwise enclosure for solutions of the Sylvester equation following the fast verification method of Ref. [3]. When the data originate from an upper triangular Schur factor T, triangular_sylvester_miyajima_enclosure extracts the blocks T₁₁, T₁₂, and T₂₂, solves the associated Sylvester system T₂₂' Y₂ - Y₂ T₁₁' = T₁₂', and returns the Miyajima enclosure for the unknown block Y₂.

BallMatrix

BallMatrix is the midpoint-radius companion of the scalar Ball type. The midpoint matrix stores the approximation we would normally compute in floating-point arithmetic, whereas the radius matrix captures all sources of uncertainty (input radii, floating-point error, subnormal padding, …). Each method documented below updates both components so the result remains a rigorous enclosure.

Constructors and accessors

The constructors delegate to the underlying BallArray to perform shape and type validation. Working through them in order provides a tour of how the storage is organised:

Arithmetic

Binary operations follow a common pattern: operate on the midpoint data as if the values were exact, then grow the radius using outward rounding. The comments inside src/types/matrix.jl walk through the steps in more detail.

julia> using BallArithmetic
julia> A = ones((2, 2))2×2 Matrix{Float64}:
 1.0  1.0
 1.0  1.0
julia> bA = BallMatrix(A, A/128)2×2 BallMatrix{Float64, Float64, Ball{Float64, Float64}, Matrix{Float64}, Matrix{Float64}}:
 1.0 ± 0.0078125  …  1.0 ± 0.0078125
 1.0 ± 0.0078125     1.0 ± 0.0078125
julia> bA^22×2 BallMatrix{Float64, Float64, Ball{Float64, Float64}, Matrix{Float64}, Matrix{Float64}}:
 2.0 ± 0.0313721  …  2.0 ± 0.0313721
 2.0 ± 0.0313721     2.0 ± 0.0313721

Rounding-mode controlled products

Some matrix enclosures benefit from explicitly steering the floating-point rounding mode. The wrapper BallArithmetic.oishi_MMul implements the Oishi–Rump product, which evaluates the real and imaginary parts of F*G with downward and upward rounding and returns the result as a BallMatrix. The routine is particularly useful when replicating the eigenvalue and singular value enclosures described in Ref. [4].

Internally we also expose the auxiliary kernels from Ref. [5]. The helpers _ccrprod, _cr, _iprod, and _ciprod implement Algorithms 4–7 and propagate rectangular or ball bounds through matrix products. They are available for advanced workflows that need direct access to the underlying interval data.

using BallArithmetic
setprecision(BigFloat, 128) do
    F = Complex{BigFloat}[1 + im 2; 3 - im 4]
    G = Complex{BigFloat}[2 - im 1; -1 3 + im]
    B = BallArithmetic.oishi_MMul(F, G)
    (mid(B), rad(B))
end

(Complex{BigFloat}[1.0 + 1.0im 7.0 + 3.0im; 1.0 - 5.0im 15.0 + 3.0im], BigFloat[0.0 0.0; 0.0 0.0])