GitHub

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

source
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.

source
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

source
BallArithmetic.epsilon_inflation — Method
epsilon_inflation(A::BallMatrix{T}, b::BallVector{T};
                  r=0.1, ϵ=1e-20, iter_max=20) where {T<:AbstractFloat}

Gives an enclosure of the solution of the square linear system $Ax=b$ using the ϵ-inflation algorithm, see algorithm 10.7 of

Input

  • A – square matrix of size n × n
  • b – vector of length n or matrix of size n × m
  • r – relative inflation, default 10%
  • ϵ – absolute inflation, default 1e-20
  • iter_max – maximum number of iterations

Output

  • x – enclosure of the solution of the linear system
  • cert – Boolean flag, if cert==true, then x is certified to contain the true

solution of the linear system, if cert==false, then the algorithm could not prove that $x$ actually contains the true solution.

Algorithm

Given the real system $Ax=b$ and an approximate solution $̃x$, we initialize $x₀ = [̃x, ̃x]$. At each iteration the algorithm computes the inflation

$y = xₖ * [1 - r, 1 + r] .+ [-ϵ, ϵ]$

and the update

$xₖ₊₁ = Z + (I - CA)y$,

where $Z = C(b - Ax₀)$ and $C$ is an approximate inverse of $A$. If the condition $xₖ₊₁ ⊂ y$ is met, then $xₖ₊₁$ is a proved enclosure of $A⁻¹b$ and cert is set to true. If the condition is not met by the maximum number of iterations, the latest computed enclosure is returned, but $cert$ is set to false, meaning the algorithm could not prove that the enclosure contains the true solution. For interval systems, $̃x$ is obtained considering the midpoint of $A$ and $b$.

Notes

Examples

source
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)
source
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)
source
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

source
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

source
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

source
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

source