API - Core Types

BallArithmetic

Tools for rigorous linear algebra with floating-point balls. The package builds on Ball numbers—midpoint/radius pairs that track rounding error—to provide array, matrix, and vector types with propagated uncertainty. High-level routines such as norm bounds, eigenvalue certification, and singular value enclosures are exposed via this module after it has been loaded.

Most functionality is available through exported symbols documented in the respective source files. Loading the module also enables auxiliary utilities such as rigorous rounding macros and promotion rules that let Ball values interact seamlessly with base Julia numeric types.

Ball{T, CT}

Closed floating-point ball with midpoint type CT and radius type T. Each value represents the set { c + δ : |δ| ≤ r } where c::CT is the stored midpoint and r::T ≥ 0 is the radius. Both real and complex midpoints are supported as long as the radius is expressed in the underlying real field. The type behaves as a number and participates in arithmetic with rigorous outward rounding.

Ball{T, CT}

Closed floating-point ball with midpoint type CT and radius type T. Each value represents the set { c + δ : |δ| ≤ r } where c::CT is the stored midpoint and r::T ≥ 0 is the radius. Both real and complex midpoints are supported as long as the radius is expressed in the underlying real field. The type behaves as a number and participates in arithmetic with rigorous outward rounding.

Ball{T, CT}

Closed floating-point ball with midpoint type CT and radius type T. Each value represents the set { c + δ : |δ| ≤ r } where c::CT is the stored midpoint and r::T ≥ 0 is the radius. Both real and complex midpoints are supported as long as the radius is expressed in the underlying real field. The type behaves as a number and participates in arithmetic with rigorous outward rounding.

mid(x)

Return the midpoint of x. For plain numbers the midpoint is the value itself, while for balls the stored center is returned.

mid(v::AbstractVector{<:Ball})

Return a vector of midpoints for a collection of balls.

mid(M::AbstractMatrix{<:Ball})

Return a matrix of midpoints for a collection of balls.

mid(A::AbstractArray)

Fallback definition that treats ordinary arrays as their own midpoint representation. Specialisations for BallArray overload this method to return the stored midpoint data.

mid(A::BallArray)

Return the stored midpoint array.

mid(A::AbstractMatrix)

Return the midpoint matrix associated with A. For plain matrices the midpoint equals the matrix itself; for BallMatrix values this method is extended elsewhere to extract the stored midpoint data.

mid(v::AbstractVector)

Treat ordinary vectors as their own midpoints.

rad(x)

Return the radius associated with x. Numbers default to a zero radius, and balls return their stored uncertainty.

rad(v::AbstractVector{<:Ball})

Return a vector of radii for a collection of balls.

rad(M::AbstractMatrix{<:Ball})

Return a matrix of radii for a collection of balls.

rad(A::AbstractArray)

Return a zero array of matching size that serves as the default radius for non-ball arrays. The storage layout of A is preserved; see _zero_radius.

rad(A::BallArray)

Return the stored radius array.

rad(A::AbstractMatrix{T})

Return a matrix of radii matching the size of A. For non-ball matrices this defaults to a zero matrix, while for BallMatrix values the method is overloaded to provide the stored uncertainty information.

The storage layout of A is preserved, so a Diagonal, triangular or sparse midpoint matrix yields a radius matrix of the same type; see _zero_radius.

rad(A::AbstractMatrix{Complex{T}})

Return a matrix of real radii matching the size of the complex matrix A. Even for complex entries the radius is measured over the underlying real field, hence the resulting matrix has element type T.

rad(v::AbstractVector)

Default radius for non-ball vectors: a zero vector of the appropriate floating-point type. The storage layout of A is preserved, so a sparse vector yields a sparse radius; see _zero_radius.

mag(x::Ball)

An upper bound of the largest absolute value of the points of x.

mig(x::Ball)

A lower bound of the smallest absolute value of the points of x, zero when x contains zero.

midtype(::Ball)

Return the type used to store the midpoint component of a Ball. This is useful for allocating arrays that mirror the internal layout of a ball or a collection of balls.

radtype(x)

Return the floating-point type used to store radii for x. The helper accepts either a ball instance or the associated type, mirroring the behaviour of midtype.

inf(x::Ball)

Return the infimum (lower endpoint) of the set represented by x by evaluating mid(x) - rad(x) with downward rounding.

sup(x::Ball)

Return the supremum (upper endpoint) of the set represented by x by evaluating mid(x) + rad(x) with outward rounding.

ball_hull(a::Ball, b::Ball)

A ball containing both a and b.

For real midpoints the result encloses [min(inf(a), inf(b)), max(sup(a), sup(b))], the interval hull of the two sets; the enclosure is the one _ball_from_bounds states.

For complex midpoints the midpoint of the result is that of a or that of b, whichever gives the smaller radius, and the radius is max(rad(x), add_up(dist_up(mid(a), mid(b)), rad(y))) with x the argument whose midpoint was kept and y the other one. The disc of x is contained because the radius is at least rad(x), and the disc of y because every one of its points is within |mid(a) − mid(b)| + rad(y) of the midpoint that was kept. The result is not the smallest enclosing disc, whose centre lies between the two midpoints.

Note that dist_up bounds the distance through abs_up(::Complex), which squares the components of the difference, so the radius is Inf once a component passes about sqrt(floatmax(T)).

intersect_ball(a::Ball, b::Ball)

Return the intersection of the real balls a and b. When the balls do not overlap the function returns nothing to indicate that the intersection is empty.

in0(B₁::Ball, B₂::Ball)

