Methods
This section describes systems methods implemented in IntervalMatrices.jl
.
Common functions
IntervalArithmetic.inf
— Functioninf(A::IntervalMatrix{T}) where {T}
Return the infimum of an interval matrix A
, which corresponds to taking the element-wise infimum of A
.
Input
A
– interval matrix
Output
A scalar matrix whose coefficients are the infima of each element in A
.
IntervalArithmetic.sup
— Functionsup(A::IntervalMatrix{T}) where {T}
Return the supremum of an interval matrix A
, which corresponds to taking the element-wise supremum of A
.
Input
A
– interval matrix
Output
A scalar matrix whose coefficients are the suprema of each element in A
.
IntervalArithmetic.mid
— Functionmid(A::IntervalMatrix{T}) where {T}
Return the midpoint of an interval matrix A
, which corresponds to taking the element-wise midpoint of A
.
Input
A
– interval matrix
Output
A scalar matrix whose coefficients are the midpoints of each element in A
.
IntervalArithmetic.diam
— Functiondiam(A::IntervalMatrix{T}) where {T}
Return a matrix whose entries describe the diameters of the intervals.
Input
A
– interval matrix
Output
A matrix B
of the same shape as A
such that B[i, j] == diam(A[i, j])
for each i
and j
.
IntervalArithmetic.radius
— Functionradius(A::IntervalMatrix{T}) where {T}
Return the radius of an interval matrix A
, which corresponds to taking the element-wise radius of A
.
Input
A
– interval matrix
Output
A scalar matrix whose coefficients are the radii of each element in A
.
IntervalArithmetic.midpoint_radius
— Functionmidpoint_radius(A::IntervalMatrix{T}) where {T}
Split an interval matrix $A$ into two scalar matrices $C$ and $S$ such that $A = C + [-S, S]$.
Input
A
– interval matrix
Output
A pair (C, S)
such that the entries of C
are the central points and the entries of S
are the (nonnegative) radii of the intervals in A
.
Base.rand
— Functionrand(::Type{IntervalMatrix}, m::Int=2, [n]::Int=m;
N=Float64, rng::AbstractRNG=GLOBAL_RNG)
Return a random interval matrix of the given size and numeric type.
Input
IntervalMatrix
– type, used for dispatchm
– (optional, default:2
) number of rowsn
– (optional, default:m
) number of columnsrng
– (optional, default:GLOBAL_RNG
) random-number generator
Output
An interval matrix of size $m × n$ whose coefficients are normally-distributed intervals of type N
with mean 0
and standard deviation 1
.
Notes
If this function is called with only one argument, it creates a square matrix, because the number of columns defaults to the number of rows.
IntervalMatrices.sample
— Functionsample(A::IntervalMatrix{T}; rng::AbstractRNG=GLOBAL_RNG) where {T}
Return a sample of the given random interval matrix.
Input
A
– interval matrixm
– (optional, default:2
) number of rowsn
– (optional, default:2
) number of columnsrng
– (optional, default:GLOBAL_RNG
) random-number generator
Output
An interval matrix of size $m × n$ whose coefficients are normally-distributed intervals of type N
with mean 0
and standard deviation 1
.
Base.:∈
— Function∈(M::AbstractMatrix, A::AbstractIntervalMatrix)
Check whether a concrete matrix is an instance of an interval matrix.
Input
M
– concrete matrixA
– interval matrix
Output
true
iff M
is an instance of A
Algorithm
We check for each entry in M
whether it belongs to the corresponding interval in A
.
IntervalArithmetic.:±
— Function±(C::MT, S::MT) where {T, MT<:AbstractMatrix{T}}
Return an interval matrix such that the center and radius of the intervals is given by the matrices C
and S
respectively.
Input
C
– center matrixS
– radii matrix
Output
An interval matrix M
such that M[i, j]
corresponds to the interval whose center is C[i, j]
and whose radius is S[i, j]
, for each i
and j
. That is, $M = C + [-S, S]$.
Notes
The radii matrix should be nonnegative, i.e. S[i, j] ≥ 0
for each i
and j
.
Examples
julia> [1 2; 3 4] ± [1 2; 4 5]
2×2 IntervalMatrix{Float64, Interval{Float64}, Matrix{Interval{Float64}}}:
[0.0, 2.0] [0.0, 4.0]
[-1.0, 7.0] [-1.0, 9.0]
Base.:⊆
— Function⊆(A::AbstractIntervalMatrix, B::AbstractIntervalMatrix)
Check whether an interval matrix is contained in another interval matrix.
Input
A
– interval matrixB
– interval matrix
Output
true
iff A[i, j] ⊆ B[i, j]
for all i, j
.
Base.:∩
— Function∩(A::IntervalMatrix, B::IntervalMatrix)
Intersect two interval matrices.
Input
A
– interval matrixB
– interval matrix (of the same shape asA
)
Output
A new matrix C
of the same shape as A
such that C[i, j] = A[i, j] ∩ B[i, j]
for each i
and j
.
Base.:∪
— Function∪(A::IntervalMatrix, B::IntervalMatrix)
Finds the interval union (hull) of two interval matrices. This is equivalent to hull
.
Input
A
– interval matrixB
– interval matrix (of the same shape asA
)
Output
A new matrix C
of the same shape as A
such that C[i, j] = A[i, j] ∪ B[i, j]
for each i
and j
.
IntervalArithmetic.hull
— Functionhull(A::IntervalMatrix, B::IntervalMatrix)
Finds the interval hull of two interval matrices. This is equivalent to ∪
.
Input
A
– interval matrixB
– interval matrix (of the same shape asA
)
Output
A new matrix C
of the same shape as A
such that C[i, j] = hull(A[i, j], B[i, j])
for each i
and j
.
Arithmetic
IntervalMatrices.square
— Functionsquare(A::IntervalMatrix)
Compute the square of an interval matrix.
Input
A
– interval matrix
Output
An interval matrix equivalent to A * A
.
Algorithm
We follow [1, Section 6].
[1] Kosheleva, Kreinovich, Mayer, Nguyen. Computing the cube of an interval matrix is NP-hard. SAC 2005.
IntervalMatrices.scale
— Functionscale(A::IntervalMatrix{T}, α::T) where {T}
Return a new interval matrix whose entries are scaled by the given factor.
Input
A
– interval matrixα
– scaling factor
Output
A new matrix B
of the same shape as A
such that B[i, j] = α*A[i, j]
for each i
and j
.
Notes
See scale!
for the in-place version of this function.
IntervalMatrices.scale!
— Functionscale!(A::IntervalMatrix{T}, α::T) where {T}
Modifies the given interval matrix, scaling its entries by the given factor.
Input
A
– interval matrixα
– scaling factor
Output
The matrix A
such that for each i
and j
, the new value of A[i, j]
is α*A[i, j]
.
Notes
This is the in-place version of scale
.
IntervalMatrices.set_multiplication_mode
— Functionset_multiplication_mode(multype)
Sets the algorithm used to perform matrix multiplication with interval matrices.
Input
multype
– symbol describing the algorithm used:slow
– uses traditional matrix multiplication algorithm.:fast
– computes an enclosure of the matrix product using the midpoint-radius notation of the matrix [RUM10].
:fast
support was removed in IntervalArithmetic
v0.22.
Notes
- By default,
:slow
is used. - Using
fast
is generally significantly faster, but it may return larger intervals, especially if midpoint and radius have the same order of magnitude (50% overestimate at most) [RUM99].
Matrix power
IntervalMatrices.increment!
— Functionincrement!(pow::IntervalMatrixPower; [algorithm=default_algorithm])
Increment a matrix power in-place (i.e., storing the result in pow
).
Input
pow
– wrapper of a matrix power (modified in this function)algorithm
– (optional; default:default_algorithm
) algorithm to compute the matrix power; available options:"multiply"
– fast computation using*
from the previous result"power"
– recomputation using^
"decompose_binary"
– decomposek = 2a + b
"intersect"
– combination of"multiply"
/"power"
/"decompose_binary"
Output
The next matrix power, reflected in the modified wrapper.
Notes
Independent of "algorithm"
, if the index is a power of two, we compute the exact result using squaring.
IntervalMatrices.increment
— Functionincrement(pow::IntervalMatrixPower; [algorithm=default_algorithm])
Increment a matrix power without modifying pow
.
Input
pow
– wrapper of a matrix poweralgorithm
– (optional; default:default_algorithm
) algorithm to compute the matrix power; seeincrement!
for available options
Output
The next matrix power.
IntervalMatrices.matrix
— Functionmatrix(pow::IntervalMatrixPower)
Return the matrix represented by a wrapper of a matrix power.
Input
pow
– wrapper of a matrix power
Output
The matrix power represented by the wrapper.
IntervalMatrices.base
— Functionbase(pow::IntervalMatrixPower)
Return the original matrix represented by a wrapper of a matrix power.
Input
pow
– wrapper of a matrix power
Output
The matrix $M$ being the basis of the matrix power $M^k$ represented by the wrapper.
IntervalMatrices.index
— Functionindex(pow::IntervalMatrixPower)
Return the current index of the wrapper of a matrix power.
Input
pow
– wrapper of a matrix power
Output
The index k
of the wrapper representing $M^k$.
Matrix exponential
Algorithms
IntervalMatrices.Horner
— TypeHorner <: AbstractExponentiationMethod
Matrix exponential using Horner's method.
Fields
K
– number of expansions in the Horner scheme
IntervalMatrices.ScaleAndSquare
— TypeScaleAndSquare <: AbstractExponentiationMethod
Fields
l
– scaling-and-squaring orderp
– order of the approximation
IntervalMatrices.TaylorOverapproximation
— TypeTaylorOverapproximation <: AbstractExponentiationMethod
Matrix exponential overapproximation using a truncated Taylor series.
Fields
p
– order of the approximation
IntervalMatrices.TaylorUnderapproximation
— TypeTaylorUnderapproximation <: AbstractExponentiationMethod
Matrix exponential underapproximation using a truncated Taylor series.
Fields
p
– order of the approximation
Implementations
IntervalMatrices.exp_overapproximation
— Functionexp_overapproximation(A::IntervalMatrix{T}, t, p)
Overapproximation of the exponential of an interval matrix, exp(A*t)
, using a truncated Taylor series.
Input
A
– interval matrixt
– exponentiation factorp
– order of the approximation
Output
A matrix enclosure of exp(A*t)
, i.e. an interval matrix M = (m_{ij})
such that [exp(A*t)]_{ij} ⊆ m_{ij}
.
Algorithm
See Theorem 1 in Reachability Analysis of Linear Systems with Uncertain Parameters and Inputs by M. Althoff, O. Stursberg, M. Buss.
IntervalMatrices.horner
— Functionhorner(A::IntervalMatrix{T}, K::Integer; [validate]::Bool=true)
Compute the matrix exponential using the Horner scheme.
Input
A
– interval matrixK
– number of expansions in the Horner schemevalidate
– (optional; default:true
) option to validate the precondition of the algorithm
Algorithm
We use the algorithm in [1, Section 4.2].
[1] Goldsztejn, Alexandre, Arnold Neumaier. "On the exponentiation of interval matrices". Reliable Computing. 2014.
IntervalMatrices.scale_and_square
— Functionscale_and_square(A::IntervalMatrix{T}, l::Integer, t, p;
[validate]::Bool=true)
Compute the matrix exponential using scaling and squaring.
Input
A
– interval matrixl
– scaling-and-squaring ordert
– non-negative time valuep
– order of the approximationvalidate
– (optional; default:true
) option to validate the precondition of the algorithm
Algorithm
We use the algorithm in [1, Section 4.3], which first scales A
by factor $2^{-l}$, computes the matrix exponential for the scaled matrix, and then squares the result $l$ times.
\[ \exp(A * 2^{-l})^{2^l}\]
[1] Goldsztejn, Alexandre, Arnold Neumaier. "On the exponentiation of interval matrices". Reliable Computing. 2014.
IntervalMatrices.exp_underapproximation
— Functionexp_underapproximation(A::IntervalMatrix{T}, t, p) where {T}
Underapproximation of the exponential of an interval matrix, exp(A*t)
, using a truncated Taylor series expansion.
Input
A
– interval matrixt
– exponentiation factorp
– order of the approximation
Output
An underapproximation of exp(A*t)
, i.e. an interval matrix M = (m_{ij})
such that m_{ij} ⊆ [exp(A*t)]_{ij}
.
Algorithm
See Theorem 2 in Reachability Analysis of Linear Systems with Uncertain Parameters and Inputs by M. Althoff, O. Stursberg, M. Buss.
Finite expansions
IntervalMatrices.quadratic_expansion
— Functionquadratic_expansion(A::IntervalMatrix, α::Real, β::Real)
Compute the quadratic expansion of an interval matrix, $αA + βA^2$, using interval arithmetic.
Input
A
– interval matrixα
– linear coefficientβ
– quadratic coefficient
Output
An interval matrix that encloses $B := αA + βA^2$.
Algorithm
This a variation of the algorithm in [1, Section 6]. If $A = (aᵢⱼ)$ and $B := αA + βA^2 = (bᵢⱼ)$, the idea is to compute each $bᵢⱼ$ by factoring out repeated expressions (thus the term single-use expressions).
First, let $i = j$. In this case,
\[bⱼⱼ = β\sum_\{k, k ≠ j} a_{jk} a_{kj} + (α + βa_{jj}) a_{jj}.\]
Now consider $i ≠ j$. Then,
\[bᵢⱼ = β\sum_\{k, k ≠ i, k ≠ j} a_{ik} a_{kj} + (α + βa_{ii} + βa_{jj}) a_{ij}.\]
[1] Kosheleva, Kreinovich, Mayer, Nguyen. Computing the cube of an interval matrix is NP-hard. SAC 2005.
Correction terms
IntervalMatrices.correction_hull
— Functioncorrection_hull(A::IntervalMatrix{T}, t, p) where {T}
Compute the correction term for the convex hull of a point and its linear map with an interval matrix in order to contain all trajectories of a linear system.
Input
A
– interval matrixt
– non-negative time valuep
– order of the approximation
Output
An interval matrix representing the correction term.
Algorithm
See Theorem 3 in [1].
[1] M. Althoff, O. Stursberg, M. Buss. Reachability Analysis of Linear Systems with Uncertain Parameters and Inputs. CDC 2007.
IntervalMatrices.input_correction
— Functioninput_correction(A::IntervalMatrix{T}, t, p) where {T}
Compute the input correction matrix for discretizing an inhomogeneous affine dynamical system with an interval matrix and an input domain not containing the origin.
Input
A
– interval matrixt
– non-negative time valuep
– order of the Taylor approximation
Output
An interval matrix representing the correction matrix.
Algorithm
See Proposition 3.4 in [1].
[1] M. Althoff. Reachability analysis and its application to the safety assessment of autonomous cars. 2010.
Norms
LinearAlgebra.opnorm
— Functionopnorm(A::IntervalMatrix, p::Real=Inf)
The matrix norm of an interval matrix.
Input
A
– interval matrixp
– (optional, default:Inf
) the class ofp
-norm
Notes
The matrix $p$-norm of an interval matrix $A$ is defined as
\[ ‖A‖_p := ‖\max(|\text{inf}(A)|, |\text{sup}(A)|)‖_p\]
where $\max$ and $|·|$ are taken elementwise.
IntervalMatrices.diam_norm
— Functiondiam_norm(A::IntervalMatrix, p=Inf)
Return the diameter norm of the interval matrix.
Input
A
– interval matrixp
– (optional, default:Inf
) thep
-norm used; valid options are:1
,2
,Inf
Output
The operator norm, in the p
-norm, of the scalar matrix obtained by taking the element-wise diam
function, where diam(x) := sup(x) - inf(x)
for an interval x
.
Notes
This function gives a measure of the width of the interval matrix.