added algorithms to bound eigenvalues of interval matrices (#77)

* added algorithms to bound eigenvalues of interval matrices

* updated function signature

* added docstring to documentation

Author: lucaferranti

docs/make.jl

@@ -27,10 +27,11 @@ makedocs(;
         ],
         "API" => [
             "Interval matrices classification" => "api/classify.md",
-            "solver interface" => "api/solve.md",
+            "Solver interface" => "api/solve.md",
             "Interval linear systems" => "api/algorithms.md",
             "Preconditioners" => "api/precondition.md",
             "Verified real linear systems" => "api/epsilon_inflation.md",
+            "Eigenvalues" => "api/eigenvalues.md",
             "Miscellaneous" => "api/misc.md"
         ],
         "References" => "references.md",

docs/src/api/eigenvalues.md

@@ -0,0 +1,10 @@
+```@index
+Pages = ["eigenvalues.md"]
+```
+
+```@autodocs
+Modules=[IntervalLinearAlgebra]
+Pages=["interval_eigenvalues.jl"]
+Private=false
+```
+

docs/src/references.md

@@ -1,5 +1,35 @@
 # [References](@id all_ref)
 
+#### [HLA13] 
+
+```@raw html
+<ul><li>
+```
+M. Hladík. Bounds on eigenvalues of real and complex interval matrices. Appl. Math. Comput., 219(10):5584–5591, 2013.
+```@raw html
+<li style="list-style: none"><details>
+<summary>bibtex</summary>
+```
+```
+@article{Hla2013a,
+ author = "Milan Hlad\'{\i}k",
+ title = "Bounds on eigenvalues of real and complex interval matrices",
+ journal = "Appl. Math. Comput.",
+ fjournal = "Applied Mathematics and Computation",
+ volume = "219",
+ number = "10",
+ pages = "5584-5591",
+ year = "2013",
+ issn = "0096-3003",
+ doi = "10.1016/j.amc.2012.11.075",
+}
+```
+```@raw html
+</details></li></ul>
+```
+---
+
+
 #### [HOR19] 
 
 ```@raw html

src/IntervalLinearAlgebra.jl

@@ -1,11 +1,11 @@
 module IntervalLinearAlgebra
 
 using StaticArrays, Requires, Reexport
+using LinearAlgebra: checksquare
 
 import Base: *
 import CommonSolve: solve
 import IntervalArithmetic: mid
-using LinearAlgebra: checksquare
 
 function  __init__()
     @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("linear_systems/oettli_nonlinear.jl")
@@ -23,7 +23,8 @@ export
     solve, enclose, epsilon_inflation,
     comparison_matrix, interval_norm, interval_isapprox, Orthants,
     is_H_matrix, is_strongly_regular, is_strictly_diagonally_dominant, is_Z_matrix, is_M_matrix,
-    rref
+    rref,
+    eigenbox
 
 
 include("linear_systems/enclosures.jl")
@@ -35,4 +36,6 @@ include("multiplication.jl")
 include("utils.jl")
 include("classify.jl")
 include("rref.jl")
+
+include("eigenvalues/interval_eigenvalues.jl")
 end

src/eigenvalues/interval_eigenvalues.jl

@@ -0,0 +1,70 @@
+"""
+    eigenbox(A)
+
+Returns an enclosure of all the eigenvalues of `A`. If `A` is symmetric, than the
+output is a real interval, otherwise it is a complex interval.
+
+### Algorithm
+
+The algorithms used by the function are described in [[HLA13]](@ref).
+
+### Notes
+
+The enclosure is not rigorous, meaning that the real eigenvalue problems solved internally
+utilize normal floating point computations.
+
+### Examples
+
+```jldoctest
+julia> A = [0 -1 -1;2 -1.399.. -0.001 0;1 0.5 -1]
+3×3 Matrix{Interval{Float64}}:
+ [0, 0]  [-1, -1]                       [-1, -1]
+ [2, 2]       [-1.39901, -0.000999999]    [0, 0]
+ [1, 1]        [0.5, 0.5]               [-1, -1]
+
+julia> eigenbox(A)
+[-1.90679, 0.970154] + [-2.51903, 2.51903]im
+```
+"""
+function eigenbox(A::Symmetric{Interval{T}, Matrix{Interval{T}}}) where {T}
+
+    AΔ = Symmetric(radius.(A))
+    Ac = Symmetric(mid.(A))
+
+    ρ = eigmax(AΔ)
+    λmax = eigmax(Ac)
+    λmin = eigmin(Ac)
+    return Interval(λmin - ρ, λmax + ρ)
+
+end
+
+function eigenbox(A::AbstractMatrix{Interval{T}}) where {T}
+
+    λ = eigenbox(Symmetric(0.5*(A + A')))
+
+    n = checksquare(A)
+    μ = eigenbox(Symmetric([zeros(n, n) 0.5*(A - A');
+                            0.5*(A' - A) zeros(n, n)]))
+
+    return λ + μ*im
+end
+
+function eigenbox(M::AbstractMatrix{Complex{Interval{T}}}) where {T}
+    A = real.(M)
+    B = imag.(M)
+    λ = eigenbox(Symmetric(0.5*[A+A' B'-B;
+                                B-B' A+A']))
+
+    μ = eigenbox(Symmetric(0.5*[B+B' A-A';
+                                A'-A B+B']))
+
+    return λ + μ*im
+end
+
+
+function eigenbox(M::Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}}) where T
+    A = real(M)
+    B = imag(M)
+    return eigenbox(Symmetric([A B';B A]))
+
+end

test/runtests.jl

@@ -7,6 +7,9 @@ include("test_classify.jl")
 include("test_multiplication.jl")
 include("test_utils.jl")
 
+
+include("test_eigenvalues/test_interval_eigenvalues.jl")
+
 include("test_solvers/test_enclosures.jl")
 include("test_solvers/test_epsilon_inflation.jl")
 include("test_solvers/test_oettli_prager.jl")

test/test_eigenvalues/test_interval_eigenvalues.jl

@@ -0,0 +1,36 @@
+@testset "Eigenvalues of interval matrices" begin
+
+    # symmetrix matrix
+    A = Symmetric([-1 0 -1..1;
+         0 -1 -1..1;
+         -1..1 -1..1 0.1])
+
+    ev = eigenbox(A)
+    @test interval_isapprox(ev, -2.4143..1.5143; atol=1e-3)
+
+    # real matrix
+    A = [-3.. -2 4..5 4..6 -1..1.5;
+         -4.. -3 -4.. -3 -4.. -3 1..2;
+         -5.. -4 2..3 -5.. -4 -1..0;
+         -1..0.1 0..1 1..2 -4..2.5]
+
+    ev = eigenbox(A)
+    @test interval_isapprox(real(ev), -8.8221..3.4408; atol=1e-3)
+    @test interval_isapprox(imag(ev), -10.7497..10.7497; atol=1e-3)
+
+
+    # hermitian matrix
+    A = Hermitian([1..2 (5..9)+(2..5)*im (3..5)+(2..4)im;
+        (5..9)+(-5.. -2)*im 2..3 (7..8)+(6..10)im;
+        (3..5)+(-4.. -2)*im (7..8)+(-10.. -6)*im 3..4])
+
+    ev = eigenbox(A)
+    @test interval_isapprox(ev, -15.4447..24.3359; atol=1e-3)
+
+    # complex matrix
+    A = [(1..2)+(3..4)*im 3..4;1+(2..3)*im 4..5]
+
+    ev = eigenbox(A)
+    @test interval_isapprox(real(ev), -1.28812..7.28812; atol=1e-3)
+    @test interval_isapprox(imag(ev), -2.04649..5.54649; atol=1e-3)
+end