Return true when B₁ is contained in the interior of B₂, which for balls is |mid(B₁) − mid(B₂)| + rad(B₁) < rad(B₂). The left-hand side is bounded above, the distance of the centres by dist_up and the sum by add_up, and compared against the stored radius, which is exact, so a true answer is a proof of the containment. (The modulus of a complex difference is not obtained from abs under an upward rounding mode: abs is hypot, which does not honour it, and the difference itself rounds a negative component toward zero.)

This is the predicate Rump writes in0 and defines in footnote 19 of

S. M. Rump, Verification methods: rigorous results using floating-point arithmetic, Acta Numerica 19 (2010) 287-449, doi 10.1017/S096249291000005X,

as "checks X ⊂ int(Y) componentwise". Every self-mapping test in the Krawczyk-Moore-Rump line needs the interior, not the closure: a fixed point on the boundary satisfies X ⊆ Y without giving the contraction the argument rests on, so in, which is not strict, is not a substitute.

in0(A₁::BallArray, A₂::BallArray)

Componentwise interior containment: true when every entry of A₁ lies in the interior of the corresponding entry of A₂, which is the predicate Rump writes in0 (see the scalar in0 for the reference). Returns false at the first entry that fails.

This is the test every self-mapping argument in the Krawczyk-Moore-Rump line needs; in, which is not strict, will not serve, since a fixed point on the boundary satisfies X ⊆ Y without giving the contraction.

in0(A₁::AbstractArray{<:Ball}, A₂::AbstractArray{<:Ball})

Componentwise interior containment for plain arrays whose entries are balls, which is what ball arithmetic on a broadcast expression returns. Same test as the BallArray method, entry by entry.

±(c, r)

Shorthand constructor for Ball(c, r). The operator mirrors the common mathematical notation c ± r for centered intervals.

_ball_from_bounds(lo, hi)

The ball enclosing the interval [lo, hi], for floating-point endpoints with lo ≤ hi. inv, sqrt and abs of a real ball, and ball_hull and intersect_ball, all build their result here.

Writing B for the result, what is guaranteed is that in exact arithmetic

B.c − B.r ≤ lo    and    B.c + B.r ≥ hi,

so that [lo, hi] ⊆ [B.c − B.r, B.c + B.r], together with B.r ≥ 0. Each of the two inequalities follows from one directed subtraction on its own: the radius is at least sub_up(B.c, lo) ≥ B.c − lo, which gives the first, and at least sub_up(hi, B.c) ≥ hi − B.c, which gives the second, so neither appeals to the direction in which the midpoint was rounded. A radius measured from one end alone, or from the half-width, reaches the other end only when the midpoint lies on the far side of (lo + hi)/2, and that is where issue #73 came from: the ball_hull of v0.3.0 rounded the midpoint to nearest and took the radius to be the half-width rounded up, so on ball_hull(Ball(-1.0, 0.5), Ball(-0.6, 0.75)) the upper end of the result fell one unit in the last place below the upper end of its second argument. The version that preceded this one rounded the midpoint up and measured the radius from lo, which encloses, but only because of the direction of that rounding.

The midpoint is lo + (hi − lo)/2 rounded up where hi − lo is representable, and lo/2 + hi/2 rounded up where it is not, so an interval wider than floatmax(T) still has a finite midpoint. Both forms give B.c ≥ (lo + hi)/2 exactly, which is what keeps the radius near the half-width rather than near the full width; the enclosure above does not depend on it.

An endpoint that is infinite returns Ball(zero(T), T(Inf)), the whole line being the only enclosure of an unbounded interval that a ball with a finite midpoint admits. Endpoints out of order, and an endpoint that is NaN, raise an ArgumentError: neither describes a set to enclose.

abs_up(x::BigFloat)

Compute the absolute value of a BigFloat with guaranteed upper bound. For positive x, this is just x. For negative x, this rounds -x upward.

abs_up(x) / abs_down(x)

Upper and lower bounds of |x|. For a real x the modulus is exact. For a complex z it is √(re² + im²) with every operation rounded in the stated direction, and not hypot, whose result is not correctly rounded.

BallArray{T, N, NT, BT, CA, RA}

Multi-dimensional array whose entries are Ball values. The type stores midpoint data c::CA and radius data r::RA separately while presenting an AbstractArray{BT, N} interface that behaves like an array of balls. The parameters mirror the element type and storage layout and are inferred automatically from the provided midpoint and radius containers.

BallMatrix{T, NT, BT, CM, RM}

Alias for the two-dimensional BallArray used to represent matrices with rigorous error control. The type parameters mirror those of BallArray and describe the base floating-point type T, the element type NT, the Ball container BT, and the concrete matrix types used for the midpoints (CM) and radii (RM).

Users typically construct instances through the BallMatrix constructors below rather than specifying these parameters explicitly.

What set of matrices a ball matrix is

A ball matrix is the set of the matrices whose stored entries lie in the stored balls: an entry with midpoint m and radius r ranges over |x − m| ≤ r, and an entry that is not stored has no tolerance of its own. The storage types CM and RM therefore say what the set is, and the routines dispatch on them.

• Midpoint and radius stored as plain matrices: every entry is stored, and the set is the matrices X with |X − M| ≤ R entrywise. With positive radii it is an open set of matrices, whatever M is: it contains matrices that are not Hermitian, not symmetric, not triangular. Enclosures of its eigenvalues are discs of the complex plane.
• Midpoint and radius stored as Diagonal, or as triangular: the entries outside the structure are zero and are not stored, so every matrix of the set is diagonal, or triangular.
• Midpoint stored as Hermitian (or real Symmetric) and radius stored as Symmetric: one triangle is stored, with one tolerance for each pair {i, j}, and the other triangle is determined by it. Every matrix of the set is Hermitian: it is the Hermitian matrices X with |X − M| ≤ R, a HermitianBallMatrix. Enclosures of its eigenvalues are real intervals, and the statements of the routines are for these Hermitian matrices.
• Midpoint stored as Hermitian and radius stored as a plain matrix: every entry has its own tolerance, so the perturbation has no structure and the set is the full entrywise ball around a Hermitian midpoint. The enclosures are discs; what the storage of the midpoint changes is the approximation that is certified, which is that of a Hermitian matrix.

