# Methods

This section describes systems methods implemented in `IntervalMatrices.jl`

.

## Common functions

`IntervalArithmetic.inf`

— Function`inf(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`

— Function`sup(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`

— Function`mid(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`

— Function`diam(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`

— Function`radius(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`

— Function`midpoint_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 dispatch`m`

– (optional, default:`2`

) number of rows`n`

– (optional, default:`m`

) number of columns`rng`

– (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 matrix`m`

– (optional, default:`2`

) number of rows`n`

– (optional, default:`2`

) number of columns`rng`

– (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 matrix`A`

– 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 matrix`S`

– 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, 2] [0, 4]
[-1, 7] [-1, 9]
```

`Base.:⊆`

— Function`⊆(A::AbstractIntervalMatrix, B::AbstractIntervalMatrix)`

Check whether an interval matrix is contained in another interval matrix.

**Input**

`A`

– interval matrix`B`

– 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 matrix`B`

– interval matrix (of the same shape as`A`

)

**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 matrix`B`

– interval matrix (of the same shape as`A`

)

**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`

— Function`hull(A::IntervalMatrix, B::IntervalMatrix)`

Finds the interval hull of two interval matrices. This is equivalent to `∪`

.

**Input**

`A`

– interval matrix`B`

– interval matrix (of the same shape as`A`

)

**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`

— Function`square(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`

— Function`scale(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!`

— Function`scale(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`

— Function`set_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].

**Notes**

- By default,
`:fast`

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!`

— Function`increment!(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"`

– decompose`k = 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`

— Function`increment(pow::IntervalMatrixPower; [algorithm=default_algorithm])`

Increment a matrix power without modifying `pow`

.

**Input**

`pow`

– wrapper of a matrix power`algorithm`

– (optional; default:`default_algorithm`

) algorithm to compute the matrix power; see`increment!`

for available options

**Output**

The next matrix power.

`Base.get`

— Function`get(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`

— Function`base(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`

— Function`index(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`

— Type`Horner <: AbstractExponentiationMethod`

Matrix exponential using Horner's method.

**Fields**

`K`

– number of expansions in the Horner scheme

`IntervalMatrices.ScaleAndSquare`

— Type`ScaleAndSquare <: AbstractExponentiationMethod`

**Fields**

`l`

– scaling-and-squaring order`p`

– order of the approximation

`IntervalMatrices.TaylorOverapproximation`

— Type`TaylorOverapproximation <: AbstractExponentiationMethod`

Matrix exponential overapproximation using a truncated Taylor series.

**Fields**

`p`

– order of the approximation

`IntervalMatrices.TaylorUnderapproximation`

— Type`TaylorUnderapproximation <: AbstractExponentiationMethod`

Matrix exponential underapproximation using a truncated Taylor series.

**Fields**

`p`

– order of the approximation

### Implementations

`IntervalMatrices.exp_overapproximation`

— Function`exp_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 matrix`t`

– exponentiation factor`p`

– 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`

— Function`horner(A::IntervalMatrix{T}, K::Integer; [validate]::Bool=true)`

Compute the matrix exponential using the Horner scheme.

**Input**

`A`

– interval matrix`K`

– number of expansions in the Horner scheme`validate`

– (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`

— Function```
scale_and_square(A::IntervalMatrix{T}, l::Integer, t, p;
[validate]::Bool=true)
```

Compute the matrix exponential using scaling and squaring.

**Input**

`A`

– interval matrix`l`

– scaling-and-squaring order`t`

– non-negative time value`p`

– order of the approximation`validate`

– (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`

— Function`exp_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 matrix`t`

– exponentiation factor`p`

– 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`

— Function`quadratic_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`

— Function`correction_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 matrix`t`

– non-negative time value`p`

– 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`

— Function`input_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 matrix`t`

– non-negative time value`p`

– 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 automonous cars. 2010.

## Norms

`LinearAlgebra.opnorm`

— Function`opnorm(A::IntervalMatrix, p::Real=Inf)`

The matrix norm of an interval matrix.

**Input**

`A`

– interval matrix`p`

– (optional, default:`Inf`

) the class of`p`

-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`

— Function`diam_norm(A::IntervalMatrix, p=Inf)`

Return the diameter norm of the interval matrix.

**Input**

`A`

– interval matrix`p`

– (optional, default:`Inf`

) the`p`

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