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., a couple containing a center matrix and a radius matrix.

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
BallArithmetic._compute_exclusion_circle_level_set_priori — Method
_compute_exclusion_circle_level_set_priori(T,
λ,
ϵ;
rel_pearl_size,
max_initial_newton)

This method bounds the resolvent on a circle centered at λ that intersects in at least one point z0 the ϵ level set. This intersection is found by a Newton step, and fixes the radius of the circle,

The value of rel_pearl_size gives us the relative radius of the pearls with respect to the radius of the circle

Some rule of thumbs for the number of SVD computations: if relpearlsize is 1/32, we are going to compute and certify 160 svds, if relpearlsize is 1/64 we are going to compute and certify 320 svds. In other words, the time of the computation scales linearly with the quality of the pearl necklace

BallArithmetic.collatz_upper_bound_L2_opnorm — Method
collatz_upper_bound_L2_opnorm(A::BallMatrix; iterates=10)

Give a rigorous upper bound on the ℓ² norm of the matrix A by using the Collatz theorem.

We use Perron theory here: if for two matrices with B positive |A| < B we have ρ(A)<=ρ(B) by Wielandt's theorem Wielandt's theorem

The keyword argument iterates is used to establish how many times we are iterating the vector of ones before we use Collatz's estimate.

BallArithmetic.compute_enclosure — Method
compute_enclosure(A::BallMatrix, r1, r2, ϵ; max_initial_newton = 30,
        max_steps = Int64(ceil(256 * π)), rel_steps = 16)

Given a BallMatrix `A`, this method follows the level lines of level `ϵ`

around the eigenvalues with modulus bound between r1 and r2.

The keyword arguments - maxinitialnewton: maximum number of newton steps to reach the level lines - maxsteps: maximum number of steps following the contour - relsteps: relative integration step for the Euler method

The method outputs an array of truples: - the first element is the eigenvalue we are enclosing (in the case of the excluding circles, it is 0.0 or the maximum modulus of the eigenvalues) - the second element is an upper bound on the resolvent norm - the third element is a list of points on the enclosing line; the resolvent is rigorously bound on circles centered at each point and of radius 5/8 the distance to the previous point

BallArithmetic.evbox — Method
evbox(A::BallMatrix{T})

Compute rigorous enclosure of each eigenvalue following Ref. [3] TODO: Using Miyajima's algorithm is overkill, may be worth using

References

  • [3] Miyajima, JCAM 246, 9 (2012)
BallArithmetic.fft — Method
fft

Computes the FFT of a BallVector using the a priori error bound in Ref. [4]

References

  • [4] Higham, Siam (1996)
BallArithmetic.gevbox — Method
gevbox(A::BallMatrix{T}, B::BallMatrix{T})

Compute rigorous enclosure of each eigenvalue in generalized eigenvalue problem following Ref. [3]

References

  • [3] Miyajima, JCAM 246, 9 (2012)
BallArithmetic.svd_bound_L2_opnorm — Method
svd_bound_L2_opnorm(A::BallMatrix{T})

Returns a rigorous upper bound on the ℓ²-norm of the ball matrix A using the rigorous enclosure for the singular values implemented in svd/svd.jl

BallArithmetic.svd_bound_L2_opnorm_inverse — Method
svd_bound_L2_opnorm_inverse(A::BallMatrix)

Returns a rigorous upper bound on the ℓ²-norm of the inverse of the ball matrix A using the rigorous enclosure for the singular values implemented in svd/svd.jl

BallArithmetic.svd_bound_L2_resolvent — Method
svd_bound_L2_resolvent(A::BallMatrix, lam::Ball)

Returns a rigorous upper bound on the ℓ²-norm of the resolvent of A at λ, i.e., ||(A-λ)^{-1}||_{ℓ²}

BallArithmetic.upper_abs — Method
upper_abs(A)

Return a floating point matrix B whose entries are bigger or equal (componentwise) any of the entries of A

BallArithmetic.upper_bound_L2_opnorm — Method
upper_bound_L_inf_opnorm(A::BallMatrix{T})

Returns a rigorous upper bound on the ℓ²-norm of the ball matrix A using the best between the Collatz bound and the interpolation bound

BallArithmetic.upper_bound_norm — Function
upper_bound_norm(A::BallMatrix, p::Real = 2)

Compute a rigorous upper bound for the Frobenius p-norm of a BallMatrix

BallArithmetic.NumericalTest.rounding_test — Method
rounding_test(n, k)

Let u=fill(2^(-53), k-1) and let A be the matrix [I u; 0 2^(-53)]

This test checks the result of A*A' in different rounding modes, running BLAS on n threads