A routine that writes a structured ball matrix out as a dense one computes on a larger set for the Hermitian kind (the two triangles then vary independently) and on the same set for the diagonal and triangular kinds; the statement it returns is the one of the kind it was given.

Sums, differences and real multiples of Hermitian ball matrices are Hermitian ball matrices; products are not, and are returned with plain storage.

BallVector{T, NT, BT, CM, RM}

Alias for the one-dimensional BallArray, representing vectors of balls.

HermitianBallMatrix{T}

The BallMatrix whose midpoint is stored as Hermitian (or real Symmetric) and whose radius is stored as Symmetric: the set of the Hermitian matrices X with |X − M| ≤ R entrywise. Building one states the hypothesis that the matrix is Hermitian; the routines that receive it return real intervals for its eigenvalues, and their statements are for the Hermitian matrices within the tolerances. See BallMatrix for the definition of the set.

HermitianBallMatrix(M, R; uplo = :U)
HermitianBallMatrix(M; uplo = :U)
HermitianBallMatrix(A::BallMatrix; uplo = :U)

BallMatrix(Hermitian(M, uplo), Symmetric(R, uplo)), with radius zero in the second form and M = mid(A), R = rad(A) in the third. Only the triangle uplo of M and of R is read, as Hermitian and Symmetric do; the other triangle is not compared with it. On the diagonal the midpoint is the real part of that of M, with the same radius.

DenseBallMatrix{T, NT}

The BallMatrix whose midpoints and radii are stored in a Matrix. Unlike BallMatrix{T, NT}, which leaves the storage types free, it is a concrete type; the fields of the result structs are declared with it, and a ball matrix with other storage is converted to it when the struct is built.

DenseBallVector{T, NT}

The BallVector whose midpoints and radii are stored in a Vector: a concrete type, for the fields of the result structs. See DenseBallMatrix.

_zero_radius(A, ::Type{T})

Return an all-zero, T-valued array laid out like A. Structured and sparse midpoints keep their storage type, so BallMatrix(Diagonal(...)) gets a Diagonal radius rather than a dense one — for a sparse midpoint the dense fallback costs orders of magnitude more memory than the midpoints themselves.

This is sound because every entry a structured type stores implicitly is an exact zero (or, for the unit triangular types, an exact one), so the corresponding radius is exactly zero and needs no storage.

fill!(similar(A, T), zero(T)) is used rather than zero(A), which would keep a complex element type, or zero(similar(A, T)), which reads the undefined references similar leaves behind for BigFloat.

For Adjoint, Transpose, SubArray and ranges similar already yields a plain dense array, matching the previous behaviour.

_check_axes(c, r)

Reject midpoint and radius containers that do not line up. Mismatched shapes used to be accepted silently: the resulting object reported size(c) while indexing or norm bounds later failed with a confusing BoundsError or DimensionMismatch far from the construction site.

Comparing axes rather than size keeps offset-indexed arrays working. The check is O(1), so it costs nothing on the hot paths.

isvalid_enclosure(A::BallArray) -> Bool

Report whether every radius of A is a genuine enclosure radius, that is non-negative and not NaN. A negative or NaN radius describes no set at all, so an array failing this test carries no rigorous meaning.

Validating the radii is the caller's responsibility. The package deliberately does not check them on construction: the radii it produces internally are non-negative by construction (accumulated under RoundUp from non-negative quantities), and scanning every entry would cost up to 76% of a BallMatrix addition at n = 100. When radii come from outside — read from a file, supplied by a user, converted from another package — call this (or check_enclosure) yourself before relying on the enclosure.

The axes of c and r are checked on construction, since that test is O(1).

check_enclosure(A::BallArray) -> A

Throw an ArgumentError unless isvalid_enclosure holds, otherwise return A unchanged so the call can be chained.

+(x::Ball)

Return x unchanged. Unary plus exists for completeness so that generic numeric code can treat balls like other scalar types.

Base.:+(x::Ball, y::Ball)

Combine two balls using addition and enlarge the radius so that the result remains a rigorous enclosure. The midpoint is the rounded sum and the radius accounts for both operands plus floating-point roundoff.

/(x::Ball, y::Ball)

Divide x by y by multiplying with the inverse of y. The operation inherits the same domain restrictions as inv.

inv(x::Ball)

Return the multiplicative inverse of a real ball. The method throws an ArgumentError when the interval straddles zero, since no rigorous inverse can be produced in that case.

inv(x::Ball{T, Complex{T}})

Return the multiplicative inverse of a complex ball. The method throws an ArgumentError when the ball is not proved to exclude zero. The image of the disc B(m, r) under 1/z is the disc B(conj(m)/D, r/D) with D = |m|² − r². D is bracketed by D_lo ≤ D ≤ D_hi with directed operations and D_lo > 0 is required; the centre is computed with D_lo, and the radius is r/D_lo plus the distance |m|(1/D_lo − 1/D_hi) between that centre and the exact one, plus the rounding of its two divisions.

sqrt(x::Ball)

Principal square root of a non-negative real ball. The method verifies that the enclosure stays within the domain of the square root and then propagates rounding errors to produce a rigorous result.

abs(x::Ball)

