Linear systems
This tutorial will show you how to solve linear systems rigorously using IntervalLinearAlgebra.jl
.
Solve interval linear systems
An interval linear system $\mathbf{Ax}=\mathbf{b}$ is a linear system where $\mathbf{A}$ and $\mathbf{b}$ contain intervals. In general, the solution set $\mathbf{x}$ can have a complex non-convex shape and can thus be hard to characterize exactly (see this article for more details). Hence we are interested in finding an interval box containing $\mathbf{x}$. In IntervalLinearAlgebra.jl
, this is achieved through the solve
function, which gives a handy interface to choose the algorithm and preconditioning mechanism. The syntax to call solve is
solve(A, b, method, precondition)
- $A$ is an interval matrix
- $b$ is an interval vector
method
is an optional parameter to choose the algorithm used to solve the interval linear system, see below for more detailsprecondition
is an optional parameter to choose the preconditioning for the problem. More details about preconditoining can be found here
Methods
The supported methods are
Direct solvers
GaussianElimination
HansenBliekRohn
LinearOettliPrager
(requires importing LazySets.jl)
Iterative solvers
LinearKrawczyk
Jacobi
GaussSeidel
NonLinearOettliPrager
(requires importing IntervalConstraintProgramming.jl)
LinearOettliPrager
and NonLinearOettliPrager
are "special" in the sense that they try to exactly characterize the solution set using Oettli-Präger and are not considered in this tutorial. More information about them can be found here. The other solvers return a vector of intervals, representing an interval enclosure of the solution set. If the method is not specified, Gaussian elimination is used by default.
Preconditioning
The supported preconditioning mechanisms are
If preconditioning is not specified, then an heuristic strategy based on the type of matrix and solver is used to choose the preconditioning. The strategy is discussed at the end of the preconditioning tutorial.
Examples
We now demonstrate a few examples using the solve function, these examples are taken from [HOR19].
using IntervalLinearAlgebra
A = [4..6 -1..1 -1..1 -1..1;-1..1 -6.. -4 -1..1 -1..1;-1..1 -1..1 9..11 -1..1;-1..1 -1..1 -1..1 -11.. -9]
4×4 Matrix{Interval{Float64}}:
[4, 6] [-1, 1] [-1, 1] [-1, 1]
[-1, 1] [-6, -4] [-1, 1] [-1, 1]
[-1, 1] [-1, 1] [9, 11] [-1, 1]
[-1, 1] [-1, 1] [-1, 1] [-11, -9]
b = [-2..4, 1..8, -4..10, 2..12]
4-element Vector{Interval{Float64}}:
[-2, 4]
[1, 8]
[-4, 10]
[2, 12]
solve(A, b, HansenBliekRohn())
4-element Vector{Interval{Float64}}:
[-2.50001, 3.10001]
[-3.90001, 1.20001]
[-1.40001, 2.15001]
[-2.35001, 0.600001]
solve(A, b, GaussianElimination())
4-element Vector{Interval{Float64}}:
[-2.60001, 3.10001]
[-3.90001, 1.50001]
[-1.43001, 2.15001]
[-2.35001, 0.600001]
solve(A, b, GaussSeidel())
4-element Vector{Interval{Float64}}:
[-2.60001, 3.10001]
[-3.90001, 1.65001]
[-1.48334, 2.15001]
[-2.35001, 0.794445]
For iterative methods, an additional optional parameter X0
representing an initial guess for the solution's enclosure can be given. If not given, a rough initial enclosure is computed using the enclose
function.
X0 = fill(-5..5, 4)
solve(A, b, GaussSeidel(), InverseMidpoint(), X0)
4-element Vector{Interval{Float64}}:
[-2.60001, 3.10001]
[-3.90001, 1.65001]
[-1.48334, 2.15001]
[-2.35001, 0.794445]
Verify real linear systems
IntervalLinearAlgebra.jl
also offers functionalities to solve real linear systems rigorously. It is of course possible to just convert the real system to an interval system and use the methods described above. In this situation, however, the system will have the property where the diameters of the intervals will be very small (zero or a few floating point units). To solve these kind of systems, it can be more efficient to use the epsilon inflation method [RUM10], especially for bigger matrices. Here is an example
A = [1.0 2;3 4]
2×2 Matrix{Float64}:
1.0 2.0
3.0 4.0
b = [3, 7]
2-element Vector{Int64}:
3
7
the real linear system $Ax=b$ can now be solved rigorously using the epsilon_inflation
function.
x, cert = epsilon_inflation(A, b)
@show cert
x
2-element Vector{Interval{Float64}}:
[0.999999, 1.00001]
[0.999999, 1.00001]
This function returns two values: an interval vector x
and a boolean certificate cert
. If cert==true
then x
is guaranteed to be an enclosure of the real linear system Ax=b
. If cert == false
then the algorithm could not verify that the enclosure is rigorous, i.e. it may or may not contain the true solution.
In the following example the epsilon inflation method returns a non-rigorous bound
A1 = [1..1+1e-16 2;3 4]
x1, cert = epsilon_inflation(A1, b)
@show cert
x1
2-element Vector{Interval{Float64}}:
[0.999999, 1.00001]
[0.999999, 1.00001]
Since the matrix A1
is non-regular (it contains the matrix $\begin{bmatrix}1&2\\3&4\end{bmatrix}$ which is singluar), the solution set is unbounded, hence the algorithm could not prove (rightly) that x1
is an enclosure of the true solution.