Return a real ball that encloses |z| for every z in x, that is the interval [max(0, |c| − r), |c| + r]. For a real ball that does not contain zero this is the ball of centre |c| and the same radius, exactly. For a complex ball the modulus of the centre is bounded below and above, since hypot is not correctly rounded.

conj(x::Ball)

Complex conjugate of a ball. The midpoint is conjugated while the radius remains unchanged.

in(x::Number, B::Ball)

Return true when the scalar x is proved to lie in the ball B: an upper bound of |c − x| is compared with the radius. A number that is not a float is first enclosed by Ball(x).

in(B₁::Ball{T}, B₂::Ball{T})

Check whether the enclosure B₁ is fully contained in B₂. The test expands the endpoints using outward rounding to ensure a rigorous decision.

in(B₁::Ball{T, Complex{T}}, B₂::Ball{T, Complex{T}})

Containment test for complex balls. The check reduces the problem to the real case by comparing the distance between midpoints with the radii of the two enclosures.

isless(a::Ball, b::Ball)

Compare two balls by their midpoints. This provides a total ordering for sorting and comparison operations. For rigorous "certainly less than" semantics, use sup(a) < inf(b).

isless(a::Ball, b::Number)

Compare a ball with a number by comparing the ball's midpoint.

<(a::Ball, b::Ball)

Compare two balls. Returns true if the ball midpoints satisfy a < b.

<=(a::Ball, b::Ball)

Compare two balls. Returns true if the ball midpoints satisfy a <= b.

>(a::Ball, b::Ball)

Compare two balls by midpoints.

>=(a::Ball, b::Ball)

Compare two balls by midpoints.

Base.convert(::Type{Ball{T, CT}}, x::Ball)

Convert a ball to the same enclosure expressed with alternative midpoint and radius types. This is typically used when promoting collections of balls to a common numeric representation.

Base.convert(::Type{Ball{T, CT}}, c::Number)

Embed a plain number into a ball whose midpoint is c converted to CT. The radius is zero when that conversion is exact and a bound on its rounding otherwise.

size(A::BallArray)

Forward the size of the underlying midpoint storage.

length(A::BallArray)

Total number of elements stored in the array, matching the midpoint container.

getindex(A::BallArray, inds...)

Indexing a BallArray returns either a single Ball or another BallArray depending on the provided indices. Midpoints and radii are looked up independently so that the enclosure remains rigorous.

setindex!(A::BallArray, x, inds...)

Assign a value x to the given indices by storing its midpoint and radius separately.

copy(A::BallArray)

Create a fresh BallArray with copies of both midpoint and radius storage.

real(A::BallArray)

Extract the real part of a BallArray. For purely real storage the result is returned unchanged, while complex arrays drop the imaginary part of the midpoint.

imag(A::BallArray)

Return the imaginary part of a BallArray. Real arrays produce a zero enclosure, while complex arrays keep the stored radii and extract the imaginary midpoints.

zeros(::Type{Ball}, dims)

Allocate a zero BallArray of the requested dimensions, using the element type's midpoint and radius types to choose the storage format.

ones(::Type{Ball}, dims)

Return a BallArray whose midpoints are filled with ones and whose radii are identically zero.

fill(x::Ball, dims...)

Create a BallArray where every element equals the ball x.

adjoint(M::BallMatrix)

Return the conjugate transpose of M, preserving rigorous enclosures by transposing both the midpoint and the radius matrices.

transpose(M::BallMatrix)

The transpose of M: midpoint and radius transposed. Without this method the generic one returns a lazy wrapper around M, which is not a BallMatrix.

Base.:+(A::BallMatrix, B::BallMatrix)

Combine two BallMatrix values elementwise using addition while tracking floating-point and enclosure errors. The midpoint matrices are added directly and the radii are enlarged using outward rounding to maintain a rigorous enclosure.

Base.:+(A::BallMatrix, B::AbstractMatrix)

Add a plain matrix to a BallMatrix by combining the midpoints directly and inflating the radii to account for the existing uncertainty in A and the floating-point error introduced by the addition.

Base.:+(B::AbstractMatrix, A::BallMatrix)

Commutative counterpart of +(::BallMatrix, ::AbstractMatrix), allowing the plain matrix to appear on the left-hand side.

+(A::BallMatrix, J::UniformScaling)

Add a UniformScaling operator (such as I) to a square BallMatrix. The diagonal of the midpoint is shifted by J and the corresponding radii are enlarged to maintain a valid enclosure.

+(J::UniformScaling, A::BallMatrix)

Symmetric counterpart to +(::BallMatrix, ::UniformScaling), allowing the uniform scaling to appear on the left-hand side.

+(A::BallMatrix, J::UniformScaling{Ball})

Add a ball-valued UniformScaling to a BallMatrix, incorporating both the midpoint and the radius of the scaling into the resulting enclosure.

Base.:+(A::BallVector, B::BallVector)

Combine two ball vectors elementwise using addition, enlarging the radius to include roundoff and the uncertainties of both operands.

-(x::Ball)

Negate the midpoint of x while keeping the radius unchanged. The result encloses the additive inverse of the represented set.

Base.:-(x::Ball, y::Ball)

Combine two balls using subtraction and enlarge the radius so that the result remains a rigorous enclosure. The midpoint is the rounded difference and the radius accounts for both operands plus floating-point roundoff.

Base.:-(A::BallMatrix, B::BallMatrix)

Combine two BallMatrix values elementwise using subtraction while tracking floating-point and enclosure errors. The midpoint matrices are subtracted directly and the radii are enlarged using outward rounding to maintain a rigorous enclosure.

Base.:-(A::BallMatrix, B::AbstractMatrix)

Subtract a plain matrix from a BallMatrix. The midpoint subtraction is performed elementwise while the radius is enlarged to remain enclosure safe.

Base.:-(B::AbstractMatrix, A::BallMatrix)

Subtract a BallMatrix from a plain matrix, reusing the implementation of -(::BallMatrix, ::AbstractMatrix) by swapping the argument order.

-(A::BallMatrix, J::UniformScaling)

Subtract a UniformScaling operator from a BallMatrix while updating the stored radii so that the result remains enclosure-safe.

-(A::BallMatrix, J::UniformScaling{Ball})

Subtract a ball-valued UniformScaling from a BallMatrix, accounting for both the midpoint and radius of the scaling in the resulting enclosure.

-(J::UniformScaling, A::BallMatrix)

Subtract a BallMatrix from a uniform scaling operator, useful for forming expressions such as λ * I - A with rigorous error bounds.

-(J::UniformScaling{Ball}, A::BallMatrix)

Subtract a BallMatrix from a ball-valued uniform scaling operator, expanding the resulting radii by both the uncertainty in J and the rounding performed during the subtraction.

Base.:-(A::BallVector, B::BallVector)

Combine two ball vectors elementwise using subtraction, enlarging the radius to include roundoff and the uncertainties of both operands.

*(x::Ball, y::Ball)

Multiply two balls and return the enclosure of the product. The midpoint is the product of the midpoints, whereas the radius collects propagated uncertainty from both operands and the intrinsic rounding error of the operation.

*(λ::Number, A::BallMatrix)

Scale the BallMatrix A by the scalar λ. Both the midpoint and the radius are scaled, and an outward-rounded padding proportional to floating-point error is added so the result remains a rigorous enclosure.

*(A::BallMatrix, λ::Number)

Right-multiplication of BallMatrix by a scalar. Delegates to left-multiplication.

*(λ::Ball, A::BallMatrix)

Scale A by a scalar Ball. The midpoint is scaled by the midpoint of λ and the radius is enlarged to account for both the uncertainty in λ and floating-point rounding.

*(A::BallMatrix, λ::Ball)

Right-multiplication of BallMatrix by a scalar Ball. Delegates to left-multiplication.

*(A::BallMatrix, B::BallMatrix)

Multiply two compatible BallMatrix values using the rigorously rounded matrix multiplication kernels provided in the MMul submodule.

*(A::BallMatrix, B::AbstractMatrix)

Multiply a BallMatrix by an arbitrary dense or structured matrix, propagating the midpoint product and rigorously bounding the resulting radii.

*(A::AbstractMatrix, B::BallMatrix)

Multiply an arbitrary dense or structured matrix with a BallMatrix. The enclosure of the result is obtained using the same rigorous matrix multiplication kernel as in the purely ball-valued case.

*(A::BallMatrix{T, Complex{T}}, B::BallMatrix{T, T})

Multiply a complex BallMatrix by a real one by dispatching to the real and imaginary parts separately.

*(A::BallMatrix{T, T}, B::BallMatrix{T, Complex{T}})

Multiply a real BallMatrix by a complex one, computing the product of midpoint matrices componentwise and rigorously propagating the radii.

*(A::BallMatrix{T, Complex{T}}, B::BallMatrix{T, Complex{T}})

Multiply two complex BallMatrix values. The product is formed by combining the real and imaginary parts via the previously defined real matrix multiplications.

*(A::AbstractMatrix{T}, B::BallMatrix{T, Complex{T}})

Multiply a real floating-point matrix, dense or sparse, by a complex BallMatrix, through the real kernel on the real and imaginary parts of B, so a sparse A stays sparse.

*(A::AbstractMatrix{Complex{T}}, B::BallMatrix{T, T})

Multiply a complex floating-point matrix, dense or sparse, by a real BallMatrix, through the real kernel on the real and imaginary parts of A.

*(A::AbstractMatrix{Complex{T}}, B::BallMatrix{T, Complex{T}})

Multiply a complex floating-point matrix, dense or sparse, by a complex BallMatrix, through the real-times-complex products of the real and imaginary parts of A; the generic plain-times-ball method reaches MMul4, which has real methods only.

*(λ::Number, v::BallVector)

Scale a ball vector by a scalar. The midpoint is scaled directly while the radius accounts for propagated uncertainty and roundoff.

*(λ::Ball, v::BallVector)

Scale a ball vector by a ball-valued scalar, combining the uncertainty in both arguments.

*(A::BallMatrix, v::BallVector)

Multiply a ball matrix with a ball vector. The vector is reshaped into a column matrix so that the existing BallMatrix multiplication handles the enclosure bookkeeping.

*(A::AbstractMatrix, v::BallVector)

Promote a plain matrix to a BallMatrix before multiplying it with a ball vector.

*(A::BallMatrix, v::AbstractVector)

Multiply a ball matrix with a plain vector by promoting the vector to a column BallMatrix and reusing the matrix-matrix multiplication kernel.

Typed on AbstractVector rather than Vector so that views, ranges and sparse vectors take the rigorous kernel too; previously they fell through to generic element-wise Ball multiplication, which is far slower and returns a Vector{Ball} instead of a BallVector.

iscertified(result) -> Bool

Whether result proves what it computes: true when the verification behind it succeeded, so that the statement in the docstring of typeof(result) holds for every matrix and every vector of the input balls, and false when some step of that verification did not apply, in which case the enclosures carried by result prove nothing and usually have infinite radii.

This is the supported way to ask the question. The field it reads is called certified in version 0.4; in 0.3.0 it was called success on the verified decompositions, on the generalised eigenvalue results and on most of the linear system solvers, verified on KrawczykResult, NKEigenpairResult, Oishi2023Result, RumpOishi2024Result and VerifiedLinearSystemResult, is_regular on RegularityResult and precomputation_success on SylvesterResolventResult. Code that calls iscertified reads the flag in either version.

A type that carries no such flag raises an ArgumentError naming it, in place of returning false: a false would read as "not proved", which is a claim about a computation that the type never made. The flags that mean something else, and so are not read here, are converged of the iterative methods and of the refinements, which is about an iteration reaching a tolerance; spectrum_covered of VerifyEigAllResult, which asserts in addition that the union of the discs holds the whole spectrum; empty of IterativeResult and of ShavingResult; and success of PreconditionerResult, whose preconditioner is an approximate inverse.

Example

B = randn(100, 100); A = B' * B + I
iscertified(verified_cholesky(Float64, A))     # true: A is positive definite
iscertified(verified_cholesky(Float64, [1.0 2.0; 2.0 1.0]))     # false: it is not

upper_bound_norm(A::BallMatrix, p::Real = 2)

An upper bound of the entrywise p-norm of every matrix in the ball A, for p = 1, 2 or Inf: the norm of |mid(A)| + rad(A), every operation rounded upwards. p = 2 is the Frobenius norm. It is not an operator norm; for those see upper_bound_L2_opnorm, upper_bound_L1_opnorm and upper_bound_L_inf_opnorm.

upper_bound_norm(v::BallVector, p::Real = 2)

An upper bound of the p-norm of every vector in the ball v, for p = 1, 2 or Inf: the norm of |mid(v)| + rad(v), every operation rounded upwards.

upper_bound_L1_opnorm(A::BallMatrix{T})

Rigorous upper bound on the ℓ¹ operator norm (largest column sum) of every matrix in the ball A: the largest column sum of upper_abs(A), rounded upward.

upper_bound_L2_opnorm(A::BallMatrix; svd = false)

Rigorous upper bound on the spectral norm of every matrix in the ball A: the least of the Collatz bound, the interpolation bound √(‖·‖₁‖·‖∞), the Frobenius bound and, when svd = true, the upper end of the largest singular value enclosed by svdbox (an O(n³) factorisation, see svd_bound_L2_opnorm). The first three are computed on upper_abs(A).

upper_bound_L_inf_opnorm(A::BallMatrix{T})

Rigorous upper bound on the ℓ^∞ operator norm (largest row sum) of every matrix in the ball A: the largest row sum of upper_abs(A), rounded upward.

collatz_upper_bound_L2_opnorm(A::BallMatrix; iterates = 5)

A rigorous upper bound on the spectral norm of every matrix in the ball A, by Collatz's inclusion applied to |A|ᵀ|A|.

Reference

S. M. Rump, Verified bounds for singular values, in particular for the spectral norm of a matrix and its inverse, BIT Numer. Math. 51 (2011) 367-384, doi 10.1007/s10543-010-0294-0, Section 3: for real or complex A,

‖A‖₂ ≤ ‖ |A| ‖₂                                                      (3.2)

and, for nonnegative B, since ‖B‖₂² = ρ(BᵀB) and BᵀB is nonnegative, "Collatz' Theorem in Perron-Frobenius Theory yields"

‖B‖₂² ≤ max_i (Bᵀ(Bx))_i / x_i      for every positive vector x.     (3.3)

"Few power set iterations starting with the vector of 1's" give the vector. Here B is upper_abs(A), an entrywise upper bound of |M| for every M in the ball, which is enough because ρ is monotone on nonnegative matrices.

iterates is the number of power iterations used to choose x; any value, including zero, gives a valid bound (with zero it is √‖BᵀB‖_∞). The result is Inf when the quotient is not finite.

svd_bound_L2_opnorm(A::BallMatrix{T})

Returns a rigorous upper bound on the ℓ²-norm (spectral norm) of the ball matrix A using the rigorous SVD enclosure from rigorous_svd.

This method computes the largest certified singular value, providing the tightest possible upper bound on ‖A‖₂. The bound is essentially exact (typically <0.01% overestimation), but requires O(n³) computation for the SVD.

Comparison with other L2 norm methods

*Collatz performs well on structured matrices (tridiagonal: ~0%, Hilbert: ~0%, diagonally dominant: ~26%) but poorly on random matrices (~200-500%).

Example

A = BallMatrix(randn(50, 50))
bound = svd_bound_L2_opnorm(A)  # Tight bound, slower
fast_bound = upper_bound_L2_opnorm(A)  # Looser bound, faster

See also: upper_bound_L2_opnorm, collatz_upper_bound_L2_opnorm, rigorous_svd

svd_bound_L2_opnorm_inverse(A::BallMatrix)

An upper bound of the spectral norm of M⁻¹ that holds for every matrix M of the ball A: the reciprocal of the lower end of the smallest singular value enclosed by svdbox, rounded up. The enclosure covers the singular values of every matrix of the ball, and σ_min(M)⁻¹ = ‖M⁻¹‖₂, so the bound is uniform over the ball; a matrix of the ball being nonsingular is part of what is proved.

Inf is returned when that lower end is not positive, which happens when the ball may contain a singular matrix and is not a statement that it does. For an empty matrix the value is zero.

triangular_inverse_two_norm_bound(U)

An upper bound of the spectral norm ‖U⁻¹‖₂, by ‖M‖₂ ≤ √(‖M‖₁‖M‖_∞) applied to the two bounds above, and Inf when either of them is. U may be a matrix with exact entries or an upper triangular BallMatrix whose balls below the diagonal are the point zero, and in the second case the bound holds for every matrix of the ball.

rump_oishi_2024_triangular_bound(T::BallMatrix, k::Int) -> upper bound of ‖T⁻¹‖₂

Theorem 1.3 of Rump and Oishi for an upper triangular T, split as T = [A B; 0 D] with A of order k, and D = D_d + D_f into its diagonal and the rest. The lower left block is zero, so the theorem reads: if ‖D_d⁻¹ D_f‖₂ < 1, then T is invertible and

‖T⁻¹‖₂ ≤ max{ ‖A⁻¹‖₂, ‖D_d⁻¹‖₂ / (1 − ‖D_d⁻¹ D_f‖₂) } · ψ(‖A⁻¹B‖₂),

with ψ(μ) the norm of [I −N; 0 I] for ‖N‖₂ = μ (Lemma 1.2, psi_schur_factor). The bound holds for every matrix of the ball T. Inf is returned where the hypothesis is not established or a diagonal entry of D may vanish. For k equal to the order the result is svd_bound_L2_opnorm_inverse of T.

How each quantity is bounded: ‖A⁻¹‖₂ by svd_bound_L2_opnorm_inverse; A⁻¹B is enclosed by backward_substitution and its norm, like that of D_d⁻¹ D_f (formed by ball divisions), bounded by collatz_upper_bound_L2_opnorm; ‖D_d⁻¹‖₂ is the largest 1/mig(d_ii); the final operations are rounded upwards.

The note states the theorem for real matrices. Complex triangular matrices are accepted here: the block factorisation, the Neumann series and Lemma 1.2 use only the norm.

Reference

S. M. Rump and S. Oishi, A note on Oishi's lower bound for the smallest singular value of linearized Galerkin equations, Japan J. Indust. Appl. Math. 41(2) (2024), 1097-1104, doi 10.1007/s13160-024-00645-7 ([8]), Lemma 1.2 and Theorem 1.3.

upper_abs(A)

Entrywise upper bound of |C| + R for a ball matrix with midpoint C and radius R, so that |X| ≤ upper_abs(A) entrywise for every X in the ball; for a plain matrix, of |A|.

Every norm bound below goes through this matrix: the spectral, Frobenius, ℓ¹ and ℓ^∞ norms are monotone in the moduli of the entries (Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., Lemma 6.6, for the spectral norm), so a bound on the norm of |C| + R bounds the norm of every matrix in the ball, and it is never larger than the norm of C plus the norm of R.

For Float32 and Float64 the moduli are formed with abs_up, as √(re² + im²) rounded upward operation by operation, because the complex abs is hypot, which is not correctly rounded and so is not a bound under RoundUp.

Oishi2023Result{T}

Result from the Oishi 2023 Schur complement bound computation.

Fields

• sigma_min_lower: Rigorous lower bound on the minimum singular value
• G_inv_upper: Upper bound on ‖G⁻¹‖₂ (= 1/σ_min)
• A_inv_norm: Upper bound on ‖A⁻¹‖₂
• A_inv_B_norm: Upper bound on ‖A⁻¹B‖₂
• C_A_inv_norm: Upper bound on ‖CA⁻¹‖₂
• schur_contraction: Upper bound on ‖Dd⁻¹(Df - CA⁻¹B)‖₂
• certified: Whether all conditions of Theorem 1 are satisfied
• block_size: The block size m used

The flag is called certified since version 0.4 and was called verified in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

oishi_2023_sigma_min_bound(G::BallMatrix{T}, m::Int) where {T}

Compute a rigorous lower bound on the minimum singular value of G using the Schur complement method from Oishi (2023).

Arguments

• G::BallMatrix: The matrix to analyze
• m::Int: Block size for the partition (A is m×m)

Method (Theorem 1 from Oishi 2023)

Partition G as:

G = [A  B]
    [C  D]

where A ∈ Mm, B ∈ M{m,n-m}, C ∈ M{n-m,m}, D ∈ M{n-m,n-m}.

Let Dd and Df be the diagonal and off-diagonal parts of D.

If the following conditions hold:

1. ‖A⁻¹B‖₂ < 1
2. ‖CA⁻¹‖₂ < 1
3. ‖Dd⁻¹(Df - CA⁻¹B)‖₂ < 1

Then G is invertible and:

‖G⁻¹‖₂ ≤ max{‖A⁻¹‖₂, ‖D_d⁻¹‖₂/(1 - ‖D_d⁻¹(D_f - CA⁻¹B)‖₂)}
         / ((1 - ‖A⁻¹B‖₂)(1 - ‖CA⁻¹‖₂))

Since σmin(G) = 1/‖G⁻¹‖₂, this gives a lower bound on σmin.

Returns

Oishi2023Result containing the bounds and verification status.

Reference

S. Oishi, Lower bounds for the smallest singular values of generalized asymptotic diagonal dominant matrices, Japan J. Indust. Appl. Math. 40 (2023), 1569-1585, doi 10.1007/s13160-023-00596-5 ([9]), Theorem 1.

oishi_2023_optimal_block_size(G::BallMatrix{T}; max_m::Int=100) where {T}

oishi_2023_sigma_min_bound for every block size m = 1, …, min(max_m, n − 1), and the m whose certified lower bound of σ_min(G) is the largest, as (best_m, best_result). Each certified result is a valid bound, so the largest is; when none is certified, best_result is the uncertified one for m = 1. The statement and its reference are those of oishi_2023_sigma_min_bound.

RumpOishi2024Result{T}

Result from the Rump-Oishi 2024 improved Schur complement bound computation.

This improves upon Oishi 2023 by:

1. Removing the conditions ‖A⁻¹B‖ < 1 and ‖CA⁻¹‖ < 1
2. Using the exact ψ(N) formula instead of 1/(1-‖N‖)
3. Optionally using a fast γ bound to avoid expensive matrix products

Fields

• sigma_min_lower: Rigorous lower bound on the minimum singular value
• G_inv_upper: Upper bound on ‖G⁻¹‖₂ (= 1/σ_min)
• A_inv_norm: Upper bound on ‖A⁻¹‖₂
• psi_A_inv_B: ψ(‖A⁻¹B‖₂) factor
• psi_C_A_inv: ψ(‖CA⁻¹‖₂) factor
• schur_contraction: Upper bound on ‖Dd⁻¹(Df - CA⁻¹B)‖₂
• used_fast_gamma: Whether the fast γ bound was used
• certified: Whether the method succeeded (only requires schur_contraction < 1)
• block_size: The block size m used

Reference

S. M. Rump and S. Oishi, Japan J. Indust. Appl. Math. 41(2) (2024), 1097-1104 ([8]), Lemma 1.2 and Theorem 1.3.

The flag is called certified since version 0.4 and was called verified in 0.3.0; the old name still resolves, with a deprecation warning, and iscertified reads it in either version.

rump_oishi_2024_sigma_min_bound(G::BallMatrix{T}, m::Int;
                                 try_fast_gamma::Bool=true) where {T}

Compute a rigorous lower bound on the minimum singular value of G using the improved Schur complement method from Rump-Oishi (2024).

Arguments

• G::BallMatrix: The matrix to analyze
• m::Int: Block size for the partition (A is m×m)
• try_fast_gamma::Bool=true: Whether to try the fast γ bound first

Improvements over Oishi 2023

1. Removes conditions 1 & 2: No longer requires ‖A⁻¹B‖ < 1 and ‖CA⁻¹‖ < 1. Only condition 3 (‖Dd⁻¹(Df - CA⁻¹B)‖ < 1) is needed.

2. Uses exact ψ(N) formula: Instead of 1/(1-‖N‖), uses ψ(μ) = √(1 + 2αμ√(1-α²) + α²μ²) which is tighter and works for any μ ≥ 0.

3. Fast γ bound (optional): Uses π(N) = √(‖N‖₁·‖N‖∞) to quickly check if the method will succeed, avoiding expensive matrix products.

Method (Theorem 1.3 from Rump-Oishi 2024)

Partition G as:

G = [A  B]
    [C  D]

If ‖Dd⁻¹(Df - CA⁻¹B)‖ < 1, then:

‖G⁻¹‖ ≤ max{‖A⁻¹‖, ‖Dd⁻¹‖/(1 - ‖Dd⁻¹(Df - CA⁻¹B)‖)} · ψ(‖A⁻¹B‖) · ψ(‖CA⁻¹‖)

Returns

RumpOishi2024Result containing the bounds and verification status.

Reference

S. M. Rump and S. Oishi, A note on Oishi's lower bound for the smallest singular value of linearized Galerkin equations, Japan J. Indust. Appl. Math. 41(2) (2024), 1097-1104, doi 10.1007/s13160-024-00645-7 ([8]), Lemma 1.2 and Theorem 1.3.

rump_oishi_2024_optimal_block_size(G::BallMatrix{T}; max_m::Int=100,
                                    try_fast_gamma::Bool=true) where {T}

rump_oishi_2024_sigma_min_bound for every block size m = 1, …, min(max_m, n − 1), and the m whose certified lower bound of σ_min(G) is the largest, as (best_m, best_result). Each certified result is a valid bound, so the largest is; when none is certified, best_result is the uncertified one for m = 1. The statement and its reference are those of rump_oishi_2024_sigma_min_bound.

psi_schur_factor(μ::T) where {T}

Compute ψ(μ) from Lemma 1.2 of [8].

For a block triangular matrix H = [I, -N; 0, I] with ‖N‖ = μ, this computes the exact value of ‖H‖ = ‖H⁻¹‖.

The paper's formula (5) is

α = √(½(1 + 1/√(1 + 4/μ²))),      ψ(μ) = √(1 + 2αμ√(1-α²) + α²μ²).

It is evaluated here in the equal form ψ(μ) = (μ + √(μ² + 4))/2: β = α² solves β² − β + 1/(4 + μ²) = 0 (the proof of Lemma 1.2), so α√(1-α²) = 1/√(4 + μ²) and, with t = μ/√(μ² + 4), β = (1 + t)/2, whence ψ² = 1 + 2t + μ²(1 + t)/2 = ((μ + √(μ² + 4))/2)². Formula (5) is the maximum over α of the quantity (6), so an α carrying a rounding error gives a value below ψ; the closed form has no such step.

This replaces the weaker bound 1/(1-μ) used in Oishi 2023, and works for any μ ≥ 0. By (10) of the paper an upper bound of ‖N‖ may be passed for μ.

Returns

An upper bound of ψ(μ).

psi_squared(μ)

An upper bound of ψ(μ)² = 1 + μ²/2 + (μ/2)√(μ² + 4), rounded up.

For S(X) = [I X; 0 I] and μ = ‖X‖₂, ‖S(X)‖₂ = ‖S(X)⁻¹‖₂ = ψ(μ), so ψ(μ)² is the condition number of S(X) in the spectral norm: with the singular value decomposition of X the matrix S(X)*S(X) splits into blocks [1 σ; σ 1 + σ²], whose largest eigenvalue is 1 + σ²/2 + (σ/2)√(σ² + 4), increasing in σ; and S(X)⁻¹ = S(−X). Being increasing, the function may be evaluated at an upper bound of ‖X‖₂.

pi_norm(M::BallMatrix{T}) where {T}

An upper bound of the ℓ² operator norm of every matrix of the ball M, π(M) = √(‖M‖₁‖M‖_∞), from the largest column sum and the largest row sum of the moduli of its entries. It is the bound written π in [8] and used there for the cheap γ of equation (9); the same inequality appears among the immediate bounds of Section 3 of [10].